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 Eigen { 
00029 
00030 namespace internal{
00031 template<typename MatrixType, int UpLo> struct LLT_Traits;
00032 }
00033 
00061  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00062   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00063   * the strict lower part does not have to store correct values.
00064   */
00065 template<typename _MatrixType, int _UpLo> class LLT
00066 {
00067   public:
00068     typedef _MatrixType MatrixType;
00069     enum {
00070       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00071       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00072       Options = MatrixType::Options,
00073       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00074     };
00075     typedef typename MatrixType::Scalar Scalar;
00076     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00077     typedef typename MatrixType::Index Index;
00078 
00079     enum {
00080       PacketSize = internal::packet_traits<Scalar>::size,
00081       AlignmentMask = int(PacketSize)-1,
00082       UpLo = _UpLo
00083     };
00084 
00085     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00086 
00093     LLT() : m_matrix(), m_isInitialized(false) {}
00094 
00101     LLT(Index size) : m_matrix(size, size),
00102                     m_isInitialized(false) {}
00103 
00104     LLT(const MatrixType& matrix)
00105       : m_matrix(matrix.rows(), matrix.cols()),
00106         m_isInitialized(false)
00107     {
00108       compute(matrix);
00109     }
00110 
00112     inline typename Traits::MatrixU matrixU() const
00113     {
00114       eigen_assert(m_isInitialized && "LLT is not initialized.");
00115       return Traits::getU(m_matrix);
00116     }
00117 
00119     inline typename Traits::MatrixL matrixL() const
00120     {
00121       eigen_assert(m_isInitialized && "LLT is not initialized.");
00122       return Traits::getL(m_matrix);
00123     }
00124 
00135     template<typename Rhs>
00136     inline const internal::solve_retval<LLT, Rhs>
00137     solve(const MatrixBase<Rhs>& b) const
00138     {
00139       eigen_assert(m_isInitialized && "LLT is not initialized.");
00140       eigen_assert(m_matrix.rows()==b.rows()
00141                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00142       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
00143     }
00144 
00145     #ifdef EIGEN2_SUPPORT
00146     template<typename OtherDerived, typename ResultType>
00147     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
00148     {
00149       *result = this->solve(b);
00150       return true;
00151     }
00152     
00153     bool isPositiveDefinite() const { return true; }
00154     #endif
00155 
00156     template<typename Derived>
00157     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00158 
00159     LLT& compute(const MatrixType& matrix);
00160 
00165     inline const MatrixType& matrixLLT() const
00166     {
00167       eigen_assert(m_isInitialized && "LLT is not initialized.");
00168       return m_matrix;
00169     }
00170 
00171     MatrixType reconstructedMatrix() const;
00172 
00173 
00179     ComputationInfo info() const
00180     {
00181       eigen_assert(m_isInitialized && "LLT is not initialized.");
00182       return m_info;
00183     }
00184 
00185     inline Index rows() const { return m_matrix.rows(); }
00186     inline Index cols() const { return m_matrix.cols(); }
00187 
00188     template<typename VectorType>
00189     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
00190 
00191   protected:
00196     MatrixType m_matrix;
00197     bool m_isInitialized;
00198     ComputationInfo m_info;
00199 };
00200 
00201 namespace internal {
00202 
00203 template<typename Scalar, int UpLo> struct llt_inplace;
00204 
00205 template<typename MatrixType, typename VectorType>
00206 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
00207 {
00208   typedef typename MatrixType::Scalar Scalar;
00209   typedef typename MatrixType::RealScalar RealScalar;
00210   typedef typename MatrixType::Index Index;
00211   typedef typename MatrixType::ColXpr ColXpr;
00212   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00213   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00214   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00215   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00216 
00217   int n = mat.cols();
00218   eigen_assert(mat.rows()==n && vec.size()==n);
00219 
00220   TempVectorType temp;
00221 
00222   if(sigma>0)
00223   {
00224     // This version is based on Givens rotations.
00225     // It is faster than the other one below, but only works for updates,
00226     // i.e., for sigma > 0
00227     temp = sqrt(sigma) * vec;
00228 
00229     for(int i=0; i<n; ++i)
00230     {
00231       JacobiRotation<Scalar> g;
00232       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00233 
00234       int rs = n-i-1;
00235       if(rs>0)
00236       {
00237         ColXprSegment x(mat.col(i).tail(rs));
00238         TempVecSegment y(temp.tail(rs));
00239         apply_rotation_in_the_plane(x, y, g);
00240       }
00241     }
00242   }
00243   else
00244   {
00245     temp = vec;
00246     RealScalar beta = 1;
00247     for(int j=0; j<n; ++j)
00248     {
00249       RealScalar Ljj = real(mat.coeff(j,j));
00250       RealScalar dj = abs2(Ljj);
00251       Scalar wj = temp.coeff(j);
00252       RealScalar swj2 = sigma*abs2(wj);
00253       RealScalar gamma = dj*beta + swj2;
00254 
00255       RealScalar x = dj + swj2/beta;
00256       if (x<=RealScalar(0))
00257         return j;
00258       RealScalar nLjj = sqrt(x);
00259       mat.coeffRef(j,j) = nLjj;
00260       beta += swj2/dj;
00261 
00262       // Update the terms of L
00263       Index rs = n-j-1;
00264       if(rs)
00265       {
00266         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00267         if(gamma != 0)
00268           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
00269       }
00270     }
00271   }
00272   return -1;
00273 }
00274 
00275 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00276 {
00277   typedef typename NumTraits<Scalar>::Real RealScalar;
00278   template<typename MatrixType>
00279   static typename MatrixType::Index unblocked(MatrixType& mat)
00280   {
00281     typedef typename MatrixType::Index Index;
00282     
00283     eigen_assert(mat.rows()==mat.cols());
00284     const Index size = mat.rows();
00285     for(Index k = 0; k < size; ++k)
00286     {
00287       Index rs = size-k-1; // remaining size
00288 
00289       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00290       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00291       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00292 
00293       RealScalar x = real(mat.coeff(k,k));
00294       if (k>0) x -= A10.squaredNorm();
00295       if (x<=RealScalar(0))
00296         return k;
00297       mat.coeffRef(k,k) = x = sqrt(x);
00298       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00299       if (rs>0) A21 *= RealScalar(1)/x;
00300     }
00301     return -1;
00302   }
00303 
00304   template<typename MatrixType>
00305   static typename MatrixType::Index blocked(MatrixType& m)
00306   {
00307     typedef typename MatrixType::Index Index;
00308     eigen_assert(m.rows()==m.cols());
00309     Index size = m.rows();
00310     if(size<32)
00311       return unblocked(m);
00312 
00313     Index blockSize = size/8;
00314     blockSize = (blockSize/16)*16;
00315     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00316 
00317     for (Index k=0; k<size; k+=blockSize)
00318     {
00319       // partition the matrix:
00320       //       A00 |  -  |  -
00321       // lu  = A10 | A11 |  -
00322       //       A20 | A21 | A22
00323       Index bs = (std::min)(blockSize, size-k);
00324       Index rs = size - k - bs;
00325       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00326       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00327       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00328 
00329       Index ret;
00330       if((ret=unblocked(A11))>=0) return k+ret;
00331       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00332       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
00333     }
00334     return -1;
00335   }
00336 
00337   template<typename MatrixType, typename VectorType>
00338   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00339   {
00340     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00341   }
00342 };
00343   
00344 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00345 {
00346   typedef typename NumTraits<Scalar>::Real RealScalar;
00347 
00348   template<typename MatrixType>
00349   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
00350   {
00351     Transpose<MatrixType> matt(mat);
00352     return llt_inplace<Scalar, Lower>::unblocked(matt);
00353   }
00354   template<typename MatrixType>
00355   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
00356   {
00357     Transpose<MatrixType> matt(mat);
00358     return llt_inplace<Scalar, Lower>::blocked(matt);
00359   }
00360   template<typename MatrixType, typename VectorType>
00361   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00362   {
00363     Transpose<MatrixType> matt(mat);
00364     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00365   }
00366 };
00367 
00368 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00369 {
00370   typedef const TriangularView<const MatrixType, Lower> MatrixL;
00371   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00372   static inline MatrixL getL(const MatrixType& m) { return m; }
00373   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00374   static bool inplace_decomposition(MatrixType& m)
00375   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00376 };
00377 
00378 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00379 {
00380   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00381   typedef const TriangularView<const MatrixType, Upper> MatrixU;
00382   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
00383   static inline MatrixU getU(const MatrixType& m) { return m; }
00384   static bool inplace_decomposition(MatrixType& m)
00385   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00386 };
00387 
00388 } // end namespace internal
00389 
00397 template<typename MatrixType, int _UpLo>
00398 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
00399 {
00400   eigen_assert(a.rows()==a.cols());
00401   const Index size = a.rows();
00402   m_matrix.resize(size, size);
00403   m_matrix = a;
00404 
00405   m_isInitialized = true;
00406   bool ok = Traits::inplace_decomposition(m_matrix);
00407   m_info = ok ? Success : NumericalIssue;
00408 
00409   return *this;
00410 }
00411 
00417 template<typename _MatrixType, int _UpLo>
00418 template<typename VectorType>
00419 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00420 {
00421   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00422   eigen_assert(v.size()==m_matrix.cols());
00423   eigen_assert(m_isInitialized);
00424   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00425     m_info = NumericalIssue;
00426   else
00427     m_info = Success;
00428 
00429   return *this;
00430 }
00431     
00432 namespace internal {
00433 template<typename _MatrixType, int UpLo, typename Rhs>
00434 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
00435   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
00436 {
00437   typedef LLT<_MatrixType,UpLo> LLTType;
00438   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
00439 
00440   template<typename Dest> void evalTo(Dest& dst) const
00441   {
00442     dst = rhs();
00443     dec().solveInPlace(dst);
00444   }
00445 };
00446 }
00447 
00461 template<typename MatrixType, int _UpLo>
00462 template<typename Derived>
00463 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00464 {
00465   eigen_assert(m_isInitialized && "LLT is not initialized.");
00466   eigen_assert(m_matrix.rows()==bAndX.rows());
00467   matrixL().solveInPlace(bAndX);
00468   matrixU().solveInPlace(bAndX);
00469 }
00470 
00474 template<typename MatrixType, int _UpLo>
00475 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00476 {
00477   eigen_assert(m_isInitialized && "LLT is not initialized.");
00478   return matrixL() * matrixL().adjoint().toDenseMatrix();
00479 }
00480 
00484 template<typename Derived>
00485 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00486 MatrixBase<Derived>::llt() const
00487 {
00488   return LLT<PlainObject>(derived());
00489 }
00490 
00494 template<typename MatrixType, unsigned int UpLo>
00495 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00496 SelfAdjointView<MatrixType, UpLo>::llt() const
00497 {
00498   return LLT<PlainObject,UpLo>(m_matrix);
00499 }
00500 
00501 } // end namespace Eigen
00502 
00503 #endif // EIGEN_LLT_H