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

[MATH-235] Added random data test for eigen decomposition, improved error handling. 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1363105 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Thomas Neidhart 2012-07-18 20:49:43 +00:00
parent 7dd77c4bdf
commit 1213c8a2dd
3 changed files with 98 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
View File

@ -18,6 +18,8 @@
package org.apache.commons.math3.linear;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
@ -102,9 +104,13 @@ public class EigenDecomposition {
/**
* Calculates the eigen decomposition of the given real matrix.
* <p>
* Supports decomposition of a general matrix since 3.1.
*
* @param matrix Matrix to decompose.
* @throws MaxCountExceededException if the algorithm fails to converge.
* @throws MathArithmeticException if the decomposition of a general matrix
* results in a matrix with zero norm
*/
public EigenDecomposition(final RealMatrix matrix) {
if (isSymmetric(matrix, false)) {
@ -218,7 +224,6 @@ public class EigenDecomposition {
}
// return the cached matrix
return cachedV;
}
/**
@ -233,6 +238,7 @@ public class EigenDecomposition {
* @see #getImagEigenvalues()
*/
public RealMatrix getD() {
if (cachedD == null) {
// cache the matrix for subsequent calls
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
@ -359,10 +365,20 @@ public class EigenDecomposition {
/**
* Gets a solver for finding the A &times; X = B solution in exact
* linear sense.
* <p>
* Since 3.1, eigen decomposition of a general matrix is supported,
* but the {@link DecompositionSolver} only supports real eigenvalues.
*
* @return a solver.
* @return a solver
* @throws MathUnsupportedOperationException if the decomposition resulted in
* complex eigenvalues
*/
public DecompositionSolver getSolver() {
for (int i = 0; i < imagEigenvalues.length; i++) {
if (imagEigenvalues[i] != 0.0) {
throw new MathUnsupportedOperationException();
}
}
return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
}
@ -726,6 +742,7 @@ public class EigenDecomposition {
* Find eigenvectors from a matrix transformed to Schur form.
*
* @param schur the schur transformation of the matrix
* @throws MathArithmeticException if the Schur form has a norm of zero
*/
private void findEigenVectorsFromSchur(final SchurTransformer schur) {
final double[][] matrixT = schur.getT().getData();
@ -741,9 +758,9 @@ public class EigenDecomposition {
}
}
if (Precision.equals(norm, 0.0)) {
// TODO: we can not handle a zero matrix, what exception to throw?
return;
// we can not handle a matrix with zero norm
if (Precision.equals(norm, 0.0, epsilon)) {
throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
}
// Backsubstitute to find vectors of upper triangular form

10
src/main/java/org/apache/commons/math3/linear/SchurTransformer.java
View File

@ -436,13 +436,17 @@ class SchurTransformer {
* Contains variable names as present in the original JAMA code.
*/
private static class ShiftInfo {
/** TODO: describe */
// CHECKSTYLE: stop all
/** x shift info */
double x;
/** TODO: describe */
/** y shift info */
double y;
/** TODO: describe */
/** w shift info */
double w;
/** Indicates an exceptional shift. */
double exShift;
// CHECKSTYLE: resume all
}
}

103
src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java
View File

@ -21,6 +21,7 @@ import java.util.Arrays;
import java.util.Random;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.After;
@ -38,7 +39,7 @@ public class EigenDecompositionTest {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
}
@ -50,7 +51,7 @@ public class EigenDecompositionTest {
{ 12.0, 66.0 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
Assert.assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
}
@ -64,7 +65,7 @@ public class EigenDecompositionTest {
{ -16560.0, 7920.0, 17300.0 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
Assert.assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
Assert.assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
@ -79,7 +80,7 @@ public class EigenDecompositionTest {
{ 15, 30, 45 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(70.0, ed.getRealEigenvalue(0), 3.0e-11);
Assert.assertEquals(0.0, ed.getRealEigenvalue(1), 3.0e-11);
Assert.assertEquals(0.0, ed.getRealEigenvalue(2), 3.0e-11);
@ -95,7 +96,7 @@ public class EigenDecompositionTest {
{ 0.000, 0.000, -0.048, 0.136 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
@ -112,7 +113,7 @@ public class EigenDecompositionTest {
{ -0.2976, 0.1152, -0.1344, 0.3872 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
@ -144,9 +145,7 @@ public class EigenDecompositionTest {
};
EigenDecomposition decomposition;
decomposition = new EigenDecomposition(mainTridiagonal,
secondaryTridiagonal,
Precision.SAFE_MIN);
decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
double[] eigenValues = decomposition.getRealEigenvalues();
for (int i = 0; i < refEigenValues.length; ++i) {
@ -189,9 +188,7 @@ public class EigenDecompositionTest {
// the following line triggers the exception
EigenDecomposition decomposition;
decomposition = new EigenDecomposition(mainTridiagonal,
secondaryTridiagonal,
Precision.SAFE_MIN);
decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
double[] eigenValues = decomposition.getRealEigenvalues();
for (int i = 0; i < refEigenValues.length; ++i) {
@ -236,9 +233,7 @@ public class EigenDecompositionTest {
// the following line triggers the exception
EigenDecomposition decomposition;
decomposition = new EigenDecomposition(mainTridiagonal,
secondaryTridiagonal,
Precision.SAFE_MIN);
decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
double[] eigenValues = decomposition.getRealEigenvalues();
for (int i = 0; i < refEigenValues.length; ++i) {
@ -268,9 +263,7 @@ public class EigenDecompositionTest {
TriDiagonalTransformer t =
new TriDiagonalTransformer(createTestMatrix(r, ref));
EigenDecomposition ed;
ed = new EigenDecomposition(t.getMainDiagonalRef(),
t.getSecondaryDiagonalRef(),
Precision.SAFE_MIN);
ed = new EigenDecomposition(t.getMainDiagonalRef(), t.getSecondaryDiagonalRef());
double[] eigenValues = ed.getRealEigenvalues();
Assert.assertEquals(ref.length, eigenValues.length);
for (int i = 0; i < ref.length; ++i) {
@ -284,7 +277,7 @@ public class EigenDecompositionTest {
public void testDimensions() {
final int m = matrix.getRowDimension();
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
Assert.assertEquals(m, ed.getV().getRowDimension());
Assert.assertEquals(m, ed.getV().getColumnDimension());
Assert.assertEquals(m, ed.getD().getColumnDimension());
@ -297,7 +290,7 @@ public class EigenDecompositionTest {
@Test
public void testEigenvalues() {
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
double[] eigenValues = ed.getRealEigenvalues();
Assert.assertEquals(refValues.length, eigenValues.length);
for (int i = 0; i < refValues.length; ++i) {
@ -315,8 +308,7 @@ public class EigenDecompositionTest {
}
Arrays.sort(bigValues);
EigenDecomposition ed;
ed = new EigenDecomposition(createTestMatrix(r, bigValues),
Precision.SAFE_MIN);
ed = new EigenDecomposition(createTestMatrix(r, bigValues));
double[] eigenValues = ed.getRealEigenvalues();
Assert.assertEquals(bigValues.length, eigenValues.length);
for (int i = 0; i < bigValues.length; ++i) {
@ -333,7 +325,7 @@ public class EigenDecompositionTest {
});
EigenDecomposition ed;
ed = new EigenDecomposition(symmetric, Precision.SAFE_MIN);
ed = new EigenDecomposition(symmetric);
RealMatrix d = ed.getD();
RealMatrix v = ed.getV();
@ -341,8 +333,6 @@ public class EigenDecompositionTest {
double norm = v.multiply(d).multiply(vT).subtract(symmetric).getNorm();
Assert.assertEquals(0, norm, 6.0e-13);
// check(A.times(V),V.times(D));
}
@Test
@ -374,8 +364,53 @@ public class EigenDecompositionTest {
checkUnsymmetricMatrix(MatrixUtils.createRealMatrix(randData2));
}
@Test
public void testRandomUnsymmetricMatrix() {
for (int run = 0; run < 100; run++) {
Random r = new Random(System.currentTimeMillis());
// matrix size
int size = r.nextInt(20) + 4;
double[][] data = new double[size][size];
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
data[i][j] = r.nextInt(100);
}
}
RealMatrix m = MatrixUtils.createRealMatrix(data);
checkUnsymmetricMatrix(m);
}
}
@Test
public void testNormalDistributionUnsymmetricMatrix() {
for (int run = 0; run < 100; run++) {
Random r = new Random(System.currentTimeMillis());
NormalDistribution dist = new NormalDistribution(0.0, r.nextDouble() * 5);
// matrix size
int size = r.nextInt(20) + 4;
double[][] data = new double[size][size];
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
data[i][j] = dist.sample();
}
}
RealMatrix m = MatrixUtils.createRealMatrix(data);
checkUnsymmetricMatrix(m);
}
}
/**
* Checks that the eigen decomposition of a general (unsymmetric) matrix is valid by
* checking: A*V = V*D
*/
private void checkUnsymmetricMatrix(final RealMatrix m) {
EigenDecomposition ed = new EigenDecomposition(m, Precision.SAFE_MIN);
EigenDecomposition ed = new EigenDecomposition(m);
RealMatrix d = ed.getD();
RealMatrix v = ed.getV();
@ -389,14 +424,14 @@ public class EigenDecompositionTest {
RealMatrix invV = new LUDecomposition(v).getSolver().getInverse();
double norm = v.multiply(d).multiply(invV).subtract(m).getNorm();
Assert.assertEquals(0.0, norm, 6.0e-13);
Assert.assertEquals(0.0, norm, 1.0e-10);
}
/** test eigenvectors */
@Test
public void testEigenvectors() {
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
for (int i = 0; i < matrix.getRowDimension(); ++i) {
double lambda = ed.getRealEigenvalue(i);
RealVector v = ed.getEigenvector(i);
@ -409,7 +444,7 @@ public class EigenDecompositionTest {
@Test
public void testAEqualVDVt() {
EigenDecomposition ed;
ed = new EigenDecomposition(matrix, Precision.SAFE_MIN);
ed = new EigenDecomposition(matrix);
RealMatrix v = ed.getV();
RealMatrix d = ed.getD();
RealMatrix vT = ed.getVT();
@ -420,7 +455,7 @@ public class EigenDecompositionTest {
/** test that V is orthogonal */
@Test
public void testVOrthogonal() {
RealMatrix v = new EigenDecomposition(matrix, Precision.SAFE_MIN).getV();
RealMatrix v = new EigenDecomposition(matrix).getV();
RealMatrix vTv = v.transpose().multiply(v);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
Assert.assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
@ -432,7 +467,7 @@ public class EigenDecompositionTest {
double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
EigenDecomposition ed;
ed = new EigenDecomposition(m, Precision.SAFE_MIN);
ed = new EigenDecomposition(m);
Assert.assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
Assert.assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
Assert.assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
@ -450,7 +485,7 @@ public class EigenDecompositionTest {
{4, 2, 3}
});
EigenDecomposition ed;
ed = new EigenDecomposition(repeated, Precision.SAFE_MIN);
ed = new EigenDecomposition(repeated);
checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
}
@ -466,7 +501,7 @@ public class EigenDecompositionTest {
{-4, -4, 8}
});
EigenDecomposition ed;
ed = new EigenDecomposition(distinct, Precision.SAFE_MIN);
ed = new EigenDecomposition(distinct);
checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
@ -484,7 +519,7 @@ public class EigenDecompositionTest {
{ -1.0,0.0, 1.0 }
});
EigenDecomposition ed;
ed = new EigenDecomposition(indefinite, Precision.SAFE_MIN);
ed = new EigenDecomposition(indefinite);
checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12);
double isqrt3 = 1/FastMath.sqrt(3.0);
checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12);