You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2012/07/18 22:49:44 UTC
svn commit: r1363105 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/linear/
test/java/org/apache/commons/math3/linear/
Author: tn
Date: Wed Jul 18 20:49:43 2012
New Revision: 1363105
URL: http://svn.apache.org/viewvc?rev=1363105&view=rev
Log:
[MATH-235] Added random data test for eigen decomposition, improved error handling.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java?rev=1363105&r1=1363104&r2=1363105&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java Wed Jul 18 20:49:43 2012
@@ -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 × X = B solution in exact
* linear sense.
- *
- * @return a solver.
+ * <p>
+ * Since 3.1, eigen decomposition of a general matrix is supported,
+ * but the {@link DecompositionSolver} only supports real eigenvalues.
+ *
+ * @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
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java?rev=1363105&r1=1363104&r2=1363105&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java Wed Jul 18 20:49:43 2012
@@ -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
}
}
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java?rev=1363105&r1=1363104&r2=1363105&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java Wed Jul 18 20:49:43 2012
@@ -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);