commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
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 GMT
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.
  * <p>The Singular Value Decomposition of matrix A is a set of three matrices:
  * U, &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>.
- * Let A be an m &times; n matrix, then U is an m &times; n orthogonal matrix,
- * &Sigma; is a n &times; n diagonal matrix with positive diagonal elements,
- * and V is an n &times; n orthogonal matrix.</p>
+ * Let A be a m &times; n matrix, then U is a m &times; p orthogonal matrix,
+ * &Sigma; is a p &times; p diagonal matrix with positive diagonal elements,
+ * V is a n &times; p orthogonal matrix (hence V<sup>T</sup> is a p &times;
n
+ * orthogonal matrix). The size p depends on the chosen algorithm:
+ * <ul>
+ *   <li>for full SVD, p is n,</li>
+ *   <li>for compact SVD, p is the rank r of the matrix
+ *       (i. e. the number of positive singular values),</li>
+ *   <li>for truncated SVD p is min(r, t) where t is user-specified.</li>
+ * </ul>
+ * </p>
  * <p>This interface is similar to the class with similar name from the
  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
  * following changes:</p>

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.
  * <p>The Singular Value Decomposition of matrix A is a set of three matrices:
  * U, &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>.
- * Let A be an m &times; n matrix, then U is an m &times; n orthogonal matrix,
- * &Sigma; is a n &times; n diagonal matrix with positive diagonal elements,
- * and V is an n &times; n orthogonal matrix.</p>
+ * Let A be a m &times; n matrix, then U is a m &times; p orthogonal matrix,
+ * &Sigma; is a p &times; p diagonal matrix with positive diagonal elements,
+ * V is a n &times; p orthogonal matrix (hence V<sup>T</sup> is a p &times;
n
+ * orthogonal matrix). The size p depends on the chosen algorithm:
+ * <ul>
+ *   <li>for full SVD, p would be n, but this is not supported by this implementation,</li>
+ *   <li>for compact SVD, p is the rank r of the matrix
+ *       (i. e. the number of positive singular values),</li>
+ *   <li>for truncated SVD p is min(r, t) where t is user-specified.</li>
+ * </ul>
+ * </p>
+ * <p>
+ * Note that since this class computes only the compact or truncated SVD and not
+ * the full SVD, the singular values computed are always positive.
+ * </p>
  *
  * @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;
-
-        /** U<sup>T</sup> 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 &times; X = B in least square sense.
@@ -356,27 +393,10 @@
          * @param b right-hand side of the equation A &times; X = B
          * @return a vector X that minimizes the two norm of A &times; 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 &times; X = B in least square sense.
@@ -385,27 +405,10 @@
          * @param b right-hand side of the equation A &times; X = B
          * @return a vector X that minimizes the two norm of A &times; 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 &times; X = B in least square sense.
@@ -414,31 +417,10 @@
          * @param b right-hand side of the equation A &times; X = B
          * @return a matrix X that minimizes the two norm of A &times; 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));
+            }
+        }
     }
 
 }



Mime
View raw message