Eigen  3.2.10
FullPivHouseholderQR.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_FULLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
18 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
19 
20 template<typename MatrixType>
21 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
22 {
23  typedef typename MatrixType::PlainObject ReturnType;
24 };
25 
26 }
27 
49 template<typename _MatrixType> class FullPivHouseholderQR
50 {
51  public:
52 
53  typedef _MatrixType MatrixType;
54  enum {
55  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
57  Options = MatrixType::Options,
58  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
59  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
60  };
61  typedef typename MatrixType::Scalar Scalar;
62  typedef typename MatrixType::RealScalar RealScalar;
63  typedef typename MatrixType::Index Index;
64  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
65  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
66  typedef Matrix<Index, 1,
67  EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
68  EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
69  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
70  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
71  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
72 
79  : m_qr(),
80  m_hCoeffs(),
81  m_rows_transpositions(),
82  m_cols_transpositions(),
83  m_cols_permutation(),
84  m_temp(),
85  m_isInitialized(false),
86  m_usePrescribedThreshold(false) {}
87 
94  FullPivHouseholderQR(Index rows, Index cols)
95  : m_qr(rows, cols),
96  m_hCoeffs((std::min)(rows,cols)),
97  m_rows_transpositions((std::min)(rows,cols)),
98  m_cols_transpositions((std::min)(rows,cols)),
99  m_cols_permutation(cols),
100  m_temp(cols),
101  m_isInitialized(false),
102  m_usePrescribedThreshold(false) {}
103 
116  FullPivHouseholderQR(const MatrixType& matrix)
117  : m_qr(matrix.rows(), matrix.cols()),
118  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
119  m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
120  m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
121  m_cols_permutation(matrix.cols()),
122  m_temp(matrix.cols()),
123  m_isInitialized(false),
124  m_usePrescribedThreshold(false)
125  {
126  compute(matrix);
127  }
128 
144  template<typename Rhs>
145  inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
146  solve(const MatrixBase<Rhs>& b) const
147  {
148  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
149  return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
150  }
151 
154  MatrixQReturnType matrixQ(void) const;
155 
158  const MatrixType& matrixQR() const
159  {
160  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
161  return m_qr;
162  }
163 
164  FullPivHouseholderQR& compute(const MatrixType& matrix);
165 
167  const PermutationType& colsPermutation() const
168  {
169  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
170  return m_cols_permutation;
171  }
172 
174  const IntDiagSizeVectorType& rowsTranspositions() const
175  {
176  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
177  return m_rows_transpositions;
178  }
179 
193  typename MatrixType::RealScalar absDeterminant() const;
194 
207  typename MatrixType::RealScalar logAbsDeterminant() const;
208 
215  inline Index rank() const
216  {
217  using std::abs;
218  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
219  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
220  Index result = 0;
221  for(Index i = 0; i < m_nonzero_pivots; ++i)
222  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
223  return result;
224  }
225 
232  inline Index dimensionOfKernel() const
233  {
234  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
235  return cols() - rank();
236  }
237 
245  inline bool isInjective() const
246  {
247  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
248  return rank() == cols();
249  }
250 
258  inline bool isSurjective() const
259  {
260  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
261  return rank() == rows();
262  }
263 
270  inline bool isInvertible() const
271  {
272  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
273  return isInjective() && isSurjective();
274  }
275  inline const
281  internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
282  inverse() const
283  {
284  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
285  return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
286  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
287  }
288 
289  inline Index rows() const { return m_qr.rows(); }
290  inline Index cols() const { return m_qr.cols(); }
291 
296  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
297 
315  FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
316  {
317  m_usePrescribedThreshold = true;
318  m_prescribedThreshold = threshold;
319  return *this;
320  }
321 
331  {
332  m_usePrescribedThreshold = false;
333  return *this;
334  }
335 
340  RealScalar threshold() const
341  {
342  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
343  return m_usePrescribedThreshold ? m_prescribedThreshold
344  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
345  // and turns out to be identical to Higham's formula used already in LDLt.
346  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
347  }
348 
356  inline Index nonzeroPivots() const
357  {
358  eigen_assert(m_isInitialized && "LU is not initialized.");
359  return m_nonzero_pivots;
360  }
361 
365  RealScalar maxPivot() const { return m_maxpivot; }
366 
367  protected:
368 
369  static void check_template_parameters()
370  {
371  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
372  }
373 
374  MatrixType m_qr;
375  HCoeffsType m_hCoeffs;
376  IntDiagSizeVectorType m_rows_transpositions;
377  IntDiagSizeVectorType m_cols_transpositions;
378  PermutationType m_cols_permutation;
379  RowVectorType m_temp;
380  bool m_isInitialized, m_usePrescribedThreshold;
381  RealScalar m_prescribedThreshold, m_maxpivot;
382  Index m_nonzero_pivots;
383  RealScalar m_precision;
384  Index m_det_pq;
385 };
386 
387 template<typename MatrixType>
388 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
389 {
390  using std::abs;
391  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
392  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
393  return abs(m_qr.diagonal().prod());
394 }
395 
396 template<typename MatrixType>
397 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
398 {
399  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
400  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
401  return m_qr.diagonal().cwiseAbs().array().log().sum();
402 }
403 
410 template<typename MatrixType>
412 {
413  check_template_parameters();
414 
415  using std::abs;
416  Index rows = matrix.rows();
417  Index cols = matrix.cols();
418  Index size = (std::min)(rows,cols);
419 
420  m_qr = matrix;
421  m_hCoeffs.resize(size);
422 
423  m_temp.resize(cols);
424 
425  m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
426 
427  m_rows_transpositions.resize(size);
428  m_cols_transpositions.resize(size);
429  Index number_of_transpositions = 0;
430 
431  RealScalar biggest(0);
432 
433  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
434  m_maxpivot = RealScalar(0);
435 
436  for (Index k = 0; k < size; ++k)
437  {
438  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
439  RealScalar biggest_in_corner;
440 
441  biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
442  .cwiseAbs()
443  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
444  row_of_biggest_in_corner += k;
445  col_of_biggest_in_corner += k;
446  if(k==0) biggest = biggest_in_corner;
447 
448  // if the corner is negligible, then we have less than full rank, and we can finish early
449  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
450  {
451  m_nonzero_pivots = k;
452  for(Index i = k; i < size; i++)
453  {
454  m_rows_transpositions.coeffRef(i) = i;
455  m_cols_transpositions.coeffRef(i) = i;
456  m_hCoeffs.coeffRef(i) = Scalar(0);
457  }
458  break;
459  }
460 
461  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
462  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
463  if(k != row_of_biggest_in_corner) {
464  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
465  ++number_of_transpositions;
466  }
467  if(k != col_of_biggest_in_corner) {
468  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
469  ++number_of_transpositions;
470  }
471 
472  RealScalar beta;
473  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
474  m_qr.coeffRef(k,k) = beta;
475 
476  // remember the maximum absolute value of diagonal coefficients
477  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
478 
479  m_qr.bottomRightCorner(rows-k, cols-k-1)
480  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
481  }
482 
483  m_cols_permutation.setIdentity(cols);
484  for(Index k = 0; k < size; ++k)
485  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
486 
487  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
488  m_isInitialized = true;
489 
490  return *this;
491 }
492 
493 namespace internal {
494 
495 template<typename _MatrixType, typename Rhs>
496 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
497  : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
498 {
499  EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
500 
501  template<typename Dest> void evalTo(Dest& dst) const
502  {
503  const Index rows = dec().rows(), cols = dec().cols();
504  eigen_assert(rhs().rows() == rows);
505 
506  // FIXME introduce nonzeroPivots() and use it here. and more generally,
507  // make the same improvements in this dec as in FullPivLU.
508  if(dec().rank()==0)
509  {
510  dst.setZero();
511  return;
512  }
513 
514  typename Rhs::PlainObject c(rhs());
515 
516  Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
517  for (Index k = 0; k < dec().rank(); ++k)
518  {
519  Index remainingSize = rows-k;
520  c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
521  c.bottomRightCorner(remainingSize, rhs().cols())
522  .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
523  dec().hCoeffs().coeff(k), &temp.coeffRef(0));
524  }
525 
526  dec().matrixQR()
527  .topLeftCorner(dec().rank(), dec().rank())
528  .template triangularView<Upper>()
529  .solveInPlace(c.topRows(dec().rank()));
530 
531  for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
532  for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
533  }
534 };
535 
542 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
543  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
544 {
545 public:
546  typedef typename MatrixType::Index Index;
548  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
549  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
550  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
551 
552  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
553  const HCoeffsType& hCoeffs,
554  const IntDiagSizeVectorType& rowsTranspositions)
555  : m_qr(qr),
556  m_hCoeffs(hCoeffs),
557  m_rowsTranspositions(rowsTranspositions)
558  {}
559 
560  template <typename ResultType>
561  void evalTo(ResultType& result) const
562  {
563  const Index rows = m_qr.rows();
564  WorkVectorType workspace(rows);
565  evalTo(result, workspace);
566  }
567 
568  template <typename ResultType>
569  void evalTo(ResultType& result, WorkVectorType& workspace) const
570  {
571  using numext::conj;
572  // compute the product H'_0 H'_1 ... H'_n-1,
573  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
574  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
575  const Index rows = m_qr.rows();
576  const Index cols = m_qr.cols();
577  const Index size = (std::min)(rows, cols);
578  workspace.resize(rows);
579  result.setIdentity(rows, rows);
580  for (Index k = size-1; k >= 0; k--)
581  {
582  result.block(k, k, rows-k, rows-k)
583  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
584  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
585  }
586  }
587 
588  Index rows() const { return m_qr.rows(); }
589  Index cols() const { return m_qr.rows(); }
590 
591 protected:
592  typename MatrixType::Nested m_qr;
593  typename HCoeffsType::Nested m_hCoeffs;
594  typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
595 };
596 
597 } // end namespace internal
598 
599 template<typename MatrixType>
600 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
601 {
602  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
603  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
604 }
605 
610 template<typename Derived>
613 {
614  return FullPivHouseholderQR<PlainObject>(eval());
615 }
616 
617 } // end namespace Eigen
618 
619 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:340
bool isInjective() const
Definition: FullPivHouseholderQR.h:245
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:232
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: ForwardDeclarations.h:223
Definition: LDLT.h:16
Definition: StdDeque.h:50
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:296
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:167
Definition: Constants.h:266
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:94
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:600
const internal::solve_retval< FullPivHouseholderQR, typename MatrixType::IdentityReturnType > inverse() const
Definition: FullPivHouseholderQR.h:282
FullPivHouseholderQR & compute(const MatrixType &matrix)
Definition: FullPivHouseholderQR.h:411
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:78
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:356
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:174
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:365
bool isSurjective() const
Definition: FullPivHouseholderQR.h:258
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:330
FullPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:116
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:612
Definition: Eigen_Colamd.h:50
const internal::solve_retval< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivHouseholderQR.h:146
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:397
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:388
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:158
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:315
bool isInvertible() const
Definition: FullPivHouseholderQR.h:270
Index rank() const
Definition: FullPivHouseholderQR.h:215