diff --git a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java index 2c5da29c8..53a40f08c 100644 --- a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java +++ b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java @@ -1699,19 +1699,21 @@ public class EigenDecompositionImpl implements EigenDecomposition { // perform an initial non-shifted LDLt decomposition final double[] d = new double[m]; final double[] l = new double[m - 1]; - double di = main[0]; + // avoid zero divide on indefinite matrix + final double mu = realEigenvalues[m-1] <= 0 && realEigenvalues[0] > 0 ? 0.5-realEigenvalues[m-1] : 0; + double di = main[0]+mu; d[0] = di; for (int i = 1; i < m; ++i) { final double eiM1 = secondary[i - 1]; final double ratio = eiM1 / di; - di = main[i] - eiM1 * ratio; + di = main[i] - eiM1 * ratio + mu; l[i - 1] = ratio; d[i] = di; } // compute eigenvectors for (int i = 0; i < m; ++i) { - eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l); + eigenvectors[i] = findEigenvector(realEigenvalues[i]+mu, d, l); } } diff --git a/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java b/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java index de595e5db..9d7563921 100644 --- a/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java +++ b/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java @@ -243,6 +243,24 @@ public class EigenDecompositionImplTest extends TestCase { checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12); } + /** + * Verifies operation on indefinite matrix + */ + public void testZeroDivide() { + RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] { + { 0.0, 1.0, -1.0 }, + { 1.0, 1.0, 0.0 }, + { -1.0,0.0, 1.0 } + }); + EigenDecomposition ed = new EigenDecompositionImpl(indefinite, MathUtils.SAFE_MIN); + checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12); + double isqrt3 = 1/Math.sqrt(3.0); + checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12); + double isqrt2 = 1/Math.sqrt(2.0); + checkEigenVector((new double[] {0.0,-isqrt2,-isqrt2}), ed, 1E-12); + double isqrt6 = 1/Math.sqrt(6.0); + checkEigenVector((new double[] {2*isqrt6,-isqrt6,isqrt6}), ed, 1E-12); + } /** * Verifies that the given EigenDecomposition has eigenvalues equivalent to * the targetValues, ignoring the order of the values and allowing @@ -257,6 +275,7 @@ public class EigenDecompositionImplTest extends TestCase { } } + /** * Returns true iff there is an entry within tolerance of value in * searchArray.