[ https://issues.apache.org/jira/browse/MATH1424?page=com.atlassian.jira.plugin.system.issuetabpanels:commenttabpanel&focusedCommentId=16076370#comment16076370
]
Gilles commented on MATH1424:

Fine, since it's unlikely that we'll be in a position to consider a release of "Commons Math"
before we finalize the initial release of ["Commons Numbers"http://commons.apache.org/proper/commonsnumbers/modules.html].
> Wrong eigen values computed by EigenDecomposition when the input matrix has large values
> 
>
> Key: MATH1424
> URL: https://issues.apache.org/jira/browse/MATH1424
> Project: Commons Math
> Issue Type: Bug
> Affects Versions: 3.6.1
> Environment: JDK 7.51 64 bits on Windows 7.
> Reporter: Jerome
> Labels: easyfix
>
> 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 (E10) 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)
