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 struct LLT_Traits; @@ -40,7 +40,7 @@ template 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 class LLT template LLT& compute(const EigenBase& 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, true >::rcond(m_l1_norm, *this); + } + /** \returns the LLT decomposition matrix * * TODO: document the storage layout @@ -164,7 +174,7 @@ template class LLT template LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); - + #ifndef EIGEN_PARSED_BY_DOXYGEN template EIGEN_DEVICE_FUNC @@ -172,17 +182,18 @@ template 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 struct llt_inplace 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 struct llt_inplace return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); } }; - + template struct llt_inplace { typedef typename NumTraits::Real RealScalar; @@ -387,12 +398,32 @@ template LLT& LLT::compute(const EigenBase& 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 template @@ -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::llt() const return LLT(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 +template struct EstimateInverseMatrixL1NormImpl {}; } // namespace internal -template +template class ConditionEstimator { public: typedef typename Decomposition::MatrixType MatrixType; @@ -101,7 +101,8 @@ class ConditionEstimator { return 0; } return internal::EstimateInverseMatrixL1NormImpl< - Decomposition, NumTraits::IsComplex>::compute(dec); + Decomposition, IsSelfAdjoint, + NumTraits::IsComplex != 0>::compute(dec); } /** @@ -116,9 +117,27 @@ class ConditionEstimator { namespace internal { +template +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 +struct solve_helper { + static inline Vector solve_adjoint(const Decomposition& dec, + const Vector& v) { + return dec.solve(v); + } +}; + + // Partial specialization for real matrices. -template -struct EstimateInverseMatrixL1NormImpl { +template +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename internal::plain_col_type::type Vector; @@ -152,7 +171,7 @@ struct EstimateInverseMatrixL1NormImpl { 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::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 { }; // Partial specialization for complex matrices. -template -struct EstimateInverseMatrixL1NormImpl { +template +struct EstimateInverseMatrixL1NormImpl { typedef typename Decomposition::MatrixType MatrixType; typedef typename internal::traits::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -238,7 +257,7 @@ struct EstimateInverseMatrixL1NormImpl { 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::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 #include +template +typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { + MatrixType symm = m.template selfadjointView(); + return symm.cwiseAbs().colwise().sum().maxCoeff(); +} + template class CholType> void test_chol_update(const MatrixType& symm) { typedef typename MatrixType::Scalar Scalar; @@ -77,7 +83,7 @@ template void cholesky(const MatrixType& m) { SquareMatrixType symmUp = symm.template triangularView(); SquareMatrixType symmLo = symm.template triangularView(); - + LLT chollo(symmLo); VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); vecX = chollo.solve(vecB); @@ -85,6 +91,14 @@ template 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(symmLo)) / + matrix_l1_norm(symmLo_inverse); + RealScalar rcond_est = chollo.rcond(); + VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); + // test the upper mode LLT cholup(symmUp); VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); @@ -93,6 +107,15 @@ template 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(symmUp)) / + matrix_l1_norm(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 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 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 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 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 void cholesky(const MatrixType& m) { RealScalar large_tol = std::sqrt(test_precision()); VERIFY((A * vecX).isApprox(vecB, large_tol)); - + ++g_test_level; VERIFY_IS_APPROX(A * vecX,vecB); --g_test_level; @@ -314,14 +337,14 @@ template 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 void cholesky_definiteness(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); MatrixType mat; LDLT 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(1,EIGEN_TEST_MAX_SIZE); + + s = internal::random(1,EIGEN_TEST_MAX_SIZE); CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) ); TEST_SET_BUT_UNUSED_VARIABLE(s) - + s = internal::random(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(10) ); CALL_SUBTEST_9( LDLT(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 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(m3, m2); @@ -199,7 +200,8 @@ template 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(m3, m2);