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 
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 
78  LLT() : m_matrix(), m_isInitialized(false) {}
79 
86  LLT(Index size) : m_matrix(size, size),
87  m_isInitialized(false) {}
88 
89  LLT(const MatrixType& matrix)
90  : m_matrix(matrix.rows(), matrix.cols()),
91  m_isInitialized(false)
92  {
93  compute(matrix);
94  }
95 
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 
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 
120  template<typename Rhs>
121  inline const internal::solve_retval<LLT, Rhs>
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>
132  bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
133  {
134  *result = this->solve(b);
135  return true;
136  }
137 
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 
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 
165  {
166  eigen_assert(m_isInitialized && "LLT is not initialized.");
167  return m_info;
168  }
169 
170  inline Index rows() const { return m_matrix.rows(); }
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:
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>
191 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
192 {
193  typedef typename MatrixType::Scalar Scalar;
194  typedef typename MatrixType::RealScalar RealScalar;
195  typedef typename MatrixType::Index Index;
196  typedef typename MatrixType::ColXpr ColXpr;
197  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
198  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
199  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
200  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
201 
202  int n = mat.cols();
203  eigen_assert(mat.rows()==n && vec.size()==n);
204 
205  TempVectorType temp;
206 
207  if(sigma>0)
208  {
209  // This version is based on Givens rotations.
210  // It is faster than the other one below, but only works for updates,
211  // i.e., for sigma > 0
212  temp = sqrt(sigma) * vec;
213 
214  for(int i=0; i<n; ++i)
215  {
216  JacobiRotation<Scalar> g;
217  g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
218 
219  int rs = n-i-1;
220  if(rs>0)
221  {
222  ColXprSegment x(mat.col(i).tail(rs));
223  TempVecSegment y(temp.tail(rs));
224  apply_rotation_in_the_plane(x, y, g);
225  }
226  }
227  }
228  else
229  {
230  temp = vec;
231  RealScalar beta = 1;
232  for(int j=0; j<n; ++j)
233  {
234  RealScalar Ljj = real(mat.coeff(j,j));
235  RealScalar dj = abs2(Ljj);
236  Scalar wj = temp.coeff(j);
237  RealScalar swj2 = sigma*abs2(wj);
238  RealScalar gamma = dj*beta + swj2;
239 
240  RealScalar x = dj + swj2/beta;
241  if (x<=RealScalar(0))
242  return j;
243  RealScalar nLjj = sqrt(x);
244  mat.coeffRef(j,j) = nLjj;
245  beta += swj2/dj;
246 
247  // Update the terms of L
248  Index rs = n-j-1;
249  if(rs)
250  {
251  temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
252  if(gamma != 0)
253  mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
254  }
255  }
256  }
257  return -1;
258 }
259 
260 template<typename Scalar> struct llt_inplace<Scalar, Lower>
261 {
262  typedef typename NumTraits<Scalar>::Real RealScalar;
263  template<typename MatrixType>
264  static typename MatrixType::Index unblocked(MatrixType& mat)
265  {
266  typedef typename MatrixType::Index Index;
267 
268  eigen_assert(mat.rows()==mat.cols());
269  const Index size = mat.rows();
270  for(Index k = 0; k < size; ++k)
271  {
272  Index rs = size-k-1; // remaining size
273 
274  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
275  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
276  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
277 
278  RealScalar x = real(mat.coeff(k,k));
279  if (k>0) x -= A10.squaredNorm();
280  if (x<=RealScalar(0))
281  return k;
282  mat.coeffRef(k,k) = x = sqrt(x);
283  if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
284  if (rs>0) A21 *= RealScalar(1)/x;
285  }
286  return -1;
287  }
288 
289  template<typename MatrixType>
290  static typename MatrixType::Index blocked(MatrixType& m)
291  {
292  typedef typename MatrixType::Index Index;
293  eigen_assert(m.rows()==m.cols());
294  Index size = m.rows();
295  if(size<32)
296  return unblocked(m);
297 
298  Index blockSize = size/8;
299  blockSize = (blockSize/16)*16;
300  blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
301 
302  for (Index k=0; k<size; k+=blockSize)
303  {
304  // partition the matrix:
305  // A00 | - | -
306  // lu = A10 | A11 | -
307  // A20 | A21 | A22
308  Index bs = (std::min)(blockSize, size-k);
309  Index rs = size - k - bs;
310  Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
311  Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
312  Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
313 
314  Index ret;
315  if((ret=unblocked(A11))>=0) return k+ret;
316  if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
317  if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
318  }
319  return -1;
320  }
321 
322  template<typename MatrixType, typename VectorType>
323  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
324  {
325  return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
326  }
327 };
328 
329 template<typename Scalar> struct llt_inplace<Scalar, Upper>
330 {
331  typedef typename NumTraits<Scalar>::Real RealScalar;
332 
333  template<typename MatrixType>
334  static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
335  {
336  Transpose<MatrixType> matt(mat);
337  return llt_inplace<Scalar, Lower>::unblocked(matt);
338  }
339  template<typename MatrixType>
340  static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
341  {
342  Transpose<MatrixType> matt(mat);
343  return llt_inplace<Scalar, Lower>::blocked(matt);
344  }
345  template<typename MatrixType, typename VectorType>
346  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
347  {
348  Transpose<MatrixType> matt(mat);
349  return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
350  }
351 };
352 
353 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
354 {
355  typedef const TriangularView<const MatrixType, Lower> MatrixL;
356  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
357  static inline MatrixL getL(const MatrixType& m) { return m; }
358  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
359  static bool inplace_decomposition(MatrixType& m)
360  { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
361 };
362 
363 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
364 {
365  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
366  typedef const TriangularView<const MatrixType, Upper> MatrixU;
367  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
368  static inline MatrixU getU(const MatrixType& m) { return m; }
369  static bool inplace_decomposition(MatrixType& m)
370  { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
371 };
372 
373 } // end namespace internal
374 
382 template<typename MatrixType, int _UpLo>
384 {
385  eigen_assert(a.rows()==a.cols());
386  const Index size = a.rows();
387  m_matrix.resize(size, size);
388  m_matrix = a;
389 
390  m_isInitialized = true;
391  bool ok = Traits::inplace_decomposition(m_matrix);
392  m_info = ok ? Success : NumericalIssue;
393 
394  return *this;
395 }
396 
402 template<typename _MatrixType, int _UpLo>
403 template<typename VectorType>
404 LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
405 {
406  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
407  eigen_assert(v.size()==m_matrix.cols());
408  eigen_assert(m_isInitialized);
409  if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
410  m_info = NumericalIssue;
411  else
412  m_info = Success;
413 
414  return *this;
415 }
416 
417 namespace internal {
418 template<typename _MatrixType, int UpLo, typename Rhs>
419 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
420  : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
421 {
422  typedef LLT<_MatrixType,UpLo> LLTType;
423  EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
424 
425  template<typename Dest> void evalTo(Dest& dst) const
426  {
427  dst = rhs();
428  dec().solveInPlace(dst);
429  }
430 };
431 }
432 
446 template<typename MatrixType, int _UpLo>
447 template<typename Derived>
448 void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
449 {
450  eigen_assert(m_isInitialized && "LLT is not initialized.");
451  eigen_assert(m_matrix.rows()==bAndX.rows());
452  matrixL().solveInPlace(bAndX);
453  matrixU().solveInPlace(bAndX);
454 }
455 
459 template<typename MatrixType, int _UpLo>
461 {
462  eigen_assert(m_isInitialized && "LLT is not initialized.");
463  return matrixL() * matrixL().adjoint().toDenseMatrix();
464 }
465 
469 template<typename Derived>
472 {
473  return LLT<PlainObject>(derived());
474 }
475 
479 template<typename MatrixType, unsigned int UpLo>
482 {
483  return LLT<PlainObject,UpLo>(m_matrix);
484 }
485 
486 } // end namespace Eigen
487 
488 #endif // EIGEN_LLT_H