implement JacobiSVD::solve() and expand the unit test

Eigen/src/SVD/JacobiSVD.h

@@ -31,16 +31,6 @@ template<typename MatrixType, int QRPreconditioner,
          bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
 struct ei_svd_precondition_2x2_block_to_be_real {};
 
-template<typename MatrixType, int QRPreconditioner,
-         bool PossiblyMoreRowsThanCols = (MatrixType::RowsAtCompileTime == Dynamic)
-                                         || (MatrixType::RowsAtCompileTime > MatrixType::ColsAtCompileTime) >
-struct ei_svd_precondition_if_more_rows_than_cols;
-
-template<typename MatrixType, int QRPreconditioner,
-         bool PossiblyMoreColsThanRows = (MatrixType::ColsAtCompileTime == Dynamic)
-                                         || (MatrixType::ColsAtCompileTime > MatrixType::RowsAtCompileTime) >
-struct ei_svd_precondition_if_more_cols_than_rows;
-
 
 /*** QR preconditioners (R-SVD) ***/
 
@@ -81,8 +71,6 @@ struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner,
   {
     if(matrix.rows() > matrix.cols())
     {
-      ei_assert(!svd.m_computeThinU && "JacobiSVD: can't compute a thin U with the FullPivHouseholderQR preconditioner. "
-                                       "Use the ColPivHouseholderQR preconditioner instead.");
       FullPivHouseholderQR<MatrixType> qr(matrix);
       svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
       if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ();
@@ -100,8 +88,6 @@ struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner,
   {
     if(matrix.cols() > matrix.rows())
     {
-      ei_assert(!svd.m_computeThinV && "JacobiSVD: can't compute a thin V with the FullPivHouseholderQR preconditioner. "
-                                       "Use the ColPivHouseholderQR preconditioner instead.");
       typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
                      MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
               TransposeTypeWithSameStorageOrder;
@@ -233,9 +219,11 @@ struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, Precon
   *
   * \sa MatrixBase::jacobiSvd()
   */
-template<typename MatrixType, int QRPreconditioner> class JacobiSVD
+template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
 {
-  private:
+  public:
+
+    typedef _MatrixType MatrixType;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef typename MatrixType::Index Index;
@@ -262,8 +250,6 @@ template<typename MatrixType, int QRPreconditioner> class JacobiSVD
                    MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
             WorkMatrixType;
 
-  public:
-
     /** \brief Default Constructor.
       *
       * The default constructor is useful in cases in which the user intends to
@@ -352,6 +338,33 @@ template<typename MatrixType, int QRPreconditioner> class JacobiSVD
     inline bool computeU() const { return m_computeFullU || m_computeThinU; }
     inline bool computeV() const { return m_computeFullV || m_computeThinV; }
 
+    /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
+      *
+      * \param b the right-hand-side of the equation to solve.
+      *
+      * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V,
+      * 
+      * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
+      * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
+      */
+    template<typename Rhs>
+    inline const ei_solve_retval<JacobiSVD, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
+      ei_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
+      return ei_solve_retval<JacobiSVD, Rhs>(*this, b.derived());
+    }
+
+    Index nonzeroSingularValues() const
+    {
+      ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
+      return m_nonzeroSingularValues;
+    }
+
+    inline Index rows() const { return m_rows; }
+    inline Index cols() const { return m_cols; }
+
   protected:
     MatrixUType m_matrixU;
     MatrixVType m_matrixV;
@@ -360,10 +373,11 @@ template<typename MatrixType, int QRPreconditioner> class JacobiSVD
     bool m_isInitialized;
     bool m_computeFullU, m_computeThinU;
     bool m_computeFullV, m_computeThinV;
+    Index m_nonzeroSingularValues, m_rows, m_cols;
 
-    template<typename _MatrixType, int _QRPreconditioner, bool _IsComplex>
+    template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
     friend struct ei_svd_precondition_2x2_block_to_be_real;
-    template<typename _MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
+    template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
     friend struct ei_qr_preconditioner_impl;
 };
 
@@ -457,9 +471,15 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
   ei_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
   ei_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
             "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
-  Index rows = matrix.rows();
+  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
-  Index cols = matrix.cols();
+  {
-  Index diagSize = std::min(rows, cols);
+    ei_assert(!(m_computeThinU || m_computeThinV) &&
+              "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
+              "Use the ColPivHouseholderQR preconditioner instead.");
+  }
+  m_rows = matrix.rows();
+  m_cols = matrix.cols();
+  Index diagSize = std::min(m_rows, m_cols);
   m_singularValues.resize(diagSize);
   const RealScalar precision = 2 * NumTraits<Scalar>::epsilon();
 
@@ -467,10 +487,10 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
   && !ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, PreconditionIfMoreRowsThanCols>::run(*this, matrix))
   {
     m_workMatrix = matrix.block(0,0,diagSize,diagSize);
-    if(m_computeFullU) m_matrixU.setIdentity(rows,rows);
+    if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
-    if(m_computeThinU) m_matrixU.setIdentity(rows,diagSize);
+    if(m_computeThinU) m_matrixU.setIdentity(m_rows,diagSize);
-    if(m_computeFullV) m_matrixV.setIdentity(cols,cols);
+    if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
-    if(m_computeThinV) m_matrixV.setIdentity(diagSize,cols);
+    if(m_computeThinV) m_matrixV.setIdentity(diagSize,m_cols);
   }
 
   bool finished = false;
@@ -507,10 +527,17 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
     if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
   }
 
+  m_nonzeroSingularValues = diagSize;
+
   for(Index i = 0; i < diagSize; i++)
   {
     Index pos;
-    m_singularValues.tail(diagSize-i).maxCoeff(&pos);
+    RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize-i).maxCoeff(&pos);
+    if(maxRemainingSingularValue == RealScalar(0))
+    {
+      m_nonzeroSingularValues = i;
+      break;
+    }
     if(pos)
     {
       pos += i;
@@ -523,4 +550,33 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
   m_isInitialized = true;
   return *this;
 }
+
+template<typename _MatrixType, int QRPreconditioner, typename Rhs>
+struct ei_solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
+  : ei_solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
+{
+  typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
+  EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    ei_assert(rhs().rows() == dec().rows());
+
+    // A = U S V^*
+    // So A^{-1} = V S^{-1} U^*
+
+    Index diagSize = std::min(dec().rows(), dec().cols());
+    typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
+
+    Index nonzeroSingVals = dec().nonzeroSingularValues();
+    invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
+    invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
+
+    dst = dec().matrixV().leftCols(diagSize)
+        * invertedSingVals.asDiagonal()
+        * dec().matrixU().leftCols(diagSize).adjoint()
+        * rhs();
+  }
+};
+
 #endif // EIGEN_JACOBISVD_H

test/jacobisvd.cpp

@@ -79,6 +79,29 @@ void jacobisvd_compare_to_full(const MatrixType& m,
     VERIFY_IS_EQUAL(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
 }
 
+template<typename MatrixType, int QRPreconditioner>
+void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::Index Index;
+  Index rows = m.rows();
+  Index cols = m.cols();
+
+  enum {
+    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = MatrixType::ColsAtCompileTime
+  };
+
+  typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
+  typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
+
+  RhsType rhs = RhsType::Random(rows, ei_random<Index>(1, cols));
+  JacobiSVD<MatrixType, QRPreconditioner> svd(m, 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);
+}
+
 template<typename MatrixType, int QRPreconditioner>
 void jacobisvd_test_all_computation_options(const MatrixType& m)
 {
@@ -87,6 +110,7 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
   JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV);
 
   jacobisvd_check_full(m, fullSvd);
+  jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV);
 
   if(QRPreconditioner == FullPivHouseholderQRPreconditioner)
     return;
@@ -102,6 +126,9 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
     jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd);
     jacobisvd_compare_to_full(m, ComputeThinU             , fullSvd);
     jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd);
+    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV);
+    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV);
+    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV);
   }
 }
 
@@ -109,29 +136,74 @@ template<typename MatrixType>
 void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
 {
   MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a;
+
   jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m);
   jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m);
   jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m);
   jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m);
 }
 
-template<typename MatrixType> void jacobisvd_verify_assert()
+template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
 {
-  MatrixType tmp;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::Index Index;
+  Index rows = m.rows();
+  Index cols = m.cols();
 
-  JacobiSVD<MatrixType> svd;
+  enum {
-  //VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp))
+    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = MatrixType::ColsAtCompileTime
+  };
+
+  typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
+
+  RhsType rhs(rows);
+
+  JacobiSVD<MatrixType, HouseholderQRPreconditioner> svd;
   VERIFY_RAISES_ASSERT(svd.matrixU())
   VERIFY_RAISES_ASSERT(svd.singularValues())
   VERIFY_RAISES_ASSERT(svd.matrixV())
-  /*VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp))
+  VERIFY_RAISES_ASSERT(svd.solve(rhs))
-  VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp))
+
-  VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp))
+  MatrixType a = MatrixType::Zero(rows, cols);
-  VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))*/
+  a.setZero();
+  svd.compute(a, 0);
+  VERIFY_RAISES_ASSERT(svd.matrixU())
+  VERIFY_RAISES_ASSERT(svd.matrixV())
+  svd.singularValues();
+  VERIFY_RAISES_ASSERT(svd.solve(rhs))
+
+  if (ColsAtCompileTime == Dynamic)
+  {
+    svd.compute(a, ComputeThinU);
+    svd.matrixU();
+    VERIFY_RAISES_ASSERT(svd.matrixV())
+    VERIFY_RAISES_ASSERT(svd.solve(rhs))
+
+    svd.compute(a, ComputeThinV);
+    svd.matrixV();
+    VERIFY_RAISES_ASSERT(svd.matrixU())
+    VERIFY_RAISES_ASSERT(svd.solve(rhs))
+
+    JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr;
+    VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV))
+    VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV))
+    VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV))
+  }
+  else
+  {
+    VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
+    VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
+  }
 }
 
 void test_jacobisvd()
 {
+  CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
+  CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) ));
+  CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) ));
+  CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) ));
+
   for(int i = 0; i < g_repeat; i++) {
     Matrix2cd m;
     m << 0, 1,
@@ -149,17 +221,14 @@ void test_jacobisvd()
     CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() ));
     CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) ));
 
-    CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(50,50)) ));
+    int r = ei_random<int>(1, 30),
-    CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(14,7)) ));
+        c = ei_random<int>(1, 30);
+    CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
+    CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
+    (void) r;
+    (void) c;
   }
-  CALL_SUBTEST_9(( jacobisvd<MatrixXf>(MatrixXf(300,200)) ));
-  CALL_SUBTEST_10(( jacobisvd<MatrixXcd>(MatrixXcd(100,150)) ));
 
-  CALL_SUBTEST_3(( jacobisvd_verify_assert<Matrix3f>() ));
+  CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(ei_random<int>(200, 300), ei_random<int>(200, 300))) ));
-  CALL_SUBTEST_3(( jacobisvd_verify_assert<Matrix3d>() ));
+  CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(ei_random<int>(100, 120), ei_random<int>(100, 120))) ));
-  CALL_SUBTEST_9(( jacobisvd_verify_assert<MatrixXf>() ));
-  CALL_SUBTEST_11(( jacobisvd_verify_assert<MatrixXd>() ));
-
-  // Test problem size constructors
-  CALL_SUBTEST_12( JacobiSVD<MatrixXf>(10, 20) );
 }