GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add condition estimation to Cholesky (LLT) factorization.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2016-04-01 16:19:45 -07:00
parent fb8dccc23e
commit f54137606e
4 changed files with 111 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
Eigen/src/Cholesky/LLT.h
View File

@ -10,7 +10,7 @@
#ifndef EIGEN_LLT_H #ifndef EIGEN_LLT_H
#define EIGEN_LLT_H #define EIGEN_LLT_H
namespace Eigen { namespace Eigen {
namespace internal{ namespace internal{
template<typename MatrixType, int UpLo> struct LLT_Traits; 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 * Example: \include LLT_example.cpp
* Output: \verbinclude LLT_example.out * Output: \verbinclude LLT_example.out
* *
* \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT * \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) /* 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> template<typename InputType>
LLT& compute(const EigenBase<InputType>& matrix); 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 /** \returns the LLT decomposition matrix
* *
* TODO: document the storage layout * TODO: document the storage layout
@ -164,7 +174,7 @@ template<typename _MatrixType, int _UpLo> class LLT
template<typename VectorType> template<typename VectorType>
LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType> template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
@ -172,17 +182,18 @@ template<typename _MatrixType, int _UpLo> class LLT
#endif #endif
protected: protected:
static void check_template_parameters() static void check_template_parameters()
{ {
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
} }
/** \internal /** \internal
* Used to compute and store L * Used to compute and store L
* The strict upper part is not used and even not initialized. * The strict upper part is not used and even not initialized.
*/ */
MatrixType m_matrix; MatrixType m_matrix;
RealScalar m_l1_norm;
bool m_isInitialized; bool m_isInitialized;
ComputationInfo m_info; ComputationInfo m_info;
}; };
@ -268,7 +279,7 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
static Index unblocked(MatrixType& mat) static Index unblocked(MatrixType& mat)
{ {
using std::sqrt; using std::sqrt;
eigen_assert(mat.rows()==mat.cols()); eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows(); const Index size = mat.rows();
for(Index k = 0; k < size; ++k) 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); return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
} }
}; };
template<typename Scalar> struct llt_inplace<Scalar, Upper> template<typename Scalar> struct llt_inplace<Scalar, Upper>
{ {
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
@ -387,12 +398,32 @@ template<typename InputType>
LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a) LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
{ {
check_template_parameters(); check_template_parameters();
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
const Index size = a.rows(); const Index size = a.rows();
m_matrix.resize(size, size); m_matrix.resize(size, size);
m_matrix = a.derived(); 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; m_isInitialized = true;
bool ok = Traits::inplace_decomposition(m_matrix); bool ok = Traits::inplace_decomposition(m_matrix);
m_info = ok ? Success : NumericalIssue; m_info = ok ? Success : NumericalIssue;
@ -419,7 +450,7 @@ LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, c
return *this; return *this;
} }
#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType,int _UpLo> template<typename _MatrixType,int _UpLo>
template<typename RhsType, typename DstType> template<typename RhsType, typename DstType>
@ -431,7 +462,7 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
#endif #endif
/** \internal use x = llt_object.solve(x); /** \internal use x = llt_object.solve(x);
* *
* This is the \em in-place version of solve(). * This is the \em in-place version of solve().
* *
* \param bAndX represents both the right-hand side matrix b and result x. * \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); return LLT<PlainObject,UpLo>(m_matrix);
} }
#endif // __CUDACC__ #endif // __CUDACC__
} // end namespace Eigen } // end namespace Eigen
#endif // EIGEN_LLT_H #endif // EIGEN_LLT_H

37
Eigen/src/Core/ConditionEstimator.h
View File

@ -13,11 +13,11 @@
namespace Eigen { namespace Eigen {
namespace internal { namespace internal {
template <typename Decomposition, bool IsComplex> template <typename Decomposition, bool IsSelfAdjoint, bool IsComplex>
struct EstimateInverseMatrixL1NormImpl {}; struct EstimateInverseMatrixL1NormImpl {};
} // namespace internal } // namespace internal
template <typename Decomposition> template <typename Decomposition, bool IsSelfAdjoint = false>
class ConditionEstimator { class ConditionEstimator {
public: public:
typedef typename Decomposition::MatrixType MatrixType; typedef typename Decomposition::MatrixType MatrixType;
@ -101,7 +101,8 @@ class ConditionEstimator {
return 0; return 0;
} }
return internal::EstimateInverseMatrixL1NormImpl< 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 { 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. // Partial specialization for real matrices.
template <typename Decomposition> template <typename Decomposition, bool IsSelfAdjoint>
struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> { struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, false> {
typedef typename Decomposition::MatrixType MatrixType; typedef typename Decomposition::MatrixType MatrixType;
typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename internal::plain_col_type<MatrixType>::type Vector; 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; int old_v_max_abs_index = v_max_abs_index;
for (int k = 0; k < 4; ++k) { for (int k = 0; k < 4; ++k) {
// argmax |inv(matrix)^T * sign_vector| // 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); v.cwiseAbs().maxCoeff(&v_max_abs_index);
if (v_max_abs_index == old_v_max_abs_index) { if (v_max_abs_index == old_v_max_abs_index) {
// Break if the solution stagnated. // Break if the solution stagnated.
@ -200,8 +219,8 @@ struct EstimateInverseMatrixL1NormImpl<Decomposition, 0> {
}; };
// Partial specialization for complex matrices. // Partial specialization for complex matrices.
template <typename Decomposition> template <typename Decomposition, bool IsSelfAdjoint>
struct EstimateInverseMatrixL1NormImpl<Decomposition, 1> { struct EstimateInverseMatrixL1NormImpl<Decomposition, IsSelfAdjoint, true> {
typedef typename Decomposition::MatrixType MatrixType; typedef typename Decomposition::MatrixType MatrixType;
typedef typename internal::traits<MatrixType>::Scalar Scalar; typedef typename internal::traits<MatrixType>::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
@ -238,7 +257,7 @@ struct EstimateInverseMatrixL1NormImpl<Decomposition, 1> {
RealVector abs_v = v.cwiseAbs(); RealVector abs_v = v.cwiseAbs();
const Vector psi = const Vector psi =
(abs_v.array() == 0).select(v.cwiseQuotient(abs_v), ones); (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(); const RealVector z = v.real();
z.cwiseAbs().maxCoeff(&v_max_abs_index); z.cwiseAbs().maxCoeff(&v_max_abs_index);
if (v_max_abs_index == old_v_max_abs_index) { if (v_max_abs_index == old_v_max_abs_index) {

47
test/cholesky.cpp
View File

@ -17,6 +17,12 @@
#include <Eigen/Cholesky> #include <Eigen/Cholesky>
#include <Eigen/QR> #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) template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
@ -77,7 +83,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
{ {
SquareMatrixType symmUp = symm.template triangularView<Upper>(); SquareMatrixType symmUp = symm.template triangularView<Upper>();
SquareMatrixType symmLo = symm.template triangularView<Lower>(); SquareMatrixType symmLo = symm.template triangularView<Lower>();
LLT<SquareMatrixType,Lower> chollo(symmLo); LLT<SquareMatrixType,Lower> chollo(symmLo);
VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix()); VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
vecX = chollo.solve(vecB); vecX = chollo.solve(vecB);
@ -85,6 +91,14 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
matX = chollo.solve(matB); matX = chollo.solve(matB);
VERIFY_IS_APPROX(symm * matX, 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 // test the upper mode
LLT<SquareMatrixType,Upper> cholup(symmUp); LLT<SquareMatrixType,Upper> cholup(symmUp);
VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix()); VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
@ -93,6 +107,15 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
matX = cholup.solve(matB); matX = cholup.solve(matB);
VERIFY_IS_APPROX(symm * matX, 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; MatrixType neg = -symmLo;
chollo.compute(neg); chollo.compute(neg);
VERIFY(chollo.info()==NumericalIssue); 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(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU())); VERIFY_IS_APPROX(MatrixType(cholup.matrixL().transpose().conjugate()), MatrixType(cholup.matrixU()));
VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL())); VERIFY_IS_APPROX(MatrixType(cholup.matrixU().transpose().conjugate()), MatrixType(cholup.matrixL()));
// test some special use cases of SelfCwiseBinaryOp: // test some special use cases of SelfCwiseBinaryOp:
MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols); MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
m2 = m1; m2 = m1;
@ -167,7 +190,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
// restore // restore
if(sign == -1) if(sign == -1)
symm = -symm; symm = -symm;
// check matrices coming from linear constraints with Lagrange multipliers // check matrices coming from linear constraints with Lagrange multipliers
if(rows>=3) if(rows>=3)
{ {
@ -183,7 +206,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
vecX = ldltlo.solve(vecB); vecX = ldltlo.solve(vecB);
VERIFY_IS_APPROX(A * vecX, vecB); VERIFY_IS_APPROX(A * vecX, vecB);
} }
// check non-full rank matrices // check non-full rank matrices
if(rows>=3) if(rows>=3)
{ {
@ -199,7 +222,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
vecX = ldltlo.solve(vecB); vecX = ldltlo.solve(vecB);
VERIFY_IS_APPROX(A * vecX, vecB); VERIFY_IS_APPROX(A * vecX, vecB);
} }
// check matrices with a wide spectrum // check matrices with a wide spectrum
if(rows>=3) if(rows>=3)
{ {
@ -225,7 +248,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
{ {
RealScalar large_tol = std::sqrt(test_precision<RealScalar>()); RealScalar large_tol = std::sqrt(test_precision<RealScalar>());
VERIFY((A * vecX).isApprox(vecB, large_tol)); VERIFY((A * vecX).isApprox(vecB, large_tol));
++g_test_level; ++g_test_level;
VERIFY_IS_APPROX(A * vecX,vecB); VERIFY_IS_APPROX(A * vecX,vecB);
--g_test_level; --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. // 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 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
{ {
eigen_assert(m.rows() == 2 && m.cols() == 2); eigen_assert(m.rows() == 2 && m.cols() == 2);
MatrixType mat; MatrixType mat;
LDLT<MatrixType> ldlt(2); LDLT<MatrixType> ldlt(2);
{ {
mat << 1, 0, 0, -1; mat << 1, 0, 0, -1;
ldlt.compute(mat); ldlt.compute(mat);
@ -384,11 +407,11 @@ void test_cholesky()
CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) ); CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) ); 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)) ); CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) );
TEST_SET_BUT_UNUSED_VARIABLE(s) TEST_SET_BUT_UNUSED_VARIABLE(s)
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
TEST_SET_BUT_UNUSED_VARIABLE(s) TEST_SET_BUT_UNUSED_VARIABLE(s)
@ -402,6 +425,6 @@ void test_cholesky()
// Test problem size constructors // Test problem size constructors
CALL_SUBTEST_9( LLT<MatrixXf>(10) ); CALL_SUBTEST_9( LLT<MatrixXf>(10) );
CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries)
} }

10
test/lu.cpp
View File

@ -151,10 +151,11 @@ template<typename MatrixType> void lu_invertible()
MatrixType m1_inverse = lu.inverse(); MatrixType m1_inverse = lu.inverse();
VERIFY_IS_APPROX(m2, m1_inverse*m3); 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); 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. const RealScalar rcond_est = lu.rcond();
VERIFY(lu.rcond() > rcond / 10 && lu.rcond() < rcond * 10); VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10);
// test solve with transposed // test solve with transposed
lu.template _solve_impl_transposed<false>(m3, m2); lu.template _solve_impl_transposed<false>(m3, m2);
@ -199,7 +200,8 @@ template<typename MatrixType> void lu_partial_piv()
// Test condition number estimation. // Test condition number estimation.
RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); 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 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 // test solve with transposed
plu.template _solve_impl_transposed<false>(m3, m2); plu.template _solve_impl_transposed<false>(m3, m2);