[jira] [Created] (MATH-1424) Wrong eigen values computed by EigenDecomposition when the input matrix has large values

["Jerome (JIRA)" <ji...@apache.org>]

Jerome created MATH-1424:
----------------------------

             Summary: Wrong eigen values computed by EigenDecomposition when the input matrix has large values
                 Key: MATH-1424
                 URL: https://issues.apache.org/jira/browse/MATH-1424
             Project: Commons Math
          Issue Type: Bug
    Affects Versions: 3.6.1
         Environment: JDK 7.51 64 bits on Windows 7.
            Reporter: Jerome


The following code gives a wrong result:
RealMatrix m = [[10_000_000.0, -1_000_000.0],[-1_000_000.1, 20_000_000.0]]; // pseudo code
EigenDecomposition ed = new EigenDecomposition(m);
double[] eigenValues = ed.getRealEigenvalues();

Computed values: [1.57E13, 1.57E13].
Expected values: [1.0E7, 2.0E7]

The problem lies in method EigenDecomposition.transformToSchur(RealMatrix).
At line 758, the value matT[i+1][i] is checked against 0.0 within an EPSILON margin.
If the precision of the computation were perfect, matT[i+1][i] == 0.0 means that matT[i][i] is a solution of the characteristic polynomial of m. In the other case there are 2 complex solutions.
But due to imprecisions, this value can be different from 0.0 while m has only real solutions.
The else part assume that the solutions are complex, which is wrong in the provided example.
To correct it, you should resolve the 2 degree polynomial without assuming the solutions are complex (that is: test whether p*p + matT[i+1][i] * matT[i][i+1] is negative for 2 complex solutions, or positive or null for 2 real solutions).
You should also avoid testing values against something within epsilon margin, because this method is almost always wrong in some cases. At least, check with a margin that depends on the amplitude of the value (ex: margin = highest absolute value of the matrix * EPSILON); this is still wrong but problems will occur less often.

The problem occurs when the input matrix has large values because matT has values of magnitude E7. MatT[1, 0] is really low (E-10) and you can not expect a better precision due to the large values on the diagonal.
The test within EPSILON margin fails, which does not occurs when the input matrix has lowest values.
Testing the code with m2 = m / pow(2, 20) will work, because matT[1, 0] is now low enough.



--
This message was sent by Atlassian JIRA
(v6.4.14#64029)