diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index c9f49e896..7edac3859 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -34,19 +34,19 @@ * * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition * - * Perform a robust Cholesky decomposition of a symmetric positive semidefinite + * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite * matrix \f$ A \f$ such that \f$ A = P^TLDL^*P \f$, where P is a permutation matrix, L * is lower triangular with a unit diagonal and D is a diagonal matrix. * - * The decomposition uses pivoting to ensure stability, such that if A is - * positive semidefinite (i.e. eigenvalues are non-negative), then L will have + * The decomposition uses pivoting to ensure stability, so that L will have * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root * on D also stabilizes the computation. * + * \sa MatrixBase::ldlt(), class LLT + */ + /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE * 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. - * - * \sa MatrixBase::ldlt(), class LLT */ template class LDLT { @@ -81,8 +81,20 @@ template class LDLT /** \returns the coefficients of the diagonal matrix D */ inline DiagonalCoeffs vectorD(void) const { return m_matrix.diagonal(); } + /** \returns true if the matrix is positive (semidefinite) */ + inline bool isPositive(void) const { return m_sign == 1; } + + /** \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 m_rank == m_matrix.rows(); } + 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. * @@ -112,15 +124,15 @@ template class LDLT MatrixType m_matrix; IntColVectorType m_p; IntColVectorType m_transpositions; - int m_rank; + int m_rank, m_sign; }; -/** Compute / recompute the LLT decomposition A = L D L^* = U^* D U of \a matrix +/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix */ template void LDLT::compute(const MatrixType& a) { - assert(a.rows()==a.cols()); + ei_assert(a.rows()==a.cols()); const int size = a.rows(); m_rank = size; @@ -129,6 +141,7 @@ void LDLT::compute(const MatrixType& a) if (size <= 1) { m_p.setZero(); m_transpositions.setZero(); + m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1; return; } @@ -141,35 +154,38 @@ void LDLT::compute(const MatrixType& a) for (int j = 0; j < size; ++j) { - // Find largest element diagonal and pivot it upward for processing next. - int row_of_biggest_in_corner, col_of_biggest_in_corner; + // Find largest diagonal element + int index_of_biggest_in_corner; biggest_in_corner = m_matrix.diagonal().end(size-j).cwise().abs() - .maxCoeff(&row_of_biggest_in_corner, - &col_of_biggest_in_corner); + .maxCoeff(&index_of_biggest_in_corner); + index_of_biggest_in_corner += j; - // The biggest overall is the point of reference to which further diagonals - // are compared; if any diagonal is negligible to machine epsilon 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 208, also by Higham. - if(j == 0) cutoff = ei_abs(precision() * size * biggest_in_corner); + if(j == 0) + { + // 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 208, also by Higham. + cutoff = ei_abs(precision() * size * biggest_in_corner); + + m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; + } // Finish early if the matrix is not full rank. if(biggest_in_corner < cutoff) { for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - m_matrix.block(j, j, size-j, size-j).fill(0); // Zero unreliable data. m_rank = j; break; } - row_of_biggest_in_corner += j; - m_transpositions.coeffRef(j) = row_of_biggest_in_corner; - if(j != row_of_biggest_in_corner) + m_transpositions.coeffRef(j) = index_of_biggest_in_corner; + if(j != index_of_biggest_in_corner) { - m_matrix.row(j).swap(m_matrix.row(row_of_biggest_in_corner)); - m_matrix.col(j).swap(m_matrix.col(row_of_biggest_in_corner)); + m_matrix.row(j).swap(m_matrix.row(index_of_biggest_in_corner)); + m_matrix.col(j).swap(m_matrix.col(index_of_biggest_in_corner)); } if (j == 0) { @@ -182,13 +198,11 @@ void LDLT::compute(const MatrixType& a) * m_matrix.col(j).start(j).conjugate()).coeff(0,0)); m_matrix.coeffRef(j,j) = Djj; - // Finish early if the matrix is not full rank or is indefinite. This - // check is deliberately not against eps, so that the decomposition works - // regardless of overall matrix scale. - if(Djj <= 0) + // 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. { for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - m_matrix.block(j, j, size-j, size-j).fill(0); // Zero unreliable data. m_rank = j; break; } @@ -228,7 +242,7 @@ bool LDLT ::solve(const MatrixBase &b, MatrixBase *result) const { const int size = m_matrix.rows(); - ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b"); + ei_assert(size==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); *result = b; return solveInPlace(*result); } diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 66b86a204..b9d004e97 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -41,10 +41,11 @@ * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other * situations like generalised eigen problems with hermitian matrices. * + * \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. - * - * \sa MatrixBase::llt(), class LDLT */ template class LLT { @@ -72,6 +73,9 @@ template class LLT /** \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 bool solve(const MatrixBase &b, MatrixBase *result) const; @@ -97,11 +101,10 @@ void LLT::compute(const MatrixType& a) assert(a.rows()==a.cols()); const int size = a.rows(); m_matrix.resize(size, size); - const RealScalar eps = precision(); - + const RealScalar reference = 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_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; @@ -110,7 +113,7 @@ void LLT::compute(const MatrixType& a) { 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)))) + if (ei_isMuchSmallerThan(x, reference) || (!ei_isMuchSmallerThan(ei_imag(tmp), reference))) { m_isPositiveDefinite = false; return; diff --git a/test/cholesky.cpp b/test/cholesky.cpp index b3e0df438..37a344029 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -25,7 +25,7 @@ #define EIGEN_NO_ASSERTION_CHECKING #include "main.h" #include -#include +#include #ifdef HAS_GSL #include "gsl_helper.h" @@ -52,6 +52,10 @@ template void cholesky(const MatrixType& m) MatrixType a1 = MatrixType::Random(rows,cols); symm += a1 * a1.adjoint(); + // to test if really Cholesky only uses the upper triangular part, uncomment the following + // FIXME: currently that fails !! + //symm.template part().setZero(); + #ifdef HAS_GSL if (ei_is_same_type::ret) { @@ -80,17 +84,6 @@ template void cholesky(const MatrixType& m) } #endif - { - LDLT ldlt(symm); - VERIFY(ldlt.isPositiveDefinite()); - // 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); - VERIFY_IS_APPROX(symm * vecX, vecB); - ldlt.solve(matB, &matX); - VERIFY_IS_APPROX(symm * matX, matB); - } - { LLT chol(symm); VERIFY(chol.isPositiveDefinite()); @@ -101,6 +94,35 @@ template void cholesky(const MatrixType& m) VERIFY_IS_APPROX(symm * matX, matB); } + int sign = ei_random()%2 ? 1 : -1; + + if(sign == -1) + { + symm = -symm; // test a negative matrix + } + + { + LDLT 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); + VERIFY_IS_APPROX(symm * vecX, vecB); + ldlt.solve(matB, &matX); + VERIFY_IS_APPROX(symm * matX, matB); + } + // test isPositiveDefinite on non definite matrix if (rows>4) { @@ -112,6 +134,52 @@ template void cholesky(const MatrixType& m) } } +template +void doSomeRankPreservingOperations(Eigen::MatrixBase& m) +{ + typedef typename Derived::RealScalar RealScalar; + for(int a = 0; a < 3*(m.rows()+m.cols()); a++) + { + RealScalar d = Eigen::ei_random(-1,1); + int i = Eigen::ei_random(0,m.rows()-1); // i is a random row number + int j; + do { + j = Eigen::ei_random(0,m.rows()-1); + } while (i==j); // j is another one (must be different) + m.row(i) += d * m.row(j); + + i = Eigen::ei_random(0,m.cols()-1); // i is a random column number + do { + j = Eigen::ei_random(0,m.cols()-1); + } while (i==j); // j is another one (must be different) + m.col(i) += d * m.col(j); + } +} + +template 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(50,200); + int rank = ei_random(40, rows-1); + + + // generate a random positive matrix a of given rank + MatrixType m = MatrixType::Random(rows,rows); + QR qr(m); + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix 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 ldlt(a); + + VERIFY( ei_abs(ldlt.rank() - rank) <= rank / 20 ); // yes, LDLT::rank is a bit inaccurate... +} + + void test_cholesky() { for(int i = 0; i < g_repeat; i++) { @@ -120,7 +188,12 @@ void test_cholesky() CALL_SUBTEST( cholesky(Matrix3f()) ); CALL_SUBTEST( cholesky(Matrix4d()) ); CALL_SUBTEST( cholesky(MatrixXcd(7,7)) ); - CALL_SUBTEST( cholesky(MatrixXf(17,17)) ); - CALL_SUBTEST( cholesky(MatrixXd(33,33)) ); + 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() ); + CALL_SUBTEST( ldlt_rank() ); + CALL_SUBTEST( ldlt_rank() ); } }