mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-13 18:37:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix LDLT::solve() if matrix singular but solution exists (bug #241).

Browse Source 
Clarify this in docs and add regression test.
This commit is contained in:
Jitse Niesen 2011-09-11 06:30:53 +01:00
parent 59b83c14fd
commit 6b006772f1
2 changed files with 42 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -158,10 +158,19 @@ template<typename _MatrixType, int _UpLo> class LDLT
}
/** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> .
*
* \note_about_checking_solutions
*
* \sa solveInPlace(), MatrixBase::ldlt()
* More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$
* by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$,
* \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then
* \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the
* least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
* computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular.
*
* \sa MatrixBase::ldlt()
*/
template<typename Rhs>
inline const internal::solve_retval<LDLT, Rhs>
@ -376,7 +385,21 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
dec().matrixL().solveInPlace(dst);
// dst = D^-1 (L^-1 P b)
dst = dec().vectorD().asDiagonal().inverse() * dst;
// more precisely, use pseudo-inverse of D (see bug 241)
using std::abs;
using std::max;
typedef typename LDLTType::MatrixType MatrixType;
typedef typename LDLTType::Scalar Scalar;
typedef typename LDLTType::RealScalar RealScalar;
const Diagonal<const MatrixType> vectorD = dec().vectorD();
RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
for (Index i = 0; i < vectorD.size(); ++i) {
if(abs(vectorD(i)) > tolerance)
dst.row(i) /= vectorD(i);
else
dst.row(i).setZero();
}
// dst = L^-T (D^-1 L^-1 P b)
dec().matrixU().solveInPlace(dst);

17
test/cholesky.cpp
View File

@ -266,6 +266,22 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
}
}
// regression test for bug 241
template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
{
eigen_assert(m.rows() == 2 && m.cols() == 2);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
MatrixType matA;
matA << 1, 1, 1, 1;
VectorType vecB;
vecB << 1, 1;
VectorType vecX = matA.ldlt().solve(vecB);
VERIFY_IS_APPROX(matA * vecX, vecB);
}
template<typename MatrixType> void cholesky_verify_assert()
{
MatrixType tmp;
@ -292,6 +308,7 @@ void test_cholesky()
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
CALL_SUBTEST_3( cholesky(Matrix2d()) );
CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) );
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);