LLT.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_LLT_H
26 #define EIGEN_LLT_H
27 
28 namespace Eigen {
29 
30 namespace internal{
31 template<typename MatrixType, int UpLo> struct LLT_Traits;
32 }
33 
61  /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
62  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
63  * the strict lower part does not have to store correct values.
64  */
65 template<typename _MatrixType, int _UpLo> class LLT
66 {
67  public:
68  typedef _MatrixType MatrixType;
69  enum {
73  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
74  };
75  typedef typename MatrixType::Scalar Scalar;
77  typedef typename MatrixType::Index Index;
78 
79  enum {
80  PacketSize = internal::packet_traits<Scalar>::size,
82  UpLo = _UpLo
83  };
84 
85  typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
86 
93  LLT() : m_matrix(), m_isInitialized(false) {}
94 
101  LLT(Index size) : m_matrix(size, size),
102  m_isInitialized(false) {}
103 
104  LLT(const MatrixType& matrix)
105  : m_matrix(matrix.rows(), matrix.cols()),
106  m_isInitialized(false)
107  {
108  compute(matrix);
109  }
110 
112  inline typename Traits::MatrixU matrixU() const
113  {
114  eigen_assert(m_isInitialized && "LLT is not initialized.");
115  return Traits::getU(m_matrix);
116  }
117 
119  inline typename Traits::MatrixL matrixL() const
120  {
121  eigen_assert(m_isInitialized && "LLT is not initialized.");
122  return Traits::getL(m_matrix);
123  }
124 
135  template<typename Rhs>
136  inline const internal::solve_retval<LLT, Rhs>
137  solve(const MatrixBase<Rhs>& b) const
138  {
139  eigen_assert(m_isInitialized && "LLT is not initialized.");
140  eigen_assert(m_matrix.rows()==b.rows()
141  && "LLT::solve(): invalid number of rows of the right hand side matrix b");
142  return internal::solve_retval<LLT, Rhs>(*this, b.derived());
143  }
144 
145  #ifdef EIGEN2_SUPPORT
146  template<typename OtherDerived, typename ResultType>
147  bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const
148  {
149  *result = this->solve(b);
150  return true;
151  }
152 
153  bool isPositiveDefinite() const { return true; }
154  #endif
155 
156  template<typename Derived>
157  void solveInPlace(MatrixBase<Derived> &bAndX) const;
158 
159  LLT& compute(const MatrixType& matrix);
160 
165  inline const MatrixType& matrixLLT() const
166  {
167  eigen_assert(m_isInitialized && "LLT is not initialized.");
168  return m_matrix;
169  }
170 
172 
173 
180  {
181  eigen_assert(m_isInitialized && "LLT is not initialized.");
182  return m_info;
183  }
184 
185  inline Index rows() const { return m_matrix.rows(); }
186  inline Index cols() const { return m_matrix.cols(); }
187 
188  template<typename VectorType>
189  LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
190 
191  protected:
199 };
200 
201 namespace internal {
202 
203 template<typename Scalar, int UpLo> struct llt_inplace;
204 
205 template<typename MatrixType, typename VectorType>
206 static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
207 {
208  typedef typename MatrixType::Scalar Scalar;
209  typedef typename MatrixType::RealScalar RealScalar;
210  typedef typename MatrixType::Index Index;
211  typedef typename MatrixType::ColXpr ColXpr;
212  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
213  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
214  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
215  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
216 
217  int n = mat.cols();
218  eigen_assert(mat.rows()==n && vec.size()==n);
219 
220  TempVectorType temp;
221 
222  if(sigma>0)
223  {
224  // This version is based on Givens rotations.
225  // It is faster than the other one below, but only works for updates,
226  // i.e., for sigma > 0
227  temp = sqrt(sigma) * vec;
228 
229  for(int i=0; i<n; ++i)
230  {
232  g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
233 
234  int rs = n-i-1;
235  if(rs>0)
236  {
237  ColXprSegment x(mat.col(i).tail(rs));
238  TempVecSegment y(temp.tail(rs));
240  }
241  }
242  }
243  else
244  {
245  temp = vec;
246  RealScalar beta = 1;
247  for(int j=0; j<n; ++j)
248  {
249  RealScalar Ljj = real(mat.coeff(j,j));
250  RealScalar dj = abs2(Ljj);
251  Scalar wj = temp.coeff(j);
252  RealScalar swj2 = sigma*abs2(wj);
253  RealScalar gamma = dj*beta + swj2;
254 
255  RealScalar x = dj + swj2/beta;
256  if (x<=RealScalar(0))
257  return j;
258  RealScalar nLjj = sqrt(x);
259  mat.coeffRef(j,j) = nLjj;
260  beta += swj2/dj;
261 
262  // Update the terms of L
263  Index rs = n-j-1;
264  if(rs)
265  {
266  temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
267  if(gamma != 0)
268  mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
269  }
270  }
271  }
272  return -1;
273 }
274 
275 template<typename Scalar> struct llt_inplace<Scalar, Lower>
276 {
277  typedef typename NumTraits<Scalar>::Real RealScalar;
278  template<typename MatrixType>
279  static typename MatrixType::Index unblocked(MatrixType& mat)
280  {
281  typedef typename MatrixType::Index Index;
282 
283  eigen_assert(mat.rows()==mat.cols());
284  const Index size = mat.rows();
285  for(Index k = 0; k < size; ++k)
286  {
287  Index rs = size-k-1; // remaining size
288 
289  Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
290  Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
291  Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
292 
293  RealScalar x = real(mat.coeff(k,k));
294  if (k>0) x -= A10.squaredNorm();
295  if (x<=RealScalar(0))
296  return k;
297  mat.coeffRef(k,k) = x = sqrt(x);
298  if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
299  if (rs>0) A21 *= RealScalar(1)/x;
300  }
301  return -1;
302  }
303 
304  template<typename MatrixType>
305  static typename MatrixType::Index blocked(MatrixType& m)
306  {
307  typedef typename MatrixType::Index Index;
308  eigen_assert(m.rows()==m.cols());
309  Index size = m.rows();
310  if(size<32)
311  return unblocked(m);
312 
313  Index blockSize = size/8;
314  blockSize = (blockSize/16)*16;
315  blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
316 
317  for (Index k=0; k<size; k+=blockSize)
318  {
319  // partition the matrix:
320  // A00 | - | -
321  // lu = A10 | A11 | -
322  // A20 | A21 | A22
323  Index bs = (std::min)(blockSize, size-k);
324  Index rs = size - k - bs;
325  Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
326  Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
327  Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
328 
329  Index ret;
330  if((ret=unblocked(A11))>=0) return k+ret;
331  if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
332  if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
333  }
334  return -1;
335  }
336 
337  template<typename MatrixType, typename VectorType>
338  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
339  {
340  return Eigen::internal::llt_rank_update_lower(mat, vec, sigma);
341  }
342 };
343 
344 template<typename Scalar> struct llt_inplace<Scalar, Upper>
345 {
346  typedef typename NumTraits<Scalar>::Real RealScalar;
347 
348  template<typename MatrixType>
349  static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
350  {
351  Transpose<MatrixType> matt(mat);
352  return llt_inplace<Scalar, Lower>::unblocked(matt);
353  }
354  template<typename MatrixType>
355  static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
356  {
357  Transpose<MatrixType> matt(mat);
358  return llt_inplace<Scalar, Lower>::blocked(matt);
359  }
360  template<typename MatrixType, typename VectorType>
361  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
362  {
363  Transpose<MatrixType> matt(mat);
364  return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
365  }
366 };
367 
368 template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
369 {
370  typedef const TriangularView<const MatrixType, Lower> MatrixL;
371  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU;
372  static inline MatrixL getL(const MatrixType& m) { return m; }
373  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
374  static bool inplace_decomposition(MatrixType& m)
375  { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; }
376 };
377 
378 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
379 {
380  typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL;
381  typedef const TriangularView<const MatrixType, Upper> MatrixU;
382  static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); }
383  static inline MatrixU getU(const MatrixType& m) { return m; }
384  static bool inplace_decomposition(MatrixType& m)
385  { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; }
386 };
387 
388 } // end namespace internal
389 
397 template<typename MatrixType, int _UpLo>
399 {
400  eigen_assert(a.rows()==a.cols());
401  const Index size = a.rows();
402  m_matrix.resize(size, size);
403  m_matrix = a;
404 
405  m_isInitialized = true;
406  bool ok = Traits::inplace_decomposition(m_matrix);
407  m_info = ok ? Success : NumericalIssue;
408 
409  return *this;
410 }
411 
417 template<typename _MatrixType, int _UpLo>
418 template<typename VectorType>
420 {
422  eigen_assert(v.size()==m_matrix.cols());
423  eigen_assert(m_isInitialized);
424  if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
425  m_info = NumericalIssue;
426  else
427  m_info = Success;
428 
429  return *this;
430 }
431 
432 namespace internal {
433 template<typename _MatrixType, int UpLo, typename Rhs>
434 struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
435  : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
436 {
437  typedef LLT<_MatrixType,UpLo> LLTType;
438  EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
439 
440  template<typename Dest> void evalTo(Dest& dst) const
441  {
442  dst = rhs();
443  dec().solveInPlace(dst);
444  }
445 };
446 }
447 
461 template<typename MatrixType, int _UpLo>
462 template<typename Derived>
464 {
465  eigen_assert(m_isInitialized && "LLT is not initialized.");
466  eigen_assert(m_matrix.rows()==bAndX.rows());
467  matrixL().solveInPlace(bAndX);
468  matrixU().solveInPlace(bAndX);
469 }
470 
474 template<typename MatrixType, int _UpLo>
476 {
477  eigen_assert(m_isInitialized && "LLT is not initialized.");
478  return matrixL() * matrixL().adjoint().toDenseMatrix();
479 }
480 
484 template<typename Derived>
487 {
488  return LLT<PlainObject>(derived());
489 }
490 
494 template<typename MatrixType, unsigned int UpLo>
497 {
498  return LLT<PlainObject,UpLo>(m_matrix);
499 }
500 
501 } // end namespace Eigen
502 
503 #endif // EIGEN_LLT_H