commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r894908 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java site/xdoc/changes.xml test/java/org/apache/commons/math/linear/SingularValueSolverTest.java
Date Thu, 31 Dec 2009 17:52:16 GMT
Author: luc
Date: Thu Dec 31 17:52:16 2009
New Revision: 894908

URL: http://svn.apache.org/viewvc?rev=894908&view=rev
Log:
Singular Value Decomposition now computes either the compact SVD (using only positive singular
values) or truncated SVD (using a user-specified maximal number of singular values).
Fixed Singular Value Decomposition solving of singular systems.
JIRA: MATH-320, MATH-321

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    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/SingularValueDecompositionImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java?rev=894908&r1=894907&r2=894908&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
Thu Dec 31 17:52:16 2009
@@ -159,27 +159,28 @@
             if (m >= n) {
                 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal
                 final RealMatrix e =
-                    eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
+                    eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1);
                 final double[][] eData = e.getData();
                 final double[][] wData = new double[m][p];
                 double[] ei1 = eData[0];
-                for (int i = 0; i < p - 1; ++i) {
+                for (int i = 0; i < p; ++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];
-                    for (int j = 0; j < p; ++j) {
-                        wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
+                    if (i < n - 1) {
+                        ei1 = eData[i + 1];
+                        final double si = secondaryBidiagonal[i];
+                        for (int j = 0; j < p; ++j) {
+                            wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
+                        }
+                    } else {
+                        for (int j = 0; j < p; ++j) {
+                            wi[j] = mi * ei0[j] / singularValues[j];
+                        }
                     }
                 }
-                // last row
-                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 = p; i < m; ++i) {
                     wData[i] = new double[p];
                 }
@@ -247,26 +248,26 @@
                 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal
                 // compute W = Bt.E.S^(-1) where E is the eigenvectors matrix
                 final RealMatrix e =
-                    eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
+                    eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
                 final double[][] eData = e.getData();
                 final double[][] wData = new double[n][p];
                 double[] ei1 = eData[0];
-                for (int i = 0; i < p - 1; ++i) {
+                for (int i = 0; i < p; ++i) {
                     final double mi = mainBidiagonal[i];
-                    final double si = secondaryBidiagonal[i];
                     final double[] ei0 = ei1;
                     final double[] wi  = wData[i];
-                    ei1 = eData[i + 1];
-                    for (int j = 0; j < p; ++j) {
-                        wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
+                    if (i < m - 1) {
+                        ei1 = eData[i + 1];
+                        final double si = secondaryBidiagonal[i];
+                        for (int j = 0; j < p; ++j) {
+                            wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
+                        }
+                    } else {
+                        for (int j = 0; j < p; ++j) {
+                            wi[j] = mi * ei0[j] / singularValues[j];
+                        }
                     }
                 }
-                // last row
-                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 = p; i < n; ++i) {
                     wData[i] = new double[p];
                 }

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=894908&r1=894907&r2=894908&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Thu Dec 31 17:52:16 2009
@@ -39,10 +39,18 @@
   </properties>
   <body>
     <release version="2.1" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-321" >
+        Singular Value Decomposition now computes either the compact SVD (using only
+        positive singular values) or truncated SVD (using a user-specified maximal
+        number of singular values).
+      </action>
+      <action dev="luc" type="fix" issue="MATH-320" >
+        Fixed Singular Value Decomposition solving of singular systems.
+      </action>
       <action dev="psteitz" type="update" issue="MATH-239" due-to="Christian Semrau">
         Added MathUtils methods to compute gcd and lcm for long arguments.
-      </actions>
-      <action dev="psteitz" tyoe="update" issue="MATH-287" due-to="Matthew Rowles">
+      </action>
+      <action dev="psteitz" type="update" issue="MATH-287" due-to="Matthew Rowles">
         Added support for weighted univariate statistics.
       </action>
       <action dev="luc" type="fix" issue="MATH-326" due-to="Jake Mannix">

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=894908&r1=894907&r2=894908&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
Thu Dec 31 17:52:16 2009
@@ -139,6 +139,32 @@
     }
 
     @Test
+    public void testTruncated() {
+
+        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))
+        };
+        SingularValueDecomposition svd =
+            new SingularValueDecompositionImpl(rm, 1);
+
+        // check we get the expected theoretical singular values
+        double[] singularValues = svd.getSingularValues();
+        Assert.assertEquals(reference.length, singularValues.length);
+        for (int i = 0; i < reference.length; ++i) {
+            Assert.assertEquals(reference[i], singularValues[i], 4.0e-13);
+        }
+
+        // check the truncated decomposition DON'T allows to recover the original matrix
+        RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
+        Assert.assertTrue(recomposed.subtract(rm).getNorm() > 1.4);
+
+    }
+
+    @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 }
@@ -149,44 +175,38 @@
         };
         SingularValueDecomposition svd =
             new SingularValueDecompositionImpl(rm);
+
+        // check we get the expected theoretical singular values
         double[] singularValues = svd.getSingularValues();
+        Assert.assertEquals(reference.length, singularValues.length);
         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);
-    }
 
+        // check the decomposition allows to recover the original matrix
+        RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
+        Assert.assertEquals(0.0, recomposed.subtract(rm).getNorm(), 5.0e-13);
+
+        // check we can solve a singular system
+        double[] b = new double[] { 5.0, 6.0, 7.0 };
+        double[] resSVD = svd.getSolver().solve(b);
+        Assert.assertEquals(rm.getColumnDimension(), resSVD.length);
+
+        // check the solution really minimizes the residuals
+        double svdMinResidual = residual(rm, resSVD, b);
+        double epsilon = 2 * Math.ulp(svdMinResidual);
+        double h = 0.1;
+        int    k = 3;
+        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[] x = new double[] { resSVD[0] + d0, resSVD[1] + d1, resSVD[2]
+ d2 };
+                    Assert.assertTrue((residual(rm, x, b) - svdMinResidual) > -epsilon);
+                }
+            }
+        }
+
+    }
 
     @Test
     public void testMath320B() {
@@ -195,18 +215,17 @@
         });
         SingularValueDecomposition svd =
             new SingularValueDecompositionImpl(rm);
-        regularElements(svd.getU());
-        regularElements(svd.getVT());
+        RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
+        Assert.assertEquals(0.0, recomposed.subtract(rm).getNorm(), 2.0e-15);
     }
 
-    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));
-            }
+    private double residual(RealMatrix a, double[] x, double[] b) {
+        double[] ax = a.operate(x);
+        double sum = 0;
+        for (int i = 0; i < ax.length; ++i) {
+            sum += (ax[i] - b[i]) * (ax[i] - b[i]);
         }
+        return Math.sqrt(sum);
     }
 
 }



Mime
View raw message