JacobiSVD: move from Lapack to Matlab strategy for the default threshold

1 files changed, 53 insertions, 0 deletions

diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index 76157c30f..e378de477 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -84,6 +84,59 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
   SolutionType x = svd.solve(rhs);
   // evaluate normal equation which works also for least-squares solutions
   VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
+  
+  // check minimal norm solutions
+  {
+    // generate a full-rank m x n problem with m<n
+    enum {
+      RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1,
+      RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1
+    };
+    typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2;
+    typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2;
+    typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T;
+    Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2);
+    MatrixType2 m2(rank,cols);
+    int guard = 0;
+    do {
+      m2.setRandom();
+    } while(m2.jacobiSvd().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());
+    Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2);
+    tmp.conservativeResize(cols);
+    tmp.tail(cols-rank).setZero();
+    SolutionType x21 = qr.householderQ() * tmp;
+    // now check with SVD
+    JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions);
+    SolutionType x22 = svd2.solve(rhs2);
+    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;
+    Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3);
+    Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank);
+    MatrixType3 m3 = C * m2;
+    RhsType3 rhs3 = C * rhs2;
+    JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions);
+    SolutionType x3 = svd3.solve(rhs3);
+    if(svd3.rank()!=rank) {
+      std::cout << m3 << "\n\n";
+      std::cout << svd3.singularValues().transpose() << "\n";
+    std::cout << svd3.rank() << " == " << rank << "\n";
+    std::cout << x21.norm() << " == " << x3.norm() << "\n";
+    }
+//     VERIFY_IS_APPROX(m3*x3, rhs3);
+    VERIFY_IS_APPROX(m3*x21, rhs3);
+    VERIFY_IS_APPROX(m2*x3, rhs2);
+    
+    VERIFY_IS_APPROX(x21, x3);
+  }
 }
 
 template<typename MatrixType, int QRPreconditioner>