IncompleteCholesky.h (15036B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 5 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H 12 #define EIGEN_INCOMPLETE_CHOlESKY_H 13 14 #include <vector> 15 #include <list> 16 17 namespace Eigen { 18 /** 19 * \brief Modified Incomplete Cholesky with dual threshold 20 * 21 * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with 22 * Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999 23 * 24 * \tparam Scalar the scalar type of the input matrices 25 * \tparam _UpLo The triangular part that will be used for the computations. It can be Lower 26 * or Upper. Default is Lower. 27 * \tparam _OrderingType The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<int>, 28 * unless EIGEN_MPL2_ONLY is defined, in which case the default is NaturalOrdering<int>. 29 * 30 * \implsparsesolverconcept 31 * 32 * It performs the following incomplete factorization: \f$ S P A P' S \approx L L' \f$ 33 * where L is a lower triangular factor, S is a diagonal scaling matrix, and P is a 34 * fill-in reducing permutation as computed by the ordering method. 35 * 36 * \b Shifting \b strategy: Let \f$ B = S P A P' S \f$ be the scaled matrix on which the factorization is carried out, 37 * and \f$ \beta \f$ be the minimum value of the diagonal. If \f$ \beta > 0 \f$ then, the factorization is directly performed 38 * on the matrix B. Otherwise, the factorization is performed on the shifted matrix \f$ B + (\sigma+|\beta| I \f$ where 39 * \f$ \sigma \f$ is the initial shift value as returned and set by setInitialShift() method. The default value is \f$ \sigma = 10^{-3} \f$. 40 * If the factorization fails, then the shift in doubled until it succeed or a maximum of ten attempts. If it still fails, as returned by 41 * the info() method, then you can either increase the initial shift, or better use another preconditioning technique. 42 * 43 */ 44 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> > 45 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > 46 { 47 protected: 48 typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base; 49 using Base::m_isInitialized; 50 public: 51 typedef typename NumTraits<Scalar>::Real RealScalar; 52 typedef _OrderingType OrderingType; 53 typedef typename OrderingType::PermutationType PermutationType; 54 typedef typename PermutationType::StorageIndex StorageIndex; 55 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType; 56 typedef Matrix<Scalar,Dynamic,1> VectorSx; 57 typedef Matrix<RealScalar,Dynamic,1> VectorRx; 58 typedef Matrix<StorageIndex,Dynamic, 1> VectorIx; 59 typedef std::vector<std::list<StorageIndex> > VectorList; 60 enum { UpLo = _UpLo }; 61 enum { 62 ColsAtCompileTime = Dynamic, 63 MaxColsAtCompileTime = Dynamic 64 }; 65 public: 66 67 /** Default constructor leaving the object in a partly non-initialized stage. 68 * 69 * You must call compute() or the pair analyzePattern()/factorize() to make it valid. 70 * 71 * \sa IncompleteCholesky(const MatrixType&) 72 */ 73 IncompleteCholesky() : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) {} 74 75 /** Constructor computing the incomplete factorization for the given matrix \a matrix. 76 */ 77 template<typename MatrixType> 78 IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) 79 { 80 compute(matrix); 81 } 82 83 /** \returns number of rows of the factored matrix */ 84 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_L.rows(); } 85 86 /** \returns number of columns of the factored matrix */ 87 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_L.cols(); } 88 89 90 /** \brief Reports whether previous computation was successful. 91 * 92 * It triggers an assertion if \c *this has not been initialized through the respective constructor, 93 * or a call to compute() or analyzePattern(). 94 * 95 * \returns \c Success if computation was successful, 96 * \c NumericalIssue if the matrix appears to be negative. 97 */ 98 ComputationInfo info() const 99 { 100 eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized."); 101 return m_info; 102 } 103 104 /** \brief Set the initial shift parameter \f$ \sigma \f$. 105 */ 106 void setInitialShift(RealScalar shift) { m_initialShift = shift; } 107 108 /** \brief Computes the fill reducing permutation vector using the sparsity pattern of \a mat 109 */ 110 template<typename MatrixType> 111 void analyzePattern(const MatrixType& mat) 112 { 113 OrderingType ord; 114 PermutationType pinv; 115 ord(mat.template selfadjointView<UpLo>(), pinv); 116 if(pinv.size()>0) m_perm = pinv.inverse(); 117 else m_perm.resize(0); 118 m_L.resize(mat.rows(), mat.cols()); 119 m_analysisIsOk = true; 120 m_isInitialized = true; 121 m_info = Success; 122 } 123 124 /** \brief Performs the numerical factorization of the input matrix \a mat 125 * 126 * The method analyzePattern() or compute() must have been called beforehand 127 * with a matrix having the same pattern. 128 * 129 * \sa compute(), analyzePattern() 130 */ 131 template<typename MatrixType> 132 void factorize(const MatrixType& mat); 133 134 /** Computes or re-computes the incomplete Cholesky factorization of the input matrix \a mat 135 * 136 * It is a shortcut for a sequential call to the analyzePattern() and factorize() methods. 137 * 138 * \sa analyzePattern(), factorize() 139 */ 140 template<typename MatrixType> 141 void compute(const MatrixType& mat) 142 { 143 analyzePattern(mat); 144 factorize(mat); 145 } 146 147 // internal 148 template<typename Rhs, typename Dest> 149 void _solve_impl(const Rhs& b, Dest& x) const 150 { 151 eigen_assert(m_factorizationIsOk && "factorize() should be called first"); 152 if (m_perm.rows() == b.rows()) x = m_perm * b; 153 else x = b; 154 x = m_scale.asDiagonal() * x; 155 x = m_L.template triangularView<Lower>().solve(x); 156 x = m_L.adjoint().template triangularView<Upper>().solve(x); 157 x = m_scale.asDiagonal() * x; 158 if (m_perm.rows() == b.rows()) 159 x = m_perm.inverse() * x; 160 } 161 162 /** \returns the sparse lower triangular factor L */ 163 const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; } 164 165 /** \returns a vector representing the scaling factor S */ 166 const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; } 167 168 /** \returns the fill-in reducing permutation P (can be empty for a natural ordering) */ 169 const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; } 170 171 protected: 172 FactorType m_L; // The lower part stored in CSC 173 VectorRx m_scale; // The vector for scaling the matrix 174 RealScalar m_initialShift; // The initial shift parameter 175 bool m_analysisIsOk; 176 bool m_factorizationIsOk; 177 ComputationInfo m_info; 178 PermutationType m_perm; 179 180 private: 181 inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); 182 }; 183 184 // Based on the following paper: 185 // C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with 186 // Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999 187 // http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf 188 template<typename Scalar, int _UpLo, typename OrderingType> 189 template<typename _MatrixType> 190 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) 191 { 192 using std::sqrt; 193 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 194 195 // 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 196 197 // Apply the fill-reducing permutation computed in analyzePattern() 198 if (m_perm.rows() == mat.rows() ) // To detect the null permutation 199 { 200 // The temporary is needed to make sure that the diagonal entry is properly sorted 201 FactorType tmp(mat.rows(), mat.cols()); 202 tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm); 203 m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>(); 204 } 205 else 206 { 207 m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); 208 } 209 210 Index n = m_L.cols(); 211 Index nnz = m_L.nonZeros(); 212 Map<VectorSx> vals(m_L.valuePtr(), nnz); //values 213 Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices 214 Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row 215 VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization 216 VectorList listCol(n); // listCol(j) is a linked list of columns to update column j 217 VectorSx col_vals(n); // Store a nonzero values in each column 218 VectorIx col_irow(n); // Row indices of nonzero elements in each column 219 VectorIx col_pattern(n); 220 col_pattern.fill(-1); 221 StorageIndex col_nnz; 222 223 224 // Computes the scaling factors 225 m_scale.resize(n); 226 m_scale.setZero(); 227 for (Index j = 0; j < n; j++) 228 for (Index k = colPtr[j]; k < colPtr[j+1]; k++) 229 { 230 m_scale(j) += numext::abs2(vals(k)); 231 if(rowIdx[k]!=j) 232 m_scale(rowIdx[k]) += numext::abs2(vals(k)); 233 } 234 235 m_scale = m_scale.cwiseSqrt().cwiseSqrt(); 236 237 for (Index j = 0; j < n; ++j) 238 if(m_scale(j)>(std::numeric_limits<RealScalar>::min)()) 239 m_scale(j) = RealScalar(1)/m_scale(j); 240 else 241 m_scale(j) = 1; 242 243 // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster) 244 245 // Scale and compute the shift for the matrix 246 RealScalar mindiag = NumTraits<RealScalar>::highest(); 247 for (Index j = 0; j < n; j++) 248 { 249 for (Index k = colPtr[j]; k < colPtr[j+1]; k++) 250 vals[k] *= (m_scale(j)*m_scale(rowIdx[k])); 251 eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored"); 252 mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag); 253 } 254 255 FactorType L_save = m_L; 256 257 RealScalar shift = 0; 258 if(mindiag <= RealScalar(0.)) 259 shift = m_initialShift - mindiag; 260 261 m_info = NumericalIssue; 262 263 // Try to perform the incomplete factorization using the current shift 264 int iter = 0; 265 do 266 { 267 // Apply the shift to the diagonal elements of the matrix 268 for (Index j = 0; j < n; j++) 269 vals[colPtr[j]] += shift; 270 271 // jki version of the Cholesky factorization 272 Index j=0; 273 for (; j < n; ++j) 274 { 275 // Left-looking factorization of the j-th column 276 // First, load the j-th column into col_vals 277 Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored 278 col_nnz = 0; 279 for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++) 280 { 281 StorageIndex l = rowIdx[i]; 282 col_vals(col_nnz) = vals[i]; 283 col_irow(col_nnz) = l; 284 col_pattern(l) = col_nnz; 285 col_nnz++; 286 } 287 { 288 typename std::list<StorageIndex>::iterator k; 289 // Browse all previous columns that will update column j 290 for(k = listCol[j].begin(); k != listCol[j].end(); k++) 291 { 292 Index jk = firstElt(*k); // First element to use in the column 293 eigen_internal_assert(rowIdx[jk]==j); 294 Scalar v_j_jk = numext::conj(vals[jk]); 295 296 jk += 1; 297 for (Index i = jk; i < colPtr[*k+1]; i++) 298 { 299 StorageIndex l = rowIdx[i]; 300 if(col_pattern[l]<0) 301 { 302 col_vals(col_nnz) = vals[i] * v_j_jk; 303 col_irow[col_nnz] = l; 304 col_pattern(l) = col_nnz; 305 col_nnz++; 306 } 307 else 308 col_vals(col_pattern[l]) -= vals[i] * v_j_jk; 309 } 310 updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol); 311 } 312 } 313 314 // Scale the current column 315 if(numext::real(diag) <= 0) 316 { 317 if(++iter>=10) 318 return; 319 320 // increase shift 321 shift = numext::maxi(m_initialShift,RealScalar(2)*shift); 322 // restore m_L, col_pattern, and listCol 323 vals = Map<const VectorSx>(L_save.valuePtr(), nnz); 324 rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz); 325 colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1); 326 col_pattern.fill(-1); 327 for(Index i=0; i<n; ++i) 328 listCol[i].clear(); 329 330 break; 331 } 332 333 RealScalar rdiag = sqrt(numext::real(diag)); 334 vals[colPtr[j]] = rdiag; 335 for (Index k = 0; k<col_nnz; ++k) 336 { 337 Index i = col_irow[k]; 338 //Scale 339 col_vals(k) /= rdiag; 340 //Update the remaining diagonals with col_vals 341 vals[colPtr[i]] -= numext::abs2(col_vals(k)); 342 } 343 // Select the largest p elements 344 // p is the original number of elements in the column (without the diagonal) 345 Index p = colPtr[j+1] - colPtr[j] - 1 ; 346 Ref<VectorSx> cvals = col_vals.head(col_nnz); 347 Ref<VectorIx> cirow = col_irow.head(col_nnz); 348 internal::QuickSplit(cvals,cirow, p); 349 // Insert the largest p elements in the matrix 350 Index cpt = 0; 351 for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++) 352 { 353 vals[i] = col_vals(cpt); 354 rowIdx[i] = col_irow(cpt); 355 // restore col_pattern: 356 col_pattern(col_irow(cpt)) = -1; 357 cpt++; 358 } 359 // Get the first smallest row index and put it after the diagonal element 360 Index jk = colPtr(j)+1; 361 updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol); 362 } 363 364 if(j==n) 365 { 366 m_factorizationIsOk = true; 367 m_info = Success; 368 } 369 } while(m_info!=Success); 370 } 371 372 template<typename Scalar, int _UpLo, typename OrderingType> 373 inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol) 374 { 375 if (jk < colPtr(col+1) ) 376 { 377 Index p = colPtr(col+1) - jk; 378 Index minpos; 379 rowIdx.segment(jk,p).minCoeff(&minpos); 380 minpos += jk; 381 if (rowIdx(minpos) != rowIdx(jk)) 382 { 383 //Swap 384 std::swap(rowIdx(jk),rowIdx(minpos)); 385 std::swap(vals(jk),vals(minpos)); 386 } 387 firstElt(col) = internal::convert_index<StorageIndex,Index>(jk); 388 listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col)); 389 } 390 } 391 392 } // end namespace Eigen 393 394 #endif