Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fix for possible zero divide on an indefinite matrix 

Fix for: MATH-297


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@826550 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
William Barker 2009-10-18 23:29:37 +00:00
parent ddc5554353
commit becc8d9704
2 changed files with 24 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

@ -1699,19 +1699,21 @@ public class EigenDecompositionImpl implements EigenDecomposition {
// perform an initial non-shifted LDLt decomposition // perform an initial non-shifted LDLt decomposition
final double[] d = new double[m]; final double[] d = new double[m];
final double[] l = new double[m - 1]; 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; d[0] = di;
for (int i = 1; i < m; ++i) { for (int i = 1; i < m; ++i) {
final double eiM1 = secondary[i - 1]; final double eiM1 = secondary[i - 1];
final double ratio = eiM1 / di; final double ratio = eiM1 / di;
di = main[i] - eiM1 * ratio; di = main[i] - eiM1 * ratio + mu;
l[i - 1] = ratio; l[i - 1] = ratio;
d[i] = di; d[i] = di;
} }
// compute eigenvectors // compute eigenvectors
for (int i = 0; i < m; ++i) { for (int i = 0; i < m; ++i) {
eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l); eigenvectors[i] = findEigenvector(realEigenvalues[i]+mu, d, l);
} }
} }

19
src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
View File

@ -243,6 +243,24 @@ public class EigenDecompositionImplTest extends TestCase {
checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12); 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 * Verifies that the given EigenDecomposition has eigenvalues equivalent to
* the targetValues, ignoring the order of the values and allowing * 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 * Returns true iff there is an entry within tolerance of value in
* searchArray. * searchArray.