LLT.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_LLT_H
00011 #define EIGEN_LLT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal{
00016 template<typename MatrixType, int UpLo> struct LLT_Traits;
00017 }
00018 
00019 /** \ingroup Cholesky_Module
00020   *
00021   * \class LLT
00022   *
00023   * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
00024   *
00025   * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
00026   * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
00027   *             The other triangular part won't be read.
00028   *
00029   * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
00030   * matrix A such that A = LL^* = U^*U, where L is lower triangular.
00031   *
00032   * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
00033   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
00034   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
00035   * situations like generalised eigen problems with hermitian matrices.
00036   *
00037   * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
00038   * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
00039   * has a solution.
00040   *
00041   * Example: \include LLT_example.cpp
00042   * Output: \verbinclude LLT_example.out
00043   *    
00044   * \sa MatrixBase::llt(), class LDLT
00045   */
00046  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00047   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00048   * the strict lower part does not have to store correct values.
00049   */
00050 template<typename _MatrixType, int _UpLo> class LLT
00051 {
00052   public:
00053     typedef _MatrixType MatrixType;
00054     enum {
00055       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00056       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00057       Options = MatrixType::Options,
00058       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00059     };
00060     typedef typename MatrixType::Scalar Scalar;
00061     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00062     typedef typename MatrixType::Index Index;
00063 
00064     enum {
00065       PacketSize = internal::packet_traits<Scalar>::size,
00066       AlignmentMask = int(PacketSize)-1,
00067       UpLo = _UpLo
00068     };
00069 
00070     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00071 
00072     /**
00073       * \brief Default Constructor.
00074       *
00075       * The default constructor is useful in cases in which the user intends to
00076       * perform decompositions via LLT::compute(const MatrixType&).
00077       */
00078     LLT() : m_matrix(), m_isInitialized(false) {}
00079 
00080     /** \brief Default Constructor with memory preallocation
00081       *
00082       * Like the default constructor but with preallocation of the internal data
00083       * according to the specified problem \a size.
00084       * \sa LLT()
00085       */
00086     LLT(Index size) : m_matrix(size, size),
00087                     m_isInitialized(false) {}
00088 
00089     LLT(const MatrixType& matrix)
00090       : m_matrix(matrix.rows(), matrix.cols()),
00091         m_isInitialized(false)
00092     {
00093       compute(matrix);
00094     }
00095 
00096     /** \returns a view of the upper triangular matrix U */
00097     inline typename Traits::MatrixU matrixU () const
00098     {
00099       eigen_assert(m_isInitialized && "LLT is not initialized.");
00100       return Traits::getU(m_matrix);
00101     }
00102 
00103     /** \returns a view of the lower triangular matrix L */
00104     inline typename Traits::MatrixL matrixL () const
00105     {
00106       eigen_assert(m_isInitialized && "LLT is not initialized.");
00107       return Traits::getL(m_matrix);
00108     }
00109 
00110     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
00111       *
00112       * Since this LLT class assumes anyway that the matrix A is invertible, the solution
00113       * theoretically exists and is unique regardless of b.
00114       *
00115       * Example: \include LLT_solve.cpp
00116       * Output: \verbinclude LLT_solve.out
00117       *
00118       * \sa solveInPlace(), MatrixBase::llt()
00119       */
00120     template<typename Rhs>
00121     inline const internal::solve_retval<LLT, Rhs>
00122     solve (const MatrixBase<Rhs>& b) const
00123     {
00124       eigen_assert(m_isInitialized && "LLT is not initialized.");
00125       eigen_assert(m_matrix.rows()==b.rows()
00126                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00127       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00128     }
00129 
00130     #ifdef EIGEN2_SUPPORT
00131     template<typename OtherDerived, typename ResultType>
00132     bool solve (const MatrixBase<OtherDerived>& b, ResultType *result) const
00133     {
00134       *result = this->solve (b);
00135       return true;
00136     }
00137     
00138     bool isPositiveDefinite() const { return true; }
00139     #endif
00140 
00141     template<typename Derived>
00142     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00143 
00144     LLT& compute(const MatrixType& matrix);
00145 
00146     /** \returns the LLT decomposition matrix
00147       *
00148       * TODO: document the storage layout
00149       */
00150     inline const MatrixType& matrixLLT () const
00151     {
00152       eigen_assert(m_isInitialized && "LLT is not initialized.");
00153       return m_matrix;
00154     }
00155 
00156     MatrixType reconstructedMatrix () const;
00157 
00158 
00159     /** \brief Reports whether previous computation was successful.
00160       *
00161       * \returns \c Success if computation was succesful,
00162       *          \c NumericalIssue if the matrix.appears to be negative.
00163       */
00164     ComputationInfo info() const
00165     {
00166       eigen_assert(m_isInitialized && "LLT is not initialized.");
00167       return m_info;
00168     }
00169 
00170     inline Index rows() const { return m_matrix.rows(); }
00171     inline Index cols() const { return m_matrix.cols(); }
00172 
00173     template<typename VectorType>
00174     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
00175 
00176   protected:
00177     
00178     static void check_template_parameters()
00179     {
00180       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00181     }
00182     
00183     /** \internal
00184       * Used to compute and store L
00185       * The strict upper part is not used and even not initialized.
00186       */
00187     MatrixType m_matrix;
00188     bool m_isInitialized;
00189     ComputationInfo m_info;
00190 };
00191 
00192 namespace internal {
00193 
00194 template<typename Scalar, int UpLo> struct llt_inplace;
00195 
00196 template<typename MatrixType, typename VectorType>
00197 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
00198 {
00199   using std::sqrt;
00200   typedef typename MatrixType::Scalar Scalar;
00201   typedef typename MatrixType::RealScalar RealScalar;
00202   typedef typename MatrixType::Index Index;
00203   typedef typename MatrixType::ColXpr ColXpr;
00204   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00205   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00206   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00207   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00208 
00209   Index n = mat.cols();
00210   eigen_assert(mat.rows()==n && vec.size()==n);
00211 
00212   TempVectorType temp;
00213 
00214   if(sigma>0)
00215   {
00216     // This version is based on Givens rotations.
00217     // It is faster than the other one below, but only works for updates,
00218     // i.e., for sigma > 0
00219     temp = sqrt(sigma) * vec;
00220 
00221     for(Index i=0; i<n; ++i)
00222     {
00223       JacobiRotation<Scalar> g;
00224       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00225 
00226       Index rs = n-i-1;
00227       if(rs>0)
00228       {
00229         ColXprSegment x(mat.col(i).tail(rs));
00230         TempVecSegment y(temp.tail(rs));
00231         apply_rotation_in_the_plane(x, y, g);
00232       }
00233     }
00234   }
00235   else
00236   {
00237     temp = vec;
00238     RealScalar beta = 1;
00239     for(Index j=0; j<n; ++j)
00240     {
00241       RealScalar Ljj = numext::real(mat.coeff(j,j));
00242       RealScalar dj = numext::abs2(Ljj);
00243       Scalar wj = temp.coeff(j);
00244       RealScalar swj2 = sigma*numext::abs2(wj);
00245       RealScalar gamma = dj*beta + swj2;
00246 
00247       RealScalar x = dj + swj2/beta;
00248       if (x<=RealScalar(0))
00249         return j;
00250       RealScalar nLjj = sqrt(x);
00251       mat.coeffRef(j,j) = nLjj;
00252       beta += swj2/dj;
00253 
00254       // Update the terms of L
00255       Index rs = n-j-1;
00256       if(rs)
00257       {
00258         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00259         if(gamma != 0)
00260           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
00261       }
00262     }
00263   }
00264   return -1;
00265 }
00266 
00267 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00268 {
00269   typedef typename NumTraits<Scalar>::Real RealScalar;
00270   template<typename MatrixType>
00271   static typename MatrixType::Index unblocked(MatrixType& mat)
00272   {
00273     using std::sqrt;
00274     typedef typename MatrixType::Index Index;
00275     
00276     eigen_assert(mat.rows()==mat.cols());
00277     const Index size = mat.rows();
00278     for(Index k = 0; k < size; ++k)
00279     {
00280       Index rs = size-k-1; // remaining size
00281 
00282       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00283       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00284       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00285 
00286       RealScalar x = numext::real(mat.coeff(k,k));
00287       if (k>0) x -= A10.squaredNorm();
00288       if (x<=RealScalar(0))
00289         return k;
00290       mat.coeffRef(k,k) = x = sqrt(x);
00291       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00292       if (rs>0) A21 /= x;
00293     }
00294     return -1;
00295   }
00296 
00297   template<typename MatrixType>
00298   static typename MatrixType::Index blocked(MatrixType& m)
00299   {
00300     typedef typename MatrixType::Index Index;
00301     eigen_assert(m.rows()==m.cols());
00302     Index size = m.rows();
00303     if(size<32)
00304       return unblocked(m);
00305 
00306     Index blockSize = size/8;
00307     blockSize = (blockSize/16)*16;
00308     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00309 
00310     for (Index k=0; k<size; k+=blockSize)
00311     {
00312       // partition the matrix:
00313       //       A00 |  -  |  -
00314       // lu  = A10 | A11 |  -
00315       //       A20 | A21 | A22
00316       Index bs = (std::min)(blockSize, size-k);
00317       Index rs = size - k - bs;
00318       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00319       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00320       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00321 
00322       Index ret;
00323       if((ret=unblocked(A11))>=0) return k+ret;
00324       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00325       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00326     }
00327     return -1;
00328   }
00329 
00330   template<typename MatrixType, typename VectorType>
00331   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00332   {
00333     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00334   }
00335 };
00336   
00337 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00338 {
00339   typedef typename NumTraits<Scalar>::Real RealScalar;
00340 
00341   template<typename MatrixType>
00342   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00343   {
00344     Transpose<MatrixType> matt(mat);
00345     return llt_inplace<Scalar, Lower>::unblocked(matt);
00346   }
00347   template<typename MatrixType>
00348   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00349   {
00350     Transpose<MatrixType> matt(mat);
00351     return llt_inplace<Scalar, Lower>::blocked(matt);
00352   }
00353   template<typename MatrixType, typename VectorType>
00354   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00355   {
00356     Transpose<MatrixType> matt(mat);
00357     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00358   }
00359 };
00360 
00361 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00362 {
00363   typedef const TriangularView<const MatrixType, Lower> MatrixL;
00364   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00365   static inline MatrixL getL(const MatrixType& m) { return m; }
00366   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00367   static bool inplace_decomposition(MatrixType& m)
00368   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00369 };
00370 
00371 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00372 {
00373   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00374   typedef const TriangularView<const MatrixType, Upper> MatrixU;
00375   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00376   static inline MatrixU getU(const MatrixType& m) { return m; }
00377   static bool inplace_decomposition(MatrixType& m)
00378   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00379 };
00380 
00381 } // end namespace internal
00382 
00383 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
00384   *
00385   * \returns a reference to *this
00386   *
00387   * Example: \include TutorialLinAlgComputeTwice.cpp
00388   * Output: \verbinclude TutorialLinAlgComputeTwice.out
00389   */
00390 template<typename MatrixType, int _UpLo>
00391 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00392 {
00393   check_template_parameters();
00394   
00395   eigen_assert(a.rows()==a.cols());
00396   const Index size = a.rows();
00397   m_matrix.resize(size, size);
00398   m_matrix = a;
00399 
00400   m_isInitialized = true;
00401   bool ok = Traits::inplace_decomposition(m_matrix);
00402   m_info = ok ? Success : NumericalIssue;
00403 
00404   return *this;
00405 }
00406 
00407 /** Performs a rank one update (or dowdate) of the current decomposition.
00408   * If A = LL^* before the rank one update,
00409   * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
00410   * of same dimension.
00411   */
00412 template<typename _MatrixType, int _UpLo>
00413 template<typename VectorType>
00414 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00415 {
00416   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00417   eigen_assert(v.size()==m_matrix.cols());
00418   eigen_assert(m_isInitialized);
00419   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00420     m_info = NumericalIssue;
00421   else
00422     m_info = Success;
00423 
00424   return *this;
00425 }
00426     
00427 namespace internal {
00428 template<typename _MatrixType, int UpLo, typename Rhs>
00429 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00430   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00431 {
00432   typedef LLT<_MatrixType,UpLo> LLTType;
00433   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00434 
00435   template<typename Dest> void evalTo(Dest& dst) const
00436   {
00437     dst = rhs();
00438     dec().solveInPlace(dst);
00439   }
00440 };
00441 }
00442 
00443 /** \internal use x = llt_object.solve(x);
00444   * 
00445   * This is the \em in-place version of solve().
00446   *
00447   * \param bAndX represents both the right-hand side matrix b and result x.
00448   *
00449   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
00450   *
00451   * This version avoids a copy when the right hand side matrix b is not
00452   * needed anymore.
00453   *
00454   * \sa LLT::solve(), MatrixBase::llt()
00455   */
00456 template<typename MatrixType, int _UpLo>
00457 template<typename Derived>
00458 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00459 {
00460   eigen_assert(m_isInitialized && "LLT is not initialized.");
00461   eigen_assert(m_matrix.rows()==bAndX.rows());
00462   matrixL().solveInPlace(bAndX);
00463   matrixU().solveInPlace(bAndX);
00464 }
00465 
00466 /** \returns the matrix represented by the decomposition,
00467  * i.e., it returns the product: L L^*.
00468  * This function is provided for debug purpose. */
00469 template<typename MatrixType, int _UpLo>
00470 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix () const
00471 {
00472   eigen_assert(m_isInitialized && "LLT is not initialized.");
00473   return matrixL() * matrixL().adjoint().toDenseMatrix();
00474 }
00475 
00476 /** \cholesky_module
00477   * \returns the LLT decomposition of \c *this
00478   */
00479 template<typename Derived>
00480 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00481 MatrixBase<Derived>::llt () const
00482 {
00483   return LLT<PlainObject>(derived());
00484 }
00485 
00486 /** \cholesky_module
00487   * \returns the LLT decomposition of \c *this
00488   */
00489 template<typename MatrixType, unsigned int UpLo>
00490 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00491 SelfAdjointView<MatrixType, UpLo>::llt () const
00492 {
00493   return LLT<PlainObject,UpLo>(m_matrix);
00494 }
00495 
00496 } // end namespace Eigen
00497 
00498 #endif // EIGEN_LLT_H