diff options
author | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-04-01 16:19:45 -0700 |
---|---|---|
committer | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-04-01 16:19:45 -0700 |
commit | f54137606eb6d68cbafd10d90013e254b26137ed (patch) | |
tree | bfb7752ec5bd7d6334c7fb139a898e148b203327 /Eigen | |
parent | fb8dccc23e5f717319c230c2701a5fbf1d3c3975 (diff) |
Add condition estimation to Cholesky (LLT) factorization.
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 53 | ||||
-rw-r--r-- | Eigen/src/Core/ConditionEstimator.h | 37 |
2 files changed, 70 insertions, 20 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) { |