Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.3.71
LLT.h
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 
56 template<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;
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 
84  LLT() : m_matrix(), m_isInitialized(false) {}
85 
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 
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 
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 
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 
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 
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 
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 
187  {
188  eigen_assert(m_isInitialized && "LLT is not initialized.");
189  return m_info;
190  }
191 
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 
222  MatrixType m_matrix;
223  RealScalar m_l1_norm;
224  bool m_isInitialized;
225  ComputationInfo m_info;
226 };
227 
228 namespace internal {
229 
230 template<typename Scalar, int UpLo> struct llt_inplace;
231 
232 template<typename MatrixType, typename VectorType>
233 static 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 
302 template<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 
370 template<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 
394 template<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 
404 template<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 
423 template<typename MatrixType, int _UpLo>
424 template<typename InputType>
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 
460 template<typename _MatrixType, int _UpLo>
461 template<typename VectorType>
462 LLT<_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
476 template<typename _MatrixType,int _UpLo>
477 template<typename RhsType, typename DstType>
478 void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
479 {
480  dst = rhs;
481  solveInPlace(dst);
482 }
483 #endif
484 
498 template<typename MatrixType, int _UpLo>
499 template<typename Derived>
500 void 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 
511 template<typename MatrixType, int _UpLo>
513 {
514  eigen_assert(m_isInitialized && "LLT is not initialized.");
515  return matrixL() * matrixL().adjoint().toDenseMatrix();
516 }
517 
522 template<typename Derived>
525 {
526  return LLT<PlainObject>(derived());
527 }
528 
533 template<typename MatrixType, unsigned int UpLo>
536 {
537  return LLT<PlainObject,UpLo>(m_matrix);
538 }
539 
540 } // end namespace Eigen
541 
542 #endif // EIGEN_LLT_H
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:434
Eigen::LLT::adjoint
const LLT & adjoint() const
Definition: LLT.h:197
Eigen::EigenBase::cols
Index cols() const
Definition: EigenBase.h:62
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:309
Eigen::LLT::rcond
RealScalar rcond() const
Definition: LLT.h:161
Eigen::MatrixBase::llt
const LLT< PlainObject > llt() const
Definition: LLT.h:524
Eigen::sqrt
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Eigen::EigenBase
Definition: EigenBase.h:29
Eigen::DenseCoeffsBase< Derived, DirectWriteAccessors >::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::SelfAdjointView::llt
const LLT< PlainObject, UpLo > llt() const
Definition: LLT.h:535
Eigen::Upper
@ Upper
Definition: Constants.h:206
Eigen::Success
@ Success
Definition: Constants.h:432
Eigen::LLT::matrixL
Traits::MatrixL matrixL() const
Definition: LLT.h:126
Eigen::LLT::reconstructedMatrix
MatrixType reconstructedMatrix() const
Definition: LLT.h:512
Eigen::Lower
@ Lower
Definition: Constants.h:204
Eigen::LLT::Index
Eigen::Index Index
Definition: LLT.h:67
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::LLT
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:56
Eigen::Solve
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Eigen::LLT::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: LLT.h:186
Eigen::LLT::LLT
LLT(Index size)
Default Constructor with memory preallocation.
Definition: LLT.h:92
Eigen::LLT::LLT
LLT(EigenBase< InputType > &matrix)
Constructs a LDLT factorization from a given matrix.
Definition: LLT.h:111
Eigen::LLT::matrixU
Traits::MatrixU matrixU() const
Definition: LLT.h:119
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::EigenBase::rows
Index rows() const
Definition: EigenBase.h:59
Eigen::LLT::solve
const Solve< LLT, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: LLT.h:144
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:430
Eigen::LLT::matrixLLT
const MatrixType & matrixLLT() const
Definition: LLT.h:172
Eigen::LLT::LLT
LLT()
Default Constructor.
Definition: LLT.h:84
Eigen::DenseCoeffsBase< Derived, DirectWriteAccessors >::rows
Index rows() const
Definition: EigenBase.h:59
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:150
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33