Wrong eigenvalues order in OrderedEigenDecomposition

[afossa]

If the input matrix has null eigenvalues, the order in which they are returned by OrderedEigenDecomposition might be wrong. As an example, with the input matrix

Array2DRowRealMatrix{{81.0,63.0,55.0,49.0,0.0,0.0},{63.0,82.0,80.0,69.0,0.0,0.0},{55.0,80.0,92.0,75.0,0.0,0.0},{49.0,69.0,75.0,73.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0,0.0}}

extracting the main diagonal of OrderedEigenDecomposition#getD() returns

[0.0, 4.793253233672134, 7.429392849462756, 36.43053556404571, 279.34681835281935, 0.0]

while it should return

[0.0, 0.0, 4.793253233672134, 7.429392849462756, 36.43053556404571, 279.34681835281935]

I acknowledge that if more than one eigenvalue is null, then there is an ambiguity in the order of the corresponding eigenvectors.

It would also be nice to expose the epsilon parameter for internal checks in the class constructor, and have the possibility to choose among ascending/descending order (or is a descending order guaranteed by the "standard" EigenDecomposition?).

I would also expect the arrays realEigenvalues and imagEigenvalues to be sorted in the same way D is.

[axkr]

I would also expect the arrays realEigenvalues and imagEigenvalues to be sorted in the same way D is.

You can define your own comparator.
See #248 (and #249 ?)

[afossa]

This is true for OrderedComplexEigenDecomposition, but not for the "real" equivalent OrderedEigenDecomposition, which is the one I am interested in.

[maisonobe]

This is much more difficult than I first thought.
The problem comes for the handling of non-symmetric matrices. In this case, there are conjugate eigenvalues and conjugate eigenvectors, and they are organized in pairs. The matrix $D$ is not diagonal anymore but has 2x2 blocks holding the real and imaginary parts of the conjugate eigenvalues pairs. The matrices $D$ and $V$ in this form are still both real matrices and the defining expression $A V = V D$ still holds and if $V$ is inversible then $A=V D V^{-1}$. However, the columns of $V$ do not represent the eigenvectors individually, one should combine two columns to get one complex eigenvector (two different combinations to get the two conjugate vectors from the two real columns).

This becomes difficult when reordering. Look for example at the matrix (extracted from one of the test cases), written as a complex decomposition:

$$\begin{bmatrix}3&-2&0\\4&-1&0\\1&1&1\end{bmatrix} = \begin{bmatrix}-2+4i&-2-4i&0\\2+6i&2-6i&0\\5&5&1\end{bmatrix} \begin{bmatrix}1+2i&0&0\\0&1-2i&0\\0&0&1\end{bmatrix} \begin{bmatrix}\frac{-3-i}{20}&\frac{2-i}{20}&0\\\frac{-3+i}{20}&\frac{2+i}{20}&0\\\frac{3}{2}&-1&1\end{bmatrix}$$

One clearly sees that the first two columns of the $V$ matrix are conjugate vectors, that the first two eigenvalues $\lambda_1 = 1+2i$ and $\lambda_2=1-2i$ are conjugate to each other and that the first two rows of the third matrix (which is $V^{-1}$) are also conjugate to each other.

The decomposition using real matrices (which is what EigenDecomposition does) and a non-diagonal $D$ matrix is:

$$\begin{bmatrix}3&-2&0\\4&-1&0\\1&1&1\end{bmatrix} = \begin{bmatrix}0&\sqrt{\frac{8}{17}}&0\\\sqrt{\frac{8}{17}}&\sqrt{\frac{8}{17}}&0\\\sqrt{\frac{8}{17}}&-\sqrt{\frac{2}{17}}&\frac{\sqrt{17}}{2}\end{bmatrix} \begin{bmatrix}1&2&0\\-2&1&0\\0&0&1\end{bmatrix} \begin{bmatrix}-\sqrt{\frac{17}{8}}&\sqrt{\frac{17}{8}}&0\\\sqrt{\frac{17}{8}}&0&0\\\frac{3}{\sqrt{17}}&-\frac{2}{\sqrt{17}}&\frac{2}{\sqrt{17}}\end{bmatrix}$$

We see that the two conjugate complex eigenvalues appear in the real decomposition as a 2x2 block in the upper left of the $D$ matrix. We also see that the $V$ matrix does not really contain the eigenvectors themselves, but something different.

If we were to sort the eigenvalues using the ComplexComparator that uses real part as primary order and imaginary part as secondary order, then the real eigenvalue $\lambda_3=1$ would be inserted between the other two values, since with this ordering we have $\lambda_2<\lambda_3<\lambda_1$. This would lead to some weird ordering of the columns and rebuilding the complex eigenvectors would imply keeping track of that.

Is this what we want?

[maisonobe]

I have fixed this in master branch.
It appeared to me that OrderedEigenDecomposition was deeply flawed and that EigenDecomposition itself was not really good. Indeed, the symmetric versus non-symmetric cases are really different:

• they are not handled the same way right from the beginning (using Tridiagonal transformations in the symmetric case and Schur transformation in the non-symmetric case),
• they do not produce similar results (real eigenvalues and eigenvectors in the symmetric case, complex eigenvalues and eigenvectors in the non-symmetric case.
• they do not produce similar middle matrix (diagonal matrix in the symmetric case, block-diagonal matrix in the non-symmetric case)
• the columns of the left matrix have slightly different semantics (eigen vectors in the symmetric case, real part and plus or minus the imaginary part in the non-symmetric case)

So I split these two different mathematical problems in two different classes: EigenDecompositionSymmetric and EigenDecompositionNonSymmetric. Of course the non-symmetric case can be fed with a symmetric matrix, and the complex values generated will have zero imaginary parts, but they will be handled as complex eigenvalues/eigenvectors.

With the split performed, ordering becomes natural in the symmetric case where all eigenvalues are real numbers, so I just added a boolean so users can select increasing/decreasing order for eigenvalues.

I did not provide any ordering capabilities in the non-symmetric case as it is really awkward. In fact, as per the example above, one can insert one real eigenvalue between two conjugate eigenvalues and in this case the D matrix is even farther from a diagonal matrix, having non-zero elements up to the corners, and the real and imaginary parts of the eigenvectors would be separated. This means that the real matrices V and D are not really useful anymore. So in this case, I would recommend just converting the initial matrix to complex form and use OrderedComplexEigenDecomposition.