Eigen  3.2.92
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 traits<FullPivHouseholderQR<_MatrixType> >
19  : traits<_MatrixType>
20 {
21  enum { Flags = 0 };
22 };
23 
24 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
25 
26 template<typename MatrixType>
27 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
28 {
29  typedef typename MatrixType::PlainObject ReturnType;
30 };
31 
32 } // end namespace internal
33 
55 template<typename _MatrixType> class FullPivHouseholderQR
56 {
57  public:
58 
59  typedef _MatrixType MatrixType;
60  enum {
61  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
62  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
63  Options = MatrixType::Options,
64  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
66  };
67  typedef typename MatrixType::Scalar Scalar;
68  typedef typename MatrixType::RealScalar RealScalar;
69  // FIXME should be int
70  typedef typename MatrixType::StorageIndex StorageIndex;
71  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
72  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
73  typedef Matrix<StorageIndex, 1,
74  EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
75  EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
76  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
77  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
78  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
79  typedef typename MatrixType::PlainObject PlainObject;
80 
87  : m_qr(),
88  m_hCoeffs(),
89  m_rows_transpositions(),
90  m_cols_transpositions(),
91  m_cols_permutation(),
92  m_temp(),
93  m_isInitialized(false),
94  m_usePrescribedThreshold(false) {}
95 
102  FullPivHouseholderQR(Index rows, Index cols)
103  : m_qr(rows, cols),
104  m_hCoeffs((std::min)(rows,cols)),
105  m_rows_transpositions((std::min)(rows,cols)),
106  m_cols_transpositions((std::min)(rows,cols)),
107  m_cols_permutation(cols),
108  m_temp(cols),
109  m_isInitialized(false),
110  m_usePrescribedThreshold(false) {}
111 
124  template<typename InputType>
126  : m_qr(matrix.rows(), matrix.cols()),
127  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
128  m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
129  m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
130  m_cols_permutation(matrix.cols()),
131  m_temp(matrix.cols()),
132  m_isInitialized(false),
133  m_usePrescribedThreshold(false)
134  {
135  compute(matrix.derived());
136  }
137 
156  template<typename Rhs>
158  solve(const MatrixBase<Rhs>& b) const
159  {
160  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
161  return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
162  }
163 
166  MatrixQReturnType matrixQ(void) const;
167 
170  const MatrixType& matrixQR() const
171  {
172  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
173  return m_qr;
174  }
175 
176  template<typename InputType>
177  FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
178 
180  const PermutationType& colsPermutation() const
181  {
182  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
183  return m_cols_permutation;
184  }
185 
187  const IntDiagSizeVectorType& rowsTranspositions() const
188  {
189  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
190  return m_rows_transpositions;
191  }
192 
206  typename MatrixType::RealScalar absDeterminant() const;
207 
220  typename MatrixType::RealScalar logAbsDeterminant() const;
221 
228  inline Index rank() const
229  {
230  using std::abs;
231  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
232  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
233  Index result = 0;
234  for(Index i = 0; i < m_nonzero_pivots; ++i)
235  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
236  return result;
237  }
238 
245  inline Index dimensionOfKernel() const
246  {
247  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
248  return cols() - rank();
249  }
250 
258  inline bool isInjective() const
259  {
260  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
261  return rank() == cols();
262  }
263 
271  inline bool isSurjective() const
272  {
273  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
274  return rank() == rows();
275  }
276 
283  inline bool isInvertible() const
284  {
285  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
286  return isInjective() && isSurjective();
287  }
288 
295  {
296  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
297  return Inverse<FullPivHouseholderQR>(*this);
298  }
299 
300  inline Index rows() const { return m_qr.rows(); }
301  inline Index cols() const { return m_qr.cols(); }
302 
307  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
308 
326  FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
327  {
328  m_usePrescribedThreshold = true;
329  m_prescribedThreshold = threshold;
330  return *this;
331  }
332 
342  {
343  m_usePrescribedThreshold = false;
344  return *this;
345  }
346 
351  RealScalar threshold() const
352  {
353  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
354  return m_usePrescribedThreshold ? m_prescribedThreshold
355  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
356  // and turns out to be identical to Higham's formula used already in LDLt.
357  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
358  }
359 
367  inline Index nonzeroPivots() const
368  {
369  eigen_assert(m_isInitialized && "LU is not initialized.");
370  return m_nonzero_pivots;
371  }
372 
376  RealScalar maxPivot() const { return m_maxpivot; }
377 
378  #ifndef EIGEN_PARSED_BY_DOXYGEN
379  template<typename RhsType, typename DstType>
380  EIGEN_DEVICE_FUNC
381  void _solve_impl(const RhsType &rhs, DstType &dst) const;
382  #endif
383 
384  protected:
385 
386  static void check_template_parameters()
387  {
388  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
389  }
390 
391  void computeInPlace();
392 
393  MatrixType m_qr;
394  HCoeffsType m_hCoeffs;
395  IntDiagSizeVectorType m_rows_transpositions;
396  IntDiagSizeVectorType m_cols_transpositions;
397  PermutationType m_cols_permutation;
398  RowVectorType m_temp;
399  bool m_isInitialized, m_usePrescribedThreshold;
400  RealScalar m_prescribedThreshold, m_maxpivot;
401  Index m_nonzero_pivots;
402  RealScalar m_precision;
403  Index m_det_pq;
404 };
405 
406 template<typename MatrixType>
407 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
408 {
409  using std::abs;
410  eigen_assert(m_isInitialized && "FullPivHouseholderQR 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 abs(m_qr.diagonal().prod());
413 }
414 
415 template<typename MatrixType>
416 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
417 {
418  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
419  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
420  return m_qr.diagonal().cwiseAbs().array().log().sum();
421 }
422 
429 template<typename MatrixType>
430 template<typename InputType>
432 {
433  check_template_parameters();
434 
435  m_qr = matrix.derived();
436 
437  computeInPlace();
438 
439  return *this;
440 }
441 
442 template<typename MatrixType>
444 {
445  using std::abs;
446  Index rows = m_qr.rows();
447  Index cols = m_qr.cols();
448  Index size = (std::min)(rows,cols);
449 
450 
451  m_hCoeffs.resize(size);
452 
453  m_temp.resize(cols);
454 
455  m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
456 
457  m_rows_transpositions.resize(size);
458  m_cols_transpositions.resize(size);
459  Index number_of_transpositions = 0;
460 
461  RealScalar biggest(0);
462 
463  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
464  m_maxpivot = RealScalar(0);
465 
466  for (Index k = 0; k < size; ++k)
467  {
468  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
469  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
470  typedef typename Scoring::result_type Score;
471 
472  Score score = m_qr.bottomRightCorner(rows-k, cols-k)
473  .unaryExpr(Scoring())
474  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
475  row_of_biggest_in_corner += k;
476  col_of_biggest_in_corner += k;
477  RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
478  if(k==0) biggest = biggest_in_corner;
479 
480  // if the corner is negligible, then we have less than full rank, and we can finish early
481  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
482  {
483  m_nonzero_pivots = k;
484  for(Index i = k; i < size; i++)
485  {
486  m_rows_transpositions.coeffRef(i) = i;
487  m_cols_transpositions.coeffRef(i) = i;
488  m_hCoeffs.coeffRef(i) = Scalar(0);
489  }
490  break;
491  }
492 
493  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
494  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
495  if(k != row_of_biggest_in_corner) {
496  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
497  ++number_of_transpositions;
498  }
499  if(k != col_of_biggest_in_corner) {
500  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
501  ++number_of_transpositions;
502  }
503 
504  RealScalar beta;
505  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
506  m_qr.coeffRef(k,k) = beta;
507 
508  // remember the maximum absolute value of diagonal coefficients
509  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
510 
511  m_qr.bottomRightCorner(rows-k, cols-k-1)
512  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
513  }
514 
515  m_cols_permutation.setIdentity(cols);
516  for(Index k = 0; k < size; ++k)
517  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
518 
519  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
520  m_isInitialized = true;
521 }
522 
523 #ifndef EIGEN_PARSED_BY_DOXYGEN
524 template<typename _MatrixType>
525 template<typename RhsType, typename DstType>
526 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
527 {
528  eigen_assert(rhs.rows() == rows());
529  const Index l_rank = rank();
530 
531  // FIXME introduce nonzeroPivots() and use it here. and more generally,
532  // make the same improvements in this dec as in FullPivLU.
533  if(l_rank==0)
534  {
535  dst.setZero();
536  return;
537  }
538 
539  typename RhsType::PlainObject c(rhs);
540 
542  for (Index k = 0; k < l_rank; ++k)
543  {
544  Index remainingSize = rows()-k;
545  c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
546  c.bottomRightCorner(remainingSize, rhs.cols())
547  .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
548  m_hCoeffs.coeff(k), &temp.coeffRef(0));
549  }
550 
551  m_qr.topLeftCorner(l_rank, l_rank)
552  .template triangularView<Upper>()
553  .solveInPlace(c.topRows(l_rank));
554 
555  for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
556  for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
557 }
558 #endif
559 
560 namespace internal {
561 
562 template<typename DstXprType, typename MatrixType, typename Scalar>
563 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
564 {
565  typedef FullPivHouseholderQR<MatrixType> QrType;
566  typedef Inverse<QrType> SrcXprType;
567  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
568  {
569  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
570  }
571 };
572 
579 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
580  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
581 {
582 public:
584  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
585  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
586  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
587 
588  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
589  const HCoeffsType& hCoeffs,
590  const IntDiagSizeVectorType& rowsTranspositions)
591  : m_qr(qr),
592  m_hCoeffs(hCoeffs),
593  m_rowsTranspositions(rowsTranspositions)
594  {}
595 
596  template <typename ResultType>
597  void evalTo(ResultType& result) const
598  {
599  const Index rows = m_qr.rows();
600  WorkVectorType workspace(rows);
601  evalTo(result, workspace);
602  }
603 
604  template <typename ResultType>
605  void evalTo(ResultType& result, WorkVectorType& workspace) const
606  {
607  using numext::conj;
608  // compute the product H'_0 H'_1 ... H'_n-1,
609  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
610  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
611  const Index rows = m_qr.rows();
612  const Index cols = m_qr.cols();
613  const Index size = (std::min)(rows, cols);
614  workspace.resize(rows);
615  result.setIdentity(rows, rows);
616  for (Index k = size-1; k >= 0; k--)
617  {
618  result.block(k, k, rows-k, rows-k)
619  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
620  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
621  }
622  }
623 
624  Index rows() const { return m_qr.rows(); }
625  Index cols() const { return m_qr.rows(); }
626 
627 protected:
628  typename MatrixType::Nested m_qr;
629  typename HCoeffsType::Nested m_hCoeffs;
630  typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
631 };
632 
633 // template<typename MatrixType>
634 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
635 // : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
636 // {};
637 
638 } // end namespace internal
639 
640 template<typename MatrixType>
641 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
642 {
643  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
644  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
645 }
646 
647 #ifndef __CUDACC__
648 
652 template<typename Derived>
655 {
656  return FullPivHouseholderQR<PlainObject>(eval());
657 }
658 #endif // __CUDACC__
659 
660 } // end namespace Eigen
661 
662 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:245
bool isInvertible() const
Definition: FullPivHouseholderQR.h:283
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:341
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: ForwardDeclarations.h:254
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:641
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:376
Definition: LDLT.h:16
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:187
bool isInjective() const
Definition: FullPivHouseholderQR.h:258
Definition: StdDeque.h:58
bool isSurjective() const
Definition: FullPivHouseholderQR.h:271
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
RowXpr row(Index i)
Definition: DenseBase.h:797
Derived & derived()
Definition: EigenBase.h:44
Definition: EigenBase.h:28
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:102
Expression of the inverse of another expression.
Definition: Inverse.h:43
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivHouseholderQR.h:158
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:307
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:86
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:180
Definition: Constants.h:322
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:326
Definition: Eigen_Colamd.h:54
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:125
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:416
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:351
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:407
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:170
Index rank() const
Definition: FullPivHouseholderQR.h:228
Pseudo expression representing a solving operation.
Definition: Solve.h:62
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:654
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:294
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:367