svn commit: r885279 - /commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java

[lu...@apache.org]

Author: luc
Date: Sun Nov 29 21:53:36 2009
New Revision: 885279

URL: http://svn.apache.org/viewvc?rev=885279&view=rev
Log:
Prevent NaN to occur for singular matrices
Numerical inaccuracies in the underlying eigendecomposition could induce
very small negative eigenvalues, so the square root produced NaNs. The
eigenvalues really cannot be negative, so it is safe to replace the negative
ones by 0.
There are remaining problems with singular matrices:
 - the singular vectors also contain NaNs
 - the solver does not really work in least square sense and
   complain about singular matrices
JIRA: MATH-320 

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.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=885279&r1=885278&r2=885279&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 Sun Nov 29 21:53:36 2009
@@ -115,7 +115,8 @@
                                        MathUtils.SAFE_MIN);
         singularValues = eigenDecomposition.getRealEigenvalues();
         for (int i = 0; i < singularValues.length; ++i) {
-            singularValues[i] = Math.sqrt(singularValues[i]);
+            final double si = singularValues[i];
+            singularValues[i] = (si < 0) ? 0.0 : Math.sqrt(si);
         }
 
     }
@@ -133,14 +134,15 @@
                 double[] ei1 = eData[0];
                 iData[0] = ei1;
                 for (int i = 0; i < n - 1; ++i) {
-                    // compute Bt.E.S^(-1) where E is the eigenvectors matrix
+                    // compute B.E.S^(-1) where E is the eigenvectors matrix
                     // we reuse the array from matrix E to store the result
+                    final double mi = mainBidiagonal[i];
+                    final double si = secondaryBidiagonal[i];
                     final double[] ei0 = ei1;
                     ei1 = eData[i + 1];
                     iData[i + 1] = ei1;
                     for (int j = 0; j < n; ++j) {
-                        ei0[j] = (mainBidiagonal[i] * ei0[j] +
-                                  secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
+                        ei0[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
                     }
                 }
                 // last row
@@ -215,12 +217,13 @@
                 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
+                    final double mi = mainBidiagonal[i];
+                    final double si = secondaryBidiagonal[i];
                     final double[] ei0 = ei1;
                     ei1 = eData[i + 1];
                     iData[i + 1] = ei1;
                     for (int j = 0; j < m; ++j) {
-                        ei0[j] = (mainBidiagonal[i] * ei0[j] +
-                                  secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
+                        ei0[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
                     }
                 }
                 // last row