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

bug #1150: make IncompleteCholesky more robust by iteratively increase the shift until the factorization succeed (with at most 10 attempts).

Browse Source
This commit is contained in:
Gael Guennebaud 2016-01-23 22:13:54 +01:00
parent cb4e53ff7f
commit 0caa4b1531
2 changed files with 135 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h
View File

@ -240,7 +240,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
else else
m_scale(j) = 1; m_scale(j) = 1;
// FIXME disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster) // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
// Scale and compute the shift for the matrix // Scale and compute the shift for the matrix
RealScalar mindiag = NumTraits<RealScalar>::highest(); RealScalar mindiag = NumTraits<RealScalar>::highest();
@ -252,16 +252,25 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag); mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
} }
FactorType L_save = m_L;
RealScalar shift = 0; RealScalar shift = 0;
if(mindiag <= RealScalar(0.)) if(mindiag <= RealScalar(0.))
shift = m_initialShift - mindiag; shift = m_initialShift - mindiag;
m_info = NumericalIssue;
// Try to perform the incomplete factorization using the current shift
int iter = 0;
do
{
// Apply the shift to the diagonal elements of the matrix // Apply the shift to the diagonal elements of the matrix
for (Index j = 0; j < n; j++) for (Index j = 0; j < n; j++)
vals[colPtr[j]] += shift; vals[colPtr[j]] += shift;
// jki version of the Cholesky factorization // jki version of the Cholesky factorization
for (Index j=0; j < n; ++j) Index j=0;
for (; j < n; ++j)
{ {
// Left-looking factorization of the j-th column // Left-looking factorization of the j-th column
// First, load the j-th column into col_vals // First, load the j-th column into col_vals
@ -305,8 +314,20 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
// Scale the current column // Scale the current column
if(numext::real(diag) <= 0) if(numext::real(diag) <= 0)
{ {
m_info = NumericalIssue; if(++iter>=10)
return; return;
// increase shift
shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
// restore m_L, col_pattern, and listCol
vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
col_pattern.fill(-1);
for(Index i=0; i<n; ++i)
listCol[i].clear();
break;
} }
RealScalar rdiag = sqrt(numext::real(diag)); RealScalar rdiag = sqrt(numext::real(diag));
@ -339,8 +360,13 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
Index jk = colPtr(j)+1; Index jk = colPtr(j)+1;
updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
} }
if(j==n)
{
m_factorizationIsOk = true; m_factorizationIsOk = true;
m_info = Success; m_info = Success;
}
} while(m_info!=Success);
} }
template<typename Scalar, int _UpLo, typename OrderingType> template<typename Scalar, int _UpLo, typename OrderingType>

30
test/incomplete_cholesky.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2015-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
// //
// This Source Code Form is subject to the terms of the Mozilla // This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
@ -34,4 +34,32 @@ void test_incomplete_cholesky()
CALL_SUBTEST_1(( test_incomplete_cholesky_T<double,int>() )); CALL_SUBTEST_1(( test_incomplete_cholesky_T<double,int>() ));
CALL_SUBTEST_2(( test_incomplete_cholesky_T<std::complex<double>, int>() )); CALL_SUBTEST_2(( test_incomplete_cholesky_T<std::complex<double>, int>() ));
CALL_SUBTEST_3(( test_incomplete_cholesky_T<double,long int>() )); CALL_SUBTEST_3(( test_incomplete_cholesky_T<double,long int>() ));
#ifdef EIGEN_TEST_PART_1
// regression for bug 1150
for(int N = 1; N<20; ++N)
{
Eigen::MatrixXd b( N, N );
b.setOnes();
Eigen::SparseMatrix<double> m( N, N );
m.reserve(Eigen::VectorXi::Constant(N,4));
for( int i = 0; i < N; ++i )
{
m.insert( i, i ) = 1;
m.coeffRef( i, i / 2 ) = 2;
m.coeffRef( i, i / 3 ) = 2;
m.coeffRef( i, i / 4 ) = 2;
}
Eigen::SparseMatrix<double> A;
A = m * m.transpose();
Eigen::ConjugateGradient<Eigen::SparseMatrix<double>,
Eigen::Lower | Eigen::Upper,
Eigen::IncompleteCholesky<double> > solver( A );
VERIFY(solver.preconditioner().info() == Eigen::Success);
VERIFY(solver.info() == Eigen::Success);
}
#endif
} }