Eigen  3.2.1
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LDLT.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
6 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
7 // Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com >
8 //
9 // This Source Code Form is subject to the terms of the Mozilla
10 // Public License v. 2.0. If a copy of the MPL was not distributed
11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
12 
13 #ifndef EIGEN_LDLT_H
14 #define EIGEN_LDLT_H
15 
16 namespace Eigen {
17 
18 namespace internal {
19  template<typename MatrixType, int UpLo> struct LDLT_Traits;
20 
21  // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
22  enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
23 }
24 
48 template<typename _MatrixType, int _UpLo> class LDLT
49 {
50  public:
51  typedef _MatrixType MatrixType;
52  enum {
53  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55  Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
56  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
58  UpLo = _UpLo
59  };
60  typedef typename MatrixType::Scalar Scalar;
61  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
62  typedef typename MatrixType::Index Index;
64 
67 
68  typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
69 
75  LDLT()
76  : m_matrix(),
77  m_transpositions(),
78  m_sign(internal::ZeroSign),
79  m_isInitialized(false)
80  {}
81 
88  LDLT(Index size)
89  : m_matrix(size, size),
90  m_transpositions(size),
91  m_temporary(size),
92  m_sign(internal::ZeroSign),
93  m_isInitialized(false)
94  {}
95 
101  LDLT(const MatrixType& matrix)
102  : m_matrix(matrix.rows(), matrix.cols()),
103  m_transpositions(matrix.rows()),
104  m_temporary(matrix.rows()),
105  m_sign(internal::ZeroSign),
106  m_isInitialized(false)
107  {
108  compute(matrix);
109  }
110 
114  void setZero()
115  {
116  m_isInitialized = false;
117  }
118 
120  inline typename Traits::MatrixU matrixU() const
121  {
122  eigen_assert(m_isInitialized && "LDLT is not initialized.");
123  return Traits::getU(m_matrix);
124  }
125 
127  inline typename Traits::MatrixL matrixL() const
128  {
129  eigen_assert(m_isInitialized && "LDLT is not initialized.");
130  return Traits::getL(m_matrix);
131  }
132 
135  inline const TranspositionType& transpositionsP() const
136  {
137  eigen_assert(m_isInitialized && "LDLT is not initialized.");
138  return m_transpositions;
139  }
140 
143  {
144  eigen_assert(m_isInitialized && "LDLT is not initialized.");
145  return m_matrix.diagonal();
146  }
147 
149  inline bool isPositive() const
150  {
151  eigen_assert(m_isInitialized && "LDLT is not initialized.");
152  return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
153  }
154 
155  #ifdef EIGEN2_SUPPORT
156  inline bool isPositiveDefinite() const
157  {
158  return isPositive();
159  }
160  #endif
161 
163  inline bool isNegative(void) const
164  {
165  eigen_assert(m_isInitialized && "LDLT is not initialized.");
166  return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
167  }
168 
184  template<typename Rhs>
185  inline const internal::solve_retval<LDLT, Rhs>
186  solve(const MatrixBase<Rhs>& b) const
187  {
188  eigen_assert(m_isInitialized && "LDLT is not initialized.");
189  eigen_assert(m_matrix.rows()==b.rows()
190  && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
191  return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
192  }
193 
194  #ifdef EIGEN2_SUPPORT
195  template<typename OtherDerived, typename ResultType>
196  bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
197  {
198  *result = this->solve(b);
199  return true;
200  }
201  #endif
202 
203  template<typename Derived>
204  bool solveInPlace(MatrixBase<Derived> &bAndX) const;
205 
206  LDLT& compute(const MatrixType& matrix);
207 
208  template <typename Derived>
209  LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1);
210 
215  inline const MatrixType& matrixLDLT() const
216  {
217  eigen_assert(m_isInitialized && "LDLT is not initialized.");
218  return m_matrix;
219  }
220 
221  MatrixType reconstructedMatrix() const;
222 
223  inline Index rows() const { return m_matrix.rows(); }
224  inline Index cols() const { return m_matrix.cols(); }
225 
232  {
233  eigen_assert(m_isInitialized && "LDLT is not initialized.");
234  return Success;
235  }
236 
237  protected:
238 
245  MatrixType m_matrix;
246  TranspositionType m_transpositions;
247  TmpMatrixType m_temporary;
248  internal::SignMatrix m_sign;
249  bool m_isInitialized;
250 };
251 
252 namespace internal {
253 
254 template<int UpLo> struct ldlt_inplace;
255 
256 template<> struct ldlt_inplace<Lower>
257 {
258  template<typename MatrixType, typename TranspositionType, typename Workspace>
259  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
260  {
261  using std::abs;
262  typedef typename MatrixType::Scalar Scalar;
263  typedef typename MatrixType::RealScalar RealScalar;
264  typedef typename MatrixType::Index Index;
265  eigen_assert(mat.rows()==mat.cols());
266  const Index size = mat.rows();
267 
268  if (size <= 1)
269  {
270  transpositions.setIdentity();
271  if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef;
272  else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef;
273  else sign = ZeroSign;
274  return true;
275  }
276 
277  RealScalar cutoff(0), biggest_in_corner;
278 
279  for (Index k = 0; k < size; ++k)
280  {
281  // Find largest diagonal element
282  Index index_of_biggest_in_corner;
283  biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
284  index_of_biggest_in_corner += k;
285 
286  if(k == 0)
287  {
288  // The biggest overall is the point of reference to which further diagonals
289  // are compared; if any diagonal is negligible compared
290  // to the largest overall, the algorithm bails.
291  cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
292  }
293 
294  // Finish early if the matrix is not full rank.
295  if(biggest_in_corner < cutoff)
296  {
297  for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
298  break;
299  }
300 
301  transpositions.coeffRef(k) = index_of_biggest_in_corner;
302  if(k != index_of_biggest_in_corner)
303  {
304  // apply the transposition while taking care to consider only
305  // the lower triangular part
306  Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
307  mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
308  mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
309  std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
310  for(int i=k+1;i<index_of_biggest_in_corner;++i)
311  {
312  Scalar tmp = mat.coeffRef(i,k);
313  mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i));
314  mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp);
315  }
316  if(NumTraits<Scalar>::IsComplex)
317  mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k));
318  }
319 
320  // partition the matrix:
321  // A00 | - | -
322  // lu = A10 | A11 | -
323  // A20 | A21 | A22
324  Index rs = size - k - 1;
325  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
326  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
327  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
328 
329  if(k>0)
330  {
331  temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
332  mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
333  if(rs>0)
334  A21.noalias() -= A20 * temp.head(k);
335  }
336  if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
337  A21 /= mat.coeffRef(k,k);
338 
339  RealScalar realAkk = numext::real(mat.coeffRef(k,k));
340  if (sign == PositiveSemiDef) {
341  if (realAkk < 0) sign = Indefinite;
342  } else if (sign == NegativeSemiDef) {
343  if (realAkk > 0) sign = Indefinite;
344  } else if (sign == ZeroSign) {
345  if (realAkk > 0) sign = PositiveSemiDef;
346  else if (realAkk < 0) sign = NegativeSemiDef;
347  }
348  }
349 
350  return true;
351  }
352 
353  // Reference for the algorithm: Davis and Hager, "Multiple Rank
354  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
355  // Trivial rearrangements of their computations (Timothy E. Holy)
356  // allow their algorithm to work for rank-1 updates even if the
357  // original matrix is not of full rank.
358  // Here only rank-1 updates are implemented, to reduce the
359  // requirement for intermediate storage and improve accuracy
360  template<typename MatrixType, typename WDerived>
361  static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1)
362  {
363  using numext::isfinite;
364  typedef typename MatrixType::Scalar Scalar;
365  typedef typename MatrixType::RealScalar RealScalar;
366  typedef typename MatrixType::Index Index;
367 
368  const Index size = mat.rows();
369  eigen_assert(mat.cols() == size && w.size()==size);
370 
371  RealScalar alpha = 1;
372 
373  // Apply the update
374  for (Index j = 0; j < size; j++)
375  {
376  // Check for termination due to an original decomposition of low-rank
377  if (!(isfinite)(alpha))
378  break;
379 
380  // Update the diagonal terms
381  RealScalar dj = numext::real(mat.coeff(j,j));
382  Scalar wj = w.coeff(j);
383  RealScalar swj2 = sigma*numext::abs2(wj);
384  RealScalar gamma = dj*alpha + swj2;
385 
386  mat.coeffRef(j,j) += swj2/alpha;
387  alpha += swj2/dj;
388 
389 
390  // Update the terms of L
391  Index rs = size-j-1;
392  w.tail(rs) -= wj * mat.col(j).tail(rs);
393  if(gamma != 0)
394  mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
395  }
396  return true;
397  }
398 
399  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
400  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
401  {
402  // Apply the permutation to the input w
403  tmp = transpositions * w;
404 
405  return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
406  }
407 };
408 
409 template<> struct ldlt_inplace<Upper>
410 {
411  template<typename MatrixType, typename TranspositionType, typename Workspace>
412  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
413  {
414  Transpose<MatrixType> matt(mat);
415  return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
416  }
417 
418  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
419  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1)
420  {
421  Transpose<MatrixType> matt(mat);
422  return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
423  }
424 };
425 
426 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
427 {
428  typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
429  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
430  static inline MatrixL getL(const MatrixType& m) { return m; }
431  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
432 };
433 
434 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
435 {
436  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
437  typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
438  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
439  static inline MatrixU getU(const MatrixType& m) { return m; }
440 };
441 
442 } // end namespace internal
443 
446 template<typename MatrixType, int _UpLo>
448 {
449  eigen_assert(a.rows()==a.cols());
450  const Index size = a.rows();
451 
452  m_matrix = a;
453 
454  m_transpositions.resize(size);
455  m_isInitialized = false;
456  m_temporary.resize(size);
457 
458  internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
459 
460  m_isInitialized = true;
461  return *this;
462 }
463 
469 template<typename MatrixType, int _UpLo>
470 template<typename Derived>
472 {
473  const Index size = w.rows();
474  if (m_isInitialized)
475  {
476  eigen_assert(m_matrix.rows()==size);
477  }
478  else
479  {
480  m_matrix.resize(size,size);
481  m_matrix.setZero();
482  m_transpositions.resize(size);
483  for (Index i = 0; i < size; i++)
484  m_transpositions.coeffRef(i) = i;
485  m_temporary.resize(size);
486  m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
487  m_isInitialized = true;
488  }
489 
490  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
491 
492  return *this;
493 }
494 
495 namespace internal {
496 template<typename _MatrixType, int _UpLo, typename Rhs>
497 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
498  : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
499 {
500  typedef LDLT<_MatrixType,_UpLo> LDLTType;
501  EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
502 
503  template<typename Dest> void evalTo(Dest& dst) const
504  {
505  eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
506  // dst = P b
507  dst = dec().transpositionsP() * rhs();
508 
509  // dst = L^-1 (P b)
510  dec().matrixL().solveInPlace(dst);
511 
512  // dst = D^-1 (L^-1 P b)
513  // more precisely, use pseudo-inverse of D (see bug 241)
514  using std::abs;
515  using std::max;
516  typedef typename LDLTType::MatrixType MatrixType;
517  typedef typename LDLTType::Scalar Scalar;
518  typedef typename LDLTType::RealScalar RealScalar;
519  const Diagonal<const MatrixType> vectorD = dec().vectorD();
520  RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
521  RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
522  for (Index i = 0; i < vectorD.size(); ++i) {
523  if(abs(vectorD(i)) > tolerance)
524  dst.row(i) /= vectorD(i);
525  else
526  dst.row(i).setZero();
527  }
528 
529  // dst = L^-T (D^-1 L^-1 P b)
530  dec().matrixU().solveInPlace(dst);
531 
532  // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
533  dst = dec().transpositionsP().transpose() * dst;
534  }
535 };
536 }
537 
551 template<typename MatrixType,int _UpLo>
552 template<typename Derived>
553 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
554 {
555  eigen_assert(m_isInitialized && "LDLT is not initialized.");
556  eigen_assert(m_matrix.rows() == bAndX.rows());
557 
558  bAndX = this->solve(bAndX);
559 
560  return true;
561 }
562 
566 template<typename MatrixType, int _UpLo>
568 {
569  eigen_assert(m_isInitialized && "LDLT is not initialized.");
570  const Index size = m_matrix.rows();
571  MatrixType res(size,size);
572 
573  // P
574  res.setIdentity();
575  res = transpositionsP() * res;
576  // L^* P
577  res = matrixU() * res;
578  // D(L^*P)
579  res = vectorD().asDiagonal() * res;
580  // L(DL^*P)
581  res = matrixL() * res;
582  // P^T (LDL^*P)
583  res = transpositionsP().transpose() * res;
584 
585  return res;
586 }
587 
591 template<typename MatrixType, unsigned int UpLo>
594 {
595  return LDLT<PlainObject,UpLo>(m_matrix);
596 }
597 
601 template<typename Derived>
604 {
605  return LDLT<PlainObject>(derived());
606 }
607 
608 } // end namespace Eigen
609 
610 #endif // EIGEN_LDLT_H
Robust Cholesky decomposition of a matrix with pivoting.
Definition: LDLT.h:48
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LDLT.h:231
LDLT(const MatrixType &matrix)
Constructor with decomposition.
Definition: LDLT.h:101
MatrixType reconstructedMatrix() const
Definition: LDLT.h:567
const TranspositionType & transpositionsP() const
Definition: LDLT.h:135
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Traits::MatrixL matrixL() const
Definition: LDLT.h:127
const internal::solve_retval< LDLT, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: LDLT.h:186
bool isPositive() const
Definition: LDLT.h:149
const LDLT< PlainObject, UpLo > ldlt() const
Definition: LDLT.h:593
Definition: Constants.h:169
Definition: Constants.h:167
LDLT & compute(const MatrixType &matrix)
Definition: LDLT.h:447
LDLT(Index size)
Default Constructor with memory preallocation.
Definition: LDLT.h:88
void setZero()
Definition: LDLT.h:114
Diagonal< const MatrixType > vectorD() const
Definition: LDLT.h:142
LDLT()
Default Constructor.
Definition: LDLT.h:75
Definition: Constants.h:376
bool isNegative(void) const
Definition: LDLT.h:163
const unsigned int RowMajorBit
Definition: Constants.h:53
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:64
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const LDLT< PlainObject > ldlt() const
Definition: LDLT.h:603
Traits::MatrixU matrixU() const
Definition: LDLT.h:120
const MatrixType & matrixLDLT() const
Definition: LDLT.h:215