From commits-return-10081-apmail-commons-commits-archive=commons.apache.org@commons.apache.org Wed Dec 30 23:17:33 2009 Return-Path: Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: (qmail 91696 invoked from network); 30 Dec 2009 23:17:33 -0000 Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.3) by minotaur.apache.org with SMTP; 30 Dec 2009 23:17:33 -0000 Received: (qmail 21262 invoked by uid 500); 30 Dec 2009 23:17:32 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 21179 invoked by uid 500); 30 Dec 2009 23:17:32 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 21170 invoked by uid 99); 30 Dec 2009 23:17:31 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Wed, 30 Dec 2009 23:17:31 +0000 X-ASF-Spam-Status: No, hits=-2.6 required=5.0 tests=BAYES_00 X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Wed, 30 Dec 2009 23:17:22 +0000 Received: by eris.apache.org (Postfix, from userid 65534) id 235872388A1A; Wed, 30 Dec 2009 23:17:02 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r894735 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/ test/java/org/apache/commons/math/linear/ Date: Wed, 30 Dec 2009 23:17:01 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20091230231702.235872388A1A@eris.apache.org> Author: luc Date: Wed Dec 30 23:17:01 2009 New Revision: 894735 URL: http://svn.apache.org/viewvc?rev=894735&view=rev Log: changed SVD to compute either compact SVD (using only positive singular values) or truncated SVD (using only singular values up to a user-specified max number) started fix of SVD solver that did not compute a least square solution the fix is not complete yet as it seems the solution does not really gives the smallest possible residuals. See for example the commented out parts of testMath320A in SingularValueSolverTest. JIRA: MATH-320 Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecomposition.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecomposition.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecomposition.java?rev=894735&r1=894734&r2=894735&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecomposition.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecomposition.java Wed Dec 30 23:17:01 2009 @@ -24,9 +24,17 @@ * Singular Value Decomposition of a real matrix. *

The Singular Value Decomposition of matrix A is a set of three matrices: * U, Σ and V such that A = U × Σ × VT. - * Let A be an m × n matrix, then U is an m × n orthogonal matrix, - * Σ is a n × n diagonal matrix with positive diagonal elements, - * and V is an n × n orthogonal matrix.

+ * Let A be a m × n matrix, then U is a m × p orthogonal matrix, + * Σ is a p × p diagonal matrix with positive diagonal elements, + * V is a n × p orthogonal matrix (hence VT is a p × n + * orthogonal matrix). The size p depends on the chosen algorithm: + *
    + *
  • for full SVD, p is n,
    • + *
  • for compact SVD, p is the rank r of the matrix + * (i. e. the number of positive singular values),
    • + *
  • for truncated SVD p is min(r, t) where t is user-specified.
    • + *
+ *

*

This interface is similar to the class with similar name from the * JAMA library, with the * following changes:

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java?rev=894735&r1=894734&r2=894735&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java Wed Dec 30 23:17:01 2009 @@ -21,12 +21,24 @@ import org.apache.commons.math.util.MathUtils; /** - * Calculates the Singular Value Decomposition of a matrix. + * Calculates the compact or truncated Singular Value Decomposition of a matrix. *

The Singular Value Decomposition of matrix A is a set of three matrices: * U, Σ and V such that A = U × Σ × VT. - * Let A be an m × n matrix, then U is an m × n orthogonal matrix, - * Σ is a n × n diagonal matrix with positive diagonal elements, - * and V is an n × n orthogonal matrix.

+ * Let A be a m × n matrix, then U is a m × p orthogonal matrix, + * Σ is a p × p diagonal matrix with positive diagonal elements, + * V is a n × p orthogonal matrix (hence VT is a p × n + * orthogonal matrix). The size p depends on the chosen algorithm: + *
    + *
  • for full SVD, p would be n, but this is not supported by this implementation,
    • + *
  • for compact SVD, p is the rank r of the matrix + * (i. e. the number of positive singular values),
    • + *
  • for truncated SVD p is min(r, t) where t is user-specified.
    • + *
+ *

+ *

+ * Note that since this class computes only the compact or truncated SVD and not + * the full SVD, the singular values computed are always positive. + *

* * @version $Revision$ $Date$ * @since 2.0 @@ -76,12 +88,24 @@ private RealMatrix cachedVt; /** + * Calculates the compact Singular Value Decomposition of the given matrix. + * @param matrix The matrix to decompose. + * @exception InvalidMatrixException (wrapping a {@link + * org.apache.commons.math.ConvergenceException} if algorithm fails to converge + */ + public SingularValueDecompositionImpl(final RealMatrix matrix) + throws InvalidMatrixException { + this(matrix, Math.min(matrix.getRowDimension(), matrix.getColumnDimension())); + } + + /** * Calculates the Singular Value Decomposition of the given matrix. * @param matrix The matrix to decompose. + * @param max maximal number of singular values to compute * @exception InvalidMatrixException (wrapping a {@link * org.apache.commons.math.ConvergenceException} if algorithm fails to converge */ - public SingularValueDecompositionImpl(RealMatrix matrix) + public SingularValueDecompositionImpl(final RealMatrix matrix, final int max) throws InvalidMatrixException { m = matrix.getRowDimension(); @@ -113,10 +137,14 @@ eigenDecomposition = new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal, MathUtils.SAFE_MIN); - singularValues = eigenDecomposition.getRealEigenvalues(); - for (int i = 0; i < singularValues.length; ++i) { - final double si = singularValues[i]; - singularValues[i] = (si < 0) ? 0.0 : Math.sqrt(si); + final double[] eigenValues = eigenDecomposition.getRealEigenvalues(); + int p = Math.min(max, eigenValues.length); + while ((p > 0) && (eigenValues[p - 1] <= 0)) { + --p; + } + singularValues = new double[p]; + for (int i = 0; i < p; ++i) { + singularValues[i] = Math.sqrt(eigenValues[i]); } } @@ -127,37 +155,41 @@ if (cachedU == null) { + final int p = singularValues.length; if (m >= n) { // the tridiagonal matrix is Bt.B, where B is upper bidiagonal - final double[][] eData = eigenDecomposition.getV().getData(); - final double[][] iData = new double[m][]; + final RealMatrix e = + eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1); + final double[][] eData = e.getData(); + final double[][] wData = new double[m][p]; double[] ei1 = eData[0]; - iData[0] = ei1; - for (int i = 0; i < n - 1; ++i) { - // compute B.E.S^(-1) where E is the eigenvectors matrix - // we reuse the array from matrix E to store the result + for (int i = 0; i < p - 1; ++i) { + // compute W = B.E.S^(-1) where E is the eigenvectors matrix final double mi = mainBidiagonal[i]; final double si = secondaryBidiagonal[i]; final double[] ei0 = ei1; + final double[] wi = wData[i]; ei1 = eData[i + 1]; - iData[i + 1] = ei1; - for (int j = 0; j < n; ++j) { - ei0[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; + for (int j = 0; j < p; ++j) { + wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; } } // last row - final double lastMain = mainBidiagonal[n - 1]; - for (int j = 0; j < n; ++j) { - ei1[j] *= lastMain / singularValues[j]; + final double lastMain = mainBidiagonal[p - 1]; + final double[] wr1 = wData[p - 1]; + for (int j = 0; j < p; ++j) { + wr1[j] = ei1[j] * lastMain / singularValues[j]; } - for (int i = n; i < m; ++i) { - iData[i] = new double[n]; + for (int i = p; i < m; ++i) { + wData[i] = new double[p]; } cachedU = - transformer.getU().multiply(MatrixUtils.createRealMatrix(iData)); + transformer.getU().multiply(MatrixUtils.createRealMatrix(wData)); } else { // the tridiagonal matrix is B.Bt, where B is lower bidiagonal - cachedU = transformer.getU().multiply(eigenDecomposition.getV()); + final RealMatrix e = + eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); + cachedU = transformer.getU().multiply(e); } } @@ -205,37 +237,41 @@ if (cachedV == null) { + final int p = singularValues.length; if (m >= n) { // the tridiagonal matrix is Bt.B, where B is upper bidiagonal - cachedV = transformer.getV().multiply(eigenDecomposition.getV()); + final RealMatrix e = + eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1); + cachedV = transformer.getV().multiply(e); } else { // the tridiagonal matrix is B.Bt, where B is lower bidiagonal - final double[][] eData = eigenDecomposition.getV().getData(); - final double[][] iData = new double[n][]; + // compute W = Bt.E.S^(-1) where E is the eigenvectors matrix + final RealMatrix e = + eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1); + final double[][] eData = e.getData(); + final double[][] wData = new double[n][p]; double[] ei1 = eData[0]; - iData[0] = ei1; - for (int i = 0; i < m - 1; ++i) { - // compute Bt.E.S^(-1) where E is the eigenvectors matrix - // we reuse the array from matrix E to store the result + for (int i = 0; i < p - 1; ++i) { final double mi = mainBidiagonal[i]; final double si = secondaryBidiagonal[i]; final double[] ei0 = ei1; + final double[] wi = wData[i]; ei1 = eData[i + 1]; - iData[i + 1] = ei1; - for (int j = 0; j < m; ++j) { - ei0[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; + for (int j = 0; j < p; ++j) { + wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; } } // last row - final double lastMain = mainBidiagonal[m - 1]; - for (int j = 0; j < m; ++j) { - ei1[j] *= lastMain / singularValues[j]; + final double lastMain = mainBidiagonal[p - 1]; + final double[] wr1 = wData[p - 1]; + for (int j = 0; j < p; ++j) { + wr1[j] = ei1[j] * lastMain / singularValues[j]; } - for (int i = m; i < n; ++i) { - iData[i] = new double[m]; + for (int i = p; i < n; ++i) { + wData[i] = new double[p]; } cachedV = - transformer.getV().multiply(MatrixUtils.createRealMatrix(iData)); + transformer.getV().multiply(MatrixUtils.createRealMatrix(wData)); } } @@ -262,8 +298,9 @@ public RealMatrix getCovariance(final double minSingularValue) { // get the number of singular values to consider + final int p = singularValues.length; int dimension = 0; - while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) { + while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) { ++dimension; } @@ -273,14 +310,14 @@ minSingularValue, singularValues[0]); } - final double[][] data = new double[dimension][n]; + final double[][] data = new double[dimension][p]; getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { /** {@inheritDoc} */ @Override public void visit(final int row, final int column, final double value) { data[row][column] = value / singularValues[row]; } - }, 0, dimension - 1, 0, n - 1); + }, 0, dimension - 1, 0, p - 1); RealMatrix jv = new Array2DRowRealMatrix(data, false); return jv.transpose().multiply(jv); @@ -317,20 +354,14 @@ /** {@inheritDoc} */ public DecompositionSolver getSolver() { return new Solver(singularValues, getUT(), getV(), - getRank() == singularValues.length); + getRank() == Math.max(m, n)); } /** Specialized solver. */ private static class Solver implements DecompositionSolver { - /** Singular values. */ - private final double[] singularValues; - - /** UT matrix of the decomposition. */ - private final RealMatrix uT; - - /** V matrix of the decomposition. */ - private final RealMatrix v; + /** Pseudo-inverse of the initial matrix. */ + private final RealMatrix pseudoInverse; /** Singularity indicator. */ private boolean nonSingular; @@ -344,10 +375,16 @@ */ private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v, final boolean nonSingular) { - this.singularValues = singularValues; - this.uT = uT; - this.v = v; - this.nonSingular = nonSingular; + double[][] suT = uT.getData(); + for (int i = 0; i < singularValues.length; ++i) { + final double a = 1.0 / singularValues[i]; + final double[] suTi = suT[i]; + for (int j = 0; j < suTi.length; ++j) { + suTi[j] *= a; + } + } + pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); + this.nonSingular = nonSingular; } /** Solve the linear equation A × X = B in least square sense. @@ -356,27 +393,10 @@ * @param b right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular */ public double[] solve(final double[] b) - throws IllegalArgumentException, InvalidMatrixException { - - if (b.length != uT.getColumnDimension()) { - throw MathRuntimeException.createIllegalArgumentException( - "vector length mismatch: got {0} but expected {1}", - b.length, uT.getColumnDimension()); - } - - final double[] w = uT.operate(b); - for (int i = 0; i < singularValues.length; ++i) { - final double si = singularValues[i]; - if (si == 0) { - throw new SingularMatrixException(); - } - w[i] /= si; - } - return v.operate(w); - + throws IllegalArgumentException { + return pseudoInverse.operate(b); } /** Solve the linear equation A × X = B in least square sense. @@ -385,27 +405,10 @@ * @param b right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular */ public RealVector solve(final RealVector b) - throws IllegalArgumentException, InvalidMatrixException { - - if (b.getDimension() != uT.getColumnDimension()) { - throw MathRuntimeException.createIllegalArgumentException( - "vector length mismatch: got {0} but expected {1}", - b.getDimension(), uT.getColumnDimension()); - } - - final RealVector w = uT.operate(b); - for (int i = 0; i < singularValues.length; ++i) { - final double si = singularValues[i]; - if (si == 0) { - throw new SingularMatrixException(); - } - w.setEntry(i, w.getEntry(i) / si); - } - return v.operate(w); - + throws IllegalArgumentException { + return pseudoInverse.operate(b); } /** Solve the linear equation A × X = B in least square sense. @@ -414,31 +417,10 @@ * @param b right-hand side of the equation A × X = B * @return a matrix X that minimizes the two norm of A × X - B * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular */ public RealMatrix solve(final RealMatrix b) - throws IllegalArgumentException, InvalidMatrixException { - - if (b.getRowDimension() != singularValues.length) { - throw MathRuntimeException.createIllegalArgumentException( - "dimensions mismatch: got {0}x{1} but expected {2}x{3}", - b.getRowDimension(), b.getColumnDimension(), - singularValues.length, "n"); - } - - final RealMatrix w = uT.multiply(b); - for (int i = 0; i < singularValues.length; ++i) { - final double si = singularValues[i]; - if (si == 0) { - throw new SingularMatrixException(); - } - final double inv = 1.0 / si; - for (int j = 0; j < b.getColumnDimension(); ++j) { - w.multiplyEntry(i, j, inv); - } - } - return v.multiply(w); - + throws IllegalArgumentException { + return pseudoInverse.multiply(b); } /** @@ -451,17 +433,9 @@ /** Get the pseudo-inverse of the decomposed matrix. * @return inverse matrix - * @throws InvalidMatrixException if decomposed matrix is singular */ - public RealMatrix getInverse() - throws InvalidMatrixException { - - if (!isNonSingular()) { - throw new SingularMatrixException(); - } - - return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length)); - + public RealMatrix getInverse() { + return pseudoInverse; } } Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java?rev=894735&r1=894734&r2=894735&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java Wed Dec 30 23:17:01 2009 @@ -17,18 +17,10 @@ package org.apache.commons.math.linear; -import junit.framework.Test; -import junit.framework.TestCase; -import junit.framework.TestSuite; - -import org.apache.commons.math.linear.DecompositionSolver; -import org.apache.commons.math.linear.InvalidMatrixException; -import org.apache.commons.math.linear.MatrixUtils; -import org.apache.commons.math.linear.RealMatrix; -import org.apache.commons.math.linear.ArrayRealVector; -import org.apache.commons.math.linear.SingularValueDecompositionImpl; +import org.junit.Assert; +import org.junit.Test; -public class SingularValueSolverTest extends TestCase { +public class SingularValueSolverTest { private double[][] testSquare = { { 24.0 / 25.0, 43.0 / 25.0 }, @@ -37,91 +29,68 @@ private static final double normTolerance = 10e-14; - public SingularValueSolverTest(String name) { - super(name); - } - - public static Test suite() { - TestSuite suite = new TestSuite(SingularValueSolverTest.class); - suite.setName("SingularValueSolver Tests"); - return suite; - } - /** test solve dimension errors */ + @Test public void testSolveDimensionErrors() { DecompositionSolver solver = new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver(); RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]); try { solver.solve(b); - fail("an exception should have been thrown"); + Assert.fail("an exception should have been thrown"); } catch (IllegalArgumentException iae) { // expected behavior } catch (Exception e) { - fail("wrong exception caught"); + Assert.fail("wrong exception caught"); } try { solver.solve(b.getColumn(0)); - fail("an exception should have been thrown"); + Assert.fail("an exception should have been thrown"); } catch (IllegalArgumentException iae) { // expected behavior } catch (Exception e) { - fail("wrong exception caught"); + Assert.fail("wrong exception caught"); } try { solver.solve(new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(0))); - fail("an exception should have been thrown"); + Assert.fail("an exception should have been thrown"); } catch (IllegalArgumentException iae) { // expected behavior } catch (Exception e) { - fail("wrong exception caught"); + Assert.fail("wrong exception caught"); } } - /** test solve singularity errors */ - public void testSolveSingularityErrors() { + /** test least square solve */ + @Test + public void testLeastSquareSolve() { RealMatrix m = MatrixUtils.createRealMatrix(new double[][] { { 1.0, 0.0 }, { 0.0, 0.0 } }); DecompositionSolver solver = new SingularValueDecompositionImpl(m).getSolver(); - RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]); - try { - solver.solve(b); - fail("an exception should have been thrown"); - } catch (InvalidMatrixException ime) { - // expected behavior - } catch (Exception e) { - fail("wrong exception caught"); - } - try { - solver.solve(b.getColumn(0)); - fail("an exception should have been thrown"); - } catch (InvalidMatrixException ime) { - // expected behavior - } catch (Exception e) { - fail("wrong exception caught"); - } - try { - solver.solve(b.getColumnVector(0)); - fail("an exception should have been thrown"); - } catch (InvalidMatrixException ime) { - // expected behavior - } catch (Exception e) { - fail("wrong exception caught"); - } - try { - solver.solve(new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(0))); - fail("an exception should have been thrown"); - } catch (InvalidMatrixException ime) { - // expected behavior - } catch (Exception e) { - fail("wrong exception caught"); - } + RealMatrix b = MatrixUtils.createRealMatrix(new double[][] { + { 11, 12 }, { 21, 22 } + }); + RealMatrix xMatrix = solver.solve(b); + Assert.assertEquals(11, xMatrix.getEntry(0, 0), 1.0e-15); + Assert.assertEquals(12, xMatrix.getEntry(0, 1), 1.0e-15); + Assert.assertEquals(0, xMatrix.getEntry(1, 0), 1.0e-15); + Assert.assertEquals(0, xMatrix.getEntry(1, 1), 1.0e-15); + double[] xCol = solver.solve(b.getColumn(0)); + Assert.assertEquals(11, xCol[0], 1.0e-15); + Assert.assertEquals(0, xCol[1], 1.0e-15); + RealVector xColVec = solver.solve(b.getColumnVector(0)); + Assert.assertEquals(11, xColVec.getEntry(0), 1.0e-15); + Assert.assertEquals(0, xColVec.getEntry(1), 1.0e-15); + RealVector xColOtherVec = solver.solve(new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(0))); + Assert.assertEquals(11, xColOtherVec.getEntry(0), 1.0e-15); + Assert.assertEquals(0, xColOtherVec.getEntry(1), 1.0e-15); } /** test solve */ + @Test public void testSolve() { DecompositionSolver solver = new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver(); @@ -134,18 +103,18 @@ }); // using RealMatrix - assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance); + Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance); // using double[] for (int i = 0; i < b.getColumnDimension(); ++i) { - assertEquals(0, + Assert.assertEquals(0, new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(), 1.0e-13); } // using Array2DRowRealMatrix for (int i = 0; i < b.getColumnDimension(); ++i) { - assertEquals(0, + Assert.assertEquals(0, solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(), 1.0e-13); } @@ -154,7 +123,7 @@ for (int i = 0; i < b.getColumnDimension(); ++i) { ArrayRealVectorTest.RealVectorTestImpl v = new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i)); - assertEquals(0, + Assert.assertEquals(0, solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(), 1.0e-13); } @@ -162,10 +131,82 @@ } /** test condition number */ + @Test public void testConditionNumber() { SingularValueDecompositionImpl svd = new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)); - assertEquals(3.0, svd.getConditionNumber(), 1.0e-15); + Assert.assertEquals(3.0, svd.getConditionNumber(), 1.0e-15); + } + + @Test + public void testMath320A() { + RealMatrix rm = new Array2DRowRealMatrix(new double[][] { + { 1.0, 2.0, 3.0 }, { 2.0, 3.0, 4.0 }, { 3.0, 5.0, 7.0 } + }); + double s439 = Math.sqrt(439.0); + double[] reference = new double[] { + Math.sqrt(3.0 * (21.0 + s439)), Math.sqrt(3.0 * (21.0 - s439)) + }; + SingularValueDecomposition svd = + new SingularValueDecompositionImpl(rm); + double[] singularValues = svd.getSingularValues(); + for (int i = 0; i < reference.length; ++i) { + Assert.assertEquals(reference[i], singularValues[i], 4.0e-13); + } + regularElements(svd.getU()); + regularElements(svd.getVT()); +// double[] b = new double[] { 5.0, 6.0, 7.0 }; +// double[] resSVD = svd.getSolver().solve(b); +// Assert.assertEquals(rm.getColumnDimension(), resSVD.length); +// System.out.println("resSVD = " + resSVD[0] + " " + resSVD[1] + " " + resSVD[2]); +// double minResidual = Double.POSITIVE_INFINITY; +// double d0Min = Double.NaN; +// double d1Min = Double.NaN; +// double d2Min = Double.NaN; +// double h = 0.01; +// int k = 100; +// for (double d0 = -k * h; d0 <= k * h; d0 += h) { +// for (double d1 = -k * h ; d1 <= k * h; d1 += h) { +// for (double d2 = -k * h; d2 <= k * h; d2 += h) { +// double[] f = rm.operate(new double[] { resSVD[0] + d0, resSVD[1] + d1, resSVD[2] + d2 }); +// double residual = Math.sqrt((f[0] - b[0]) * (f[0] - b[0]) + +// (f[1] - b[1]) * (f[1] - b[1]) + +// (f[2] - b[2]) * (f[2] - b[2])); +// if (residual < minResidual) { +// d0Min = d0; +// d1Min = d1; +// d2Min = d2; +// minResidual = residual; +// } +// } +// } +// } +// System.out.println(d0Min + " " + d1Min + " " + d2Min + " -> " + minResidual); +// Assert.assertEquals(0, d0Min, 1.0e-15); +// Assert.assertEquals(0, d1Min, 1.0e-15); +// Assert.assertEquals(0, d2Min, 1.0e-15); + } + + + @Test + public void testMath320B() { + RealMatrix rm = new Array2DRowRealMatrix(new double[][] { + { 1.0, 2.0 }, { 1.0, 2.0 } + }); + SingularValueDecomposition svd = + new SingularValueDecompositionImpl(rm); + regularElements(svd.getU()); + regularElements(svd.getVT()); + } + + private void regularElements(RealMatrix m) { + for (int i = 0; i < m.getRowDimension(); ++i) { + for (int j = 0; j < m.getColumnDimension(); ++j) { + double mij = m.getEntry(i, j); + Assert.assertFalse(Double.isInfinite(mij)); + Assert.assertFalse(Double.isNaN(mij)); + } + } } }