+ * 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. + *
+ * 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 diff --git a/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java b/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java index f2db4fd61..15e561958 100644 --- a/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java +++ b/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java @@ -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 } } diff --git a/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java b/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java index f3ff9c673..a4f0fbd00 100644 --- a/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java +++ b/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java @@ -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);