Main MRPT website > C++ reference
MRPT logo

LLT.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 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_LLT_H
00026 #define EIGEN_LLT_H
00027 
00028 namespace internal{
00029 template<typename MatrixType, int UpLo> struct LLT_Traits;
00030 }
00031 
00032 /** \ingroup cholesky_Module
00033   *
00034   * \class LLT
00035   *
00036   * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
00037   *
00038   * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
00039   *
00040   * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
00041   * matrix A such that A = LL^* = U^*U, where L is lower triangular.
00042   *
00043   * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
00044   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
00045   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
00046   * situations like generalised eigen problems with hermitian matrices.
00047   *
00048   * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
00049   * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
00050   * has a solution.
00051   *
00052   * \sa MatrixBase::llt(), class LDLT
00053   */
00054  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00055   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00056   * the strict lower part does not have to store correct values.
00057   */
00058 template<typename _MatrixType, int _UpLo> class LLT
00059 {
00060   public:
00061     typedef _MatrixType MatrixType;
00062     enum {
00063       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065       Options = MatrixType::Options,
00066       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00067     };
00068     typedef typename MatrixType::Scalar Scalar;
00069     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00070     typedef typename MatrixType::Index Index;
00071 
00072     enum {
00073       PacketSize = internal::packet_traits<Scalar>::size,
00074       AlignmentMask = int(PacketSize)-1,
00075       UpLo = _UpLo
00076     };
00077 
00078     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00079 
00080     /**
00081       * \brief Default Constructor.
00082       *
00083       * The default constructor is useful in cases in which the user intends to
00084       * perform decompositions via LLT::compute(const MatrixType&).
00085       */
00086     LLT() : m_matrix(), m_isInitialized(false) {}
00087 
00088     /** \brief Default Constructor with memory preallocation
00089       *
00090       * Like the default constructor but with preallocation of the internal data
00091       * according to the specified problem \a size.
00092       * \sa LLT()
00093       */
00094     LLT(Index size) : m_matrix(size, size),
00095                     m_isInitialized(false) {}
00096 
00097     LLT(const MatrixType& matrix)
00098       : m_matrix(matrix.rows(), matrix.cols()),
00099         m_isInitialized(false)
00100     {
00101       compute(matrix);
00102     }
00103 
00104     /** \returns a view of the upper triangular matrix U */
00105     inline typename Traits::MatrixU matrixU() const
00106     {
00107       eigen_assert(m_isInitialized && "LLT is not initialized.");
00108       return Traits::getU(m_matrix);
00109     }
00110 
00111     /** \returns a view of the lower triangular matrix L */
00112     inline typename Traits::MatrixL matrixL() const
00113     {
00114       eigen_assert(m_isInitialized && "LLT is not initialized.");
00115       return Traits::getL(m_matrix);
00116     }
00117 
00118     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
00119       *
00120       * Since this LLT class assumes anyway that the matrix A is invertible, the solution
00121       * theoretically exists and is unique regardless of b.
00122       *
00123       * Example: \include LLT_solve.cpp
00124       * Output: \verbinclude LLT_solve.out
00125       *
00126       * \sa solveInPlace(), MatrixBase::llt()
00127       */
00128     template<typename Rhs>
00129     inline const internal::solve_retval<LLT, Rhs>
00130     solve(const MatrixBase<Rhs>& b) const
00131     {
00132       eigen_assert(m_isInitialized && "LLT is not initialized.");
00133       eigen_assert(m_matrix.rows()==b.rows()
00134                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00135       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00136     }
00137 
00138     template<typename Derived>
00139     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00140 
00141     LLT& compute(const MatrixType& matrix);
00142 
00143     /** \returns the LLT decomposition matrix
00144       *
00145       * TODO: document the storage layout
00146       */
00147     inline const MatrixType& matrixLLT() const
00148     {
00149       eigen_assert(m_isInitialized && "LLT is not initialized.");
00150       return m_matrix;
00151     }
00152 
00153     MatrixType reconstructedMatrix() const;
00154 
00155 
00156     /** \brief Reports whether previous computation was successful.
00157       *
00158       * \returns \c Success if computation was succesful,
00159       *          \c NumericalIssue if the matrix.appears to be negative.
00160       */
00161     ComputationInfo info() const
00162     {
00163       eigen_assert(m_isInitialized && "LLT is not initialized.");
00164       return m_info;
00165     }
00166 
00167     inline Index rows() const { return m_matrix.rows(); }
00168     inline Index cols() const { return m_matrix.cols(); }
00169 
00170   protected:
00171     /** \internal
00172       * Used to compute and store L
00173       * The strict upper part is not used and even not initialized.
00174       */
00175     MatrixType m_matrix;
00176     bool m_isInitialized;
00177     ComputationInfo m_info;
00178 };
00179 
00180 namespace internal {
00181 
00182 template<int UpLo> struct llt_inplace;
00183 
00184 template<> struct llt_inplace<Lower>
00185 {
00186   template<typename MatrixType>
00187   static bool unblocked(MatrixType& mat)
00188   {
00189     typedef typename MatrixType::Scalar Scalar;
00190     typedef typename MatrixType::RealScalar RealScalar;
00191     typedef typename MatrixType::Index Index;
00192     eigen_assert(mat.rows()==mat.cols());
00193     const Index size = mat.rows();
00194     for(Index k = 0; k < size; ++k)
00195     {
00196       Index rs = size-k-1; // remaining size
00197 
00198       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00199       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00200       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00201 
00202       RealScalar x = real(mat.coeff(k,k));
00203       if (k>0) x -= A10.squaredNorm();
00204       if (x<=RealScalar(0))
00205         return false;
00206       mat.coeffRef(k,k) = x = sqrt(x);
00207       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00208       if (rs>0) A21 *= RealScalar(1)/x;
00209     }
00210     return true;
00211   }
00212 
00213   template<typename MatrixType>
00214   static bool blocked(MatrixType& m)
00215   {
00216     typedef typename MatrixType::Index Index;
00217     eigen_assert(m.rows()==m.cols());
00218     Index size = m.rows();
00219     if(size<32)
00220       return unblocked(m);
00221 
00222     Index blockSize = size/8;
00223     blockSize = (blockSize/16)*16;
00224     blockSize = std::min(std::max(blockSize,Index(8)), Index(128));
00225 
00226     for (Index k=0; k<size; k+=blockSize)
00227     {
00228       // partition the matrix:
00229       //       A00 |  -  |  -
00230       // lu  = A10 | A11 |  -
00231       //       A20 | A21 | A22
00232       Index bs = std::min(blockSize, size-k);
00233       Index rs = size - k - bs;
00234       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00235       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00236       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00237 
00238       if(!unblocked(A11)) return false;
00239       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00240       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00241     }
00242     return true;
00243   }
00244 };
00245 
00246 template<> struct llt_inplace<Upper>
00247 {
00248   template<typename MatrixType>
00249   static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
00250   {
00251     Transpose<MatrixType> matt(mat);
00252     return llt_inplace<Lower>::unblocked(matt);
00253   }
00254   template<typename MatrixType>
00255   static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
00256   {
00257     Transpose<MatrixType> matt(mat);
00258     return llt_inplace<Lower>::blocked(matt);
00259   }
00260 };
00261 
00262 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00263 {
00264   typedef TriangularView<MatrixType, Lower> MatrixL;
00265   typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU;
00266   inline static MatrixL getL(const MatrixType& m) { return m; }
00267   inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00268   static bool inplace_decomposition(MatrixType& m)
00269   { return llt_inplace<Lower>::blocked(m); }
00270 };
00271 
00272 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00273 {
00274   typedef TriangularView<typename MatrixType::AdjointReturnType, Lower> MatrixL;
00275   typedef TriangularView<MatrixType, Upper> MatrixU;
00276   inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00277   inline static MatrixU getU(const MatrixType& m) { return m; }
00278   static bool inplace_decomposition(MatrixType& m)
00279   { return llt_inplace<Upper>::blocked(m); }
00280 };
00281 
00282 } // end namespace internal
00283 
00284 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
00285   *
00286   *
00287   * \returns a reference to *this
00288   */
00289 template<typename MatrixType, int _UpLo>
00290 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00291 {
00292   assert(a.rows()==a.cols());
00293   const Index size = a.rows();
00294   m_matrix.resize(size, size);
00295   m_matrix = a;
00296 
00297   m_isInitialized = true;
00298   bool ok = Traits::inplace_decomposition(m_matrix);
00299   m_info = ok ? Success : NumericalIssue;
00300 
00301   return *this;
00302 }
00303 
00304 namespace internal {
00305 template<typename _MatrixType, int UpLo, typename Rhs>
00306 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00307   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00308 {
00309   typedef LLT<_MatrixType,UpLo> LLTType;
00310   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00311 
00312   template<typename Dest> void evalTo(Dest& dst) const
00313   {
00314     dst = rhs();
00315     dec().solveInPlace(dst);
00316   }
00317 };
00318 }
00319 
00320 /** \internal use x = llt_object.solve(x);
00321   * 
00322   * This is the \em in-place version of solve().
00323   *
00324   * \param bAndX represents both the right-hand side matrix b and result x.
00325   *
00326   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
00327   *
00328   * This version avoids a copy when the right hand side matrix b is not
00329   * needed anymore.
00330   *
00331   * \sa LLT::solve(), MatrixBase::llt()
00332   */
00333 template<typename MatrixType, int _UpLo>
00334 template<typename Derived>
00335 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00336 {
00337   eigen_assert(m_isInitialized && "LLT is not initialized.");
00338   eigen_assert(m_matrix.rows()==bAndX.rows());
00339   matrixL().solveInPlace(bAndX);
00340   matrixU().solveInPlace(bAndX);
00341 }
00342 
00343 /** \returns the matrix represented by the decomposition,
00344  * i.e., it returns the product: L L^*.
00345  * This function is provided for debug purpose. */
00346 template<typename MatrixType, int _UpLo>
00347 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00348 {
00349   eigen_assert(m_isInitialized && "LLT is not initialized.");
00350   return matrixL() * matrixL().adjoint().toDenseMatrix();
00351 }
00352 
00353 /** \cholesky_module
00354   * \returns the LLT decomposition of \c *this
00355   */
00356 template<typename Derived>
00357 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00358 MatrixBase<Derived>::llt() const
00359 {
00360   return LLT<PlainObject>(derived());
00361 }
00362 
00363 /** \cholesky_module
00364   * \returns the LLT decomposition of \c *this
00365   */
00366 template<typename MatrixType, unsigned int UpLo>
00367 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00368 SelfAdjointView<MatrixType, UpLo>::llt() const
00369 {
00370   return LLT<PlainObject,UpLo>(m_matrix);
00371 }
00372 
00373 #endif // EIGEN_LLT_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011