mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #785: Make Cholesky decomposition work for empty matrices

Browse Source
This commit is contained in:
Christoph Hertzberg 2018-12-03 16:18:15 +01:00
parent 0ea7ae7213
commit 919414b9fe
3 changed files with 14 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -304,7 +304,8 @@ template<> struct ldlt_inplace<Lower>
if (size <= 1) if (size <= 1)
{ {
transpositions.setIdentity(); transpositions.setIdentity();
if (numext::real(mat.coeff(0,0)) > static_cast<RealScalar>(0) ) sign = PositiveSemiDef; if(size==0) sign = ZeroSign;
else if (numext::real(mat.coeff(0,0)) > static_cast<RealScalar>(0) ) sign = PositiveSemiDef;
else if (numext::real(mat.coeff(0,0)) < static_cast<RealScalar>(0)) sign = NegativeSemiDef; else if (numext::real(mat.coeff(0,0)) < static_cast<RealScalar>(0)) sign = NegativeSemiDef;
else sign = ZeroSign; else sign = ZeroSign;
return true; return true;

2
Eigen/src/Core/ConditionEstimator.h
View File

@ -160,7 +160,7 @@ rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Deco
{ {
typedef typename Decomposition::RealScalar RealScalar; typedef typename Decomposition::RealScalar RealScalar;
eigen_assert(dec.rows() == dec.cols()); eigen_assert(dec.rows() == dec.cols());
if (dec.rows() == 0) return RealScalar(1); if (dec.rows() == 0) return RealScalar(1)/RealScalar(0);
if (matrix_norm == RealScalar(0)) return RealScalar(0); if (matrix_norm == RealScalar(0)) return RealScalar(0);
if (dec.rows() == 1) return RealScalar(1); if (dec.rows() == 1) return RealScalar(1);
const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec); const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);

16
test/cholesky.cpp
View File

@ -19,6 +19,7 @@
template<typename MatrixType, int UpLo> template<typename MatrixType, int UpLo>
typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) {
if(m.cols()==0) return typename MatrixType::RealScalar(0);
MatrixType symm = m.template selfadjointView<UpLo>(); MatrixType symm = m.template selfadjointView<UpLo>();
return symm.cwiseAbs().colwise().sum().maxCoeff(); return symm.cwiseAbs().colwise().sum().maxCoeff();
} }
@ -96,7 +97,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
RealScalar rcond_est = chollo.rcond(); RealScalar rcond_est = chollo.rcond();
// Verify that the estimated condition number is within a factor of 10 of the // Verify that the estimated condition number is within a factor of 10 of the
// truth. // truth.
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
// test the upper mode // test the upper mode
LLT<SquareMatrixType,Upper> cholup(symmUp); LLT<SquareMatrixType,Upper> cholup(symmUp);
@ -112,12 +113,12 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) / rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
matrix_l1_norm<MatrixType, Upper>(symmUp_inverse); matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
rcond_est = cholup.rcond(); rcond_est = cholup.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
MatrixType neg = -symmLo; MatrixType neg = -symmLo;
chollo.compute(neg); chollo.compute(neg);
VERIFY(chollo.info()==NumericalIssue); VERIFY(neg.size()==0 || chollo.info()==NumericalIssue);
VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU())); VERIFY_IS_APPROX(MatrixType(chollo.matrixL().transpose().conjugate()), MatrixType(chollo.matrixU()));
VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL())); VERIFY_IS_APPROX(MatrixType(chollo.matrixU().transpose().conjugate()), MatrixType(chollo.matrixL()));
@ -166,7 +167,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
RealScalar rcond_est = ldltlo.rcond(); RealScalar rcond_est = ldltlo.rcond();
// Verify that the estimated condition number is within a factor of 10 of the // Verify that the estimated condition number is within a factor of 10 of the
// truth. // truth.
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
LDLT<SquareMatrixType,Upper> ldltup(symmUp); LDLT<SquareMatrixType,Upper> ldltup(symmUp);
@ -183,7 +184,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) / rcond = (RealScalar(1) / matrix_l1_norm<MatrixType, Upper>(symmUp)) /
matrix_l1_norm<MatrixType, Upper>(symmUp_inverse); matrix_l1_norm<MatrixType, Upper>(symmUp_inverse);
rcond_est = ldltup.rcond(); rcond_est = ldltup.rcond();
VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); VERIFY(rcond_est >= rcond / 10 && rcond_est <= rcond * 10);
VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU())); VERIFY_IS_APPROX(MatrixType(ldltlo.matrixL().transpose().conjugate()), MatrixType(ldltlo.matrixU()));
VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL())); VERIFY_IS_APPROX(MatrixType(ldltlo.matrixU().transpose().conjugate()), MatrixType(ldltlo.matrixL()));
@ -507,6 +508,11 @@ EIGEN_DECLARE_TEST(cholesky)
CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) ); CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
TEST_SET_BUT_UNUSED_VARIABLE(s) TEST_SET_BUT_UNUSED_VARIABLE(s)
} }
// empty matrix, regression test for Bug 785:
CALL_SUBTEST_2( cholesky(MatrixXd(0,0)) );
// This does not work yet:
// CALL_SUBTEST_2( cholesky(Matrix<double,0,0>()) );
CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() ); CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() ); CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );