Strange behaviour in EigenDecomposition

[QuentinRhone]

Hello all, I am using EigenDecomposition class these days and sometimes the diagonalized matrix I get in return seems wrong. A typical case is when I set this matrix : {{23473.684554963584, 4273.093076392109}, {4273.093076392048, 4462.13956661408}}

I know the matrix is not perfectly symmetric (yet close to double precision) but I got this matrix from a AMtranspose(A) computation so I can not get a more symmetric matrix. The algorithm then computes complex Eigen values matrix : {{13967.9120607888, 10422.0456317615}, {-10422.0456317615, 13967.9120607888}}

I have checked with matlab and with a home-made 2*2 matrix diagonalizer and I get with both solutions the same real Eigen values : 24389.95769255035 and 3545.86642902732 which have nothing to see with the complex ones (in contradiction with the theoretical unicity of the Eigen values)

issue hypparchus.txt .

I don't know if the symmetricity test is just too strict or if there is something else wrong or that I didn't understand but I find it suspicious ;)

Thank you in advance for your clues, All the best,

Quentin

[maisonobe]

Indeed, there is clearly a problem here. The eigen decomposition implementation has two different paths, one for symmetric matrices (using a transformation to tri-diagonal) and one for general matrices (using transformation to Shur form). The problem seems to lie in the transformation to Shur form as the matrix is not symmetric enough to use the tri-diagonal path. As in your case you know the matrix must be symmetric by construction, you can enforce it in a preprocessing step so you are sure you don't go the Shur route. This can be done inserting this code snippet:

for (int i = 0; i < m.getRowDimension(); ++i) {
    for (int j = 0; j < i; ++j) {
        final double mean = 0.5 * (m.getEntry(i, j) + m.getEntry(j, i));
        m.setEntry(i, j, mean);
        m.setEntry(j, i, mean);
    }
}

This is only a workaround, we still have to fix the Shur transformation.

[maisonobe]

The issue seems to be fixed in the master branch now. The Schur transform by itself was correct, but how we extracted eigenvalues from the quasi-triangular matrix did not work correctly on matrices with large coefficients due to numerical rounding. I still think that in your case, as you know the matrix must be symmetric, the preprocessing step to enforce this property should still be done. It should improve numerical stability because you already know non-symmetry is due to previous numerical rounding, and this step at least reduces it so the matrix is more consistent.