Can I rely on eigenvectors_sym returning orthogonal eigenvectors?

[lciti]

Summary:

The function eigenvectors_sym returns the factorisation QDQ’ for a real symmetric matrix. The documentation describes it as a generic eigendecomposition without specifying that the eigenvectors are orthogonal. This factorisation also coincides with the SVD. This is useful information that I think could be added to the manual.

Description:

Copy-pasting from this post on discourse:

I need to factorise 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 factorisation 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?

I have opened a corresponding issue for Eigen.

Current Version:

v2.18.0

[bob-carpenter]

can I rely on this behaviour or not?

Yes, we're calling the Eigen algorithm. And that's what we tested so behavior shouldn't change in the future.

It's very easy to add a paragraph to the doc if you'd like to create a PR.

[jgabry]

I guess the issue is that the eigenvectors being orthogonal isn't documented by Eigen, so we're not even sure if we can count on that. In theory Eigen could change their underlying implementation without contradicting their documentation, but then our doc would be wrong right?

[lciti]

@jgabry I agree with you. I have opened a corresponding issue for Eigen. If they are reasonably sure they will not change the underlying implementation (symmetric QR algorithm) then they may be happy to make their docs more specific (which may benefit other Eigen users as well). In that case we could propagate that information to the Stan documentation.

[bob-carpenter]

@lciti: Document away. I forgot we were talking about eigenvectors_sym. When the matrix being decomposed is real and symmetric, the eigenvectors are guaranteed to be orthogonal.

I've been working on a non-symmetric eigendecomposition which operates in the complex domain, but that won't change how our symmetric eigendecomposition function works.

[jgabry]

@lciti Thanks for checking with the Eigen folks on this.

When the matrix being decomposed is real and symmetric, the eigenvectors are guaranteed to be orthogonal.

Oh, in that case then it does sound like orthogonality is guaranteed and so yeah documenting it would be great.

[lciti]

@bob-carpenter

When the matrix being decomposed is real and symmetric, the eigenvectors are guaranteed to be orthogonal.

That's not necessarily the case. The eigendecomposition is not unique and while for real symmetric matrices an eigendecomposition with orthogonal eigenvectors always exists, there are other eigendecompositions that are not.

For example, the matrix:

can (in principle) be decomposed as:

(the first and third matrices are one the inverse of the other but not one the transpose of the other).

In practice, for real symmetric matrices, an efficient algorithm to find the eigendecomposition is the symmetric QR algorithm which does return orthogonal eigenvectors. This is the implementation (currently) used by Eigen but unfortunately their documentation does not "promise" that their function will always return a set of orthogonal eigenvectors. @bgoodri 's objection in the discourse discussion is that in principle Eigen could one day change their implementation and stop returning orthogonal eigenvectors. I think this is quite hypothetical but I see his point.

What do you recommend? I think in principle we could even: a) commit to always return orthogonal eigenvectors; b) add this to the documentation; c) add a specific test in the testsuite that verifies this is the case; d) deal with the problem of returning orthogonal eigenvectors on the Stan side in the very unlikely event that Eigen's SelfAdjointEigenSolver will one day change implementation.

[bob-carpenter]

I recommend working it out with @bgoodri.

[bgoodri]

I think if Eigen changed their algorithm to not produce orthogonal eigenvectors (and we knew about it), we would have to change the Stan function to orthogonalize them. The stat literature mostly only does eigendecompositions with symmetric matrices and I've never seen a stat application in that case want non-orthogonal eigenvectors. Still, I think it is up to Eigen to document the properties of the algorithms they are currently using..

On Wed, Jul 29, 2020 at 11:54 AM Bob Carpenter notifications@github.com wrote:

I recommend working it out with @bgoodri https://github.com/bgoodri.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/docs/issues/235#issuecomment-665748061, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZ2XKS2V3KLXWM7ATJ34DDR6BA3NANCNFSM4PG7REYQ .

[lciti]

(and we knew about it)

We could create test cases to verify that Eigen still produces orthogonal eigenvectors. If we test Eigen's results with matrices of different sizes and structures, it's unlikely the test will pass if the new algorithm stops reliably producing orthogonal eigenvectors. We will need to create test matrices with eigenvalues with multiplicity larger than one because these are the ones that could lead to non-orthogonal eigenvectors (if the eigenvalues are distinct, the corresponding eigenvectors are guaranteed to be orthogonal).

we would have to change the Stan function to orthogonalize them

That's right. We could still call Eigen and only orthogonalise the eigenvectors corresponding to eigenvalues with multiplicity larger than one (we would need to determine when two eigenvalues are close enough to be considered identical).

[lciti]

I have just found a possible test matrix for which a standard eigendecomposition algorithm may find non-orthogonal eigenvectors. This is in Python but hopefully sufficiently readable for everyone:

    In [1]: import numpy as np

    In [2]: M = np.array([[1.3, -.4, -.2], [-.4, 1.3, .2], [-.2, .2, 1.]])

    In [3]: λ, v = np.linalg.eig(M)

    In [4]: print(v.T @ v)  # print inner product of pairs of eigenvectors (should give the identity matrix)
    [[1.00000000e+00 4.16333634e-17 1.93987781e-01]
     [4.16333634e-17 1.00000000e+00 0.00000000e+00]
     [1.93987781e-01 0.00000000e+00 1.00000000e+00]]

    In [5]: λ, v = np.linalg.eigh(M)

    In [6]: print(v.T @ v)  # print inner product of pairs of eigenvectors (should give the identity matrix)
    [[ 1.00000000e+00 -5.60938093e-17 -1.09959626e-16]
     [-5.60938093e-17  1.00000000e+00 -2.63176746e-17]
     [-1.09959626e-16 -2.63176746e-17  1.00000000e+00]]

Note: np.linalg.eig is a standard eigendecomposition algorithm not assuming symmetry (uses LAPACK's _geev under the hood), while np.linalg.eigh is for real symmetric or complex Hermitian matrices (uses LAPACK's _syevd under the hood).

[lciti]

The eigen folks addressed this issue and updated their documentation (through merge request !439 (merged) and commit b8502a9d). I think it's safe to reflect this into Stan's documentation and say that eigenvectors_sym returns orthonormal eigenvectors for a real symmetric matrix. This is potentially useful because eigenvectors_sym and eigenvalues_sym used together provide a decomposition QDQ’, which is an SVD of the matrix.

[bob-carpenter]

Thanks for following up @Iciti. Feel free to make a PR---it's really easy with doc as you can just do the edit without having to pull down a copy of the repo. Otherwise, I'm sure we'll get to it relatively soon---I need to make a pass over all the doc comments in the coming months.