From 06c63da6c713adc0b17e7f7e0daa71b1079fd84a Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Tue, 30 Oct 2012 00:29:06 +0000 Subject: [PATCH] MATH-883 Added "getSquareRoot()" method. Implementation only supports symmetric, diagonalizable matrices. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1403590 13f79535-47bb-0310-9956-ffa450edef68 --- src/changes/changes.xml | 4 ++ .../math3/linear/EigenDecomposition.java | 35 +++++++++++++- .../math3/linear/EigenDecompositionTest.java | 48 +++++++++++++++++++ 3 files changed, 86 insertions(+), 1 deletion(-) diff --git a/src/changes/changes.xml b/src/changes/changes.xml index c5e603da8..e9708cd8b 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -52,6 +52,10 @@ If the output is not quite correct, check for invisible trailing spaces! + + New "getSquareRoot" method in class "EigenDecomposition" (package + "o.a.c.m.linear"). + Added "isSymmetric" and "checkSymmetric" in "MatrixUtils" (package "o.a.c.m.linear"). diff --git a/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java b/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java index cdf815083..7acde9193 100644 --- a/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java +++ b/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java @@ -100,6 +100,8 @@ public class EigenDecomposition { private RealMatrix cachedD; /** Cached value of Vt. */ private RealMatrix cachedVt; + /** Whether the matrix is symmetric. */ + private final boolean isSymmetric; /** * Calculates the eigen decomposition of the given real matrix. @@ -113,7 +115,8 @@ public class EigenDecomposition { */ public EigenDecomposition(final RealMatrix matrix) throws MathArithmeticException { - if (isSymmetric(matrix, false)) { + isSymmetric = isSymmetric(matrix, false); + if (isSymmetric) { transformToTridiagonal(matrix); findEigenVectors(transformer.getQ().getData()); } else { @@ -149,6 +152,7 @@ public class EigenDecomposition { * @throws MaxCountExceededException if the algorithm fails to converge. */ public EigenDecomposition(final double[] main, final double[] secondary) { + isSymmetric = true; this.main = main.clone(); this.secondary = secondary.clone(); transformer = null; @@ -385,6 +389,35 @@ public class EigenDecomposition { return determinant; } + /** + * Computes the square-root of the matrix. + * This implementation assumes that the matrix is symmetric and postive + * definite. + * + * @return the square-root of the matrix. + * @throws MathUnsupportedOperationException if the matrix is not + * symmetric or not positive definite. + */ + public RealMatrix getSquareRoot() { + if (!isSymmetric) { + throw new MathUnsupportedOperationException(); + } + + final double[] sqrtEigenValues = new double[realEigenvalues.length]; + for (int i = 0; i < realEigenvalues.length; i++) { + final double eigen = realEigenvalues[i]; + if (eigen <= 0) { + throw new MathUnsupportedOperationException(); + } + sqrtEigenValues[i] = FastMath.sqrt(eigen); + } + final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues); + final RealMatrix v = getV(); + final RealMatrix vT = getVT(); + + return v.multiply(sqrtEigen).multiply(vT); + } + /** * Gets a solver for finding the A × X = B solution in exact * linear sense. 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 48a8ec46d..04805acdf 100644 --- a/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java +++ b/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java @@ -24,6 +24,7 @@ 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.apache.commons.math3.exception.MathUnsupportedOperationException; import org.junit.After; import org.junit.Assert; import org.junit.Before; @@ -336,6 +337,53 @@ public class EigenDecompositionTest { Assert.assertEquals(0, norm, 6.0e-13); } + @Test + public void testSquareRoot() { + final double[][] data = { + { 33, 24, 7 }, + { 24, 57, 11 }, + { 7, 11, 9 } + }; + + final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); + final RealMatrix sqrtM = dec.getSquareRoot(); + + // Reconstruct initial matrix. + final RealMatrix m = sqrtM.multiply(sqrtM); + + final int dim = data.length; + for (int r = 0; r < dim; r++) { + for (int c = 0; c < dim; c++) { + Assert.assertEquals("m[" + r + "][" + c + "]", + data[r][c], m.getEntry(r, c), 1e-13); + } + } + } + + @Test(expected=MathUnsupportedOperationException.class) + public void testSquareRootNonSymmetric() { + final double[][] data = { + { 1, 2, 4 }, + { 2, 3, 5 }, + { 11, 5, 9 } + }; + + final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); + final RealMatrix sqrtM = dec.getSquareRoot(); + } + + @Test(expected=MathUnsupportedOperationException.class) + public void testSquareRootNonPositiveDefinite() { + final double[][] data = { + { 1, 2, 4 }, + { 2, 3, 5 }, + { 4, 5, -9 } + }; + + final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); + final RealMatrix sqrtM = dec.getSquareRoot(); + } + @Test public void testUnsymmetric() { // Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4)