diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index c4e7fdffd..3d8789d6e 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -1,5 +1,5 @@ // This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. +// for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // @@ -41,11 +41,16 @@ * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other * situations like generalised eigen problems with hermitian matrices. * - * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, - * the strict lower part does not have to store correct values. + * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices, + * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations + * has a solution. * * \sa MatrixBase::llt(), class LDLT */ + /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) + * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, + * the strict lower part does not have to store correct values. + */ template class LLT { private: @@ -60,17 +65,30 @@ template class LLT public: + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via LLT::compute(const MatrixType&). + */ + LLT() : m_matrix(), m_isInitialized(false) {} + LLT(const MatrixType& matrix) - : m_matrix(matrix.rows(), matrix.cols()) + : m_matrix(matrix.rows(), matrix.cols()), + m_isInitialized(false) { compute(matrix); } /** \returns the lower triangular matrix L */ - inline Part matrixL(void) const { return m_matrix; } - - /** \returns true if the matrix is positive definite */ - inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } + inline Part matrixL(void) const + { + ei_assert(m_isInitialized && "LLT is not initialized."); + return m_matrix; + } + + /** \deprecated */ + inline bool isPositiveDefinite(void) const { return m_isInitialized && m_isPositiveDefinite; } template bool solve(const MatrixBase &b, MatrixBase *result) const; @@ -86,6 +104,7 @@ template class LLT * The strict upper part is not used and even not initialized. */ MatrixType m_matrix; + bool m_isInitialized; bool m_isPositiveDefinite; }; @@ -95,24 +114,34 @@ template void LLT::compute(const MatrixType& a) { assert(a.rows()==a.cols()); + m_isPositiveDefinite = true; const int size = a.rows(); m_matrix.resize(size, size); - const RealScalar eps = ei_sqrt(precision()); - + // The biggest overall is the point of reference to which further diagonals + // are compared; if any diagonal is negligible compared + // to the largest overall, the algorithm bails. This cutoff is suggested + // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by + // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical + // Algorithms" page 217, also by Higham. + const RealScalar cutoff = machine_epsilon() * size * a.diagonal().cwise().abs().maxCoeff(); RealScalar x; x = ei_real(a.coeff(0,0)); - m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1)); m_matrix.coeffRef(0,0) = ei_sqrt(x); + if(size==1) + { + m_isInitialized = true; + return; + } m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); for (int j = 1; j < size; ++j) { - Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); - x = ei_real(tmp); - if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1)))) + x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); + if (x < cutoff) { m_isPositiveDefinite = false; - return; + continue; } + m_matrix.coeffRef(j,j) = x = ei_sqrt(x); int endSize = size-j-1; @@ -127,12 +156,14 @@ void LLT::compute(const MatrixType& a) - m_matrix.col(j).end(endSize) ) / x; } } + + m_isInitialized = true; } /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A. * The result is stored in \a result * - * \returns true in case of success, false otherwise. + * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. * * In other words, it computes \f$ b = A^{-1} b \f$ with * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left. @@ -146,6 +177,7 @@ template template bool LLT::solve(const MatrixBase &b, MatrixBase *result) const { + ei_assert(m_isInitialized && "LLT is not initialized."); const int size = m_matrix.rows(); ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b"); return solveInPlace((*result) = b); @@ -155,6 +187,8 @@ bool LLT::solve(const MatrixBase &b, MatrixBase template bool LLT::solveInPlace(MatrixBase &bAndX) const { + ei_assert(m_isInitialized && "LLT is not initialized."); const int size = m_matrix.rows(); ei_assert(size==bAndX.rows()); - if (!m_isPositiveDefinite) - return false; matrixL().solveTriangularInPlace(bAndX); m_matrix.adjoint().template part().solveTriangularInPlace(bAndX); return true; diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 8f56ef40e..77b9a02cf 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -26,6 +26,7 @@ #define EIGEN_MATHFUNCTIONS_H template inline typename NumTraits::Real precision(); +template inline typename NumTraits::Real machine_epsilon(); template inline T ei_random(T a, T b); template inline T ei_random(); template inline T ei_random_amplitude() @@ -49,6 +50,7 @@ template inline T ei_hypot(T x, T y) **************/ template<> inline int precision() { return 0; } +template<> inline int machine_epsilon() { return 0; } inline int ei_real(int x) { return x; } inline int ei_imag(int) { return 0; } inline int ei_conj(int x) { return x; } @@ -93,6 +95,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision()) **************/ template<> inline float precision() { return 1e-5f; } +template<> inline float machine_epsilon() { return 1.192e-07f; } inline float ei_real(float x) { return x; } inline float ei_imag(float) { return 0.f; } inline float ei_conj(float x) { return x; } @@ -138,6 +141,8 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision inline double precision() { return 1e-11; } +template<> inline double machine_epsilon() { return 2.220e-16; } + inline double ei_real(double x) { return x; } inline double ei_imag(double) { return 0.; } inline double ei_conj(double x) { return x; } @@ -183,6 +188,7 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision inline float precision >() { return precision(); } +template<> inline float machine_epsilon >() { return machine_epsilon(); } inline float ei_real(const std::complex& x) { return std::real(x); } inline float ei_imag(const std::complex& x) { return std::imag(x); } inline std::complex ei_conj(const std::complex& x) { return std::conj(x); } @@ -216,6 +222,7 @@ inline bool ei_isApprox(const std::complex& a, const std::complex& **********************/ template<> inline double precision >() { return precision(); } +template<> inline double machine_epsilon >() { return machine_epsilon(); } inline double ei_real(const std::complex& x) { return std::real(x); } inline double ei_imag(const std::complex& x) { return std::imag(x); } inline std::complex ei_conj(const std::complex& x) { return std::conj(x); } @@ -250,6 +257,7 @@ inline bool ei_isApprox(const std::complex& a, const std::complex inline long double precision() { return precision(); } +template<> inline long double machine_epsilon() { return 1.084e-19l; } inline long double ei_real(long double x) { return x; } inline long double ei_imag(long double) { return 0.; } inline long double ei_conj(long double x) { return x; }