Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/Eigen/src/Cholesky/LLT.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_LLT_H
0011 #define EIGEN_LLT_H
0012 
0013 namespace Eigen {
0014 
0015 namespace internal{
0016 
0017 template<typename _MatrixType, int _UpLo> struct traits<LLT<_MatrixType, _UpLo> >
0018  : traits<_MatrixType>
0019 {
0020   typedef MatrixXpr XprKind;
0021   typedef SolverStorage StorageKind;
0022   typedef int StorageIndex;
0023   enum { Flags = 0 };
0024 };
0025 
0026 template<typename MatrixType, int UpLo> struct LLT_Traits;
0027 }
0028 
0029 /** \ingroup Cholesky_Module
0030   *
0031   * \class LLT
0032   *
0033   * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
0034   *
0035   * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
0036   * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
0037   *               The other triangular part won't be read.
0038   *
0039   * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
0040   * matrix A such that A = LL^* = U^*U, where L is lower triangular.
0041   *
0042   * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
0043   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
0044   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
0045   * situations like generalised eigen problems with hermitian matrices.
0046   *
0047   * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
0048   * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
0049   * has a solution.
0050   *
0051   * Example: \include LLT_example.cpp
0052   * Output: \verbinclude LLT_example.out
0053   *
0054   * \b Performance: for best performance, it is recommended to use a column-major storage format
0055   * with the Lower triangular part (the default), or, equivalently, a row-major storage format
0056   * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization
0057   * step, and rank-updates can be up to 3 times slower.
0058   *
0059   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0060   *
0061   * Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered.
0062   * Therefore, the strict lower part does not have to store correct values.
0063   *
0064   * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
0065   */
0066 template<typename _MatrixType, int _UpLo> class LLT
0067         : public SolverBase<LLT<_MatrixType, _UpLo> >
0068 {
0069   public:
0070     typedef _MatrixType MatrixType;
0071     typedef SolverBase<LLT> Base;
0072     friend class SolverBase<LLT>;
0073 
0074     EIGEN_GENERIC_PUBLIC_INTERFACE(LLT)
0075     enum {
0076       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0077     };
0078 
0079     enum {
0080       PacketSize = internal::packet_traits<Scalar>::size,
0081       AlignmentMask = int(PacketSize)-1,
0082       UpLo = _UpLo
0083     };
0084 
0085     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
0086 
0087     /**
0088       * \brief Default Constructor.
0089       *
0090       * The default constructor is useful in cases in which the user intends to
0091       * perform decompositions via LLT::compute(const MatrixType&).
0092       */
0093     LLT() : m_matrix(), m_isInitialized(false) {}
0094 
0095     /** \brief Default Constructor with memory preallocation
0096       *
0097       * Like the default constructor but with preallocation of the internal data
0098       * according to the specified problem \a size.
0099       * \sa LLT()
0100       */
0101     explicit LLT(Index size) : m_matrix(size, size),
0102                     m_isInitialized(false) {}
0103 
0104     template<typename InputType>
0105     explicit LLT(const EigenBase<InputType>& matrix)
0106       : m_matrix(matrix.rows(), matrix.cols()),
0107         m_isInitialized(false)
0108     {
0109       compute(matrix.derived());
0110     }
0111 
0112     /** \brief Constructs a LLT factorization from a given matrix
0113       *
0114       * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
0115       * \c MatrixType is a Eigen::Ref.
0116       *
0117       * \sa LLT(const EigenBase&)
0118       */
0119     template<typename InputType>
0120     explicit LLT(EigenBase<InputType>& matrix)
0121       : m_matrix(matrix.derived()),
0122         m_isInitialized(false)
0123     {
0124       compute(matrix.derived());
0125     }
0126 
0127     /** \returns a view of the upper triangular matrix U */
0128     inline typename Traits::MatrixU matrixU() const
0129     {
0130       eigen_assert(m_isInitialized && "LLT is not initialized.");
0131       return Traits::getU(m_matrix);
0132     }
0133 
0134     /** \returns a view of the lower triangular matrix L */
0135     inline typename Traits::MatrixL matrixL() const
0136     {
0137       eigen_assert(m_isInitialized && "LLT is not initialized.");
0138       return Traits::getL(m_matrix);
0139     }
0140 
0141     #ifdef EIGEN_PARSED_BY_DOXYGEN
0142     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
0143       *
0144       * Since this LLT class assumes anyway that the matrix A is invertible, the solution
0145       * theoretically exists and is unique regardless of b.
0146       *
0147       * Example: \include LLT_solve.cpp
0148       * Output: \verbinclude LLT_solve.out
0149       *
0150       * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
0151       */
0152     template<typename Rhs>
0153     inline const Solve<LLT, Rhs>
0154     solve(const MatrixBase<Rhs>& b) const;
0155     #endif
0156 
0157     template<typename Derived>
0158     void solveInPlace(const MatrixBase<Derived> &bAndX) const;
0159 
0160     template<typename InputType>
0161     LLT& compute(const EigenBase<InputType>& matrix);
0162 
0163     /** \returns an estimate of the reciprocal condition number of the matrix of
0164       *  which \c *this is the Cholesky decomposition.
0165       */
0166     RealScalar rcond() const
0167     {
0168       eigen_assert(m_isInitialized && "LLT is not initialized.");
0169       eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
0170       return internal::rcond_estimate_helper(m_l1_norm, *this);
0171     }
0172 
0173     /** \returns the LLT decomposition matrix
0174       *
0175       * TODO: document the storage layout
0176       */
0177     inline const MatrixType& matrixLLT() const
0178     {
0179       eigen_assert(m_isInitialized && "LLT is not initialized.");
0180       return m_matrix;
0181     }
0182 
0183     MatrixType reconstructedMatrix() const;
0184 
0185 
0186     /** \brief Reports whether previous computation was successful.
0187       *
0188       * \returns \c Success if computation was successful,
0189       *          \c NumericalIssue if the matrix.appears not to be positive definite.
0190       */
0191     ComputationInfo info() const
0192     {
0193       eigen_assert(m_isInitialized && "LLT is not initialized.");
0194       return m_info;
0195     }
0196 
0197     /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.
0198       *
0199       * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
0200       * \code x = decomposition.adjoint().solve(b) \endcode
0201       */
0202     const LLT& adjoint() const EIGEN_NOEXCEPT { return *this; };
0203 
0204     inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
0205     inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
0206 
0207     template<typename VectorType>
0208     LLT & rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
0209 
0210     #ifndef EIGEN_PARSED_BY_DOXYGEN
0211     template<typename RhsType, typename DstType>
0212     void _solve_impl(const RhsType &rhs, DstType &dst) const;
0213 
0214     template<bool Conjugate, typename RhsType, typename DstType>
0215     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
0216     #endif
0217 
0218   protected:
0219 
0220     static void check_template_parameters()
0221     {
0222       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0223     }
0224 
0225     /** \internal
0226       * Used to compute and store L
0227       * The strict upper part is not used and even not initialized.
0228       */
0229     MatrixType m_matrix;
0230     RealScalar m_l1_norm;
0231     bool m_isInitialized;
0232     ComputationInfo m_info;
0233 };
0234 
0235 namespace internal {
0236 
0237 template<typename Scalar, int UpLo> struct llt_inplace;
0238 
0239 template<typename MatrixType, typename VectorType>
0240 static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
0241 {
0242   using std::sqrt;
0243   typedef typename MatrixType::Scalar Scalar;
0244   typedef typename MatrixType::RealScalar RealScalar;
0245   typedef typename MatrixType::ColXpr ColXpr;
0246   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
0247   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
0248   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
0249   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
0250 
0251   Index n = mat.cols();
0252   eigen_assert(mat.rows()==n && vec.size()==n);
0253 
0254   TempVectorType temp;
0255 
0256   if(sigma>0)
0257   {
0258     // This version is based on Givens rotations.
0259     // It is faster than the other one below, but only works for updates,
0260     // i.e., for sigma > 0
0261     temp = sqrt(sigma) * vec;
0262 
0263     for(Index i=0; i<n; ++i)
0264     {
0265       JacobiRotation<Scalar> g;
0266       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
0267 
0268       Index rs = n-i-1;
0269       if(rs>0)
0270       {
0271         ColXprSegment x(mat.col(i).tail(rs));
0272         TempVecSegment y(temp.tail(rs));
0273         apply_rotation_in_the_plane(x, y, g);
0274       }
0275     }
0276   }
0277   else
0278   {
0279     temp = vec;
0280     RealScalar beta = 1;
0281     for(Index j=0; j<n; ++j)
0282     {
0283       RealScalar Ljj = numext::real(mat.coeff(j,j));
0284       RealScalar dj = numext::abs2(Ljj);
0285       Scalar wj = temp.coeff(j);
0286       RealScalar swj2 = sigma*numext::abs2(wj);
0287       RealScalar gamma = dj*beta + swj2;
0288 
0289       RealScalar x = dj + swj2/beta;
0290       if (x<=RealScalar(0))
0291         return j;
0292       RealScalar nLjj = sqrt(x);
0293       mat.coeffRef(j,j) = nLjj;
0294       beta += swj2/dj;
0295 
0296       // Update the terms of L
0297       Index rs = n-j-1;
0298       if(rs)
0299       {
0300         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
0301         if(gamma != 0)
0302           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
0303       }
0304     }
0305   }
0306   return -1;
0307 }
0308 
0309 template<typename Scalar> struct llt_inplace<Scalar, Lower>
0310 {
0311   typedef typename NumTraits<Scalar>::Real RealScalar;
0312   template<typename MatrixType>
0313   static Index unblocked(MatrixType& mat)
0314   {
0315     using std::sqrt;
0316 
0317     eigen_assert(mat.rows()==mat.cols());
0318     const Index size = mat.rows();
0319     for(Index k = 0; k < size; ++k)
0320     {
0321       Index rs = size-k-1; // remaining size
0322 
0323       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
0324       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
0325       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
0326 
0327       RealScalar x = numext::real(mat.coeff(k,k));
0328       if (k>0) x -= A10.squaredNorm();
0329       if (x<=RealScalar(0))
0330         return k;
0331       mat.coeffRef(k,k) = x = sqrt(x);
0332       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
0333       if (rs>0) A21 /= x;
0334     }
0335     return -1;
0336   }
0337 
0338   template<typename MatrixType>
0339   static Index blocked(MatrixType& m)
0340   {
0341     eigen_assert(m.rows()==m.cols());
0342     Index size = m.rows();
0343     if(size<32)
0344       return unblocked(m);
0345 
0346     Index blockSize = size/8;
0347     blockSize = (blockSize/16)*16;
0348     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
0349 
0350     for (Index k=0; k<size; k+=blockSize)
0351     {
0352       // partition the matrix:
0353       //       A00 |  -  |  -
0354       // lu  = A10 | A11 |  -
0355       //       A20 | A21 | A22
0356       Index bs = (std::min)(blockSize, size-k);
0357       Index rs = size - k - bs;
0358       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
0359       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
0360       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
0361 
0362       Index ret;
0363       if((ret=unblocked(A11))>=0) return k+ret;
0364       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
0365       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck
0366     }
0367     return -1;
0368   }
0369 
0370   template<typename MatrixType, typename VectorType>
0371   static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
0372   {
0373     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
0374   }
0375 };
0376 
0377 template<typename Scalar> struct llt_inplace<Scalar, Upper>
0378 {
0379   typedef typename NumTraits<Scalar>::Real RealScalar;
0380 
0381   template<typename MatrixType>
0382   static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat)
0383   {
0384     Transpose<MatrixType> matt(mat);
0385     return llt_inplace<Scalar, Lower>::unblocked(matt);
0386   }
0387   template<typename MatrixType>
0388   static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat)
0389   {
0390     Transpose<MatrixType> matt(mat);
0391     return llt_inplace<Scalar, Lower>::blocked(matt);
0392   }
0393   template<typename MatrixType, typename VectorType>
0394   static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
0395   {
0396     Transpose<MatrixType> matt(mat);
0397     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
0398   }
0399 };
0400 
0401 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
0402 {
0403   typedef const TriangularView<const MatrixType, Lower> MatrixL;
0404   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
0405   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
0406   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
0407   static bool inplace_decomposition(MatrixType& m)
0408   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
0409 };
0410 
0411 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
0412 {
0413   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
0414   typedef const TriangularView<const MatrixType, Upper> MatrixU;
0415   static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
0416   static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
0417   static bool inplace_decomposition(MatrixType& m)
0418   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
0419 };
0420 
0421 } // end namespace internal
0422 
0423 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
0424   *
0425   * \returns a reference to *this
0426   *
0427   * Example: \include TutorialLinAlgComputeTwice.cpp
0428   * Output: \verbinclude TutorialLinAlgComputeTwice.out
0429   */
0430 template<typename MatrixType, int _UpLo>
0431 template<typename InputType>
0432 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
0433 {
0434   check_template_parameters();
0435 
0436   eigen_assert(a.rows()==a.cols());
0437   const Index size = a.rows();
0438   m_matrix.resize(size, size);
0439   if (!internal::is_same_dense(m_matrix, a.derived()))
0440     m_matrix = a.derived();
0441 
0442   // Compute matrix L1 norm = max abs column sum.
0443   m_l1_norm = RealScalar(0);
0444   // TODO move this code to SelfAdjointView
0445   for (Index col = 0; col < size; ++col) {
0446     RealScalar abs_col_sum;
0447     if (_UpLo == Lower)
0448       abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
0449     else
0450       abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
0451     if (abs_col_sum > m_l1_norm)
0452       m_l1_norm = abs_col_sum;
0453   }
0454 
0455   m_isInitialized = true;
0456   bool ok = Traits::inplace_decomposition(m_matrix);
0457   m_info = ok ? Success : NumericalIssue;
0458 
0459   return *this;
0460 }
0461 
0462 /** Performs a rank one update (or dowdate) of the current decomposition.
0463   * If A = LL^* before the rank one update,
0464   * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
0465   * of same dimension.
0466   */
0467 template<typename _MatrixType, int _UpLo>
0468 template<typename VectorType>
0469 LLT<_MatrixType,_UpLo> & LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
0470 {
0471   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
0472   eigen_assert(v.size()==m_matrix.cols());
0473   eigen_assert(m_isInitialized);
0474   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
0475     m_info = NumericalIssue;
0476   else
0477     m_info = Success;
0478 
0479   return *this;
0480 }
0481 
0482 #ifndef EIGEN_PARSED_BY_DOXYGEN
0483 template<typename _MatrixType,int _UpLo>
0484 template<typename RhsType, typename DstType>
0485 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
0486 {
0487   _solve_impl_transposed<true>(rhs, dst);
0488 }
0489 
0490 template<typename _MatrixType,int _UpLo>
0491 template<bool Conjugate, typename RhsType, typename DstType>
0492 void LLT<_MatrixType,_UpLo>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
0493 {
0494     dst = rhs;
0495 
0496     matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst);
0497     matrixU().template conjugateIf<!Conjugate>().solveInPlace(dst);
0498 }
0499 #endif
0500 
0501 /** \internal use x = llt_object.solve(x);
0502   *
0503   * This is the \em in-place version of solve().
0504   *
0505   * \param bAndX represents both the right-hand side matrix b and result x.
0506   *
0507   * This version avoids a copy when the right hand side matrix b is not needed anymore.
0508   *
0509   * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
0510   * This function will const_cast it, so constness isn't honored here.
0511   *
0512   * \sa LLT::solve(), MatrixBase::llt()
0513   */
0514 template<typename MatrixType, int _UpLo>
0515 template<typename Derived>
0516 void LLT<MatrixType,_UpLo>::solveInPlace(const MatrixBase<Derived> &bAndX) const
0517 {
0518   eigen_assert(m_isInitialized && "LLT is not initialized.");
0519   eigen_assert(m_matrix.rows()==bAndX.rows());
0520   matrixL().solveInPlace(bAndX);
0521   matrixU().solveInPlace(bAndX);
0522 }
0523 
0524 /** \returns the matrix represented by the decomposition,
0525  * i.e., it returns the product: L L^*.
0526  * This function is provided for debug purpose. */
0527 template<typename MatrixType, int _UpLo>
0528 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
0529 {
0530   eigen_assert(m_isInitialized && "LLT is not initialized.");
0531   return matrixL() * matrixL().adjoint().toDenseMatrix();
0532 }
0533 
0534 /** \cholesky_module
0535   * \returns the LLT decomposition of \c *this
0536   * \sa SelfAdjointView::llt()
0537   */
0538 template<typename Derived>
0539 inline const LLT<typename MatrixBase<Derived>::PlainObject>
0540 MatrixBase<Derived>::llt() const
0541 {
0542   return LLT<PlainObject>(derived());
0543 }
0544 
0545 /** \cholesky_module
0546   * \returns the LLT decomposition of \c *this
0547   * \sa SelfAdjointView::llt()
0548   */
0549 template<typename MatrixType, unsigned int UpLo>
0550 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
0551 SelfAdjointView<MatrixType, UpLo>::llt() const
0552 {
0553   return LLT<PlainObject,UpLo>(m_matrix);
0554 }
0555 
0556 } // end namespace Eigen
0557 
0558 #endif // EIGEN_LLT_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!