IncompleteCholesky.h
Go to the documentation of this file.
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 {
44 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
45 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
46 {
47  protected:
50  public:
52  typedef _OrderingType OrderingType;
53  typedef typename OrderingType::PermutationType PermutationType;
54  typedef typename PermutationType::StorageIndex StorageIndex;
59  typedef std::vector<std::list<StorageIndex> > VectorList;
60  enum { UpLo = _UpLo };
61  enum {
64  };
65  public:
66 
74 
77  template<typename MatrixType>
79  {
80  compute(matrix);
81  }
82 
85 
88 
89 
99  {
100  eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
101  return m_info;
102  }
103 
106  void setInitialShift(RealScalar shift) { m_initialShift = shift; }
107 
110  template<typename MatrixType>
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 
131  template<typename MatrixType>
132  void factorize(const MatrixType& mat);
133 
140  template<typename MatrixType>
141  void compute(const MatrixType& mat)
142  {
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 
163  const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }
164 
166  const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }
167 
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
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>
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
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>
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
Eigen::SparseMatrix::resize
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:626
Eigen::SparseMatrix::cols
Index cols() const
Definition: SparseMatrix.h:140
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:444
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::SparseMatrix< Scalar, ColMajor, StorageIndex >
Eigen::IncompleteCholesky::m_L
FactorType m_L
Definition: IncompleteCholesky.h:172
col
m col(1)
Eigen::IncompleteCholesky::m_initialShift
RealScalar m_initialShift
Definition: IncompleteCholesky.h:174
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Eigen::IncompleteCholesky::FactorType
SparseMatrix< Scalar, ColMajor, StorageIndex > FactorType
Definition: IncompleteCholesky.h:55
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::IncompleteCholesky::updateList
void updateList(Ref< const VectorIx > colPtr, Ref< VectorIx > rowIdx, Ref< VectorSx > vals, const Index &col, const Index &jk, VectorIx &firstElt, VectorList &listCol)
Definition: IncompleteCholesky.h:373
Eigen::IncompleteCholesky::_solve_impl
void _solve_impl(const Rhs &b, Dest &x) const
Definition: IncompleteCholesky.h:149
Eigen::IncompleteCholesky::compute
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:141
b
Scalar * b
Definition: benchVecAdd.cpp:17
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
gtsam::diag
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Eigen::IncompleteCholesky::IncompleteCholesky
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:78
EIGEN_CONSTEXPR
#define EIGEN_CONSTEXPR
Definition: Macros.h:787
Eigen::Success
@ Success
Definition: Constants.h:442
real
float real
Definition: datatypes.h:10
Eigen::IncompleteCholesky
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:45
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
eigen_internal_assert
#define eigen_internal_assert(x)
Definition: Macros.h:1043
Eigen::IncompleteCholesky::m_perm
PermutationType m_perm
Definition: IncompleteCholesky.h:178
Eigen::IncompleteCholesky::m_analysisIsOk
bool m_analysisIsOk
Definition: IncompleteCholesky.h:175
n
int n
Definition: BiCGSTAB_simple.cpp:1
Eigen::SparseMatrix::valuePtr
const Scalar * valuePtr() const
Definition: SparseMatrix.h:150
Eigen::IncompleteCholesky::Base
SparseSolverBase< IncompleteCholesky< Scalar, _UpLo, _OrderingType > > Base
Definition: IncompleteCholesky.h:48
Eigen::IncompleteCholesky::PermutationType
OrderingType::PermutationType PermutationType
Definition: IncompleteCholesky.h:53
Eigen::SparseMatrix::innerIndexPtr
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:159
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::IncompleteCholesky::VectorList
std::vector< std::list< StorageIndex > > VectorList
Definition: IncompleteCholesky.h:59
Eigen::IncompleteCholesky::m_factorizationIsOk
bool m_factorizationIsOk
Definition: IncompleteCholesky.h:176
Eigen::IncompleteCholesky::factorize
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
Eigen::IncompleteCholesky::VectorIx
Matrix< StorageIndex, Dynamic, 1 > VectorIx
Definition: IncompleteCholesky.h:58
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::numext::mini
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1085
Eigen::IncompleteCholesky::cols
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:87
l
static const Line3 l(Rot3(), 1, 1)
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
Eigen::IncompleteCholesky::OrderingType
_OrderingType OrderingType
Definition: IncompleteCholesky.h:52
Eigen::IncompleteCholesky::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: IncompleteCholesky.h:63
Eigen::IncompleteCholesky::scalingS
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:166
Eigen::IncompleteCholesky::permutationP
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:169
conj
AnnoyingScalar conj(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:104
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::SparseMatrix::outerIndexPtr
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:168
Eigen::IncompleteCholesky::IncompleteCholesky
IncompleteCholesky()
Definition: IncompleteCholesky.h:73
Eigen::IncompleteCholesky::ColsAtCompileTime
@ ColsAtCompileTime
Definition: IncompleteCholesky.h:62
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::IncompleteCholesky::RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: IncompleteCholesky.h:51
Eigen::IncompleteCholesky::matrixL
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:163
Eigen::SparseSolverBase
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Eigen::numext::abs2
EIGEN_DEVICE_FUNC bool abs2(bool x)
Definition: Eigen/src/Core/MathFunctions.h:1294
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
Eigen::IncompleteCholesky::VectorRx
Matrix< RealScalar, Dynamic, 1 > VectorRx
Definition: IncompleteCholesky.h:57
Eigen::internal::Rhs
@ Rhs
Definition: TensorContractionMapper.h:18
Eigen::IncompleteCholesky::StorageIndex
PermutationType::StorageIndex StorageIndex
Definition: IncompleteCholesky.h:54
iter
iterator iter(handle obj)
Definition: pytypes.h:2475
Eigen::IncompleteCholesky::m_info
ComputationInfo m_info
Definition: IncompleteCholesky.h:177
p
float * p
Definition: Tutorial_Map_using.cpp:9
Eigen::internal::QuickSplit
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
Definition: IncompleteLUT.h:29
Eigen::IncompleteCholesky::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:98
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::Matrix< Scalar, Dynamic, 1 >
Eigen::IncompleteCholesky::UpLo
@ UpLo
Definition: IncompleteCholesky.h:60
EIGEN_NOEXCEPT
#define EIGEN_NOEXCEPT
Definition: Macros.h:1418
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
Eigen::IncompleteCholesky::rows
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:84
Eigen::IncompleteCholesky::VectorSx
Matrix< Scalar, Dynamic, 1 > VectorSx
Definition: IncompleteCholesky.h:56
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:440
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
Eigen::SparseMatrix::rows
Index rows() const
Definition: SparseMatrix.h:138
Eigen::IncompleteCholesky::m_scale
VectorRx m_scale
Definition: IncompleteCholesky.h:173
Eigen::IncompleteCholesky::setInitialShift
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:106
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::IncompleteCholesky::analyzePattern
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:111
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Eigen::SparseSolverBase::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:27