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

Fix bug #736: LDLT isPositive returns false for a positive semidefinite matrix

Browse Source 
Add unit test covering this case.
This commit is contained in:
Jitse Niesen 2014-02-06 11:06:06 +00:00
parent 6c527bd811
commit ff8d81762d
2 changed files with 52 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -17,6 +17,9 @@ namespace Eigen {
namespace internal { namespace internal {
template<typename MatrixType, int UpLo> struct LDLT_Traits; template<typename MatrixType, int UpLo> struct LDLT_Traits;
// PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
} }
/** \ingroup Cholesky_Module /** \ingroup Cholesky_Module
@ -69,7 +72,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
* perform decompositions via LDLT::compute(const MatrixType&). * perform decompositions via LDLT::compute(const MatrixType&).
*/ */
LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {} LDLT()
: m_matrix(),
m_transpositions(),
m_sign(internal::ZeroSign),
m_isInitialized(false)
{}
/** \brief Default Constructor with memory preallocation /** \brief Default Constructor with memory preallocation
* *
@ -81,6 +89,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
: m_matrix(size, size), : m_matrix(size, size),
m_transpositions(size), m_transpositions(size),
m_temporary(size), m_temporary(size),
m_sign(internal::ZeroSign),
m_isInitialized(false) m_isInitialized(false)
{} {}
@ -93,6 +102,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
: m_matrix(matrix.rows(), matrix.cols()), : m_matrix(matrix.rows(), matrix.cols()),
m_transpositions(matrix.rows()), m_transpositions(matrix.rows()),
m_temporary(matrix.rows()), m_temporary(matrix.rows()),
m_sign(internal::ZeroSign),
m_isInitialized(false) m_isInitialized(false)
{ {
compute(matrix); compute(matrix);
@ -139,7 +149,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
inline bool isPositive() const inline bool isPositive() const
{ {
eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == 1; return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
} }
#ifdef EIGEN2_SUPPORT #ifdef EIGEN2_SUPPORT
@ -153,7 +163,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
inline bool isNegative(void) const inline bool isNegative(void) const
{ {
eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == -1; return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
} }
/** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
@ -235,7 +245,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
MatrixType m_matrix; MatrixType m_matrix;
TranspositionType m_transpositions; TranspositionType m_transpositions;
TmpMatrixType m_temporary; TmpMatrixType m_temporary;
int m_sign; internal::SignMatrix m_sign;
bool m_isInitialized; bool m_isInitialized;
}; };
@ -246,7 +256,7 @@ template<int UpLo> struct ldlt_inplace;
template<> struct ldlt_inplace<Lower> template<> struct ldlt_inplace<Lower>
{ {
template<typename MatrixType, typename TranspositionType, typename Workspace> template<typename MatrixType, typename TranspositionType, typename Workspace>
static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
{ {
using std::abs; using std::abs;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
@ -258,8 +268,9 @@ template<> struct ldlt_inplace<Lower>
if (size <= 1) if (size <= 1)
{ {
transpositions.setIdentity(); transpositions.setIdentity();
if(sign) if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef;
*sign = numext::real(mat.coeff(0,0))>0 ? 1:-1; else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef;
else sign = ZeroSign;
return true; return true;
} }
@ -284,7 +295,6 @@ template<> struct ldlt_inplace<Lower>
if(biggest_in_corner < cutoff) if(biggest_in_corner < cutoff)
{ {
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
if(sign) *sign = 0;
break; break;
} }
@ -326,14 +336,14 @@ template<> struct ldlt_inplace<Lower>
if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k); A21 /= mat.coeffRef(k,k);
if(sign) RealScalar realAkk = numext::real(mat.coeffRef(k,k));
{ if (sign == PositiveSemiDef) {
// LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right if (realAkk < 0) sign = Indefinite;
int newSign = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; } else if (sign == NegativeSemiDef) {
if(k == 0) if (realAkk > 0) sign = Indefinite;
*sign = newSign; } else if (sign == ZeroSign) {
else if(*sign != newSign) if (realAkk > 0) sign = PositiveSemiDef;
*sign = 0; else if (realAkk < 0) sign = NegativeSemiDef;
} }
} }
@ -399,7 +409,7 @@ template<> struct ldlt_inplace<Lower>
template<> struct ldlt_inplace<Upper> template<> struct ldlt_inplace<Upper>
{ {
template<typename MatrixType, typename TranspositionType, typename Workspace> template<typename MatrixType, typename TranspositionType, typename Workspace>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
{ {
Transpose<MatrixType> matt(mat); Transpose<MatrixType> matt(mat);
return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
@ -445,7 +455,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
m_isInitialized = false; m_isInitialized = false;
m_temporary.resize(size); m_temporary.resize(size);
internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign); internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;
@ -473,7 +483,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri
for (Index i = 0; i < size; i++) for (Index i = 0; i < size; i++)
m_transpositions.coeffRef(i) = i; m_transpositions.coeffRef(i) = i;
m_temporary.resize(size); m_temporary.resize(size);
m_sign = sigma>=0 ? 1 : -1; m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
m_isInitialized = true; m_isInitialized = true;
} }

24
test/cholesky.cpp
View File

@ -263,8 +263,8 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
// This test checks that LDLT reports correctly that matrix is indefinite. // This test checks that LDLT reports correctly that matrix is indefinite.
// See http://forum.kde.org/viewtopic.php?f=74&t=106942 // See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
template<typename MatrixType> void cholesky_indefinite(const MatrixType& m) template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
{ {
eigen_assert(m.rows() == 2 && m.cols() == 2); eigen_assert(m.rows() == 2 && m.cols() == 2);
MatrixType mat; MatrixType mat;
@ -280,6 +280,24 @@ template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isNegative());
VERIFY(!ldlt.isPositive()); VERIFY(!ldlt.isPositive());
} }
{
mat << 0, 0, 0, 0;
LDLT<MatrixType> ldlt(mat);
VERIFY(ldlt.isNegative());
VERIFY(ldlt.isPositive());
}
{
mat << 0, 0, 0, 1;
LDLT<MatrixType> ldlt(mat);
VERIFY(!ldlt.isNegative());
VERIFY(ldlt.isPositive());
}
{
mat << -1, 0, 0, 0;
LDLT<MatrixType> ldlt(mat);
VERIFY(ldlt.isNegative());
VERIFY(!ldlt.isPositive());
}
} }
template<typename MatrixType> void cholesky_verify_assert() template<typename MatrixType> void cholesky_verify_assert()
@ -309,7 +327,7 @@ void test_cholesky()
CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) ); CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
CALL_SUBTEST_3( cholesky(Matrix2d()) ); CALL_SUBTEST_3( cholesky(Matrix2d()) );
CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) ); CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
CALL_SUBTEST_3( cholesky_indefinite(Matrix2d()) ); CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) ); CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) ); CALL_SUBTEST_5( cholesky(Matrix4d()) );
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);