Eigen  3.3.3
LDLT.h
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   // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef
00022   enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite };
00023 }
00024 
00050 template<typename _MatrixType, int _UpLo> class LDLT
00051 {
00052   public:
00053     typedef _MatrixType MatrixType;
00054     enum {
00055       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00056       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00057       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00058       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00059       UpLo = _UpLo
00060     };
00061     typedef typename MatrixType::Scalar Scalar;
00062     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00063     typedef Eigen::Index Index; 
00064     typedef typename MatrixType::StorageIndex StorageIndex;
00065     typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType;
00066 
00067     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00068     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00069 
00070     typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
00071 
00077     LDLT()
00078       : m_matrix(),
00079         m_transpositions(),
00080         m_sign(internal::ZeroSign),
00081         m_isInitialized(false)
00082     {}
00083 
00090     explicit LDLT(Index size)
00091       : m_matrix(size, size),
00092         m_transpositions(size),
00093         m_temporary(size),
00094         m_sign(internal::ZeroSign),
00095         m_isInitialized(false)
00096     {}
00097 
00104     template<typename InputType>
00105     explicit LDLT(const EigenBase<InputType>& matrix)
00106       : m_matrix(matrix.rows(), matrix.cols()),
00107         m_transpositions(matrix.rows()),
00108         m_temporary(matrix.rows()),
00109         m_sign(internal::ZeroSign),
00110         m_isInitialized(false)
00111     {
00112       compute(matrix.derived());
00113     }
00114 
00121     template<typename InputType>
00122     explicit LDLT(EigenBase<InputType>& matrix)
00123       : m_matrix(matrix.derived()),
00124         m_transpositions(matrix.rows()),
00125         m_temporary(matrix.rows()),
00126         m_sign(internal::ZeroSign),
00127         m_isInitialized(false)
00128     {
00129       compute(matrix.derived());
00130     }
00131 
00135     void setZero()
00136     {
00137       m_isInitialized = false;
00138     }
00139 
00141     inline typename Traits::MatrixU matrixU() const
00142     {
00143       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00144       return Traits::getU(m_matrix);
00145     }
00146 
00148     inline typename Traits::MatrixL matrixL() const
00149     {
00150       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00151       return Traits::getL(m_matrix);
00152     }
00153 
00156     inline const TranspositionType& transpositionsP() const
00157     {
00158       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00159       return m_transpositions;
00160     }
00161 
00163     inline Diagonal<const MatrixType> vectorD() const
00164     {
00165       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00166       return m_matrix.diagonal();
00167     }
00168 
00170     inline bool isPositive() const
00171     {
00172       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00173       return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign;
00174     }
00175 
00177     inline bool isNegative(void) const
00178     {
00179       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00180       return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign;
00181     }
00182 
00198     template<typename Rhs>
00199     inline const Solve<LDLT, Rhs>
00200     solve(const MatrixBase<Rhs>& b) const
00201     {
00202       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00203       eigen_assert(m_matrix.rows()==b.rows()
00204                 && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
00205       return Solve<LDLT, Rhs>(*this, b.derived());
00206     }
00207 
00208     template<typename Derived>
00209     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
00210 
00211     template<typename InputType>
00212     LDLT& compute(const EigenBase<InputType>& matrix);
00213 
00217     RealScalar rcond() const
00218     {
00219       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00220       return internal::rcond_estimate_helper(m_l1_norm, *this);
00221     }
00222 
00223     template <typename Derived>
00224     LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1);
00225 
00230     inline const MatrixType& matrixLDLT() const
00231     {
00232       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00233       return m_matrix;
00234     }
00235 
00236     MatrixType reconstructedMatrix() const;
00237 
00243     const LDLT& adjoint() const { return *this; };
00244 
00245     inline Index rows() const { return m_matrix.rows(); }
00246     inline Index cols() const { return m_matrix.cols(); }
00247 
00253     ComputationInfo info() const
00254     {
00255       eigen_assert(m_isInitialized && "LDLT is not initialized.");
00256       return m_info;
00257     }
00258 
00259     #ifndef EIGEN_PARSED_BY_DOXYGEN
00260     template<typename RhsType, typename DstType>
00261     EIGEN_DEVICE_FUNC
00262     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00263     #endif
00264 
00265   protected:
00266 
00267     static void check_template_parameters()
00268     {
00269       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00270     }
00271 
00278     MatrixType m_matrix;
00279     RealScalar m_l1_norm;
00280     TranspositionType m_transpositions;
00281     TmpMatrixType m_temporary;
00282     internal::SignMatrix m_sign;
00283     bool m_isInitialized;
00284     ComputationInfo m_info;
00285 };
00286 
00287 namespace internal {
00288 
00289 template<int UpLo> struct ldlt_inplace;
00290 
00291 template<> struct ldlt_inplace<Lower>
00292 {
00293   template<typename MatrixType, typename TranspositionType, typename Workspace>
00294   static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
00295   {
00296     using std::abs;
00297     typedef typename MatrixType::Scalar Scalar;
00298     typedef typename MatrixType::RealScalar RealScalar;
00299     typedef typename TranspositionType::StorageIndex IndexType;
00300     eigen_assert(mat.rows()==mat.cols());
00301     const Index size = mat.rows();
00302     bool found_zero_pivot = false;
00303     bool ret = true;
00304 
00305     if (size <= 1)
00306     {
00307       transpositions.setIdentity();
00308       if (numext::real(mat.coeff(0,0)) > static_cast<RealScalar>(0) ) sign = PositiveSemiDef;
00309       else if (numext::real(mat.coeff(0,0)) < static_cast<RealScalar>(0)) sign = NegativeSemiDef;
00310       else sign = ZeroSign;
00311       return true;
00312     }
00313 
00314     for (Index k = 0; k < size; ++k)
00315     {
00316       // Find largest diagonal element
00317       Index index_of_biggest_in_corner;
00318       mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
00319       index_of_biggest_in_corner += k;
00320 
00321       transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner);
00322       if(k != index_of_biggest_in_corner)
00323       {
00324         // apply the transposition while taking care to consider only
00325         // the lower triangular part
00326         Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
00327         mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
00328         mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
00329         std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
00330         for(Index i=k+1;i<index_of_biggest_in_corner;++i)
00331         {
00332           Scalar tmp = mat.coeffRef(i,k);
00333           mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i));
00334           mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp);
00335         }
00336         if(NumTraits<Scalar>::IsComplex)
00337           mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k));
00338       }
00339 
00340       // partition the matrix:
00341       //       A00 |  -  |  -
00342       // lu  = A10 | A11 |  -
00343       //       A20 | A21 | A22
00344       Index rs = size - k - 1;
00345       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00346       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00347       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00348 
00349       if(k>0)
00350       {
00351         temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
00352         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
00353         if(rs>0)
00354           A21.noalias() -= A20 * temp.head(k);
00355       }
00356 
00357       // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot
00358       // was smaller than the cutoff value. However, since LDLT is not rank-revealing
00359       // we should only make sure that we do not introduce INF or NaN values.
00360       // Remark that LAPACK also uses 0 as the cutoff value.
00361       RealScalar realAkk = numext::real(mat.coeffRef(k,k));
00362       bool pivot_is_valid = (abs(realAkk) > RealScalar(0));
00363 
00364       if(k==0 && !pivot_is_valid)
00365       {
00366         // The entire diagonal is zero, there is nothing more to do
00367         // except filling the transpositions, and checking whether the matrix is zero.
00368         sign = ZeroSign;
00369         for(Index j = 0; j<size; ++j)
00370         {
00371           transpositions.coeffRef(j) = IndexType(j);
00372           ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).all();
00373         }
00374         return ret;
00375       }
00376 
00377       if((rs>0) && pivot_is_valid)
00378         A21 /= realAkk;
00379 
00380       if(found_zero_pivot && pivot_is_valid) ret = false; // factorization failed
00381       else if(!pivot_is_valid) found_zero_pivot = true;
00382 
00383       if (sign == PositiveSemiDef) {
00384         if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite;
00385       } else if (sign == NegativeSemiDef) {
00386         if (realAkk > static_cast<RealScalar>(0)) sign = Indefinite;
00387       } else if (sign == ZeroSign) {
00388         if (realAkk > static_cast<RealScalar>(0)) sign = PositiveSemiDef;
00389         else if (realAkk < static_cast<RealScalar>(0)) sign = NegativeSemiDef;
00390       }
00391     }
00392 
00393     return ret;
00394   }
00395 
00396   // Reference for the algorithm: Davis and Hager, "Multiple Rank
00397   // Modifications of a Sparse Cholesky Factorization" (Algorithm 1)
00398   // Trivial rearrangements of their computations (Timothy E. Holy)
00399   // allow their algorithm to work for rank-1 updates even if the
00400   // original matrix is not of full rank.
00401   // Here only rank-1 updates are implemented, to reduce the
00402   // requirement for intermediate storage and improve accuracy
00403   template<typename MatrixType, typename WDerived>
00404   static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1)
00405   {
00406     using numext::isfinite;
00407     typedef typename MatrixType::Scalar Scalar;
00408     typedef typename MatrixType::RealScalar RealScalar;
00409 
00410     const Index size = mat.rows();
00411     eigen_assert(mat.cols() == size && w.size()==size);
00412 
00413     RealScalar alpha = 1;
00414 
00415     // Apply the update
00416     for (Index j = 0; j < size; j++)
00417     {
00418       // Check for termination due to an original decomposition of low-rank
00419       if (!(isfinite)(alpha))
00420         break;
00421 
00422       // Update the diagonal terms
00423       RealScalar dj = numext::real(mat.coeff(j,j));
00424       Scalar wj = w.coeff(j);
00425       RealScalar swj2 = sigma*numext::abs2(wj);
00426       RealScalar gamma = dj*alpha + swj2;
00427 
00428       mat.coeffRef(j,j) += swj2/alpha;
00429       alpha += swj2/dj;
00430 
00431 
00432       // Update the terms of L
00433       Index rs = size-j-1;
00434       w.tail(rs) -= wj * mat.col(j).tail(rs);
00435       if(gamma != 0)
00436         mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
00437     }
00438     return true;
00439   }
00440 
00441   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00442   static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1)
00443   {
00444     // Apply the permutation to the input w
00445     tmp = transpositions * w;
00446 
00447     return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma);
00448   }
00449 };
00450 
00451 template<> struct ldlt_inplace<Upper>
00452 {
00453   template<typename MatrixType, typename TranspositionType, typename Workspace>
00454   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign)
00455   {
00456     Transpose<MatrixType> matt(mat);
00457     return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
00458   }
00459 
00460   template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType>
00461   static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1)
00462   {
00463     Transpose<MatrixType> matt(mat);
00464     return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
00465   }
00466 };
00467 
00468 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower>
00469 {
00470   typedef const TriangularView<const MatrixType, UnitLower> MatrixL;
00471   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU;
00472   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
00473   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
00474 };
00475 
00476 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
00477 {
00478   typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL;
00479   typedef const TriangularView<const MatrixType, UnitUpper> MatrixU;
00480   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
00481   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
00482 };
00483 
00484 } // end namespace internal
00485 
00488 template<typename MatrixType, int _UpLo>
00489 template<typename InputType>
00490 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
00491 {
00492   check_template_parameters();
00493 
00494   eigen_assert(a.rows()==a.cols());
00495   const Index size = a.rows();
00496 
00497   m_matrix = a.derived();
00498 
00499   // Compute matrix L1 norm = max abs column sum.
00500   m_l1_norm = RealScalar(0);
00501   // TODO move this code to SelfAdjointView
00502   for (Index col = 0; col < size; ++col) {
00503     RealScalar abs_col_sum;
00504     if (_UpLo == Lower)
00505       abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
00506     else
00507       abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
00508     if (abs_col_sum > m_l1_norm)
00509       m_l1_norm = abs_col_sum;
00510   }
00511 
00512   m_transpositions.resize(size);
00513   m_isInitialized = false;
00514   m_temporary.resize(size);
00515   m_sign = internal::ZeroSign;
00516 
00517   m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success : NumericalIssue;
00518 
00519   m_isInitialized = true;
00520   return *this;
00521 }
00522 
00528 template<typename MatrixType, int _UpLo>
00529 template<typename Derived>
00530 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename LDLT<MatrixType,_UpLo>::RealScalar& sigma)
00531 {
00532   typedef typename TranspositionType::StorageIndex IndexType;
00533   const Index size = w.rows();
00534   if (m_isInitialized)
00535   {
00536     eigen_assert(m_matrix.rows()==size);
00537   }
00538   else
00539   {
00540     m_matrix.resize(size,size);
00541     m_matrix.setZero();
00542     m_transpositions.resize(size);
00543     for (Index i = 0; i < size; i++)
00544       m_transpositions.coeffRef(i) = IndexType(i);
00545     m_temporary.resize(size);
00546     m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef;
00547     m_isInitialized = true;
00548   }
00549 
00550   internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma);
00551 
00552   return *this;
00553 }
00554 
00555 #ifndef EIGEN_PARSED_BY_DOXYGEN
00556 template<typename _MatrixType, int _UpLo>
00557 template<typename RhsType, typename DstType>
00558 void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
00559 {
00560   eigen_assert(rhs.rows() == rows());
00561   // dst = P b
00562   dst = m_transpositions * rhs;
00563 
00564   // dst = L^-1 (P b)
00565   matrixL().solveInPlace(dst);
00566 
00567   // dst = D^-1 (L^-1 P b)
00568   // more precisely, use pseudo-inverse of D (see bug 241)
00569   using std::abs;
00570   const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD());
00571   // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
00572   // as motivated by LAPACK's xGELSS:
00573   // RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
00574   // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest
00575   // diagonal element is not well justified and leads to numerical issues in some cases.
00576   // Moreover, Lapack's xSYTRS routines use 0 for the tolerance.
00577   RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest();
00578 
00579   for (Index i = 0; i < vecD.size(); ++i)
00580   {
00581     if(abs(vecD(i)) > tolerance)
00582       dst.row(i) /= vecD(i);
00583     else
00584       dst.row(i).setZero();
00585   }
00586 
00587   // dst = L^-T (D^-1 L^-1 P b)
00588   matrixU().solveInPlace(dst);
00589 
00590   // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
00591   dst = m_transpositions.transpose() * dst;
00592 }
00593 #endif
00594 
00608 template<typename MatrixType,int _UpLo>
00609 template<typename Derived>
00610 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00611 {
00612   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00613   eigen_assert(m_matrix.rows() == bAndX.rows());
00614 
00615   bAndX = this->solve(bAndX);
00616 
00617   return true;
00618 }
00619 
00623 template<typename MatrixType, int _UpLo>
00624 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
00625 {
00626   eigen_assert(m_isInitialized && "LDLT is not initialized.");
00627   const Index size = m_matrix.rows();
00628   MatrixType res(size,size);
00629 
00630   // P
00631   res.setIdentity();
00632   res = transpositionsP() * res;
00633   // L^* P
00634   res = matrixU() * res;
00635   // D(L^*P)
00636   res = vectorD().real().asDiagonal() * res;
00637   // L(DL^*P)
00638   res = matrixL() * res;
00639   // P^T (LDL^*P)
00640   res = transpositionsP().transpose() * res;
00641 
00642   return res;
00643 }
00644 
00649 template<typename MatrixType, unsigned int UpLo>
00650 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00651 SelfAdjointView<MatrixType, UpLo>::ldlt() const
00652 {
00653   return LDLT<PlainObject,UpLo>(m_matrix);
00654 }
00655 
00660 template<typename Derived>
00661 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
00662 MatrixBase<Derived>::ldlt() const
00663 {
00664   return LDLT<PlainObject>(derived());
00665 }
00666 
00667 } // end namespace Eigen
00668 
00669 #endif // EIGEN_LDLT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends