diff options
author | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2011-09-11 06:30:53 +0100 |
---|---|---|
committer | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2011-09-11 06:30:53 +0100 |
commit | 6b006772f1c8083df609d70218f5848bce7dfc28 (patch) | |
tree | 0e2db2d2d77187ee7567cd2664d8611f9f27f927 | |
parent | 59b83c14fd2bec0b8c8afa7a2fa0357af7f0f827 (diff) |
Fix LDLT::solve() if matrix singular but solution exists (bug #241).
Clarify this in docs and add regression test.
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 27 | ||||
-rw-r--r-- | test/cholesky.cpp | 17 |
2 files changed, 42 insertions, 2 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 5e2352caa..f47b2ea56 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -159,9 +159,18 @@ 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); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index b2139563b..35f5a3b68 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -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); |