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