1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_LLT_H
11 #define EIGEN_LLT_H
12 
13 namespace Eigen {
14 
15 namespace internal{
16 template<typename MatrixType, int UpLo> struct LLT_Traits;
17 }
18 
19 /** \ingroup Cholesky_Module
20   *
21   * \class LLT
22   *
23   * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
24   *
25   * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
26   * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
27   *             The other triangular part won't be read.
28   *
29   * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
30   * matrix A such that A = LL^* = U^*U, where L is lower triangular.
31   *
32   * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
33   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
34   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
35   * situations like generalised eigen problems with hermitian matrices.
36   *
37   * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
38   * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
39   * has a solution.
40   *
41   * Example: \include LLT_example.cpp
42   * Output: \verbinclude LLT_example.out
43   *
44   * \sa MatrixBase::llt(), class LDLT
45   */
46  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
47   * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
48   * the strict lower part does not have to store correct values.
49   */
50 template<typename _MatrixType, int _UpLo> class LLT
51 {
52   public:
53     typedef _MatrixType MatrixType;
54     enum {
55       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
57       Options = MatrixType::Options,
58       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
59     };
60     typedef typename MatrixType::Scalar Scalar;
61     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
62     typedef typename MatrixType::Index Index;
63 
64     enum {
65       PacketSize = internal::packet_traits<Scalar>::size,
66       AlignmentMask = int(PacketSize)-1,
67       UpLo = _UpLo
68     };
69 
70     typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
71 
72     /**
73       * \brief Default Constructor.
74       *
75       * The default constructor is useful in cases in which the user intends to
76       * perform decompositions via LLT::compute(const MatrixType&).
77       */
LLT()78     LLT() : m_matrix(), m_isInitialized(false) {}
79 
80     /** \brief Default Constructor with memory preallocation
81       *
82       * Like the default constructor but with preallocation of the internal data
83       * according to the specified problem \a size.
84       * \sa LLT()
85       */
LLT(Index size)86     LLT(Index size) : m_matrix(size, size),
87                     m_isInitialized(false) {}
88 
LLT(const MatrixType & matrix)89     LLT(const MatrixType& matrix)
90       : m_matrix(matrix.rows(), matrix.cols()),
91         m_isInitialized(false)
92     {
93       compute(matrix);
94     }
95 
96     /** \returns a view of the upper triangular matrix U */
matrixU()97     inline typename Traits::MatrixU matrixU() const
98     {
99       eigen_assert(m_isInitialized && "LLT is not initialized.");
100       return Traits::getU(m_matrix);
101     }
102 
103     /** \returns a view of the lower triangular matrix L */
matrixL()104     inline typename Traits::MatrixL matrixL() const
105     {
106       eigen_assert(m_isInitialized && "LLT is not initialized.");
107       return Traits::getL(m_matrix);
108     }
109 
110     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
111       *
112       * Since this LLT class assumes anyway that the matrix A is invertible, the solution
113       * theoretically exists and is unique regardless of b.
114       *
115       * Example: \include LLT_solve.cpp
116       * Output: \verbinclude LLT_solve.out
117       *
118       * \sa solveInPlace(), MatrixBase::llt()
119       */
120     template<typename Rhs>
121     inline const internal::solve_retval<LLT, Rhs>
solve(const MatrixBase<Rhs> & b)122     solve(const MatrixBase<Rhs>& b) const
123     {
124       eigen_assert(m_isInitialized && "LLT is not initialized.");
125       eigen_assert(m_matrix.rows()==b.rows()
126                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
127       return internal::solve_retval<LLT, Rhs>(*this, b.derived());
128     }
129 
130     #ifdef EIGEN2_SUPPORT
131     template<typename OtherDerived, typename ResultType>
solve(const MatrixBase<OtherDerived> & b,ResultType * result)132     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
133     {
134       *result = this->solve(b);
135       return true;
136     }
137 
isPositiveDefinite()138     bool isPositiveDefinite() const { return true; }
139     #endif
140 
141     template<typename Derived>
142     void solveInPlace(MatrixBase<Derived> &bAndX) const;
143 
144     LLT& compute(const MatrixType& matrix);
145 
146     /** \returns the LLT decomposition matrix
147       *
148       * TODO: document the storage layout
149       */
matrixLLT()150     inline const MatrixType& matrixLLT() const
151     {
152       eigen_assert(m_isInitialized && "LLT is not initialized.");
153       return m_matrix;
154     }
155 
156     MatrixType reconstructedMatrix() const;
157 
158 
159     /** \brief Reports whether previous computation was successful.
160       *
161       * \returns \c Success if computation was succesful,
162       *          \c NumericalIssue if the matrix.appears to be negative.
163       */
info()164     ComputationInfo info() const
165     {
166       eigen_assert(m_isInitialized && "LLT is not initialized.");
167       return m_info;
168     }
169 
rows()170     inline Index rows() const { return m_matrix.rows(); }
cols()171     inline Index cols() const { return m_matrix.cols(); }
172 
173     template<typename VectorType>
174     LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
175 
176   protected:
177     /** \internal
178       * Used to compute and store L
179       * The strict upper part is not used and even not initialized.
180       */
181     MatrixType m_matrix;
182     bool m_isInitialized;
183     ComputationInfo m_info;
184 };
185 
186 namespace internal {
187 
188 template<typename Scalar, int UpLo> struct llt_inplace;
189 
190 template<typename MatrixType, typename VectorType>
llt_rank_update_lower(MatrixType & mat,const VectorType & vec,const typename MatrixType::RealScalar & sigma)191 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
192 {
193   using std::sqrt;
194   typedef typename MatrixType::Scalar Scalar;
195   typedef typename MatrixType::RealScalar RealScalar;
196   typedef typename MatrixType::Index Index;
197   typedef typename MatrixType::ColXpr ColXpr;
198   typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
199   typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
200   typedef Matrix<Scalar,Dynamic,1> TempVectorType;
201   typedef typename TempVectorType::SegmentReturnType TempVecSegment;
202 
203   Index n = mat.cols();
204   eigen_assert(mat.rows()==n && vec.size()==n);
205 
206   TempVectorType temp;
207 
208   if(sigma>0)
209   {
210     // This version is based on Givens rotations.
211     // It is faster than the other one below, but only works for updates,
212     // i.e., for sigma > 0
213     temp = sqrt(sigma) * vec;
214 
215     for(Index i=0; i<n; ++i)
216     {
217       JacobiRotation<Scalar> g;
218       g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
219 
220       Index rs = n-i-1;
221       if(rs>0)
222       {
223         ColXprSegment x(mat.col(i).tail(rs));
224         TempVecSegment y(temp.tail(rs));
225         apply_rotation_in_the_plane(x, y, g);
226       }
227     }
228   }
229   else
230   {
231     temp = vec;
232     RealScalar beta = 1;
233     for(Index j=0; j<n; ++j)
234     {
235       RealScalar Ljj = numext::real(mat.coeff(j,j));
236       RealScalar dj = numext::abs2(Ljj);
237       Scalar wj = temp.coeff(j);
238       RealScalar swj2 = sigma*numext::abs2(wj);
239       RealScalar gamma = dj*beta + swj2;
240 
241       RealScalar x = dj + swj2/beta;
242       if (x<=RealScalar(0))
243         return j;
244       RealScalar nLjj = sqrt(x);
245       mat.coeffRef(j,j) = nLjj;
246       beta += swj2/dj;
247 
248       // Update the terms of L
249       Index rs = n-j-1;
250       if(rs)
251       {
252         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
253         if(gamma != 0)
254           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
255       }
256     }
257   }
258   return -1;
259 }
260 
261 template<typename Scalar> struct llt_inplace<Scalar, Lower>
262 {
263   typedef typename NumTraits<Scalar>::Real RealScalar;
264   template<typename MatrixType>
265   static typename MatrixType::Index unblocked(MatrixType& mat)
266   {
267     using std::sqrt;
268     typedef typename MatrixType::Index Index;
269 
270     eigen_assert(mat.rows()==mat.cols());
271     const Index size = mat.rows();
272     for(Index k = 0; k < size; ++k)
273     {
274       Index rs = size-k-1; // remaining size
275 
276       Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
277       Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
278       Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
279 
280       RealScalar x = numext::real(mat.coeff(k,k));
281       if (k>0) x -= A10.squaredNorm();
282       if (x<=RealScalar(0))
283         return k;
284       mat.coeffRef(k,k) = x = sqrt(x);
285       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
286       if (rs>0) A21 *= RealScalar(1)/x;
287     }
288     return -1;
289   }
290 
291   template<typename MatrixType>
292   static typename MatrixType::Index blocked(MatrixType& m)
293   {
294     typedef typename MatrixType::Index Index;
295     eigen_assert(m.rows()==m.cols());
296     Index size = m.rows();
297     if(size<32)
298       return unblocked(m);
299 
300     Index blockSize = size/8;
301     blockSize = (blockSize/16)*16;
302     blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
303 
304     for (Index k=0; k<size; k+=blockSize)
305     {
306       // partition the matrix:
307       //       A00 |  -  |  -
308       // lu  = A10 | A11 |  -
309       //       A20 | A21 | A22
310       Index bs = (std::min)(blockSize, size-k);
311       Index rs = size - k - bs;
312       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
313       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
314       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
315 
316       Index ret;
317       if((ret=unblocked(A11))>=0) return k+ret;
318       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
319       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
320     }
321     return -1;
322   }
323 
324   template<typename MatrixType, typename VectorType>
325   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
326   {
327     return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
328   }
329 };
330 
331 template<typename Scalar> struct llt_inplace<Scalar, Upper>
332 {
333   typedef typename NumTraits<Scalar>::Real RealScalar;
334 
335   template<typename MatrixType>
336   static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
337   {
338     Transpose<MatrixType> matt(mat);
339     return llt_inplace<Scalar, Lower>::unblocked(matt);
340   }
341   template<typename MatrixType>
342   static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
343   {
344     Transpose<MatrixType> matt(mat);
345     return llt_inplace<Scalar, Lower>::blocked(matt);
346   }
347   template<typename MatrixType, typename VectorType>
348   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
349   {
350     Transpose<MatrixType> matt(mat);
351     return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
352   }
353 };
354 
355 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
356 {
357   typedef const TriangularView<const MatrixType, Lower> MatrixL;
358   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
359   static inline MatrixL getL(const MatrixType& m) { return m; }
360   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
361   static bool inplace_decomposition(MatrixType& m)
362   { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
363 };
364 
365 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
366 {
367   typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
368   typedef const TriangularView<const MatrixType, Upper> MatrixU;
369   static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
370   static inline MatrixU getU(const MatrixType& m) { return m; }
371   static bool inplace_decomposition(MatrixType& m)
372   { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
373 };
374 
375 } // end namespace internal
376 
377 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
378   *
379   * \returns a reference to *this
380   *
381   * Example: \include TutorialLinAlgComputeTwice.cpp
382   * Output: \verbinclude TutorialLinAlgComputeTwice.out
383   */
384 template<typename MatrixType, int _UpLo>
385 LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
386 {
387   eigen_assert(a.rows()==a.cols());
388   const Index size = a.rows();
389   m_matrix.resize(size, size);
390   m_matrix = a;
391 
392   m_isInitialized = true;
393   bool ok = Traits::inplace_decomposition(m_matrix);
394   m_info = ok ? Success : NumericalIssue;
395 
396   return *this;
397 }
398 
399 /** Performs a rank one update (or dowdate) of the current decomposition.
400   * If A = LL^* before the rank one update,
401   * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
402   * of same dimension.
403   */
404 template<typename _MatrixType, int _UpLo>
405 template<typename VectorType>
406 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
407 {
408   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
409   eigen_assert(v.size()==m_matrix.cols());
410   eigen_assert(m_isInitialized);
411   if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
412     m_info = NumericalIssue;
413   else
414     m_info = Success;
415 
416   return *this;
417 }
418 
419 namespace internal {
420 template<typename _MatrixType, int UpLo, typename Rhs>
421 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
422   : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
423 {
424   typedef LLT<_MatrixType,UpLo> LLTType;
425   EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
426 
427   template<typename Dest> void evalTo(Dest& dst) const
428   {
429     dst = rhs();
430     dec().solveInPlace(dst);
431   }
432 };
433 }
434 
435 /** \internal use x = llt_object.solve(x);
436   *
437   * This is the \em in-place version of solve().
438   *
439   * \param bAndX represents both the right-hand side matrix b and result x.
440   *
441   * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
442   *
443   * This version avoids a copy when the right hand side matrix b is not
444   * needed anymore.
445   *
446   * \sa LLT::solve(), MatrixBase::llt()
447   */
448 template<typename MatrixType, int _UpLo>
449 template<typename Derived>
450 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
451 {
452   eigen_assert(m_isInitialized && "LLT is not initialized.");
453   eigen_assert(m_matrix.rows()==bAndX.rows());
454   matrixL().solveInPlace(bAndX);
455   matrixU().solveInPlace(bAndX);
456 }
457 
458 /** \returns the matrix represented by the decomposition,
459  * i.e., it returns the product: L L^*.
460  * This function is provided for debug purpose. */
461 template<typename MatrixType, int _UpLo>
462 MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
463 {
464   eigen_assert(m_isInitialized && "LLT is not initialized.");
465   return matrixL() * matrixL().adjoint().toDenseMatrix();
466 }
467 
468 /** \cholesky_module
469   * \returns the LLT decomposition of \c *this
470   */
471 template<typename Derived>
472 inline const LLT<typename MatrixBase<Derived>::PlainObject>
473 MatrixBase<Derived>::llt() const
474 {
475   return LLT<PlainObject>(derived());
476 }
477 
478 /** \cholesky_module
479   * \returns the LLT decomposition of \c *this
480   */
481 template<typename MatrixType, unsigned int UpLo>
482 inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
483 SelfAdjointView<MatrixType, UpLo>::llt() const
484 {
485   return LLT<PlainObject,UpLo>(m_matrix);
486 }
487 
488 } // end namespace Eigen
489 
490 #endif // EIGEN_LLT_H
491