Add condition estimation to Cholesky (LLT) factorization.

4 files changed, 111 insertions, 36 deletions

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 74cf5bfe1..b55c5bebf 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -10,7 +10,7 @@
 #ifndef EIGEN_LLT_H
 #define EIGEN_LLT_H
 
-namespace Eigen { 
+namespace Eigen {
 
 namespace internal{
 template<typename MatrixType, int UpLo> struct LLT_Traits;
@@ -40,7 +40,7 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
   *
   * Example: \include LLT_example.cpp
   * Output: \verbinclude LLT_example.out
-  *    
+  *
   * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
   */
  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
@@ -135,6 +135,16 @@ template<typename _MatrixType, int _UpLo> class LLT
     template<typename InputType>
     LLT& compute(const EigenBase<InputType>& matrix);
 
+    /** \returns an estimate of the reciprocal condition number of the matrix of
+     *  which *this is the Cholesky decomposition.
+     */
+    RealScalar rcond() const
+    {
+      eigen_assert(m_isInitialized && "LLT is not initialized.");
+      eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
+      return ConditionEstimator<LLT<MatrixType, UpLo>, true >::rcond(m_l1_norm, *this);
+    }
+
     /** \returns the LLT decomposition matrix
       *
       * TODO: document the storage layout
@@ -164,7 +174,7 @@ template<typename _MatrixType, int _UpLo> class LLT
 
     template<typename VectorType>
     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
-    
+
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename RhsType, typename DstType>
     EIGEN_DEVICE_FUNC
@@ -172,17 +182,18 @@ template<typename _MatrixType, int _UpLo> class LLT
     #endif
 
   protected:
-    
+
     static void check_template_parameters()
     {
       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
     }
-    
+
     /** \internal
       * Used to compute and store L
       * The strict upper part is not used and even not initialized.
       */
     MatrixType m_matrix;
+    RealScalar m_l1_norm;
     bool m_isInitialized;
     ComputationInfo m_info;
 };
@@ -268,7 +279,7 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
   static Index unblocked(MatrixType& mat)
   {
     using std::sqrt;
-    
+
     eigen_assert(mat.rows()==mat.cols());
     const Index size = mat.rows();
     for(Index k = 0; k < size; ++k)
@@ -328,7 +339,7 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
   }
 };
-  
+
 template<typename Scalar> struct llt_inplace<Scalar, Upper>
 {
   typedef typename NumTraits<Scalar>::Real RealScalar;
@@ -387,12 +398,32 @@ template<typename InputType>
 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
 {
   check_template_parameters();
-  
+
   eigen_assert(a.rows()==a.cols());
   const Index size = a.rows();
   m_matrix.resize(size, size);
   m_matrix = a.derived();
 
+  // Compute matrix L1 norm = max abs column sum.
+  m_l1_norm = RealScalar(0);
+  if (_UpLo == Lower) {
+    for (int col = 0; col < size; ++col) {
+      const RealScalar abs_col_sum = m_matrix.col(col).tail(size - col).cwiseAbs().sum() +
+          m_matrix.row(col).tail(col).cwiseAbs().sum();
+      if (abs_col_sum > m_l1_norm) {
+        m_l1_norm = abs_col_sum;
+      }
+    }
+  } else {
+    for (int col = 0; col < a.cols(); ++col) {
+      const RealScalar abs_col_sum = m_matrix.col(col).tail(col).cwiseAbs().sum() +
+          m_matrix.row(col).tail(size - col).cwiseAbs().sum();
+      if (abs_col_sum > m_l1_norm) {
+        m_l1_norm = abs_col_sum;
+      }
+    }
+  }
+
   m_isInitialized = true;
   bool ok = Traits::inplace_decomposition(m_matrix);
   m_info = ok ? Success : NumericalIssue;
@@ -419,7 +450,7 @@ LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, c
 
   return *this;
 }
- 
+
 #ifndef EIGEN_PARSED_BY_DOXYGEN
 template<typename _MatrixType,int _UpLo>
 template<typename RhsType, typename DstType>
@@ -431,7 +462,7 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
 #endif
 
 /** \internal use x = llt_object.solve(x);
-  * 
+  *
   * This is the \em in-place version of solve().
   *
   * \param bAndX represents both the right-hand side matrix b and result x.
@@ -483,7 +514,7 @@ SelfAdjointView<MatrixType, UpLo>::llt() const
   return LLT<PlainObject,UpLo>(m_matrix);
 }
 #endif // __CUDACC__
-  
+
 } // end namespace Eigen
 
 #endif // EIGEN_LLT_H
diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h
index b65306d56..dca3da417 100644
--- a/Eigen/src/Core/ConditionEstimator.h
+++ b/Eigen/src/Core/ConditionEstimator.h
@@ -13,11 +13,11 @@
 namespace Eigen {
 
 namespace internal {
-template <typename Decomposition, bool IsComplex>
+template <typename Decomposition, bool IsSelfAdjoint, bool IsComplex>
 struct EstimateInverseMatrixL1NormImpl {};
 }  // namespace internal
 
-template <typename Decomposition>
+template <typename Decomposition, bool IsSelfAdjoint = false>
 class ConditionEstimator {
  public:
   typedef typename Decomposition::MatrixType MatrixType;
@@ -101,7 +101,8 @@ class ConditionEstimator {
       return 0;
     }
     return internal::EstimateInverseMatrixL1NormImpl<
-        Decomposition, NumTraits<Scalar>::IsComplex>::compute(dec);
+        Decomposition, IsSelfAdjoint,
+        NumTraits<Scalar>::IsComplex != 0>::compute(dec);
   }
 
   /**
@@ -116,9 +117,27 @@ class ConditionEstimator {
 
 namespace internal {
 
+template <typename Decomposition, typename Vector, bool IsSelfAdjoint = false>
+struct solve_helper {
+  static inline Vector solve_adjoint(const Decomposition& dec,
+                                     const Vector& v) {
+    return dec.adjoint().solve(v);
+  }
+};
+
+// Partial specialization for self_adjoint matrices.
+template <typename Decomposition, typename Vector>
+struct solve_helper<Decomposition, Vector, true> {
+  static inline Vector solve_adjoint(const Decomposition& dec,
+                                     const Vector& v) {
+    return dec.solve(v);
+  }
+};
+
+
 // Partial specialization for real matrices.
-template <typename Decomposition>
-struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> {
+template <typename Decomposition, bool IsSelfAdjoint>
+struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, false> {
   typedef typename Decomposition::MatrixType MatrixType;
   typedef typename internal::traits<MatrixType>::Scalar Scalar;
   typedef typename internal::plain_col_type<MatrixType>::type Vector;
@@ -152,7 +171,7 @@ struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> {
     int old_v_max_abs_index = v_max_abs_index;
     for (int k = 0; k < 4; ++k) {
       // argmax |inv(matrix)^T * sign_vector|
-      v = dec.adjoint().solve(sign_vector);
+      v = solve_helper<Decomposition, Vector, IsSelfAdjoint>::solve_adjoint(dec, sign_vector);
       v.cwiseAbs().maxCoeff(&v_max_abs_index);
       if (v_max_abs_index == old_v_max_abs_index) {
         // Break if the solution stagnated.
@@ -200,8 +219,8 @@ struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> {
 };
 
 // Partial specialization for complex matrices.
-template <typename Decomposition>
-struct EstimateInverseMatrixL1NormImpl<Decomposition, 1> {
+template <typename Decomposition, bool IsSelfAdjoint>
+struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, true> {
   typedef typename Decomposition::MatrixType MatrixType;
   typedef typename internal::traits<MatrixType>::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
@@ -238,7 +257,7 @@ struct EstimateInverseMatrixL1NormImpl<Decomposition, 1> {
       RealVector abs_v = v.cwiseAbs();
       const Vector psi =
           (abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones);
-      v = dec.adjoint().solve(psi);
+      v = solve_helper<Decomposition, Vector, IsSelfAdjoint>::solve_adjoint(dec, psi);
       const RealVector z = v.real();
       z.cwiseAbs().maxCoeff(&v_max_abs_index);
       if (v_max_abs_index == old_v_max_abs_index) {
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index d652af5bf..8a21cdbd5 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -17,6 +17,12 @@
 #include <Eigen/Cholesky>
 #include <Eigen/QR>
 
+template<typename MatrixType, int UpLo>
+typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
+  MatrixType symm = m.template selfadjointView<UpLo>();
+  return symm.cwiseAbs().colwise().sum().maxCoeff();
+}
+
 template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
 {
   typedef typename MatrixType::Scalar Scalar;
@@ -77,7 +83,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   {
     SquareMatrixType symmUp = symm.template triangularView<Upper>();
     SquareMatrixType symmLo = symm.template triangularView<Lower>();
-    
+
     LLT<SquareMatrixType,Lower> chollo(symmLo);
     VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
     vecX = chollo.solve(vecB);
@@ -85,6 +91,14 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     matX = chollo.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
 
+    // Verify that the estimated condition number is within a factor of 10 of the
+    // truth.
+    const MatrixType symmLo_inverse = chollo.solve(MatrixType::Identity(rows,cols));
+    RealScalar rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Lower>(symmLo)) /
+                             matrix_l1_norm<MatrixType, Lower>(symmLo_inverse);
+    RealScalar rcond_est = chollo.rcond();
+    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+
     // test the upper mode
     LLT<SquareMatrixType,Upper> cholup(symmUp);
     VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
@@ -93,6 +107,15 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     matX = cholup.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
 
+    // Verify that the estimated condition number is within a factor of 10 of the
+    // truth.
+    const MatrixType symmUp_inverse = cholup.solve(MatrixType::Identity(rows,cols));
+    rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
+                             matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
+    rcond_est = cholup.rcond();
+    VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
+
+
     MatrixType neg = -symmLo;
     chollo.compute(neg);
     VERIFY(chollo.info()==NumericalIssue);
@@ -101,7 +124,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
     VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU()));
     VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL()));
-    
+
     // test some special use cases of SelfCwiseBinaryOp:
     MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
     m2 = m1;
@@ -167,7 +190,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     // restore
     if(sign == -1)
       symm = -symm;
-    
+
     // check matrices coming from linear constraints with Lagrange multipliers
     if(rows>=3)
     {
@@ -183,7 +206,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       vecX = ldltlo.solve(vecB);
       VERIFY_IS_APPROX(A * vecX, vecB);
     }
-    
+
     // check non-full rank matrices
     if(rows>=3)
     {
@@ -199,7 +222,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       vecX = ldltlo.solve(vecB);
       VERIFY_IS_APPROX(A * vecX, vecB);
     }
-    
+
     // check matrices with a wide spectrum
     if(rows>=3)
     {
@@ -225,7 +248,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       {
         RealScalar large_tol =  std::sqrt(test_precision<RealScalar>());
         VERIFY((A * vecX).isApprox(vecB, large_tol));
-        
+
         ++g_test_level;
         VERIFY_IS_APPROX(A * vecX,vecB);
         --g_test_level;
@@ -314,14 +337,14 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
 }
 
 // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
-// This test checks that LDLT reports correctly that matrix is indefinite. 
+// This test checks that LDLT reports correctly that matrix is indefinite.
 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
 template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
 {
   eigen_assert(m.rows() == 2 && m.cols() == 2);
   MatrixType mat;
   LDLT<MatrixType> ldlt(2);
-  
+
   {
     mat << 1, 0, 0, -1;
     ldlt.compute(mat);
@@ -384,11 +407,11 @@ void test_cholesky()
     CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
     CALL_SUBTEST_4( cholesky(Matrix3f()) );
     CALL_SUBTEST_5( cholesky(Matrix4d()) );
-    
-    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);    
+
+    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
     CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) );
     TEST_SET_BUT_UNUSED_VARIABLE(s)
-    
+
     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
     CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
     TEST_SET_BUT_UNUSED_VARIABLE(s)
@@ -402,6 +425,6 @@ void test_cholesky()
   // Test problem size constructors
   CALL_SUBTEST_9( LLT<MatrixXf>(10) );
   CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
-  
+
   TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries)
 }
diff --git a/test/lu.cpp b/test/lu.cpp
index 9e8059f58..53b3fcee4 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -151,10 +151,11 @@ template<typename MatrixType> void lu_invertible()
   MatrixType m1_inverse = lu.inverse();
   VERIFY_IS_APPROX(m2, m1_inverse*m3);
 
-  // Test condition number estimation.
+  // Verify that the estimated condition number is within a factor of 10 of the
+  // truth.
   RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse);
-  // Verify that the estimate is within a factor of 10 of the truth.
-  VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10);
+  const RealScalar rcond_est = lu.rcond();
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
 
   // test solve with transposed
   lu.template _solve_impl_transposed<false>(m3, m2);
@@ -199,7 +200,8 @@ template<typename MatrixType> void lu_partial_piv()
   // Test condition number estimation.
   RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse);
   // Verify that the estimate is within a factor of 10 of the truth.
-  VERIFY(plu.rcond() > rcond / 10 && plu.rcond() < rcond * 10);
+  const RealScalar rcond_est = plu.rcond();
+  VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
 
   // test solve with transposed
   plu.template _solve_impl_transposed<false>(m3, m2);