From e31fc50280d412ccbc1dfedb625c48846e4060e3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 4 Aug 2015 16:16:02 +0200 Subject: Numerous fixes for IncompleteCholesky. Still have to make it fully exploit the sparse structure of the L factor, and improve robustness to illconditionned problems. --- .../src/IterativeSolvers/IncompleteCholesky.h | 177 ++++++++++++--------- 1 file changed, 103 insertions(+), 74 deletions(-) (limited to 'unsupported') diff --git a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h index 35cfa315d..a62672201 100644 --- a/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h +++ b/unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h @@ -26,25 +26,29 @@ namespace Eigen { * \tparam _OrderingType */ -template > +template > class IncompleteCholesky : public SparseSolverBase > { protected: typedef SparseSolverBase > Base; using Base::m_isInitialized; public: - typedef SparseMatrix MatrixType; + typedef typename NumTraits::Real RealScalar; typedef _OrderingType OrderingType; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - typedef PermutationMatrix PermutationType; - typedef Matrix ScalarType; - typedef Matrix IndexType; - typedef std::vector > VectorList; + typedef typename OrderingType::PermutationType PermutationType; + typedef typename PermutationType::StorageIndex StorageIndex; + typedef SparseMatrix FactorType; + typedef FactorType MatrixType; + typedef Matrix VectorSx; + typedef Matrix VectorRx; + typedef Matrix VectorIx; + typedef std::vector > VectorList; enum { UpLo = _UpLo }; public: - IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {} - IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false) + IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {} + + template + IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false) { compute(matrix); } @@ -68,7 +72,7 @@ class IncompleteCholesky : public SparseSolverBase(), m_perm); + PermutationType pinv; + ord(mat.template selfadjointView(), pinv); + if(pinv.size()>0) m_perm = pinv.inverse(); + else m_perm.resize(0); m_analysisIsOk = true; } @@ -85,7 +92,7 @@ class IncompleteCholesky : public SparseSolverBase - void compute (const MatrixType& matrix) + void compute(const MatrixType& matrix) { analyzePattern(matrix); factorize(matrix); @@ -95,30 +102,28 @@ class IncompleteCholesky : public SparseSolverBase().solve(x); + x = m_L.adjoint().template triangularView().solve(x); + x = m_scale.asDiagonal() * x; if (m_perm.rows() == b.rows()) - x = m_perm.inverse() * b; - else - x = b; - x = m_scal.asDiagonal() * x; - x = m_L.template triangularView().solve(x); - x = m_L.adjoint().template triangularView().solve(x); - if (m_perm.rows() == b.rows()) - x = m_perm * x; - x = m_scal.asDiagonal() * x; + x = m_perm.inverse() * x; + } protected: - SparseMatrix m_L; // The lower part stored in CSC - ScalarType m_scal; // The vector for scaling the matrix - Scalar m_shift; //The initial shift parameter + FactorType m_L; // The lower part stored in CSC + VectorRx m_scale; // The vector for scaling the matrix + Scalar m_initialShift; // The initial shift parameter bool m_analysisIsOk; bool m_factorizationIsOk; ComputationInfo m_info; PermutationType m_perm; private: - template - inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol); + inline void updateList(Ref colPtr, Ref rowIdx, Ref vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); }; template @@ -128,93 +133,118 @@ void IncompleteCholesky::factorize(const _MatrixType using std::sqrt; eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); - // Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added + // Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added + + m_L.resize(mat.rows(), mat.cols()); // Apply the fill-reducing permutation computed in analyzePattern() if (m_perm.rows() == mat.rows() ) // To detect the null permutation - m_L.template selfadjointView() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); + { + // The temporary is needed to make sure that the diagonal entry is properly sorted + FactorType tmp(mat.rows(), mat.cols()); + tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); + m_L.template selfadjointView() = tmp.template selfadjointView(); + } else + { m_L.template selfadjointView() = mat.template selfadjointView<_UpLo>(); + } Index n = m_L.cols(); Index nnz = m_L.nonZeros(); - Map vals(m_L.valuePtr(), nnz); //values - Map rowIdx(m_L.innerIndexPtr(), nnz); //Row indices - Map colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row - IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization - VectorList listCol(n); // listCol(j) is a linked list of columns to update column j - ScalarType curCol(n); // Store a nonzero values in each column - IndexType irow(n); // Row indices of nonzero elements in each column + Map vals(m_L.valuePtr(), nnz); //values + Map rowIdx(m_L.innerIndexPtr(), nnz); //Row indices + Map colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row + VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization + VectorList listCol(n); // listCol(j) is a linked list of columns to update column j + VectorSx curCol(n); // Store a nonzero values in each column + VectorIx irow(n); // Row indices of nonzero elements in each column // Computes the scaling factors - m_scal.resize(n); - for (int j = 0; j < n; j++) - { - m_scal(j) = m_L.col(j).norm(); - m_scal(j) = sqrt(m_scal(j)); - } + m_scale.resize(n); + m_scale.setZero(); + for (Index j = 0; j < n; j++) + for (Index k = colPtr[j]; k < colPtr[j+1]; k++) + { + m_scale(j) += numext::abs2(vals(k)); + if(rowIdx[k]!=j) + m_scale(rowIdx[k]) += numext::abs2(vals(k)); + } + + m_scale = m_scale.cwiseSqrt().cwiseSqrt(); + // Scale and compute the shift for the matrix - Scalar mindiag = vals[0]; - for (int j = 0; j < n; j++){ - for (int k = colPtr[j]; k < colPtr[j+1]; k++) - vals[k] /= (m_scal(j) * m_scal(rowIdx[k])); - mindiag = numext::mini(vals[colPtr[j]], mindiag); + RealScalar mindiag = NumTraits::highest(); + for (Index j = 0; j < n; j++) + { + for (Index k = colPtr[j]; k < colPtr[j+1]; k++) + vals[k] /= (m_scale(j)*m_scale(rowIdx[k])); + 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); } - if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag; + RealScalar shift = 0; + if(mindiag <= RealScalar(0.)) + shift = m_initialShift - mindiag; + // Apply the shift to the diagonal elements of the matrix - for (int j = 0; j < n; j++) - vals[colPtr[j]] += m_shift; + for (Index j = 0; j < n; j++) + vals[colPtr[j]] += shift; + // jki version of the Cholesky factorization - for (int j=0; j < n; ++j) + for (Index j=0; j < n; ++j) { - //Left-looking factorize the column j - // First, load the jth column into curCol + // Left-looking factorization of the j-th column + // First, load the j-th column into curCol Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored curCol.setZero(); - irow.setLinSpaced(n,0,n-1); - for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++) + irow.setLinSpaced(n,0,internal::convert_index(n-1)); + for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++) { curCol(rowIdx[i]) = vals[i]; irow(rowIdx[i]) = rowIdx[i]; } - std::list::iterator k; + typename std::list::iterator k; // Browse all previous columns that will update column j for(k = listCol[j].begin(); k != listCol[j].end(); k++) { - int jk = firstElt(*k); // First element to use in the column + 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 (int i = jk; i < colPtr[*k+1]; i++) - { - curCol(rowIdx[i]) -= vals[i] * vals[jk] ; - } + for (Index i = jk; i < colPtr[*k+1]; i++) + curCol(rowIdx[i]) -= vals[i] * v_j_jk; updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); } // Scale the current column - if(RealScalar(diag) <= 0) + if(numext::real(diag) <= 0) { - std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n"; + //std::cerr << "\nNegative diagonal during Incomplete factorization at position " << j << " (value = " << diag << ")\n"; m_info = NumericalIssue; return; } - RealScalar rdiag = sqrt(RealScalar(diag)); + + RealScalar rdiag = sqrt(numext::real(diag)); vals[colPtr[j]] = rdiag; - for (int i = j+1; i < n; i++) + // TODO, the following should iterate on the structurally non-zeros only + for (Index i = j+1; i < n; i++) { //Scale curCol(i) /= rdiag; //Update the remaining diagonals with curCol - vals[colPtr[i]] -= curCol(i) * curCol(i); + vals[colPtr[i]] -= numext::abs2(curCol(i)); } // Select the largest p elements - // p is the original number of elements in the column (without the diagonal) - int p = colPtr[j+1] - colPtr[j] - 1 ; + // p is the original number of elements in the column (without the diagonal) + // TODO, QuickSplit should operate on the structurally non zeros only. + Index p = colPtr[j+1] - colPtr[j] - 1 ; internal::QuickSplit(curCol, irow, p); // Insert the largest p elements in the matrix - int cpt = 0; - for (int i = colPtr[j]+1; i < colPtr[j+1]; i++) + Index cpt = 0; + for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++) { vals[i] = curCol(cpt); rowIdx[i] = irow(cpt); @@ -230,8 +260,7 @@ void IncompleteCholesky::factorize(const _MatrixType } template -template -inline void IncompleteCholesky::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol) +inline void IncompleteCholesky::updateList(Ref colPtr, Ref rowIdx, Ref vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol) { if (jk < colPtr(col+1) ) { @@ -245,8 +274,8 @@ inline void IncompleteCholesky::updateList(const Idx std::swap(rowIdx(jk),rowIdx(minpos)); std::swap(vals(jk),vals(minpos)); } - firstElt(col) = jk; - listCol[rowIdx(jk)].push_back(col); + firstElt(col) = internal::convert_index(jk); + listCol[rowIdx(jk)].push_back(internal::convert_index(col)); } } -- cgit v1.2.3