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

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

    [ https://issues.apache.org/jira/browse/MATH-1424?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16076328#comment-16076328 ] 

Gilles commented on MATH-1424:
------------------------------

Thanks for the report.
Could you prepare a patch, with unit test, against the current development version ("master" branch of the repository)?


> 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
>              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 (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)