LDLT.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
00006 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00007 // Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com >
00008 //
00009 // This Source Code Form is subject to the terms of the Mozilla
00010 // Public License v. 2.0. If a copy of the MPL was not distributed
00011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00012 
00013 #ifndef EIGEN_LDLT_H
00014 #define EIGEN_LDLT_H
00015 
00016 namespace Eigen { 
00017 
00018 namespace internal {
00019 template<typename MatrixType, int UpLo> struct LDLT_Traits;
00020 }
00021 
00045 template<typename _MatrixType, int _UpLo> class LDLT
00046 {
00047   public:
00048     typedef _MatrixType MatrixType;
00049     enum {
00050       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00051       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00052       Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here!
00053       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00054       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00055       UpLo = _UpLo
00056     };
00057     typedef typename MatrixType::Scalar Scalar;
00058     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00059     typedef typename MatrixType::Index Index;
00060     typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
00061 
00062     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00063     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00064 
00065     typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
00066 
00072     LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
00073 
00080     LDLT(Index size)
00081       : m_matrix(size, size),
00082         m_transpositions(size),
00083         m_temporary(size),
00084         m_isInitialized(false)
00085     {}
00086 
00092     LDLT(const MatrixType& matrix)
00093       : m_matrix(matrix.rows(), matrix.cols()),
00094         m_transpositions(matrix.rows()),
00095         m_temporary(matrix.rows()),
00096         m_isInitialized(false)
00097     {
00098       compute(matrix);
00099     }
00100 
00104     void setZero()
00105     {
00106       m_isInitialized = false;
00107     }
00108 
00110     inline typename Traits::MatrixU matrixU() const
00111     {
00112       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00113       return Traits::getU(m_matrix);
00114     }
00115 
00117     inline typename Traits::MatrixL matrixL() const
00118     {
00119       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00120       return Traits::getL(m_matrix);
00121     }
00122 
00125     inline const TranspositionType& transpositionsP() const
00126     {
00127       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00128       return m_transpositions;
00129     }
00130 
00132     inline Diagonal<const MatrixType> vectorD() const
00133     {
00134       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00135       return m_matrix.diagonal();
00136     }
00137 
00139     inline bool isPositive() const
00140     {
00141       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00142       return m_sign == 1;
00143     }
00144     
00145     #ifdef EIGEN2_SUPPORT
00146     inline bool isPositiveDefinite() const
00147     {
00148       return isPositive();
00149     }
00150     #endif
00151 
00153     inline bool isNegative(void) const
00154     {
00155       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00156       return m_sign == -1;
00157     }
00158 
00174     template<typename Rhs>
00175     inline const internal::solve_retval<LDLT, Rhs>
00176     solve(const MatrixBase<Rhs>& b) const
00177     {
00178       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00179       eigen_assert(m_matrix.rows()==b.rows()
00180                 && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
00181       return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
00182     }
00183 
00184     #ifdef EIGEN2_SUPPORT
00185     template<typename OtherDerived, typename ResultType>
00186     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00187     {
00188       *result = this->solve(b);
00189       return true;
00190     }
00191     #endif
00192 
00193     template<typename Derived>
00194     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
00195 
00196     LDLT& compute(const MatrixType& matrix);
00197 
00198     template <typename Derived>
00199     LDLT& rankUpdate(const MatrixBase<Derived>& w,RealScalar alpha=1);
00200 
00205     inline const MatrixType& matrixLDLT() const
00206     {
00207       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00208       return m_matrix;
00209     }
00210 
00211     MatrixType reconstructedMatrix() const;
00212 
00213     inline Index rows() const { return m_matrix.rows(); }
00214     inline Index cols() const { return m_matrix.cols(); }
00215 
00221     ComputationInfo info() const
00222     {
00223       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00224       return Success;
00225     }
00226 
00227   protected:
00228 
00235     MatrixType m_matrix;
00236     TranspositionType m_transpositions;
00237     TmpMatrixType m_temporary;
00238     int m_sign;
00239     bool m_isInitialized;
00240 };
00241 
00242 namespace internal {
00243 
00244 template<int UpLo> struct ldlt_inplace;
00245 
00246 template<> struct ldlt_inplace<Lower>
00247 {
00248   template<typename MatrixType, typename TranspositionType, typename Workspace>
00249   static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00250   {
00251     typedef typename MatrixType::Scalar Scalar;
00252     typedef typename MatrixType::RealScalar RealScalar;
00253     typedef typename MatrixType::Index Index;
00254     eigen_assert(mat.rows()==mat.cols());
00255     const Index size = mat.rows();
00256 
00257     if (size <= 1)
00258     {
00259       transpositions.setIdentity();
00260       if(sign)
00261         *sign = real(mat.coeff(0,0))>0 ? 1:-1;
00262       return true;
00263     }
00264 
00265     RealScalar cutoff(0), biggest_in_corner;
00266 
00267     for (Index k = 0; k < size; ++k)
00268     {
00269       // Find largest diagonal element
00270       Index index_of_biggest_in_corner;
00271       biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
00272       index_of_biggest_in_corner += k;
00273 
00274       if(k == 0)
00275       {
00276         // The biggest overall is the point of reference to which further diagonals
00277         // are compared; if any diagonal is negligible compared
00278         // to the largest overall, the algorithm bails.
00279         cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
00280 
00281         if(sign)
00282           *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
00283       }
00284 
00285       // Finish early if the matrix is not full rank.
00286       if(biggest_in_corner < cutoff)
00287       {
00288         for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
00289         break;
00290       }
00291 
00292       transpositions.coeffRef(k) = index_of_biggest_in_corner;
00293       if(k != index_of_biggest_in_corner)
00294       {
00295         // apply the transposition while taking care to consider only
00296         // the lower triangular part
00297         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
00298         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
00299         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
00300         std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
00301         for(int i=k+1;i<index_of_biggest_in_corner;++i)
00302         {
00303           Scalar tmp = mat.coeffRef(i,k);
00304           mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
00305           mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
00306         }
00307         if(NumTraits<Scalar>::IsComplex)
00308           mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
00309       }
00310 
00311       // partition the matrix:
00312       //       A00 |  -  |  -
00313       // lu  = A10 | A11 |  -
00314       //       A20 | A21 | A22
00315       Index rs = size - k - 1;
00316       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00317       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00318       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00319 
00320       if(k>0)
00321       {
00322         temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
00323         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
00324         if(rs>0)
00325           A21.noalias() -= A20 * temp.head(k);
00326       }
00327       if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
00328         A21 /= mat.coeffRef(k,k);
00329     }
00330 
00331     return true;
00332   }
00333 
00334   // Reference for the algorithm: Davis and Hager, "Multiple Rank
00335   // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
00336   // Trivial rearrangements of their computations (Timothy E. Holy)
00337   // allow their algorithm to work for rank-1 updates even if the
00338   // original matrix is not of full rank.
00339   // Here only rank-1 updates are implemented, to reduce the
00340   // requirement for intermediate storage and improve accuracy
00341   template<typename MatrixType, typename WDerived>
00342   static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, typename MatrixType::RealScalar sigma=1)
00343   {
00344     using internal::isfinite;
00345     typedef typename MatrixType::Scalar Scalar;
00346     typedef typename MatrixType::RealScalar RealScalar;
00347     typedef typename MatrixType::Index Index;
00348 
00349     const Index size = mat.rows();
00350     eigen_assert(mat.cols() == size && w.size()==size);
00351 
00352     RealScalar alpha = 1;
00353 
00354     // Apply the update
00355     for (Index j = 0; j < size; j++)
00356     {
00357       // Check for termination due to an original decomposition of low-rank
00358       if (!(isfinite)(alpha))
00359         break;
00360 
00361       // Update the diagonal terms
00362       RealScalar dj = real(mat.coeff(j,j));
00363       Scalar wj = w.coeff(j);
00364       RealScalar swj2 = sigma*abs2(wj);
00365       RealScalar gamma = dj*alpha + swj2;
00366 
00367       mat.coeffRef(j,j) += swj2/alpha;
00368       alpha += swj2/dj;
00369 
00370 
00371       // Update the terms of L
00372       Index rs = size-j-1;
00373       w.tail(rs) -= wj * mat.col(j).tail(rs);
00374       if(gamma != 0)
00375         mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
00376     }
00377     return true;
00378   }
00379 
00380   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00381   static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, typename MatrixType::RealScalar sigma=1)
00382   {
00383     // Apply the permutation to the input w
00384     tmp = transpositions * w;
00385 
00386     return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
00387   }
00388 };
00389 
00390 template<> struct ldlt_inplace<Upper>
00391 {
00392   template<typename MatrixType, typename TranspositionType, typename Workspace>
00393   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
00394   {
00395     Transpose<MatrixType> matt(mat);
00396     return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
00397   }
00398 
00399   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00400   static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1)
00401   {
00402     Transpose<MatrixType> matt(mat);
00403     return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
00404   }
00405 };
00406 
00407 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
00408 {
00409   typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
00410   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
00411   static inline MatrixL getL(const MatrixType& m) { return m; }
00412   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00413 };
00414 
00415 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
00416 {
00417   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
00418   typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
00419   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00420   static inline MatrixU getU(const MatrixType& m) { return m; }
00421 };
00422 
00423 } // end namespace internal
00424 
00427 template<typename MatrixType, int _UpLo>
00428 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00429 {
00430   eigen_assert(a.rows()==a.cols());
00431   const Index size = a.rows();
00432 
00433   m_matrix = a;
00434 
00435   m_transpositions.resize(size);
00436   m_isInitialized = false;
00437   m_temporary.resize(size);
00438 
00439   internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
00440 
00441   m_isInitialized = true;
00442   return *this;
00443 }
00444 
00450 template<typename MatrixType, int _UpLo>
00451 template<typename Derived>
00452 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w,typename NumTraits<typename MatrixType::Scalar>::Real sigma)
00453 {
00454   const Index size = w.rows();
00455   if (m_isInitialized)
00456   {
00457     eigen_assert(m_matrix.rows()==size);
00458   }
00459   else
00460   {    
00461     m_matrix.resize(size,size);
00462     m_matrix.setZero();
00463     m_transpositions.resize(size);
00464     for (Index i = 0; i < size; i++)
00465       m_transpositions.coeffRef(i) = i;
00466     m_temporary.resize(size);
00467     m_sign = sigma>=0 ? 1 : -1;
00468     m_isInitialized = true;
00469   }
00470 
00471   internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
00472 
00473   return *this;
00474 }
00475 
00476 namespace internal {
00477 template<typename _MatrixType, int _UpLo, typename Rhs>
00478 struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
00479   : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
00480 {
00481   typedef LDLT<_MatrixType,_UpLo> LDLTType;
00482   EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
00483 
00484   template<typename Dest> void evalTo(Dest& dst) const
00485   {
00486     eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
00487     // dst = P b
00488     dst = dec().transpositionsP() * rhs();
00489 
00490     // dst = L^-1 (P b)
00491     dec().matrixL().solveInPlace(dst);
00492 
00493     // dst = D^-1 (L^-1 P b)
00494     // more precisely, use pseudo-inverse of D (see bug 241)
00495     using std::abs;
00496     using std::max;
00497     typedef typename LDLTType::MatrixType MatrixType;
00498     typedef typename LDLTType::Scalar Scalar;
00499     typedef typename LDLTType::RealScalar RealScalar;
00500     const Diagonal<const MatrixType> vectorD = dec().vectorD();
00501     RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
00502                                  RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
00503     for (Index i = 0; i < vectorD.size(); ++i) {
00504       if(abs(vectorD(i)) > tolerance)
00505         dst.row(i) /= vectorD(i);
00506       else
00507         dst.row(i).setZero();
00508     }
00509 
00510     // dst = L^-T (D^-1 L^-1 P b)
00511     dec().matrixU().solveInPlace(dst);
00512 
00513     // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
00514     dst = dec().transpositionsP().transpose() * dst;
00515   }
00516 };
00517 }
00518 
00532 template<typename MatrixType,int _UpLo>
00533 template<typename Derived>
00534 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00535 {
00536   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00537   const Index size = m_matrix.rows();
00538   eigen_assert(size == bAndX.rows());
00539 
00540   bAndX = this->solve(bAndX);
00541 
00542   return true;
00543 }
00544 
00548 template<typename MatrixType, int _UpLo>
00549 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
00550 {
00551   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00552   const Index size = m_matrix.rows();
00553   MatrixType res(size,size);
00554 
00555   // P
00556   res.setIdentity();
00557   res = transpositionsP() * res;
00558   // L^* P
00559   res = matrixU() * res;
00560   // D(L^*P)
00561   res = vectorD().asDiagonal() * res;
00562   // L(DL^*P)
00563   res = matrixL() * res;
00564   // P^T (LDL^*P)
00565   res = transpositionsP().transpose() * res;
00566 
00567   return res;
00568 }
00569 
00573 template<typename MatrixType, unsigned int UpLo>
00574 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00575 SelfAdjointView<MatrixType, UpLo>::ldlt() const
00576 {
00577   return LDLT<PlainObject,UpLo>(m_matrix);
00578 }
00579 
00583 template<typename Derived>
00584 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
00585 MatrixBase<Derived>::ldlt() const
00586 {
00587   return LDLT<PlainObject>(derived());
00588 }
00589 
00590 } // end namespace Eigen
00591 
00592 #endif // EIGEN_LDLT_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:08