ColPivHouseholderQR.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-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;
59 
60  private:
61 
63 
64  public:
65 
73  : m_qr(),
74  m_hCoeffs(),
77  m_temp(),
78  m_colSqNorms(),
79  m_isInitialized(false) {}
80 
88  : m_qr(rows, cols),
89  m_hCoeffs((std::min)(rows,cols)),
90  m_colsPermutation(PermIndexType(cols)),
92  m_temp(cols),
93  m_colSqNorms(cols),
94  m_isInitialized(false),
95  m_usePrescribedThreshold(false) {}
96 
109  ColPivHouseholderQR(const MatrixType& matrix)
110  : m_qr(matrix.rows(), matrix.cols()),
111  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
112  m_colsPermutation(PermIndexType(matrix.cols())),
113  m_colsTranspositions(matrix.cols()),
114  m_temp(matrix.cols()),
115  m_colSqNorms(matrix.cols()),
116  m_isInitialized(false),
118  {
119  compute(matrix);
120  }
121 
139  template<typename Rhs>
141  solve(const MatrixBase<Rhs>& b) const
142  {
143  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
144  return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
145  }
146 
147  HouseholderSequenceType householderQ(void) const;
148  HouseholderSequenceType matrixQ(void) const
149  {
150  return householderQ();
151  }
152 
155  const MatrixType& matrixQR() const
156  {
157  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
158  return m_qr;
159  }
160 
170  const MatrixType& matrixR() const
171  {
172  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
173  return m_qr;
174  }
175 
176  ColPivHouseholderQR& compute(const MatrixType& matrix);
177 
179  const PermutationType& colsPermutation() const
180  {
181  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
182  return m_colsPermutation;
183  }
184 
198  typename MatrixType::RealScalar absDeterminant() const;
199 
212  typename MatrixType::RealScalar logAbsDeterminant() const;
213 
220  inline Index rank() const
221  {
222  using std::abs;
223  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
224  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
225  Index result = 0;
226  for(Index i = 0; i < m_nonzero_pivots; ++i)
227  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
228  return result;
229  }
230 
237  inline Index dimensionOfKernel() const
238  {
239  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
240  return cols() - rank();
241  }
242 
250  inline bool isInjective() const
251  {
252  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
253  return rank() == cols();
254  }
255 
263  inline bool isSurjective() const
264  {
265  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
266  return rank() == rows();
267  }
268 
275  inline bool isInvertible() const
276  {
277  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
278  return isInjective() && isSurjective();
279  }
280 
286  inline const
288  inverse() const
289  {
290  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
292  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
293  }
294 
295  inline Index rows() const { return m_qr.rows(); }
296  inline Index cols() const { return m_qr.cols(); }
297 
302  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
303 
322  {
325  return *this;
326  }
327 
337  {
338  m_usePrescribedThreshold = false;
339  return *this;
340  }
341 
346  RealScalar threshold() const
347  {
350  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
351  // and turns out to be identical to Higham's formula used already in LDLt.
352  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
353  }
354 
362  inline Index nonzeroPivots() const
363  {
364  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
365  return m_nonzero_pivots;
366  }
367 
371  RealScalar maxPivot() const { return m_maxpivot; }
372 
380  {
381  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
382  return Success;
383  }
384 
385  protected:
386  MatrixType m_qr;
387  HCoeffsType m_hCoeffs;
388  PermutationType m_colsPermutation;
389  IntRowVectorType m_colsTranspositions;
390  RowVectorType m_temp;
391  RealRowVectorType m_colSqNorms;
395  Index m_det_pq;
396 };
397 
398 template<typename MatrixType>
399 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
400 {
401  using std::abs;
402  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
403  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
404  return abs(m_qr.diagonal().prod());
405 }
406 
407 template<typename MatrixType>
408 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
409 {
410  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
411  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
412  return m_qr.diagonal().cwiseAbs().array().log().sum();
413 }
414 
421 template<typename MatrixType>
423 {
424  using std::abs;
425  Index rows = matrix.rows();
426  Index cols = matrix.cols();
427  Index size = matrix.diagonalSize();
428 
429  // the column permutation is stored as int indices, so just to be sure:
431 
432  m_qr = matrix;
433  m_hCoeffs.resize(size);
434 
435  m_temp.resize(cols);
436 
437  m_colsTranspositions.resize(matrix.cols());
438  Index number_of_transpositions = 0;
439 
440  m_colSqNorms.resize(cols);
441  for(Index k = 0; k < cols; ++k)
442  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
443 
444  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
445 
446  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
447  m_maxpivot = RealScalar(0);
448 
449  for(Index k = 0; k < size; ++k)
450  {
451  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
452  Index biggest_col_index;
453  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
454  biggest_col_index += k;
455 
456  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
457  // the actual squared norm of the selected column.
458  // Note that not doing so does result in solve() sometimes returning inf/nan values
459  // when running the unit test with 1000 repetitions.
460  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
461 
462  // we store that back into our table: it can't hurt to correct our table.
463  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
464 
465  // if the current biggest column is smaller than epsilon times the initial biggest column,
466  // terminate to avoid generating nan/inf values.
467  // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
468  // repetitions of the unit test, with the result of solve() filled with large values of the order
469  // of 1/(size*epsilon).
470  if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
471  {
472  m_nonzero_pivots = k;
473  m_hCoeffs.tail(size-k).setZero();
474  m_qr.bottomRightCorner(rows-k,cols-k)
475  .template triangularView<StrictlyLower>()
476  .setZero();
477  break;
478  }
479 
480  // apply the transposition to the columns
481  m_colsTranspositions.coeffRef(k) = biggest_col_index;
482  if(k != biggest_col_index) {
483  m_qr.col(k).swap(m_qr.col(biggest_col_index));
484  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
485  ++number_of_transpositions;
486  }
487 
488  // generate the householder vector, store it below the diagonal
489  RealScalar beta;
490  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
491 
492  // apply the householder transformation to the diagonal coefficient
493  m_qr.coeffRef(k,k) = beta;
494 
495  // remember the maximum absolute value of diagonal coefficients
496  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
497 
498  // apply the householder transformation
499  m_qr.bottomRightCorner(rows-k, cols-k-1)
500  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
501 
502  // update our table of squared norms of the columns
503  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
504  }
505 
507  for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)
509 
510  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
511  m_isInitialized = true;
512 
513  return *this;
514 }
515 
516 namespace internal {
517 
518 template<typename _MatrixType, typename Rhs>
519 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
520  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
521 {
523 
524  template<typename Dest> void evalTo(Dest& dst) const
525  {
526  eigen_assert(rhs().rows() == dec().rows());
527 
528  const Index cols = dec().cols(),
529  nonzero_pivots = dec().nonzeroPivots();
530 
531  if(nonzero_pivots == 0)
532  {
533  dst.setZero();
534  return;
535  }
536 
537  typename Rhs::PlainObject c(rhs());
538 
539  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
540  c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
541  .setLength(dec().nonzeroPivots())
542  .transpose()
543  );
544 
545  dec().matrixR()
546  .topLeftCorner(nonzero_pivots, nonzero_pivots)
547  .template triangularView<Upper>()
548  .solveInPlace(c.topRows(nonzero_pivots));
549 
550  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
551  for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
552  }
553 };
554 
555 } // end namespace internal
556 
558 template<typename MatrixType>
561 {
562  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
563  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
564 }
565 
570 template<typename Derived>
573 {
574  return ColPivHouseholderQR<PlainObject>(eval());
575 }
576 
577 } // end namespace Eigen
578 
579 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
PermutationMatrix< ColsAtCompileTime, MaxColsAtCompileTime > PermutationType
const PermutationType & colsPermutation() const
PermutationType::Index PermIndexType
ColPivHouseholderQR & compute(const MatrixType &matrix)
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
ColPivHouseholderQR()
Default Constructor.
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
const internal::solve_retval< ColPivHouseholderQR, typename MatrixType::IdentityReturnType > inverse() const
MatrixType::RealScalar logAbsDeterminant() const
internal::plain_diag_type< MatrixType >::type HCoeffsType
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Derived & applyTranspositionOnTheRight(Index i, Index j)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
ComputationInfo info() const
Reports whether the QR factorization was succesful.
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
Sequence of Householder reflections acting on subspaces with decreasing size.
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
const internal::solve_retval< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
ColPivHouseholderQR & setThreshold(Default_t)
HouseholderSequenceType matrixQ(void) const
internal::plain_row_type< MatrixType >::type RowVectorType
HouseholderSequenceType householderQ(void) const
const MatrixType & matrixR() const
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
const MatrixType & matrixQR() const
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
const HCoeffsType & hCoeffs() const
Matrix< Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime > MatrixQType
const ColPivHouseholderQR< PlainObject > colPivHouseholderQr() const
internal::traits< Derived >::Index Index
Definition: EigenBase.h:31
internal::plain_row_type< MatrixType, RealScalar >::type RealRowVectorType
void rhs(const real_t *x, real_t *f)
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
MatrixType::RealScalar absDeterminant() const
ColPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
internal::plain_row_type< MatrixType, Index >::type IntRowVectorType
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
IntRowVectorType m_colsTranspositions
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
#define eigen_assert(x)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
MatrixType::RealScalar RealScalar


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:30