Can I rely on eigenvectors_sym returning orthogonal eigenvectors?

I need to facotrise a real symmetric matrix A as QDQ’. Since A is symmetric, the eigendecomposition can return such factorisation (but it may in principle return a different one TDT⁻¹ such that TT’≠I). Neither the Stan manual for eigenvectors_sym nor the corresponding eigen function mention that the eigenvectors returned are orthogonal (or, equivalently, that the operation coincides with the SVD).
I looked into the eigen documentation and it looks like they use the symmetric QR algorithm that, as far as I know, does return the facotrisation I want. Testing Stan with a few random matrices seems to confirm that the eigenvectors returned are indeed orthogonal.

Since this is not documented in Stan, can I rely on this behaviour or not? If so, should it be also written in the documentation?

Thank you in advance.

The symmetric QR algorithm yields orthogonal eigenvectors. The Stan documentation tends not to go into that much detail about computations in external libraries.

Thank you for your answer.

I understand your point but I wonder if this would be valuable information (also because it’s basically also an SVD, which is otherwise missing). I wouldn’t consider it a small implementation detail but an interesting feature that could be useful for some users.

Also, I see the documentation as a contract with the user. Unless it’s written in the documentation, I feel uneasy relying on this behaviour in case one day the implementation is changed. To be on the safe side, I can add a test case in the transformed data block and throw a reject if the result if far from orthogonal.

If you want to add it to the doc, go ahead. But I wouldn’t sign a contract presented by our documentation. If anything, the code is the contract.

I think this should be documented. It’s not fair to users (developers yes, not users) to expect them to look at the underlying code. I wouldn’t sign a contract presented by our documentation, but we should be trying to make the documentation sufficient.

1 Like

@lciti Would you be willing to open an issue at https://github.com/stan-dev/docs about this?

Sure. By the way, I meant “contract” in loose terms, in the sense that if it’s written in the documentation it’s less likely to change and go unnoticed for long.

1 Like

I think there is a different level of scrutiny for Stan code vs. libraries. It is not as if eigenvectors_sym is orthogonalizing the eigenvectors. If Eigen chose to change their algorithm (which they wouldn’t) to not make the eigenvectors orthogonal anymore, we probably wouldn’t notice. I don’t think we really even know what SUNDIALS does, etc.

@bgoodri, I see your point.

I have opened an issue, I leave it to the developers to decide.

Ah ok, @bgoodri I see what you mean. If we were relying on that behavior from eigen then I think it should be documented, but in this case it sounds like we’re not relying on eigen making the eigenvectors orthogonal, that just happens to be the case.

I have opened an issue for Eigen, let’s see what they say. If they make the documentation explicit, then I think it could propagate to Stan too. In a sense, they would be signing that contract ;-)

I think Eigen should have it in their docs but I would also understand if they didn’t.

I went and got Strang out (4th international edition) for a definitive answer and section 5.5 (complex matrices), subsection “Hermitian Matrices” , Property 3 is:

Two eigenvectors of a real symmetric matrix or a Hermitian matrix, if they come from different eigenvalues, are orthogonal to one another.

The thing that doesn’t seem to be constrained is whether the entire basis is orthonormal (orthogonal, plus unit length).

@Bob_Carpenter
The unit length is easy to fix. The main issue is that the orthogonality only applies across eigenspaces corresponding to distinct eigenvalues. Within the eigenspace corresponding to a single eigenvalue with multiplicity larger than one, there is no guarantee that the eigenvectors generated by a generic eigendecomposition algorithm are orthogonal. See my example here.

1 Like

Cool. Thanks!