GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-21 09:09:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Many improvements in LLT and LDLT:

Browse Source 
* in LDLT, support the negative semidefinite case
* fix bad floating-point comparisons, improves greatly the accuracy of methods like
  isPositiveDefinite() and rank()
* simplifications
* identify (but not resolve) bug: claim that only triangular part is used, is inaccurate
* expanded unit-tests
This commit is contained in:
Benoit Jacob 2009-03-30 21:45:45 +00:00
parent df9dfa1455
commit a1ba995f05
3 changed files with 142 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
Eigen/src/Cholesky/LDLT.h
View File

@ -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<typename MatrixType> class LDLT
{
@ -81,8 +81,20 @@ template<typename MatrixType> class LDLT
/** \returns the coefficients of the diagonal matrix D */
inline DiagonalCoeffs<MatrixType> 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<typename MatrixType> 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<typename MatrixType>
void LDLT<MatrixType>::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<MatrixType>::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<MatrixType>::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<RealScalar>() * 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<RealScalar>() * 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<MatrixType>::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<MatrixType>
::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *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);
}

15
Eigen/src/Cholesky/LLT.h
View File

@ -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<typename MatrixType> class LLT
{
@ -72,6 +73,9 @@ template<typename MatrixType> 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<typename RhsDerived, typename ResDerived>
bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
@ -97,11 +101,10 @@ void LLT<MatrixType>::compute(const MatrixType& a)
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
const RealScalar eps = precision<Scalar>();
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<MatrixType>::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;

101
test/cholesky.cpp
View File

@ -25,7 +25,7 @@
#define EIGEN_NO_ASSERTION_CHECKING
#include "main.h"
#include <Eigen/Cholesky>
#include <Eigen/LU>
#include <Eigen/QR>
#ifdef HAS_GSL
#include "gsl_helper.h"
@ -52,6 +52,10 @@ template<typename MatrixType> 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<StrictlyLowerTriangular>().setZero();
#ifdef HAS_GSL
if (ei_is_same_type<RealScalar,double>::ret)
{
@ -80,17 +84,6 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
}
#endif
{
LDLT<SquareMatrixType> 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<SquareMatrixType> chol(symm);
VERIFY(chol.isPositiveDefinite());
@ -101,6 +94,35 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(symm * matX, matB);
}
int sign = ei_random<int>()%2 ? 1 : -1;
if(sign == -1)
{
symm = -symm; // test a negative matrix
}
{
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);
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<typename MatrixType> void cholesky(const MatrixType& m)
}
}
template<typename Derived>
void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
{
typedef typename Derived::RealScalar RealScalar;
for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
{
RealScalar d = Eigen::ei_random<RealScalar>(-1,1);
int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
int j;
do {
j = Eigen::ei_random<int>(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<int>(0,m.cols()-1); // i is a random column number
do {
j = Eigen::ei_random<int>(0,m.cols()-1);
} while (i==j); // j is another one (must be different)
m.col(i) += d * m.col(j);
}
}
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()
{
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<MatrixXd>() );
CALL_SUBTEST( ldlt_rank<MatrixXf>() );
CALL_SUBTEST( ldlt_rank<MatrixXcd>() );
}
}