Eigen  3.2.92
IncompleteCholesky.h
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 {
42 template <typename Scalar, int _UpLo = Lower, typename _OrderingType =
43 #ifndef EIGEN_MPL2_ONLY
44 AMDOrdering<int>
45 #else
46 NaturalOrdering<int>
47 #endif
48 >
49 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
50 {
51  protected:
53  using Base::m_isInitialized;
54  public:
55  typedef typename NumTraits<Scalar>::Real RealScalar;
56  typedef _OrderingType OrderingType;
57  typedef typename OrderingType::PermutationType PermutationType;
58  typedef typename PermutationType::StorageIndex StorageIndex;
63  typedef std::vector<std::list<StorageIndex> > VectorList;
64  enum { UpLo = _UpLo };
65  enum {
66  ColsAtCompileTime = Dynamic,
67  MaxColsAtCompileTime = Dynamic
68  };
69  public:
70 
77  IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
78 
81  template<typename MatrixType>
82  IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
83  {
84  compute(matrix);
85  }
86 
88  Index rows() const { return m_L.rows(); }
89 
91  Index cols() const { return m_L.cols(); }
92 
93 
103  {
104  eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
105  return m_info;
106  }
107 
110  void setInitialShift(RealScalar shift) { m_initialShift = shift; }
111 
114  template<typename MatrixType>
115  void analyzePattern(const MatrixType& mat)
116  {
117  OrderingType ord;
118  PermutationType pinv;
119  ord(mat.template selfadjointView<UpLo>(), pinv);
120  if(pinv.size()>0) m_perm = pinv.inverse();
121  else m_perm.resize(0);
122  m_L.resize(mat.rows(), mat.cols());
123  m_analysisIsOk = true;
124  m_isInitialized = true;
125  m_info = Success;
126  }
127 
135  template<typename MatrixType>
136  void factorize(const MatrixType& mat);
137 
144  template<typename MatrixType>
145  void compute(const MatrixType& mat)
146  {
147  analyzePattern(mat);
148  factorize(mat);
149  }
150 
151  // internal
152  template<typename Rhs, typename Dest>
153  void _solve_impl(const Rhs& b, Dest& x) const
154  {
155  eigen_assert(m_factorizationIsOk && "factorize() should be called first");
156  if (m_perm.rows() == b.rows()) x = m_perm * b;
157  else x = b;
158  x = m_scale.asDiagonal() * x;
159  x = m_L.template triangularView<Lower>().solve(x);
160  x = m_L.adjoint().template triangularView<Upper>().solve(x);
161  x = m_scale.asDiagonal() * x;
162  if (m_perm.rows() == b.rows())
163  x = m_perm.inverse() * x;
164  }
165 
167  const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }
168 
170  const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }
171 
173  const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; }
174 
175  protected:
176  FactorType m_L; // The lower part stored in CSC
177  VectorRx m_scale; // The vector for scaling the matrix
178  RealScalar m_initialShift; // The initial shift parameter
179  bool m_analysisIsOk;
180  bool m_factorizationIsOk;
181  ComputationInfo m_info;
182  PermutationType m_perm;
183 
184  private:
185  inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
186 };
187 
188 template<typename Scalar, int _UpLo, typename OrderingType>
189 template<typename _MatrixType>
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  // FIXME 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  RealScalar shift = 0;
256  if(mindiag <= RealScalar(0.))
257  shift = m_initialShift - mindiag;
258 
259  // Apply the shift to the diagonal elements of the matrix
260  for (Index j = 0; j < n; j++)
261  vals[colPtr[j]] += shift;
262 
263  // jki version of the Cholesky factorization
264  for (Index j=0; j < n; ++j)
265  {
266  // Left-looking factorization of the j-th column
267  // First, load the j-th column into col_vals
268  Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
269  col_nnz = 0;
270  for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
271  {
272  StorageIndex l = rowIdx[i];
273  col_vals(col_nnz) = vals[i];
274  col_irow(col_nnz) = l;
275  col_pattern(l) = col_nnz;
276  col_nnz++;
277  }
278  {
279  typename std::list<StorageIndex>::iterator k;
280  // Browse all previous columns that will update column j
281  for(k = listCol[j].begin(); k != listCol[j].end(); k++)
282  {
283  Index jk = firstElt(*k); // First element to use in the column
284  eigen_internal_assert(rowIdx[jk]==j);
285  Scalar v_j_jk = numext::conj(vals[jk]);
286 
287  jk += 1;
288  for (Index i = jk; i < colPtr[*k+1]; i++)
289  {
290  StorageIndex l = rowIdx[i];
291  if(col_pattern[l]<0)
292  {
293  col_vals(col_nnz) = vals[i] * v_j_jk;
294  col_irow[col_nnz] = l;
295  col_pattern(l) = col_nnz;
296  col_nnz++;
297  }
298  else
299  col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
300  }
301  updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
302  }
303  }
304 
305  // Scale the current column
306  if(numext::real(diag) <= 0)
307  {
308  m_info = NumericalIssue;
309  return;
310  }
311 
312  RealScalar rdiag = sqrt(numext::real(diag));
313  vals[colPtr[j]] = rdiag;
314  for (Index k = 0; k<col_nnz; ++k)
315  {
316  Index i = col_irow[k];
317  //Scale
318  col_vals(k) /= rdiag;
319  //Update the remaining diagonals with col_vals
320  vals[colPtr[i]] -= numext::abs2(col_vals(k));
321  }
322  // Select the largest p elements
323  // p is the original number of elements in the column (without the diagonal)
324  Index p = colPtr[j+1] - colPtr[j] - 1 ;
325  Ref<VectorSx> cvals = col_vals.head(col_nnz);
326  Ref<VectorIx> cirow = col_irow.head(col_nnz);
327  internal::QuickSplit(cvals,cirow, p);
328  // Insert the largest p elements in the matrix
329  Index cpt = 0;
330  for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
331  {
332  vals[i] = col_vals(cpt);
333  rowIdx[i] = col_irow(cpt);
334  // restore col_pattern:
335  col_pattern(col_irow(cpt)) = -1;
336  cpt++;
337  }
338  // Get the first smallest row index and put it after the diagonal element
339  Index jk = colPtr(j)+1;
340  updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
341  }
342  m_factorizationIsOk = true;
343  m_info = Success;
344 }
345 
346 template<typename Scalar, int _UpLo, typename OrderingType>
347 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)
348 {
349  if (jk < colPtr(col+1) )
350  {
351  Index p = colPtr(col+1) - jk;
352  Index minpos;
353  rowIdx.segment(jk,p).minCoeff(&minpos);
354  minpos += jk;
355  if (rowIdx(minpos) != rowIdx(jk))
356  {
357  //Swap
358  std::swap(rowIdx(jk),rowIdx(minpos));
359  std::swap(vals(jk),vals(minpos));
360  }
361  firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
362  listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
363  }
364 }
365 
366 } // end namespace Eigen
367 
368 #endif
Index cols() const
Definition: IncompleteCholesky.h:91
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
Index cols() const
Definition: SparseMatrix.h:133
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:49
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:173
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:82
const Solve< IncompleteCholesky< Scalar, _UpLo, _OrderingType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:74
Definition: LDLT.h:16
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:161
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
IncompleteCholesky()
Definition: IncompleteCholesky.h:77
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:520
Definition: Constants.h:204
Definition: Constants.h:434
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:170
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:611
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:110
Definition: Constants.h:432
Index rows() const
Definition: IncompleteCholesky.h:88
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:152
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:186
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:115
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:167
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:145
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:102
const Scalar * valuePtr() const
Definition: SparseMatrix.h:143
ComputationInfo
Definition: Constants.h:430
Index rows() const
Definition: SparseMatrix.h:131