Various numerical fixes in D&C SVD: I cannot make it fail with double, but still need to tune for single precision, and carefully test with duplicated singular values

test/svd_common.h

@@ -146,6 +146,7 @@ void svd_min_norm(const MatrixType& m, unsigned int computationOptions)
     m2.setRandom();
   } while(SVD_FOR_MIN_NORM(MatrixType2)(m2).setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10);
   VERIFY(guard<10);
+
   RhsType2 rhs2 = RhsType2::Random(rank);
   // use QR to find a reference minimal norm solution
   HouseholderQR<MatrixType2T> qr(m2.adjoint());
@@ -159,7 +160,7 @@ void svd_min_norm(const MatrixType& m, unsigned int computationOptions)
   VERIFY_IS_APPROX(m2*x21, rhs2);
   VERIFY_IS_APPROX(m2*x22, rhs2);
   VERIFY_IS_APPROX(x21, x22);
-  
+
   // Now check with a rank deficient matrix
   typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3;
   typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3;
@@ -172,7 +173,6 @@ void svd_min_norm(const MatrixType& m, unsigned int computationOptions)
   VERIFY_IS_APPROX(m3*x3, rhs3);
   VERIFY_IS_APPROX(m3*x21, rhs3);
   VERIFY_IS_APPROX(m2*x3, rhs2);
-  
   VERIFY_IS_APPROX(x21, x3);
 }
 
@@ -209,7 +209,7 @@ void svd_test_all_computation_options(const MatrixType& m, bool full_only)
     CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeThinV) ));
     CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeFullV) ));
     CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeThinV) ));
-    
+
     CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeThinV) ));
     CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeFullV) ));
     CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeThinV) ));

unsupported/Eigen/src/BDCSVD/BDCSVD.h

@@ -272,8 +272,8 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
     }
   }
 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
-  std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
+//   std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
-  std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
+//   std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
 #endif
   if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
   else              copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
@@ -612,22 +612,23 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
   Index n = col0.size();
   Index actual_n = n;
   while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
-  Index m = 0;
+//   Index m = 0;
-  Array<Index,1,Dynamic> perm(actual_n);
+//   Array<Index,1,Dynamic> perm(actual_n);
-  {
+//   {
-    for(Index k=0;k<actual_n;++k)
+//     for(Index k=0;k<actual_n;++k)
-    {
+//     {
-      if(col0(k)!=0)
+//       if(col0(k)!=0)
-        perm(m++) = k;
+//         perm(m++) = k;
-    }
+//     }
-  }
+//   }
-  perm.conservativeResize(m);
+//   perm.conservativeResize(m);
   
   for (Index k = 0; k < n; ++k)
   {
-    if (col0(k) == 0)
+    if (col0(k) == 0 || actual_n==1)
     {
-      // entry is deflated, so singular value is on diagonal
+      // if col0(k) == 0, then entry is deflated, so singular value is on diagonal
+      // if actual_n==1, then the deflated problem is already diagonalized
       singVals(k) = diag(k);
       mus(k) = 0;
       shifts(k) = diag(k);
@@ -636,7 +637,16 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
 
     // otherwise, use secular equation to find singular value
     RealScalar left = diag(k);
-    RealScalar right = (k != actual_n-1) ? diag(k+1) : (diag(actual_n-1) + col0.matrix().norm());
+    RealScalar right; // was: = (k != actual_n-1) ? diag(k+1) : (diag(actual_n-1) + col0.matrix().norm());
+    if(k==actual_n-1)
+      right = (diag(actual_n-1) + col0.matrix().norm());
+    else
+    {
+      // Skip deflated singular values
+      Index l = k+1;
+      while(col0(l)==0) { ++l; eigen_internal_assert(l<actual_n); }
+      right = diag(l);
+    }
 
     // first decide whether it's closer to the left end or the right end
     RealScalar mid = left + (right-left) / 2;
@@ -653,7 +663,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
     {
       muPrev = (right - left) * 0.1;
       if (k == actual_n-1) muCur = right - left;
-      else          muCur = (right - left) * 0.5; 
+      else                 muCur = (right - left) * 0.5; 
     }
     else
     {
@@ -671,7 +681,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
 
     // rational interpolation: fit a function of the form a / mu + b through the two previous
     // iterates and use its zero to compute the next iterate
-    bool useBisection = false;
+    bool useBisection = fPrev*fCur>0;
     while (abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
     {
       ++m_numIters;
@@ -684,32 +694,36 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
       muCur = -a / b;
       fCur = secularEq(muCur, col0, diag, diagShifted, shift, actual_n);
       
-      if (shift == left && (muCur < 0 || muCur > right - left)) useBisection = true;
+      if (shift == left  && (muCur < 0 || muCur > right - left)) useBisection = true;
       if (shift == right && (muCur < -(right - left) || muCur > 0)) useBisection = true;
     }
 
     // fall back on bisection method if rational interpolation did not work
     if (useBisection)
     {
+#ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
+      std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
+#endif
       RealScalar leftShifted, rightShifted;
       if (shift == left)
       {
-        leftShifted = 1e-30;
+        leftShifted = RealScalar(1)/NumTraits<RealScalar>::highest();
-        if (k == 0) rightShifted = right - left;
+        // I don't understand why the case k==0 would be special there:
-        else        rightShifted = (right - left) * 0.6; // theoretically we can take 0.5, but let's be safe
+        // if (k == 0) rightShifted = right - left; else 
+        rightShifted = (right - left) * 0.6; // theoretically we can take 0.5, but let's be safe
       }
       else
       {
         leftShifted = -(right - left) * 0.6;
-        rightShifted = -1e-30;
+        rightShifted = -RealScalar(1)/NumTraits<RealScalar>::highest();
       }
       
       RealScalar fLeft = secularEq(leftShifted, col0, diag, diagShifted, shift, actual_n);
       RealScalar fRight = secularEq(rightShifted, col0, diag, diagShifted, shift, actual_n);
 
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-      if(fLeft * fRight>0)
+      if(fLeft * fRight>=0)
-        std::cout << fLeft << " * " << fRight << " == " << fLeft * fRight << "  ;  " << left << " - " << right << " -> " <<  leftShifted << " " << rightShifted << "   shift=" << shift << "\n";
+        std::cout << k << " : " <<  fLeft << " * " << fRight << " == " << fLeft * fRight << "  ;  " << left << " - " << right << " -> " <<  leftShifted << " " << rightShifted << "   shift=" << shift << "\n";
 #endif
       eigen_internal_assert(fLeft * fRight < 0);
       
@@ -738,8 +752,9 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
 
     // perturb singular value slightly if it equals diagonal entry to avoid division by zero later
     // (deflation is supposed to avoid this from happening)
-    if (singVals[k] == left) singVals[k] *= 1 + NumTraits<RealScalar>::epsilon();
+    // - this does no seem to be necessary anymore -
-    if (singVals[k] == right) singVals[k] *= 1 - NumTraits<RealScalar>::epsilon();
+//     if (singVals[k] == left) singVals[k] *= 1 + NumTraits<RealScalar>::epsilon();
+//     if (singVals[k] == right) singVals[k] *= 1 - NumTraits<RealScalar>::epsilon();
   }
 }
 
@@ -989,7 +1004,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
 #endif
       deflation43(firstCol, shift, i, length);
     }
-    
+
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   assert(m_naiveU.allFinite());
   assert(m_naiveV.allFinite());
@@ -1027,7 +1042,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
       realInd[pos] = pos;
     }
     
-    for(Index i = 1; i < length - 1; i++)
+    for(Index i = 1; i < length; i++)
     {
       const Index pi = permutation[length - i];
       const Index J = realCol[pi];
@@ -1056,7 +1071,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
   
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
   for(int k=2;k<length;++k)
-    assert(diag(k-1)<diag(k) || diag(k)==0);
+    assert(diag(k-1)<=diag(k) || diag(k)==0);
 #endif
   
   //condition 4.4

unsupported/Eigen/src/BDCSVD/TODOBdcsvd.txt

@@ -1,29 +1,13 @@
 TO DO LIST
 
+- check more carefully single precision
+- check with duplicated singularvalues
+- no-malloc mode
 
+(optional optimization)
+- do all the allocations in the allocate part 
+- support static matrices
+- return a error at compilation time when using integer matrices (int, long, std::complex<int>, ...)
+- To solve the secular equation using FMM:
+http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
 
-(optional optimization) - do all the allocations in the allocate part 
-                        - support static matrices
-                        - return a error at compilation time when using integer matrices (int, long, std::complex<int>, ...)
-
-to finish the algorithm :
-			-implement the last part of the algorithm as described on the reference paper. 
-			    You may find more information on that part on this paper
-
-			-to replace the call to JacobiSVD at the end of the divide algorithm, just after the call to 
-			    deflation.
-
-(suggested step by step resolution)
-                       0) comment the call to Jacobi in the last part of the divide method and everything right after
-                               until the end of the method. What is commented can be a guideline to steps 3) 4) and 6)
-                       1) solve the secular equation (Characteristic equation) on the values that are not null (zi!=0 and di!=0), after the deflation
-                               wich should be uncommented in the divide method
-                       2) remember the values of the singular values that are already computed (zi=0)
-                       3) assign the singular values found in m_computed at the right places (with the ones found in step 2) )
-                               in decreasing order
-                       4) set the firstcol to zero (except the first element) in m_computed
-                       5) compute all the singular vectors when CompV is set to true and only the left vectors when
-                               CompV is set to false
-                       6) multiply naiveU and naiveV to the right by the matrices found, only naiveU when CompV is set to
-                               false, /!\ if CompU is false NaiveU has only 2 rows
-                       7) delete everything commented in step 0)

unsupported/Eigen/src/BDCSVD/doneInBDCSVD.txt

@@ -3,19 +3,6 @@ This unsupported package is about a divide and conquer algorithm to compute SVD.
 The implementation follows as closely as possible the following reference paper : 
 http://www.cs.yale.edu/publications/techreports/tr933.pdf
 
-The code documentation uses the same names for variables as the reference paper. The code, deflation included, is
+To solve the secular equation using FMM:
-working  but there are a few things that could be optimised as explained in the TODOBdsvd. 
-
-In the code comments were put at the line where would be the third step of the algorithm so one could simply add the call 
-of a function doing the last part of the algorithm and that would not require any knowledge of the part we implemented.
-
-In the TODOBdcsvd we explain what is the main difficulty of the last part and suggest a reference paper to help solve it.
-
-The implemented has trouble with fixed size matrices. 
-
-In the actual implementation, it returns matrices of zero when ask to do a svd on an int matrix. 
-
-
-Paper for the third part:
 http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
 

unsupported/test/bdcsvd.cpp

@@ -21,8 +21,7 @@
 
 
 #define SVD_DEFAULT(M) BDCSVD<M>
-// #define SVD_FOR_MIN_NORM(M) BDCSVD<M>
+#define SVD_FOR_MIN_NORM(M) BDCSVD<M>
-#define SVD_FOR_MIN_NORM(M) JacobiSVD<M,ColPivHouseholderQRPreconditioner>
 #include "../../test/svd_common.h"
 
 // Check all variants of JacobiSVD