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
13namespace Eigen {
14
15namespace internal{
16template<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 * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
26 * \tparam _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 * \b Performance: for best performance, it is recommended to use a column-major storage format
45 * with the Lower triangular part (the default), or, equivalently, a row-major storage format
46 * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization
47 * step, and rank-updates can be up to 3 times slower.
48 *
49 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
50 *
51 * Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered.
52 * Therefore, the strict lower part does not have to store correct values.
53 *
54 * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
55 */
56template<typename _MatrixType, int _UpLo> class LLT
57{
58 public:
59 typedef _MatrixType MatrixType;
60 enum {
61 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
62 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64 };
65 typedef typename MatrixType::Scalar Scalar;
66 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
67 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
68 typedef typename MatrixType::StorageIndex StorageIndex;
69
70 enum {
71 PacketSize = internal::packet_traits<Scalar>::size,
72 AlignmentMask = int(PacketSize)-1,
73 UpLo = _UpLo
74 };
75
76 typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
77
78 /**
79 * \brief Default Constructor.
80 *
81 * The default constructor is useful in cases in which the user intends to
82 * perform decompositions via LLT::compute(const MatrixType&).
83 */
84 LLT() : m_matrix(), m_isInitialized(false) {}
85
86 /** \brief Default Constructor with memory preallocation
87 *
88 * Like the default constructor but with preallocation of the internal data
89 * according to the specified problem \a size.
90 * \sa LLT()
91 */
92 explicit LLT(Index size) : m_matrix(size, size),
93 m_isInitialized(false) {}
94
95 template<typename InputType>
96 explicit LLT(const EigenBase<InputType>& matrix)
97 : m_matrix(matrix.rows(), matrix.cols()),
98 m_isInitialized(false)
99 {
100 compute(matrix.derived());
101 }
102
103 /** \brief Constructs a LDLT factorization from a given matrix
104 *
105 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
106 * \c MatrixType is a Eigen::Ref.
107 *
108 * \sa LLT(const EigenBase&)
109 */
110 template<typename InputType>
111 explicit LLT(EigenBase<InputType>& matrix)
112 : m_matrix(matrix.derived()),
113 m_isInitialized(false)
114 {
115 compute(matrix.derived());
116 }
117
118 /** \returns a view of the upper triangular matrix U */
119 inline typename Traits::MatrixU matrixU() const
120 {
121 eigen_assert(m_isInitialized && "LLT is not initialized.");
122 return Traits::getU(m_matrix);
123 }
124
125 /** \returns a view of the lower triangular matrix L */
126 inline typename Traits::MatrixL matrixL() const
127 {
128 eigen_assert(m_isInitialized && "LLT is not initialized.");
129 return Traits::getL(m_matrix);
130 }
131
132 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
133 *
134 * Since this LLT class assumes anyway that the matrix A is invertible, the solution
135 * theoretically exists and is unique regardless of b.
136 *
137 * Example: \include LLT_solve.cpp
138 * Output: \verbinclude LLT_solve.out
139 *
140 * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt()
141 */
142 template<typename Rhs>
143 inline const Solve<LLT, Rhs>
144 solve(const MatrixBase<Rhs>& b) const
145 {
146 eigen_assert(m_isInitialized && "LLT is not initialized.");
147 eigen_assert(m_matrix.rows()==b.rows()
148 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
149 return Solve<LLT, Rhs>(*this, b.derived());
150 }
151
152 template<typename Derived>
153 void solveInPlace(const MatrixBase<Derived> &bAndX) const;
154
155 template<typename InputType>
156 LLT& compute(const EigenBase<InputType>& matrix);
157
158 /** \returns an estimate of the reciprocal condition number of the matrix of
159 * which \c *this is the Cholesky decomposition.
160 */
161 RealScalar rcond() const
162 {
163 eigen_assert(m_isInitialized && "LLT is not initialized.");
164 eigen_assert(m_info == Success && "LLT failed because matrix appears to be negative");
165 return internal::rcond_estimate_helper(m_l1_norm, *this);
166 }
167
168 /** \returns the LLT decomposition matrix
169 *
170 * TODO: document the storage layout
171 */
172 inline const MatrixType& matrixLLT() const
173 {
174 eigen_assert(m_isInitialized && "LLT is not initialized.");
175 return m_matrix;
176 }
177
178 MatrixType reconstructedMatrix() const;
179
180
181 /** \brief Reports whether previous computation was successful.
182 *
183 * \returns \c Success if computation was succesful,
184 * \c NumericalIssue if the matrix.appears not to be positive definite.
185 */
186 ComputationInfo info() const
187 {
188 eigen_assert(m_isInitialized && "LLT is not initialized.");
189 return m_info;
190 }
191
192 /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.
193 *
194 * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as:
195 * \code x = decomposition.adjoint().solve(b) \endcode
196 */
197 const LLT& adjoint() const { return *this; };
198
199 inline Index rows() const { return m_matrix.rows(); }
200 inline Index cols() const { return m_matrix.cols(); }
201
202 template<typename VectorType>
203 LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
204
205 #ifndef EIGEN_PARSED_BY_DOXYGEN
206 template<typename RhsType, typename DstType>
207 EIGEN_DEVICE_FUNC
208 void _solve_impl(const RhsType &rhs, DstType &dst) const;
209 #endif
210
211 protected:
212
213 static void check_template_parameters()
214 {
215 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
216 }
217
218 /** \internal
219 * Used to compute and store L
220 * The strict upper part is not used and even not initialized.
221 */
222 MatrixType m_matrix;
223 RealScalar m_l1_norm;
224 bool m_isInitialized;
225 ComputationInfo m_info;
226};
227
228namespace internal {
229
230template<typename Scalar, int UpLo> struct llt_inplace;
231
232template<typename MatrixType, typename VectorType>
233static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
234{
235 using std::sqrt;
236 typedef typename MatrixType::Scalar Scalar;
237 typedef typename MatrixType::RealScalar RealScalar;
238 typedef typename MatrixType::ColXpr ColXpr;
239 typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
240 typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
241 typedef Matrix<Scalar,Dynamic,1> TempVectorType;
242 typedef typename TempVectorType::SegmentReturnType TempVecSegment;
243
244 Index n = mat.cols();
245 eigen_assert(mat.rows()==n && vec.size()==n);
246
247 TempVectorType temp;
248
249 if(sigma>0)
250 {
251 // This version is based on Givens rotations.
252 // It is faster than the other one below, but only works for updates,
253 // i.e., for sigma > 0
254 temp = sqrt(sigma) * vec;
255
256 for(Index i=0; i<n; ++i)
257 {
258 JacobiRotation<Scalar> g;
259 g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
260
261 Index rs = n-i-1;
262 if(rs>0)
263 {
264 ColXprSegment x(mat.col(i).tail(rs));
265 TempVecSegment y(temp.tail(rs));
266 apply_rotation_in_the_plane(x, y, g);
267 }
268 }
269 }
270 else
271 {
272 temp = vec;
273 RealScalar beta = 1;
274 for(Index j=0; j<n; ++j)
275 {
276 RealScalar Ljj = numext::real(mat.coeff(j,j));
277 RealScalar dj = numext::abs2(Ljj);
278 Scalar wj = temp.coeff(j);
279 RealScalar swj2 = sigma*numext::abs2(wj);
280 RealScalar gamma = dj*beta + swj2;
281
282 RealScalar x = dj + swj2/beta;
283 if (x<=RealScalar(0))
284 return j;
285 RealScalar nLjj = sqrt(x);
286 mat.coeffRef(j,j) = nLjj;
287 beta += swj2/dj;
288
289 // Update the terms of L
290 Index rs = n-j-1;
291 if(rs)
292 {
293 temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
294 if(gamma != 0)
295 mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
296 }
297 }
298 }
299 return -1;
300}
301
302template<typename Scalar> struct llt_inplace<Scalar, Lower>
303{
304 typedef typename NumTraits<Scalar>::Real RealScalar;
305 template<typename MatrixType>
306 static Index unblocked(MatrixType& mat)
307 {
308 using std::sqrt;
309
310 eigen_assert(mat.rows()==mat.cols());
311 const Index size = mat.rows();
312 for(Index k = 0; k < size; ++k)
313 {
314 Index rs = size-k-1; // remaining size
315
316 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
317 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
318 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
319
320 RealScalar x = numext::real(mat.coeff(k,k));
321 if (k>0) x -= A10.squaredNorm();
322 if (x<=RealScalar(0))
323 return k;
324 mat.coeffRef(k,k) = x = sqrt(x);
325 if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
326 if (rs>0) A21 /= x;
327 }
328 return -1;
329 }
330
331 template<typename MatrixType>
332 static Index blocked(MatrixType& m)
333 {
334 eigen_assert(m.rows()==m.cols());
335 Index size = m.rows();
336 if(size<32)
337 return unblocked(m);
338
339 Index blockSize = size/8;
340 blockSize = (blockSize/16)*16;
341 blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
342
343 for (Index k=0; k<size; k+=blockSize)
344 {
345 // partition the matrix:
346 // A00 | - | -
347 // lu = A10 | A11 | -
348 // A20 | A21 | A22
349 Index bs = (std::min)(blockSize, size-k);
350 Index rs = size - k - bs;
351 Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
352 Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
353 Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
354
355 Index ret;
356 if((ret=unblocked(A11))>=0) return k+ret;
357 if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
358 if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,typename NumTraits<RealScalar>::Literal(-1)); // bottleneck
359 }
360 return -1;
361 }
362
363 template<typename MatrixType, typename VectorType>
364 static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
365 {
366 return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
367 }
368};
369
370template<typename Scalar> struct llt_inplace<Scalar, Upper>
371{
372 typedef typename NumTraits<Scalar>::Real RealScalar;
373
374 template<typename MatrixType>
375 static EIGEN_STRONG_INLINE Index unblocked(MatrixType& mat)
376 {
377 Transpose<MatrixType> matt(mat);
378 return llt_inplace<Scalar, Lower>::unblocked(matt);
379 }
380 template<typename MatrixType>
381 static EIGEN_STRONG_INLINE Index blocked(MatrixType& mat)
382 {
383 Transpose<MatrixType> matt(mat);
384 return llt_inplace<Scalar, Lower>::blocked(matt);
385 }
386 template<typename MatrixType, typename VectorType>
387 static Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
388 {
389 Transpose<MatrixType> matt(mat);
390 return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
391 }
392};
393
394template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
395{
396 typedef const TriangularView<const MatrixType, Lower> MatrixL;
397 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
398 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
399 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
400 static bool inplace_decomposition(MatrixType& m)
401 { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
402};
403
404template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
405{
406 typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
407 typedef const TriangularView<const MatrixType, Upper> MatrixU;
408 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); }
409 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); }
410 static bool inplace_decomposition(MatrixType& m)
411 { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
412};
413
414} // end namespace internal
415
416/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
417 *
418 * \returns a reference to *this
419 *
420 * Example: \include TutorialLinAlgComputeTwice.cpp
421 * Output: \verbinclude TutorialLinAlgComputeTwice.out
422 */
423template<typename MatrixType, int _UpLo>
424template<typename InputType>
425LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a)
426{
427 check_template_parameters();
428
429 eigen_assert(a.rows()==a.cols());
430 const Index size = a.rows();
431 m_matrix.resize(size, size);
432 if (!internal::is_same_dense(m_matrix, a.derived()))
433 m_matrix = a.derived();
434
435 // Compute matrix L1 norm = max abs column sum.
436 m_l1_norm = RealScalar(0);
437 // TODO move this code to SelfAdjointView
438 for (Index col = 0; col < size; ++col) {
439 RealScalar abs_col_sum;
440 if (_UpLo == Lower)
441 abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>();
442 else
443 abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>();
444 if (abs_col_sum > m_l1_norm)
445 m_l1_norm = abs_col_sum;
446 }
447
448 m_isInitialized = true;
449 bool ok = Traits::inplace_decomposition(m_matrix);
450 m_info = ok ? Success : NumericalIssue;
451
452 return *this;
453}
454
455/** Performs a rank one update (or dowdate) of the current decomposition.
456 * If A = LL^* before the rank one update,
457 * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
458 * of same dimension.
459 */
460template<typename _MatrixType, int _UpLo>
461template<typename VectorType>
462LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
463{
464 EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
465 eigen_assert(v.size()==m_matrix.cols());
466 eigen_assert(m_isInitialized);
467 if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
468 m_info = NumericalIssue;
469 else
470 m_info = Success;
471
472 return *this;
473}
474
475#ifndef EIGEN_PARSED_BY_DOXYGEN
476template<typename _MatrixType,int _UpLo>
477template<typename RhsType, typename DstType>
478void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
479{
480 dst = rhs;
481 solveInPlace(dst);
482}
483#endif
484
485/** \internal use x = llt_object.solve(x);
486 *
487 * This is the \em in-place version of solve().
488 *
489 * \param bAndX represents both the right-hand side matrix b and result x.
490 *
491 * This version avoids a copy when the right hand side matrix b is not needed anymore.
492 *
493 * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
494 * This function will const_cast it, so constness isn't honored here.
495 *
496 * \sa LLT::solve(), MatrixBase::llt()
497 */
498template<typename MatrixType, int _UpLo>
499template<typename Derived>
500void LLT<MatrixType,_UpLo>::solveInPlace(const MatrixBase<Derived> &bAndX) const
501{
502 eigen_assert(m_isInitialized && "LLT is not initialized.");
503 eigen_assert(m_matrix.rows()==bAndX.rows());
504 matrixL().solveInPlace(bAndX);
505 matrixU().solveInPlace(bAndX);
506}
507
508/** \returns the matrix represented by the decomposition,
509 * i.e., it returns the product: L L^*.
510 * This function is provided for debug purpose. */
511template<typename MatrixType, int _UpLo>
512MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
513{
514 eigen_assert(m_isInitialized && "LLT is not initialized.");
515 return matrixL() * matrixL().adjoint().toDenseMatrix();
516}
517
518/** \cholesky_module
519 * \returns the LLT decomposition of \c *this
520 * \sa SelfAdjointView::llt()
521 */
522template<typename Derived>
523inline const LLT<typename MatrixBase<Derived>::PlainObject>
524MatrixBase<Derived>::llt() const
525{
526 return LLT<PlainObject>(derived());
527}
528
529/** \cholesky_module
530 * \returns the LLT decomposition of \c *this
531 * \sa SelfAdjointView::llt()
532 */
533template<typename MatrixType, unsigned int UpLo>
534inline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
535SelfAdjointView<MatrixType, UpLo>::llt() const
536{
537 return LLT<PlainObject,UpLo>(m_matrix);
538}
539
540} // end namespace Eigen
541
542#endif // EIGEN_LLT_H
543