[ https://issues.apache.org/jira/browse/MATH-1424?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16076370#comment-16076370 ]
Gilles commented on MATH-1424: ------------------------------ 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/commons-numbers/modules.html]. > 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)