Eigen  3.2.2
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ColPivHouseholderQR.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
13 
14 namespace Eigen {
15 
37 template<typename _MatrixType> class ColPivHouseholderQR
38 {
39  public:
40 
41  typedef _MatrixType MatrixType;
42  enum {
43  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
44  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
45  Options = MatrixType::Options,
46  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
47  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
48  };
49  typedef typename MatrixType::Scalar Scalar;
50  typedef typename MatrixType::RealScalar RealScalar;
51  typedef typename MatrixType::Index Index;
52  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
53  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
54  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
55  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
56  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
57  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
58  typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
59 
60  private:
61 
62  typedef typename PermutationType::Index PermIndexType;
63 
64  public:
65 
73  : m_qr(),
74  m_hCoeffs(),
75  m_colsPermutation(),
76  m_colsTranspositions(),
77  m_temp(),
78  m_colSqNorms(),
79  m_isInitialized(false),
80  m_usePrescribedThreshold(false) {}
81 
88  ColPivHouseholderQR(Index rows, Index cols)
89  : m_qr(rows, cols),
90  m_hCoeffs((std::min)(rows,cols)),
91  m_colsPermutation(PermIndexType(cols)),
92  m_colsTranspositions(cols),
93  m_temp(cols),
94  m_colSqNorms(cols),
95  m_isInitialized(false),
96  m_usePrescribedThreshold(false) {}
97 
110  ColPivHouseholderQR(const MatrixType& matrix)
111  : m_qr(matrix.rows(), matrix.cols()),
112  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
113  m_colsPermutation(PermIndexType(matrix.cols())),
114  m_colsTranspositions(matrix.cols()),
115  m_temp(matrix.cols()),
116  m_colSqNorms(matrix.cols()),
117  m_isInitialized(false),
118  m_usePrescribedThreshold(false)
119  {
120  compute(matrix);
121  }
122 
140  template<typename Rhs>
141  inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
142  solve(const MatrixBase<Rhs>& b) const
143  {
144  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
145  return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
146  }
147 
148  HouseholderSequenceType householderQ(void) const;
149  HouseholderSequenceType matrixQ(void) const
150  {
151  return householderQ();
152  }
153 
156  const MatrixType& matrixQR() const
157  {
158  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
159  return m_qr;
160  }
161 
171  const MatrixType& matrixR() const
172  {
173  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
174  return m_qr;
175  }
176 
177  ColPivHouseholderQR& compute(const MatrixType& matrix);
178 
180  const PermutationType& colsPermutation() const
181  {
182  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
183  return m_colsPermutation;
184  }
185 
199  typename MatrixType::RealScalar absDeterminant() const;
200 
213  typename MatrixType::RealScalar logAbsDeterminant() const;
214 
221  inline Index rank() const
222  {
223  using std::abs;
224  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
225  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
226  Index result = 0;
227  for(Index i = 0; i < m_nonzero_pivots; ++i)
228  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
229  return result;
230  }
231 
238  inline Index dimensionOfKernel() const
239  {
240  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
241  return cols() - rank();
242  }
243 
251  inline bool isInjective() const
252  {
253  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
254  return rank() == cols();
255  }
256 
264  inline bool isSurjective() const
265  {
266  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
267  return rank() == rows();
268  }
269 
276  inline bool isInvertible() const
277  {
278  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
279  return isInjective() && isSurjective();
280  }
281 
287  inline const
288  internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
289  inverse() const
290  {
291  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
292  return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
293  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
294  }
295 
296  inline Index rows() const { return m_qr.rows(); }
297  inline Index cols() const { return m_qr.cols(); }
298 
303  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
304 
323  {
324  m_usePrescribedThreshold = true;
325  m_prescribedThreshold = threshold;
326  return *this;
327  }
328 
338  {
339  m_usePrescribedThreshold = false;
340  return *this;
341  }
342 
347  RealScalar threshold() const
348  {
349  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
350  return m_usePrescribedThreshold ? m_prescribedThreshold
351  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
352  // and turns out to be identical to Higham's formula used already in LDLt.
353  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
354  }
355 
363  inline Index nonzeroPivots() const
364  {
365  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
366  return m_nonzero_pivots;
367  }
368 
372  RealScalar maxPivot() const { return m_maxpivot; }
373 
381  {
382  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
383  return Success;
384  }
385 
386  protected:
387  MatrixType m_qr;
388  HCoeffsType m_hCoeffs;
389  PermutationType m_colsPermutation;
390  IntRowVectorType m_colsTranspositions;
391  RowVectorType m_temp;
392  RealRowVectorType m_colSqNorms;
393  bool m_isInitialized, m_usePrescribedThreshold;
394  RealScalar m_prescribedThreshold, m_maxpivot;
395  Index m_nonzero_pivots;
396  Index m_det_pq;
397 };
398 
399 template<typename MatrixType>
400 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
401 {
402  using std::abs;
403  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
404  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
405  return abs(m_qr.diagonal().prod());
406 }
407 
408 template<typename MatrixType>
409 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
410 {
411  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
412  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
413  return m_qr.diagonal().cwiseAbs().array().log().sum();
414 }
415 
422 template<typename MatrixType>
424 {
425  using std::abs;
426  Index rows = matrix.rows();
427  Index cols = matrix.cols();
428  Index size = matrix.diagonalSize();
429 
430  // the column permutation is stored as int indices, so just to be sure:
431  eigen_assert(cols<=NumTraits<int>::highest());
432 
433  m_qr = matrix;
434  m_hCoeffs.resize(size);
435 
436  m_temp.resize(cols);
437 
438  m_colsTranspositions.resize(matrix.cols());
439  Index number_of_transpositions = 0;
440 
441  m_colSqNorms.resize(cols);
442  for(Index k = 0; k < cols; ++k)
443  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
444 
445  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
446 
447  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
448  m_maxpivot = RealScalar(0);
449 
450  for(Index k = 0; k < size; ++k)
451  {
452  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
453  Index biggest_col_index;
454  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
455  biggest_col_index += k;
456 
457  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
458  // the actual squared norm of the selected column.
459  // Note that not doing so does result in solve() sometimes returning inf/nan values
460  // when running the unit test with 1000 repetitions.
461  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
462 
463  // we store that back into our table: it can't hurt to correct our table.
464  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
465 
466  // if the current biggest column is smaller than epsilon times the initial biggest column,
467  // terminate to avoid generating nan/inf values.
468  // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
469  // repetitions of the unit test, with the result of solve() filled with large values of the order
470  // of 1/(size*epsilon).
471  if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
472  {
473  m_nonzero_pivots = k;
474  m_hCoeffs.tail(size-k).setZero();
475  m_qr.bottomRightCorner(rows-k,cols-k)
476  .template triangularView<StrictlyLower>()
477  .setZero();
478  break;
479  }
480 
481  // apply the transposition to the columns
482  m_colsTranspositions.coeffRef(k) = biggest_col_index;
483  if(k != biggest_col_index) {
484  m_qr.col(k).swap(m_qr.col(biggest_col_index));
485  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
486  ++number_of_transpositions;
487  }
488 
489  // generate the householder vector, store it below the diagonal
490  RealScalar beta;
491  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
492 
493  // apply the householder transformation to the diagonal coefficient
494  m_qr.coeffRef(k,k) = beta;
495 
496  // remember the maximum absolute value of diagonal coefficients
497  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
498 
499  // apply the householder transformation
500  m_qr.bottomRightCorner(rows-k, cols-k-1)
501  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
502 
503  // update our table of squared norms of the columns
504  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
505  }
506 
507  m_colsPermutation.setIdentity(PermIndexType(cols));
508  for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)
509  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
510 
511  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
512  m_isInitialized = true;
513 
514  return *this;
515 }
516 
517 namespace internal {
518 
519 template<typename _MatrixType, typename Rhs>
520 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
521  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
522 {
523  EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
524 
525  template<typename Dest> void evalTo(Dest& dst) const
526  {
527  eigen_assert(rhs().rows() == dec().rows());
528 
529  const Index cols = dec().cols(),
530  nonzero_pivots = dec().nonzeroPivots();
531 
532  if(nonzero_pivots == 0)
533  {
534  dst.setZero();
535  return;
536  }
537 
538  typename Rhs::PlainObject c(rhs());
539 
540  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
541  c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
542  .setLength(dec().nonzeroPivots())
543  .transpose()
544  );
545 
546  dec().matrixR()
547  .topLeftCorner(nonzero_pivots, nonzero_pivots)
548  .template triangularView<Upper>()
549  .solveInPlace(c.topRows(nonzero_pivots));
550 
551  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
552  for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
553  }
554 };
555 
556 } // end namespace internal
557 
559 template<typename MatrixType>
560 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
562 {
563  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
564  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
565 }
566 
571 template<typename Derived>
574 {
575  return ColPivHouseholderQR<PlainObject>(eval());
576 }
577 
578 } // end namespace Eigen
579 
580 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:156
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ForwardDeclarations.h:222
HouseholderSequenceType householderQ(void) const
Definition: ColPivHouseholderQR.h:561
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:303
bool isInjective() const
Definition: ColPivHouseholderQR.h:251
ComputationInfo info() const
Reports whether the QR factorization was succesful.
Definition: ColPivHouseholderQR.h:380
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:372
Definition: LDLT.h:16
Index rank() const
Definition: ColPivHouseholderQR.h:221
Definition: StdDeque.h:58
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:400
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:422
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:180
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: ColPivHouseholderQR.h:88
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:322
ColPivHouseholderQR & compute(const MatrixType &matrix)
Definition: ColPivHouseholderQR.h:423
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:171
bool isInvertible() const
Definition: ColPivHouseholderQR.h:276
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:347
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:409
ColPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:110
bool isSurjective() const
Definition: ColPivHouseholderQR.h:264
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:337
Definition: Eigen_Colamd.h:54
ColPivHouseholderQR()
Default Constructor.
Definition: ColPivHouseholderQR.h:72
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:363
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:238
const internal::solve_retval< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: ColPivHouseholderQR.h:142
Definition: Constants.h:376
const internal::solve_retval< ColPivHouseholderQR, typename MatrixType::IdentityReturnType > inverse() const
Definition: ColPivHouseholderQR.h:289
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const ColPivHouseholderQR< PlainObject > colPivHouseholderQr() const
Definition: ColPivHouseholderQR.h:573