diff options
Diffstat (limited to 'Eigen/src/Cholesky/LDLT.h')
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 40 |
1 files changed, 10 insertions, 30 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)); |