Eigen  3.3.3
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 
00048  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
00049   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
00050   * the strict lower part does not have to store correct values.
00051   */
00052 template<typename _MatrixType, int _UpLo> class LLT
00053 {
00054   public:
00055     typedef _MatrixType MatrixType;
00056     enum {
00057       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00058       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00059       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
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 
00066     enum {
00067       PacketSize = internal::packet_traits<Scalar>::size,
00068       AlignmentMask = int(PacketSize)-1,
00069       UpLo = _UpLo
00070     };
00071 
00072     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
00073 
00080     LLT() : m_matrix(), m_isInitialized(false) {}
00081 
00088     explicit LLT(Index size) : m_matrix(size, size),
00089                     m_isInitialized(false) {}
00090 
00091     template<typename InputType>
00092     explicit LLT(const EigenBase<InputType>& matrix)
00093       : m_matrix(matrix.rows(), matrix.cols()),
00094         m_isInitialized(false)
00095     {
00096       compute(matrix.derived());
00097     }
00098 
00106     template<typename InputType>
00107     explicit LLT(EigenBase<InputType>& matrix)
00108       : m_matrix(matrix.derived()),
00109         m_isInitialized(false)
00110     {
00111       compute(matrix.derived());
00112     }
00113 
00115     inline typename Traits::MatrixU matrixU() const
00116     {
00117       eigen_assert(m_isInitialized && "LLT is not initialized.");
00118       return Traits::getU(m_matrix);
00119     }
00120 
00122     inline typename Traits::MatrixL matrixL() const
00123     {
00124       eigen_assert(m_isInitialized && "LLT is not initialized.");
00125       return Traits::getL(m_matrix);
00126     }
00127 
00138     template<typename Rhs>
00139     inline const Solve<LLT, Rhs>
00140     solve(const MatrixBase<Rhs>& b) const
00141     {
00142       eigen_assert(m_isInitialized && "LLT is not initialized.");
00143       eigen_assert(m_matrix.rows()==b.rows()
00144                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
00145       return Solve<LLT, Rhs>(*this, b.derived());
00146     }
00147 
00148     template<typename Derived>
00149     void solveInPlace(MatrixBase<Derived> &bAndX) const;
00150 
00151     template<typename InputType>
00152     LLT& compute(const EigenBase<InputType>& matrix);
00153 
00157     RealScalar rcond() const
00158     {
00159       eigen_assert(m_isInitialized && "LLT is not initialized.");
00160       eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
00161       return internal::rcond_estimate_helper(m_l1_norm, *this);
00162     }
00163 
00168     inline const MatrixType& matrixLLT() const
00169     {
00170       eigen_assert(m_isInitialized && "LLT is not initialized.");
00171       return m_matrix;
00172     }
00173 
00174     MatrixType reconstructedMatrix() const;
00175 
00176 
00182     ComputationInfo info() const
00183     {
00184       eigen_assert(m_isInitialized && "LLT is not initialized.");
00185       return m_info;
00186     }
00187 
00193     const LLT& adjoint() const { return *this; };
00194 
00195     inline Index rows() const { return m_matrix.rows(); }
00196     inline Index cols() const { return m_matrix.cols(); }
00197 
00198     template<typename VectorType>
00199     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
00200 
00201     #ifndef EIGEN_PARSED_BY_DOXYGEN
00202     template<typename RhsType, typename DstType>
00203     EIGEN_DEVICE_FUNC
00204     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00205     #endif
00206 
00207   protected:
00208 
00209     static void check_template_parameters()
00210     {
00211       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00212     }
00213 
00218     MatrixType m_matrix;
00219     RealScalar m_l1_norm;
00220     bool m_isInitialized;
00221     ComputationInfo m_info;
00222 };
00223 
00224 namespace internal {
00225 
00226 template<typename Scalar, int UpLo> struct llt_inplace;
00227 
00228 template<typename MatrixType, typename VectorType>
00229 static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
00230 {
00231   using std::sqrt;
00232   typedef typename MatrixType::Scalar Scalar;
00233   typedef typename MatrixType::RealScalar RealScalar;
00234   typedef typename MatrixType::ColXpr ColXpr;
00235   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
00236   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
00237   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
00238   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
00239 
00240   Index n = mat.cols();
00241   eigen_assert(mat.rows()==n && vec.size()==n);
00242 
00243   TempVectorType temp;
00244 
00245   if(sigma>0)
00246   {
00247     // This version is based on Givens rotations.
00248     // It is faster than the other one below, but only works for updates,
00249     // i.e., for sigma > 0
00250     temp = sqrt(sigma) * vec;
00251 
00252     for(Index i=0; i<n; ++i)
00253     {
00254       JacobiRotation<Scalar> g;
00255       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
00256 
00257       Index rs = n-i-1;
00258       if(rs>0)
00259       {
00260         ColXprSegment x(mat.col(i).tail(rs));
00261         TempVecSegment y(temp.tail(rs));
00262         apply_rotation_in_the_plane(x, y, g);
00263       }
00264     }
00265   }
00266   else
00267   {
00268     temp = vec;
00269     RealScalar beta = 1;
00270     for(Index j=0; j<n; ++j)
00271     {
00272       RealScalar Ljj = numext::real(mat.coeff(j,j));
00273       RealScalar dj = numext::abs2(Ljj);
00274       Scalar wj = temp.coeff(j);
00275       RealScalar swj2 = sigma*numext::abs2(wj);
00276       RealScalar gamma = dj*beta + swj2;
00277 
00278       RealScalar x = dj + swj2/beta;
00279       if (x<=RealScalar(0))
00280         return j;
00281       RealScalar nLjj = sqrt(x);
00282       mat.coeffRef(j,j) = nLjj;
00283       beta += swj2/dj;
00284 
00285       // Update the terms of L
00286       Index rs = n-j-1;
00287       if(rs)
00288       {
00289         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
00290         if(gamma != 0)
00291           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
00292       }
00293     }
00294   }
00295   return -1;
00296 }
00297 
00298 template<typename Scalar> struct llt_inplace<Scalar, Lower>
00299 {
00300   typedef typename NumTraits<Scalar>::Real RealScalar;
00301   template<typename MatrixType>
00302   static Index unblocked(MatrixType& mat)
00303   {
00304     using std::sqrt;
00305 
00306     eigen_assert(mat.rows()==mat.cols());
00307     const Index size = mat.rows();
00308     for(Index k = 0; k < size; ++k)
00309     {
00310       Index rs = size-k-1; // remaining size
00311 
00312       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
00313       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
00314       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
00315 
00316       RealScalar x = numext::real(mat.coeff(k,k));
00317       if (k>0) x -= A10.squaredNorm();
00318       if (x<=RealScalar(0))
00319         return k;
00320       mat.coeffRef(k,k) = x = sqrt(x);
00321       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
00322       if (rs>0) A21 /= x;
00323     }
00324     return -1;
00325   }
00326 
00327   template<typename MatrixType>
00328   static Index blocked(MatrixType& m)
00329   {
00330     eigen_assert(m.rows()==m.cols());
00331     Index size = m.rows();
00332     if(size<32)
00333       return unblocked(m);
00334 
00335     Index blockSize = size/8;
00336     blockSize = (blockSize/16)*16;
00337     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
00338 
00339     for (Index k=0; k<size; k+=blockSize)
00340     {
00341       // partition the matrix:
00342       //       A00 |  -  |  -
00343       // lu  = A10 | A11 |  -
00344       //       A20 | A21 | A22
00345       Index bs = (std::min)(blockSize, size-k);
00346       Index rs = size - k - bs;
00347       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
00348       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
00349       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
00350 
00351       Index ret;
00352       if((ret=unblocked(A11))>=0) return k+ret;
00353       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
00354       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck
00355     }
00356     return -1;
00357   }
00358 
00359   template<typename MatrixType, typename VectorType>
00360   static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00361   {
00362     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
00363   }
00364 };
00365 
00366 template<typename Scalar> struct llt_inplace<Scalar, Upper>
00367 {
00368   typedef typename NumTraits<Scalar>::Real RealScalar;
00369 
00370   template<typename MatrixType>
00371   static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat)
00372   {
00373     Transpose<MatrixType> matt(mat);
00374     return llt_inplace<Scalar, Lower>::unblocked(matt);
00375   }
00376   template<typename MatrixType>
00377   static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat)
00378   {
00379     Transpose<MatrixType> matt(mat);
00380     return llt_inplace<Scalar, Lower>::blocked(matt);
00381   }
00382   template<typename MatrixType, typename VectorType>
00383   static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
00384   {
00385     Transpose<MatrixType> matt(mat);
00386     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
00387   }
00388 };
00389 
00390 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
00391 {
00392   typedef const TriangularView<const MatrixType, Lower> MatrixL;
00393   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
00394   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
00395   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
00396   static bool inplace_decomposition(MatrixType& m)
00397   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
00398 };
00399 
00400 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
00401 {
00402   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
00403   typedef const TriangularView<const MatrixType, Upper> MatrixU;
00404   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
00405   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
00406   static bool inplace_decomposition(MatrixType& m)
00407   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
00408 };
00409 
00410 } // end namespace internal
00411 
00419 template<typename MatrixType, int _UpLo>
00420 template<typename InputType>
00421 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
00422 {
00423   check_template_parameters();
00424 
00425   eigen_assert(a.rows()==a.cols());
00426   const Index size = a.rows();
00427   m_matrix.resize(size, size);
00428   m_matrix = a.derived();
00429 
00430   // Compute matrix L1 norm = max abs column sum.
00431   m_l1_norm = RealScalar(0);
00432   // TODO move this code to SelfAdjointView
00433   for (Index col = 0; col < size; ++col) {
00434     RealScalar abs_col_sum;
00435     if (_UpLo == Lower)
00436       abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
00437     else
00438       abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
00439     if (abs_col_sum > m_l1_norm)
00440       m_l1_norm = abs_col_sum;
00441   }
00442 
00443   m_isInitialized = true;
00444   bool ok = Traits::inplace_decomposition(m_matrix);
00445   m_info = ok ? Success : NumericalIssue;
00446 
00447   return *this;
00448 }
00449 
00455 template<typename _MatrixType, int _UpLo>
00456 template<typename VectorType>
00457 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
00458 {
00459   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
00460   eigen_assert(v.size()==m_matrix.cols());
00461   eigen_assert(m_isInitialized);
00462   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
00463     m_info = NumericalIssue;
00464   else
00465     m_info = Success;
00466 
00467   return *this;
00468 }
00469 
00470 #ifndef EIGEN_PARSED_BY_DOXYGEN
00471 template<typename _MatrixType,int _UpLo>
00472 template<typename RhsType, typename DstType>
00473 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
00474 {
00475   dst = rhs;
00476   solveInPlace(dst);
00477 }
00478 #endif
00479 
00490 template<typename MatrixType, int _UpLo>
00491 template<typename Derived>
00492 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
00493 {
00494   eigen_assert(m_isInitialized && "LLT is not initialized.");
00495   eigen_assert(m_matrix.rows()==bAndX.rows());
00496   matrixL().solveInPlace(bAndX);
00497   matrixU().solveInPlace(bAndX);
00498 }
00499 
00503 template<typename MatrixType, int _UpLo>
00504 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
00505 {
00506   eigen_assert(m_isInitialized && "LLT is not initialized.");
00507   return matrixL() * matrixL().adjoint().toDenseMatrix();
00508 }
00509 
00514 template<typename Derived>
00515 inline const LLT<typename MatrixBase<Derived>::PlainObject>
00516 MatrixBase<Derived>::llt() const
00517 {
00518   return LLT<PlainObject>(derived());
00519 }
00520 
00525 template<typename MatrixType, unsigned int UpLo>
00526 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
00527 SelfAdjointView<MatrixType, UpLo>::llt() const
00528 {
00529   return LLT<PlainObject,UpLo>(m_matrix);
00530 }
00531 
00532 } // end namespace Eigen
00533 
00534 #endif // EIGEN_LLT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends