LDLT.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) 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 // Eigen is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 3 of the License, or (at your option) any later version.
13 //
14 // Alternatively, you can redistribute it and/or
15 // modify it under the terms of the GNU General Public License as
16 // published by the Free Software Foundation; either version 2 of
17 // the License, or (at your option) any later version.
18 //
19 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
20 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU Lesser General Public
25 // License and a copy of the GNU General Public License along with
26 // Eigen. If not, see <http://www.gnu.org/licenses/>.
27 
28 #ifndef EIGEN_LDLT_H
29 #define EIGEN_LDLT_H
30 
31 namespace Eigen {
32 
33 namespace internal {
34 template<typename MatrixType, int UpLo> struct LDLT_Traits;
35 }
36 
60 template<typename _MatrixType, int _UpLo> class LDLT
61 {
62  public:
63  typedef _MatrixType MatrixType;
64  enum {
67  Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
70  UpLo = _UpLo
71  };
72  typedef typename MatrixType::Scalar Scalar;
74  typedef typename MatrixType::Index Index;
76 
79 
80  typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
81 
88 
95  LDLT(Index size)
96  : m_matrix(size, size),
97  m_transpositions(size),
98  m_temporary(size),
99  m_isInitialized(false)
100  {}
101 
107  LDLT(const MatrixType& matrix)
108  : m_matrix(matrix.rows(), matrix.cols()),
109  m_transpositions(matrix.rows()),
110  m_temporary(matrix.rows()),
111  m_isInitialized(false)
112  {
113  compute(matrix);
114  }
115 
119  void setZero()
120  {
121  m_isInitialized = false;
122  }
123 
125  inline typename Traits::MatrixU matrixU() const
126  {
127  eigen_assert(m_isInitialized && "LDLT is not initialized.");
128  return Traits::getU(m_matrix);
129  }
130 
132  inline typename Traits::MatrixL matrixL() const
133  {
134  eigen_assert(m_isInitialized && "LDLT is not initialized.");
135  return Traits::getL(m_matrix);
136  }
137 
140  inline const TranspositionType& transpositionsP() const
141  {
142  eigen_assert(m_isInitialized && "LDLT is not initialized.");
143  return m_transpositions;
144  }
145 
148  {
149  eigen_assert(m_isInitialized && "LDLT is not initialized.");
150  return m_matrix.diagonal();
151  }
152 
154  inline bool isPositive() const
155  {
156  eigen_assert(m_isInitialized && "LDLT is not initialized.");
157  return m_sign == 1;
158  }
159 
160  #ifdef EIGEN2_SUPPORT
161  inline bool isPositiveDefinite() const
162  {
163  return isPositive();
164  }
165  #endif
166 
168  inline bool isNegative(void) const
169  {
170  eigen_assert(m_isInitialized && "LDLT is not initialized.");
171  return m_sign == -1;
172  }
173 
189  template<typename Rhs>
190  inline const internal::solve_retval<LDLT, Rhs>
191  solve(const MatrixBase<Rhs>& b) const
192  {
193  eigen_assert(m_isInitialized && "LDLT is not initialized.");
194  eigen_assert(m_matrix.rows()==b.rows()
195  && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
196  return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
197  }
198 
199  #ifdef EIGEN2_SUPPORT
200  template<typename OtherDerived, typename ResultType>
201  bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
202  {
203  *result = this->solve(b);
204  return true;
205  }
206  #endif
207 
208  template<typename Derived>
209  bool solveInPlace(MatrixBase<Derived> &bAndX) const;
210 
211  LDLT& compute(const MatrixType& matrix);
212 
213  template <typename Derived>
214  LDLT& rankUpdate(const MatrixBase<Derived>& w,RealScalar alpha=1);
215 
220  inline const MatrixType& matrixLDLT() const
221  {
222  eigen_assert(m_isInitialized && "LDLT is not initialized.");
223  return m_matrix;
224  }
225 
227 
228  inline Index rows() const { return m_matrix.rows(); }
229  inline Index cols() const { return m_matrix.cols(); }
230 
237  {
238  eigen_assert(m_isInitialized && "LDLT is not initialized.");
239  return Success;
240  }
241 
242  protected:
243 
253  int m_sign;
255 };
256 
257 namespace internal {
258 
259 template<int UpLo> struct ldlt_inplace;
260 
261 template<> struct ldlt_inplace<Lower>
262 {
263  template<typename MatrixType, typename TranspositionType, typename Workspace>
264  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
265  {
266  typedef typename MatrixType::Scalar Scalar;
267  typedef typename MatrixType::RealScalar RealScalar;
268  typedef typename MatrixType::Index Index;
269  eigen_assert(mat.rows()==mat.cols());
270  const Index size = mat.rows();
271 
272  if (size <= 1)
273  {
274  transpositions.setIdentity();
275  if(sign)
276  *sign = real(mat.coeff(0,0))>0 ? 1:-1;
277  return true;
278  }
279 
280  RealScalar cutoff(0), biggest_in_corner;
281 
282  for (Index k = 0; k < size; ++k)
283  {
284  // Find largest diagonal element
285  Index index_of_biggest_in_corner;
286  biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
287  index_of_biggest_in_corner += k;
288 
289  if(k == 0)
290  {
291  // The biggest overall is the point of reference to which further diagonals
292  // are compared; if any diagonal is negligible compared
293  // to the largest overall, the algorithm bails.
294  cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
295 
296  if(sign)
297  *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
298  }
299 
300  // Finish early if the matrix is not full rank.
301  if(biggest_in_corner < cutoff)
302  {
303  for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
304  break;
305  }
306 
307  transpositions.coeffRef(k) = index_of_biggest_in_corner;
308  if(k != index_of_biggest_in_corner)
309  {
310  // apply the transposition while taking care to consider only
311  // the lower triangular part
312  Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
313  mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
314  mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
315  std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
316  for(int i=k+1;i<index_of_biggest_in_corner;++i)
317  {
318  Scalar tmp = mat.coeffRef(i,k);
319  mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
320  mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
321  }
323  mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
324  }
325 
326  // partition the matrix:
327  // A00 | - | -
328  // lu = A10 | A11 | -
329  // A20 | A21 | A22
330  Index rs = size - k - 1;
331  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
332  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
333  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
334 
335  if(k>0)
336  {
337  temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
338  mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
339  if(rs>0)
340  A21.noalias() -= A20 * temp.head(k);
341  }
342  if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
343  A21 /= mat.coeffRef(k,k);
344  }
345 
346  return true;
347  }
348 
349  // Reference for the algorithm: Davis and Hager, "Multiple Rank
350  // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
351  // Trivial rearrangements of their computations (Timothy E. Holy)
352  // allow their algorithm to work for rank-1 updates even if the
353  // original matrix is not of full rank.
354  // Here only rank-1 updates are implemented, to reduce the
355  // requirement for intermediate storage and improve accuracy
356  template<typename MatrixType, typename WDerived>
357  static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, typename MatrixType::RealScalar sigma=1)
358  {
359  using internal::isfinite;
360  typedef typename MatrixType::Scalar Scalar;
361  typedef typename MatrixType::RealScalar RealScalar;
362  typedef typename MatrixType::Index Index;
363 
364  const Index size = mat.rows();
365  eigen_assert(mat.cols() == size && w.size()==size);
366 
367  RealScalar alpha = 1;
368 
369  // Apply the update
370  for (Index j = 0; j < size; j++)
371  {
372  // Check for termination due to an original decomposition of low-rank
373  if (!isfinite(alpha))
374  break;
375 
376  // Update the diagonal terms
377  RealScalar dj = real(mat.coeff(j,j));
378  Scalar wj = w.coeff(j);
379  RealScalar swj2 = sigma*abs2(wj);
380  RealScalar gamma = dj*alpha + swj2;
381 
382  mat.coeffRef(j,j) += swj2/alpha;
383  alpha += swj2/dj;
384 
385 
386  // Update the terms of L
387  Index rs = size-j-1;
388  w.tail(rs) -= wj * mat.col(j).tail(rs);
389  if(gamma != 0)
390  mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
391  }
392  return true;
393  }
394 
395  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
396  static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, typename MatrixType::RealScalar sigma=1)
397  {
398  // Apply the permutation to the input w
399  tmp = transpositions * w;
400 
401  return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
402  }
403 };
404 
405 template<> struct ldlt_inplace<Upper>
406 {
407  template<typename MatrixType, typename TranspositionType, typename Workspace>
408  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
409  {
410  Transpose<MatrixType> matt(mat);
411  return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
412  }
413 
414  template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
415  static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1)
416  {
417  Transpose<MatrixType> matt(mat);
418  return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
419  }
420 };
421 
422 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
423 {
424  typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
425  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
426  static inline MatrixL getL(const MatrixType& m) { return m; }
427  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
428 };
429 
430 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
431 {
432  typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
433  typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
434  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
435  static inline MatrixU getU(const MatrixType& m) { return m; }
436 };
437 
438 } // end namespace internal
439 
442 template<typename MatrixType, int _UpLo>
444 {
445  eigen_assert(a.rows()==a.cols());
446  const Index size = a.rows();
447 
448  m_matrix = a;
449 
450  m_transpositions.resize(size);
451  m_isInitialized = false;
452  m_temporary.resize(size);
453 
454  internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
455 
456  m_isInitialized = true;
457  return *this;
458 }
459 
465 template<typename MatrixType, int _UpLo>
466 template<typename Derived>
468 {
469  const Index size = w.rows();
470  if (m_isInitialized)
471  {
472  eigen_assert(m_matrix.rows()==size);
473  }
474  else
475  {
476  m_matrix.resize(size,size);
477  m_matrix.setZero();
478  m_transpositions.resize(size);
479  for (Index i = 0; i < size; i++)
480  m_transpositions.coeffRef(i) = i;
481  m_temporary.resize(size);
482  m_sign = sigma>=0 ? 1 : -1;
483  m_isInitialized = true;
484  }
485 
486  internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
487 
488  return *this;
489 }
490 
491 namespace internal {
492 template<typename _MatrixType, int _UpLo, typename Rhs>
493 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
494  : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
495 {
496  typedef LDLT<_MatrixType,_UpLo> LDLTType;
497  EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
498 
499  template<typename Dest> void evalTo(Dest& dst) const
500  {
501  eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
502  // dst = P b
503  dst = dec().transpositionsP() * rhs();
504 
505  // dst = L^-1 (P b)
506  dec().matrixL().solveInPlace(dst);
507 
508  // dst = D^-1 (L^-1 P b)
509  // more precisely, use pseudo-inverse of D (see bug 241)
510  using std::abs;
511  using std::max;
512  typedef typename LDLTType::MatrixType MatrixType;
513  typedef typename LDLTType::Scalar Scalar;
514  typedef typename LDLTType::RealScalar RealScalar;
515  const Diagonal<const MatrixType> vectorD = dec().vectorD();
516  RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
517  RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
518  for (Index i = 0; i < vectorD.size(); ++i) {
519  if(abs(vectorD(i)) > tolerance)
520  dst.row(i) /= vectorD(i);
521  else
522  dst.row(i).setZero();
523  }
524 
525  // dst = L^-T (D^-1 L^-1 P b)
526  dec().matrixU().solveInPlace(dst);
527 
528  // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
529  dst = dec().transpositionsP().transpose() * dst;
530  }
531 };
532 }
533 
547 template<typename MatrixType,int _UpLo>
548 template<typename Derived>
550 {
551  eigen_assert(m_isInitialized && "LDLT is not initialized.");
552  const Index size = m_matrix.rows();
553  eigen_assert(size == bAndX.rows());
554 
555  bAndX = this->solve(bAndX);
556 
557  return true;
558 }
559 
563 template<typename MatrixType, int _UpLo>
565 {
566  eigen_assert(m_isInitialized && "LDLT is not initialized.");
567  const Index size = m_matrix.rows();
568  MatrixType res(size,size);
569 
570  // P
571  res.setIdentity();
572  res = transpositionsP() * res;
573  // L^* P
574  res = matrixU() * res;
575  // D(L^*P)
576  res = vectorD().asDiagonal() * res;
577  // L(DL^*P)
578  res = matrixL() * res;
579  // P^T (LDL^*P)
580  res = transpositionsP().transpose() * res;
581 
582  return res;
583 }
584 
588 template<typename MatrixType, unsigned int UpLo>
591 {
592  return LDLT<PlainObject,UpLo>(m_matrix);
593 }
594 
598 template<typename Derived>
601 {
602  return LDLT<PlainObject>(derived());
603 }
604 
605 } // end namespace Eigen
606 
607 #endif // EIGEN_LDLT_H