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)
