diff options
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 40 | ||||
| -rw-r--r-- | Eigen/src/Cholesky/LLT.h | 36 | ||||
| -rw-r--r-- | Eigen/src/Core/MathFunctions.h | 8 | ||||
| -rw-r--r-- | test/cholesky.cpp | 50 |
4 files changed, 35 insertions, 99 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index a32a5b8c5..06f80fd87 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -43,6 +43,9 @@ * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root * on D also stabilizes the computation. * + * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky decomposition to determine + * whether a system of equations has a solution. + * * \sa MatrixBase::ldlt(), class LLT */ /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE @@ -88,25 +91,6 @@ template<typename MatrixType> class LDLT /** \returns true if the matrix is negative (semidefinite) */ inline bool isNegative(void) const { return m_sign == -1; } - /** \returns true if the matrix is invertible */ - inline bool isInvertible(void) const { return m_rank == m_matrix.rows(); } - - /** \returns true if the matrix is positive definite */ - inline bool isPositiveDefinite(void) const { return isPositive() && isInvertible(); } - - /** \returns true if the matrix is negative definite */ - inline bool isNegativeDefinite(void) const { return isNegative() && isInvertible(); } - - /** \returns the rank of the matrix of which *this is the LDLT decomposition. - * - * \note This is computed at the time of the construction of the LDLT decomposition. This - * method does not perform any further computation. - */ - inline int rank() const - { - return m_rank; - } - template<typename RhsDerived, typename ResDerived> bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const; @@ -125,7 +109,7 @@ template<typename MatrixType> class LDLT MatrixType m_matrix; IntColVectorType m_p; IntColVectorType m_transpositions; - int m_rank, m_sign; + int m_sign; }; /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix @@ -135,7 +119,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a) { ei_assert(a.rows()==a.cols()); const int size = a.rows(); - m_rank = size; m_matrix = a; @@ -168,8 +151,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a) // 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 208, also by Higham. - cutoff = ei_abs(precision<RealScalar>() * size * biggest_in_corner); + // Algorithms" page 217, also by Higham. + cutoff = ei_abs(machine_epsilon<Scalar>() * size * biggest_in_corner); m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; } @@ -178,7 +161,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a) if(biggest_in_corner < cutoff) { for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - m_rank = j; break; } @@ -200,11 +182,9 @@ void LDLT<MatrixType>::compute(const MatrixType& a) m_matrix.coeffRef(j,j) = Djj; // Finish early if the matrix is not full rank. - if(ei_abs(Djj) < cutoff) // i made experiments, this is better than isMuchSmallerThan(biggest_in_corner), and of course - // much better than plain sign comparison as used to be done before. + if(ei_abs(Djj) < cutoff) { for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - m_rank = j; break; } @@ -230,7 +210,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a) /** 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$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left. @@ -252,6 +232,8 @@ bool LDLT<MatrixType> * * \param bAndX represents both the right-hand side matrix b and result x. * + * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. + * * This version avoids a copy when the right hand side matrix b is not * needed anymore. * @@ -264,8 +246,6 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const const int size = m_matrix.rows(); ei_assert(size == bAndX.rows()); - if (m_rank != size) return false; - // z = P b for(int i = 0; i < size; ++i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index b9d004e97..14ec51a68 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -41,6 +41,10 @@ * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other * situations like generalised eigen problems with hermitian matrices. * + * 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) @@ -70,12 +74,6 @@ template<typename MatrixType> class LLT /** \returns the lower triangular matrix L */ inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; } - /** \returns true if the matrix is positive definite */ - inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; } - - /** \returns true if the matrix is invertible, which in this context is equivalent to positive definite */ - inline bool isInvertible(void) const { return m_isPositiveDefinite; } - template<typename RhsDerived, typename ResDerived> bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const; @@ -90,7 +88,6 @@ template<typename MatrixType> class LLT * The strict upper part is not used and even not initialized. */ MatrixType m_matrix; - bool m_isPositiveDefinite; }; /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix @@ -101,23 +98,24 @@ void LLT<MatrixType>::compute(const MatrixType& a) assert(a.rows()==a.cols()); const int size = a.rows(); m_matrix.resize(size, size); - const RealScalar reference = size * a.diagonal().cwise().abs().maxCoeff(); + // 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<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff(); RealScalar x; x = ei_real(a.coeff(0,0)); - m_isPositiveDefinite = !ei_isMuchSmallerThan(x, reference) && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), reference); m_matrix.coeffRef(0,0) = ei_sqrt(x); if(size==1) 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 (ei_isMuchSmallerThan(x, reference) || (!ei_isMuchSmallerThan(ei_imag(tmp), reference))) - { - m_isPositiveDefinite = false; - return; - } + x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); + if (ei_abs(x) < cutoff) continue; + m_matrix.coeffRef(j,j) = x = ei_sqrt(x); int endSize = size-j-1; @@ -137,7 +135,7 @@ void LLT<MatrixType>::compute(const MatrixType& a) /** 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. @@ -160,6 +158,8 @@ bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDeriv * * \param bAndX represents both the right-hand side matrix b and result x. * + * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. + * * This version avoids a copy when the right hand side matrix b is not * needed anymore. * @@ -171,8 +171,6 @@ bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const { const int size = m_matrix.rows(); ei_assert(size==bAndX.rows()); - if (!m_isPositiveDefinite) - return false; matrixL().solveTriangularInPlace(bAndX); m_matrix.adjoint().template part<UpperTriangular>().solveTriangularInPlace(bAndX); return true; diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 3bc8f4d5b..e201f98b2 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -26,6 +26,7 @@ #define EIGEN_MATHFUNCTIONS_H template<typename T> inline typename NumTraits<T>::Real precision(); +template<typename T> inline typename NumTraits<T>::Real machine_epsilon(); template<typename T> inline T ei_random(T a, T b); template<typename T> inline T ei_random(); template<typename T> inline T ei_random_amplitude() @@ -50,6 +51,7 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y) **************/ template<> inline int precision<int>() { return 0; } +template<> inline int machine_epsilon<int>() { 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; } @@ -102,6 +104,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>()) **************/ template<> inline float precision<float>() { return 1e-5f; } +template<> inline float machine_epsilon<float>() { 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; } @@ -147,6 +150,8 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float **************/ template<> inline double precision<double>() { return 1e-11; } +template<> inline double machine_epsilon<double>() { 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; } @@ -192,6 +197,7 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<do *********************/ template<> inline float precision<std::complex<float> >() { return precision<float>(); } +template<> inline float machine_epsilon<std::complex<float> >() { return machine_epsilon<float>(); } inline float ei_real(const std::complex<float>& x) { return std::real(x); } inline float ei_imag(const std::complex<float>& x) { return std::imag(x); } inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); } @@ -225,6 +231,7 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& **********************/ template<> inline double precision<std::complex<double> >() { return precision<double>(); } +template<> inline double machine_epsilon<std::complex<double> >() { return machine_epsilon<double>(); } inline double ei_real(const std::complex<double>& x) { return std::real(x); } inline double ei_imag(const std::complex<double>& x) { return std::imag(x); } inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); } @@ -259,6 +266,7 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double ******************/ template<> inline long double precision<long double>() { return precision<double>(); } +template<> inline long double machine_epsilon<long double>() { 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; } diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 37a344029..b0e0dd33c 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -86,7 +86,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m) { LLT<SquareMatrixType> chol(symm); - VERIFY(chol.isPositiveDefinite()); VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint()); chol.solve(vecB, &vecX); VERIFY_IS_APPROX(symm * vecX, vecB); @@ -103,18 +102,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m) { LDLT<SquareMatrixType> ldlt(symm); - VERIFY(ldlt.isInvertible()); - if(sign == 1) - { - VERIFY(ldlt.isPositive()); - VERIFY(ldlt.isPositiveDefinite()); - } - if(sign == -1) - { - VERIFY(ldlt.isNegative()); - VERIFY(ldlt.isNegativeDefinite()); - } - // TODO(keir): This doesn't make sense now that LDLT pivots. //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint()); ldlt.solve(vecB, &vecX); @@ -123,15 +110,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m) VERIFY_IS_APPROX(symm * matX, matB); } - // test isPositiveDefinite on non definite matrix - if (rows>4) - { - SquareMatrixType symm = a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint(); - LLT<SquareMatrixType> chol(symm); - VERIFY(!chol.isPositiveDefinite()); - LDLT<SquareMatrixType> cholnosqrt(symm); - VERIFY(!cholnosqrt.isPositiveDefinite()); - } } template<typename Derived> @@ -156,29 +134,6 @@ void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m) } } -template<typename MatrixType> void ldlt_rank() -{ - // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function - int rows = ei_random<int>(50,200); - int rank = ei_random<int>(40, rows-1); - - - // generate a random positive matrix a of given rank - MatrixType m = MatrixType::Random(rows,rows); - QR<MatrixType> qr(m); - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> DiagVectorType; - DiagVectorType d(rows); - d.setZero(); - for(int i = 0; i < rank; i++) d(i)=RealScalar(1); - MatrixType a = qr.matrixQ() * d.asDiagonal() * qr.matrixQ().adjoint(); - - LDLT<MatrixType> ldlt(a); - - VERIFY( ei_abs(ldlt.rank() - rank) <= rank / 20 ); // yes, LDLT::rank is a bit inaccurate... -} - void test_cholesky() { @@ -191,9 +146,4 @@ void test_cholesky() CALL_SUBTEST( cholesky(MatrixXd(17,17)) ); CALL_SUBTEST( cholesky(MatrixXf(200,200)) ); } - for(int i = 0; i < g_repeat/3; i++) { - CALL_SUBTEST( ldlt_rank<MatrixXd>() ); - CALL_SUBTEST( ldlt_rank<MatrixXf>() ); - CALL_SUBTEST( ldlt_rank<MatrixXcd>() ); - } } |