From 0caa4b1531def27bde0ffeb942cc10f9917a47c6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 23 Jan 2016 22:13:54 +0100 Subject: [PATCH] bug #1150: make IncompleteCholesky more robust by iteratively increase the shift until the factorization succeed (with at most 10 attempts). --- .../IncompleteCholesky.h | 186 ++++++++++-------- test/incomplete_cholesky.cpp | 30 ++- 2 files changed, 135 insertions(+), 81 deletions(-) diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h index 18e9d1466..babc14d9e 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h @@ -240,7 +240,7 @@ void IncompleteCholesky::factorize(const _MatrixType else 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 RealScalar mindiag = NumTraits::highest(); @@ -251,96 +251,122 @@ void IncompleteCholesky::factorize(const _MatrixType eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored"); mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag); } + + FactorType L_save = m_L; RealScalar shift = 0; if(mindiag <= RealScalar(0.)) shift = m_initialShift - mindiag; - // Apply the shift to the diagonal elements of the matrix - for (Index j = 0; j < n; j++) - vals[colPtr[j]] += shift; - - // jki version of the Cholesky factorization - for (Index j=0; j < n; ++j) - { - // Left-looking factorization of the j-th column - // First, load the j-th column into col_vals - Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored - col_nnz = 0; - for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++) + 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 + for (Index j = 0; j < n; j++) + vals[colPtr[j]] += shift; + + // jki version of the Cholesky factorization + Index j=0; + for (; j < n; ++j) { - StorageIndex l = rowIdx[i]; - col_vals(col_nnz) = vals[i]; - col_irow(col_nnz) = l; - col_pattern(l) = col_nnz; - col_nnz++; - } - { - typename std::list::iterator k; - // Browse all previous columns that will update column j - for(k = listCol[j].begin(); k != listCol[j].end(); k++) + // Left-looking factorization of the j-th column + // First, load the j-th column into col_vals + Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored + col_nnz = 0; + for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++) { - Index jk = firstElt(*k); // First element to use in the column - eigen_internal_assert(rowIdx[jk]==j); - Scalar v_j_jk = numext::conj(vals[jk]); - - jk += 1; - for (Index i = jk; i < colPtr[*k+1]; i++) - { - StorageIndex l = rowIdx[i]; - if(col_pattern[l]<0) - { - col_vals(col_nnz) = vals[i] * v_j_jk; - col_irow[col_nnz] = l; - col_pattern(l) = col_nnz; - col_nnz++; - } - else - col_vals(col_pattern[l]) -= vals[i] * v_j_jk; - } - updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); + StorageIndex l = rowIdx[i]; + col_vals(col_nnz) = vals[i]; + col_irow(col_nnz) = l; + col_pattern(l) = col_nnz; + col_nnz++; } + { + typename std::list::iterator k; + // Browse all previous columns that will update column j + for(k = listCol[j].begin(); k != listCol[j].end(); k++) + { + Index jk = firstElt(*k); // First element to use in the column + eigen_internal_assert(rowIdx[jk]==j); + Scalar v_j_jk = numext::conj(vals[jk]); + + jk += 1; + for (Index i = jk; i < colPtr[*k+1]; i++) + { + StorageIndex l = rowIdx[i]; + if(col_pattern[l]<0) + { + col_vals(col_nnz) = vals[i] * v_j_jk; + col_irow[col_nnz] = l; + col_pattern(l) = col_nnz; + col_nnz++; + } + else + col_vals(col_pattern[l]) -= vals[i] * v_j_jk; + } + updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); + } + } + + // Scale the current column + if(numext::real(diag) <= 0) + { + if(++iter>=10) + return; + + // increase shift + shift = numext::maxi(m_initialShift,RealScalar(2)*shift); + // restore m_L, col_pattern, and listCol + vals = Map(L_save.valuePtr(), nnz); + rowIdx = Map(L_save.innerIndexPtr(), nnz); + colPtr = Map(L_save.outerIndexPtr(), n+1); + col_pattern.fill(-1); + for(Index i=0; i cvals = col_vals.head(col_nnz); + Ref cirow = col_irow.head(col_nnz); + internal::QuickSplit(cvals,cirow, p); + // Insert the largest p elements in the matrix + Index cpt = 0; + for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++) + { + vals[i] = col_vals(cpt); + rowIdx[i] = col_irow(cpt); + // restore col_pattern: + col_pattern(col_irow(cpt)) = -1; + cpt++; + } + // Get the first smallest row index and put it after the diagonal element + Index jk = colPtr(j)+1; + updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); } - - // Scale the current column - if(numext::real(diag) <= 0) + + if(j==n) { - m_info = NumericalIssue; - return; + m_factorizationIsOk = true; + m_info = Success; } - - RealScalar rdiag = sqrt(numext::real(diag)); - vals[colPtr[j]] = rdiag; - for (Index k = 0; k cvals = col_vals.head(col_nnz); - Ref cirow = col_irow.head(col_nnz); - internal::QuickSplit(cvals,cirow, p); - // Insert the largest p elements in the matrix - Index cpt = 0; - for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++) - { - vals[i] = col_vals(cpt); - rowIdx[i] = col_irow(cpt); - // restore col_pattern: - col_pattern(col_irow(cpt)) = -1; - cpt++; - } - // Get the first smallest row index and put it after the diagonal element - Index jk = colPtr(j)+1; - updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); - } - m_factorizationIsOk = true; - m_info = Success; + } while(m_info!=Success); } template diff --git a/test/incomplete_cholesky.cpp b/test/incomplete_cholesky.cpp index 7acad9872..59ffe9259 100644 --- a/test/incomplete_cholesky.cpp +++ b/test/incomplete_cholesky.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2015 Gael Guennebaud +// Copyright (C) 2015-2016 Gael Guennebaud // // 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 @@ -34,4 +34,32 @@ void test_incomplete_cholesky() CALL_SUBTEST_1(( test_incomplete_cholesky_T() )); CALL_SUBTEST_2(( test_incomplete_cholesky_T, int>() )); CALL_SUBTEST_3(( test_incomplete_cholesky_T() )); + +#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 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 A; + A = m * m.transpose(); + + Eigen::ConjugateGradient, + Eigen::Lower | Eigen::Upper, + Eigen::IncompleteCholesky > solver( A ); + VERIFY(solver.preconditioner().info() == Eigen::Success); + VERIFY(solver.info() == Eigen::Success); + } +#endif }