Eigen  3.3.3
IncompleteCholesky.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H
00012 #define EIGEN_INCOMPLETE_CHOlESKY_H
00013 
00014 #include <vector>
00015 #include <list>
00016 
00017 namespace Eigen {  
00044 template <typename Scalar, int _UpLo = Lower, typename _OrderingType =
00045 #ifndef EIGEN_MPL2_ONLY
00046 AMDOrdering<int>
00047 #else
00048 NaturalOrdering<int>
00049 #endif
00050 >
00051 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
00052 {
00053   protected:
00054     typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
00055     using Base::m_isInitialized;
00056   public:
00057     typedef typename NumTraits<Scalar>::Real RealScalar; 
00058     typedef _OrderingType OrderingType;
00059     typedef typename OrderingType::PermutationType PermutationType;
00060     typedef typename PermutationType::StorageIndex StorageIndex; 
00061     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
00062     typedef Matrix<Scalar,Dynamic,1> VectorSx;
00063     typedef Matrix<RealScalar,Dynamic,1> VectorRx;
00064     typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
00065     typedef std::vector<std::list<StorageIndex> > VectorList; 
00066     enum { UpLo = _UpLo };
00067     enum {
00068       ColsAtCompileTime = Dynamic,
00069       MaxColsAtCompileTime = Dynamic
00070     };
00071   public:
00072 
00079     IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
00080     
00083     template<typename MatrixType>
00084     IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
00085     {
00086       compute(matrix);
00087     }
00088     
00090     Index rows() const { return m_L.rows(); }
00091     
00093     Index cols() const { return m_L.cols(); }
00094     
00095 
00104     ComputationInfo info() const
00105     {
00106       eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
00107       return m_info;
00108     }
00109     
00112     void setInitialShift(RealScalar shift) { m_initialShift = shift; }
00113     
00116     template<typename MatrixType>
00117     void analyzePattern(const MatrixType& mat)
00118     {
00119       OrderingType ord; 
00120       PermutationType pinv;
00121       ord(mat.template selfadjointView<UpLo>(), pinv); 
00122       if(pinv.size()>0) m_perm = pinv.inverse();
00123       else              m_perm.resize(0);
00124       m_L.resize(mat.rows(), mat.cols());
00125       m_analysisIsOk = true;
00126       m_isInitialized = true;
00127       m_info = Success;
00128     }
00129     
00137     template<typename MatrixType>
00138     void factorize(const MatrixType& mat);
00139     
00146     template<typename MatrixType>
00147     void compute(const MatrixType& mat)
00148     {
00149       analyzePattern(mat);
00150       factorize(mat);
00151     }
00152     
00153     // internal
00154     template<typename Rhs, typename Dest>
00155     void _solve_impl(const Rhs& b, Dest& x) const
00156     {
00157       eigen_assert(m_factorizationIsOk && "factorize() should be called first");
00158       if (m_perm.rows() == b.rows())  x = m_perm * b;
00159       else                            x = b;
00160       x = m_scale.asDiagonal() * x;
00161       x = m_L.template triangularView<Lower>().solve(x);
00162       x = m_L.adjoint().template triangularView<Upper>().solve(x);
00163       x = m_scale.asDiagonal() * x;
00164       if (m_perm.rows() == b.rows())
00165         x = m_perm.inverse() * x;
00166     }
00167 
00169     const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }
00170 
00172     const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }
00173 
00175     const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; }
00176 
00177   protected:
00178     FactorType m_L;              // The lower part stored in CSC
00179     VectorRx m_scale;            // The vector for scaling the matrix 
00180     RealScalar m_initialShift;   // The initial shift parameter
00181     bool m_analysisIsOk; 
00182     bool m_factorizationIsOk; 
00183     ComputationInfo m_info;
00184     PermutationType m_perm; 
00185 
00186   private:
00187     inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); 
00188 }; 
00189 
00190 // Based on the following paper:
00191 //   C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
00192 //   Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
00193 //   http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
00194 template<typename Scalar, int _UpLo, typename OrderingType>
00195 template<typename _MatrixType>
00196 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
00197 {
00198   using std::sqrt;
00199   eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 
00200     
00201   // 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
00202   
00203   // Apply the fill-reducing permutation computed in analyzePattern()
00204   if (m_perm.rows() == mat.rows() ) // To detect the null permutation
00205   {
00206     // The temporary is needed to make sure that the diagonal entry is properly sorted
00207     FactorType tmp(mat.rows(), mat.cols());
00208     tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
00209     m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
00210   }
00211   else
00212   {
00213     m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
00214   }
00215   
00216   Index n = m_L.cols(); 
00217   Index nnz = m_L.nonZeros();
00218   Map<VectorSx> vals(m_L.valuePtr(), nnz);         //values
00219   Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);  //Row indices
00220   Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
00221   VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
00222   VectorList listCol(n);  // listCol(j) is a linked list of columns to update column j
00223   VectorSx col_vals(n);   // Store a  nonzero values in each column
00224   VectorIx col_irow(n);   // Row indices of nonzero elements in each column
00225   VectorIx col_pattern(n);
00226   col_pattern.fill(-1);
00227   StorageIndex col_nnz;
00228   
00229   
00230   // Computes the scaling factors 
00231   m_scale.resize(n);
00232   m_scale.setZero();
00233   for (Index j = 0; j < n; j++)
00234     for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
00235     {
00236       m_scale(j) += numext::abs2(vals(k));
00237       if(rowIdx[k]!=j)
00238         m_scale(rowIdx[k]) += numext::abs2(vals(k));
00239     }
00240   
00241   m_scale = m_scale.cwiseSqrt().cwiseSqrt();
00242 
00243   for (Index j = 0; j < n; ++j)
00244     if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
00245       m_scale(j) = RealScalar(1)/m_scale(j);
00246     else
00247       m_scale(j) = 1;
00248 
00249   // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
00250   
00251   // Scale and compute the shift for the matrix 
00252   RealScalar mindiag = NumTraits<RealScalar>::highest();
00253   for (Index j = 0; j < n; j++)
00254   {
00255     for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
00256       vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
00257     eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
00258     mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
00259   }
00260 
00261   FactorType L_save = m_L;
00262   
00263   RealScalar shift = 0;
00264   if(mindiag <= RealScalar(0.))
00265     shift = m_initialShift - mindiag;
00266 
00267   m_info = NumericalIssue;
00268 
00269   // Try to perform the incomplete factorization using the current shift
00270   int iter = 0;
00271   do
00272   {
00273     // Apply the shift to the diagonal elements of the matrix
00274     for (Index j = 0; j < n; j++)
00275       vals[colPtr[j]] += shift;
00276 
00277     // jki version of the Cholesky factorization
00278     Index j=0;
00279     for (; j < n; ++j)
00280     {
00281       // Left-looking factorization of the j-th column
00282       // First, load the j-th column into col_vals
00283       Scalar diag = vals[colPtr[j]];  // It is assumed that only the lower part is stored
00284       col_nnz = 0;
00285       for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
00286       {
00287         StorageIndex l = rowIdx[i];
00288         col_vals(col_nnz) = vals[i];
00289         col_irow(col_nnz) = l;
00290         col_pattern(l) = col_nnz;
00291         col_nnz++;
00292       }
00293       {
00294         typename std::list<StorageIndex>::iterator k;
00295         // Browse all previous columns that will update column j
00296         for(k = listCol[j].begin(); k != listCol[j].end(); k++)
00297         {
00298           Index jk = firstElt(*k); // First element to use in the column
00299           eigen_internal_assert(rowIdx[jk]==j);
00300           Scalar v_j_jk = numext::conj(vals[jk]);
00301 
00302           jk += 1;
00303           for (Index i = jk; i < colPtr[*k+1]; i++)
00304           {
00305             StorageIndex l = rowIdx[i];
00306             if(col_pattern[l]<0)
00307             {
00308               col_vals(col_nnz) = vals[i] * v_j_jk;
00309               col_irow[col_nnz] = l;
00310               col_pattern(l) = col_nnz;
00311               col_nnz++;
00312             }
00313             else
00314               col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
00315           }
00316           updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
00317         }
00318       }
00319 
00320       // Scale the current column
00321       if(numext::real(diag) <= 0)
00322       {
00323         if(++iter>=10)
00324           return;
00325 
00326         // increase shift
00327         shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
00328         // restore m_L, col_pattern, and listCol
00329         vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
00330         rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
00331         colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
00332         col_pattern.fill(-1);
00333         for(Index i=0; i<n; ++i)
00334           listCol[i].clear();
00335 
00336         break;
00337       }
00338 
00339       RealScalar rdiag = sqrt(numext::real(diag));
00340       vals[colPtr[j]] = rdiag;
00341       for (Index k = 0; k<col_nnz; ++k)
00342       {
00343         Index i = col_irow[k];
00344         //Scale
00345         col_vals(k) /= rdiag;
00346         //Update the remaining diagonals with col_vals
00347         vals[colPtr[i]] -= numext::abs2(col_vals(k));
00348       }
00349       // Select the largest p elements
00350       // p is the original number of elements in the column (without the diagonal)
00351       Index p = colPtr[j+1] - colPtr[j] - 1 ;
00352       Ref<VectorSx> cvals = col_vals.head(col_nnz);
00353       Ref<VectorIx> cirow = col_irow.head(col_nnz);
00354       internal::QuickSplit(cvals,cirow, p);
00355       // Insert the largest p elements in the matrix
00356       Index cpt = 0;
00357       for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
00358       {
00359         vals[i] = col_vals(cpt);
00360         rowIdx[i] = col_irow(cpt);
00361         // restore col_pattern:
00362         col_pattern(col_irow(cpt)) = -1;
00363         cpt++;
00364       }
00365       // Get the first smallest row index and put it after the diagonal element
00366       Index jk = colPtr(j)+1;
00367       updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
00368     }
00369 
00370     if(j==n)
00371     {
00372       m_factorizationIsOk = true;
00373       m_info = Success;
00374     }
00375   } while(m_info!=Success);
00376 }
00377 
00378 template<typename Scalar, int _UpLo, typename OrderingType>
00379 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)
00380 {
00381   if (jk < colPtr(col+1) )
00382   {
00383     Index p = colPtr(col+1) - jk;
00384     Index minpos; 
00385     rowIdx.segment(jk,p).minCoeff(&minpos);
00386     minpos += jk;
00387     if (rowIdx(minpos) != rowIdx(jk))
00388     {
00389       //Swap
00390       std::swap(rowIdx(jk),rowIdx(minpos));
00391       std::swap(vals(jk),vals(minpos));
00392     }
00393     firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
00394     listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
00395   }
00396 }
00397 
00398 } // end namespace Eigen 
00399 
00400 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends