Eigen  3.2.93
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 
57 template<typename _MatrixType> class FullPivHouseholderQR
58 {
59  public:
60 
61  typedef _MatrixType MatrixType;
62  enum {
63  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
64  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
65  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
66  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67  };
68  typedef typename MatrixType::Scalar Scalar;
69  typedef typename MatrixType::RealScalar RealScalar;
70  // FIXME should be int
71  typedef typename MatrixType::StorageIndex StorageIndex;
72  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
73  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74  typedef Matrix<StorageIndex, 1,
75  EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
76  EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
77  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
78  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
79  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
80  typedef typename MatrixType::PlainObject PlainObject;
81 
88  : m_qr(),
89  m_hCoeffs(),
90  m_rows_transpositions(),
91  m_cols_transpositions(),
92  m_cols_permutation(),
93  m_temp(),
94  m_isInitialized(false),
95  m_usePrescribedThreshold(false) {}
96 
104  : m_qr(rows, cols),
105  m_hCoeffs((std::min)(rows,cols)),
106  m_rows_transpositions((std::min)(rows,cols)),
107  m_cols_transpositions((std::min)(rows,cols)),
108  m_cols_permutation(cols),
109  m_temp(cols),
110  m_isInitialized(false),
111  m_usePrescribedThreshold(false) {}
112 
125  template<typename InputType>
127  : m_qr(matrix.rows(), matrix.cols()),
128  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
129  m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
130  m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
131  m_cols_permutation(matrix.cols()),
132  m_temp(matrix.cols()),
133  m_isInitialized(false),
134  m_usePrescribedThreshold(false)
135  {
136  compute(matrix.derived());
137  }
138 
145  template<typename InputType>
147  : m_qr(matrix.derived()),
148  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
149  m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
150  m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
151  m_cols_permutation(matrix.cols()),
152  m_temp(matrix.cols()),
153  m_isInitialized(false),
154  m_usePrescribedThreshold(false)
155  {
156  computeInPlace();
157  }
158 
177  template<typename Rhs>
179  solve(const MatrixBase<Rhs>& b) const
180  {
181  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
182  return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
183  }
184 
187  MatrixQReturnType matrixQ(void) const;
188 
191  const MatrixType& matrixQR() const
192  {
193  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
194  return m_qr;
195  }
196 
197  template<typename InputType>
198  FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
199 
201  const PermutationType& colsPermutation() const
202  {
203  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
204  return m_cols_permutation;
205  }
206 
208  const IntDiagSizeVectorType& rowsTranspositions() const
209  {
210  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
211  return m_rows_transpositions;
212  }
213 
227  typename MatrixType::RealScalar absDeterminant() const;
228 
241  typename MatrixType::RealScalar logAbsDeterminant() const;
242 
249  inline Index rank() const
250  {
251  using std::abs;
252  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
253  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
254  Index result = 0;
255  for(Index i = 0; i < m_nonzero_pivots; ++i)
256  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
257  return result;
258  }
259 
266  inline Index dimensionOfKernel() const
267  {
268  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
269  return cols() - rank();
270  }
271 
279  inline bool isInjective() const
280  {
281  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
282  return rank() == cols();
283  }
284 
292  inline bool isSurjective() const
293  {
294  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
295  return rank() == rows();
296  }
297 
304  inline bool isInvertible() const
305  {
306  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
307  return isInjective() && isSurjective();
308  }
309 
316  {
317  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
318  return Inverse<FullPivHouseholderQR>(*this);
319  }
320 
321  inline Index rows() const { return m_qr.rows(); }
322  inline Index cols() const { return m_qr.cols(); }
323 
328  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
329 
347  FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
348  {
349  m_usePrescribedThreshold = true;
350  m_prescribedThreshold = threshold;
351  return *this;
352  }
353 
363  {
364  m_usePrescribedThreshold = false;
365  return *this;
366  }
367 
372  RealScalar threshold() const
373  {
374  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
375  return m_usePrescribedThreshold ? m_prescribedThreshold
376  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
377  // and turns out to be identical to Higham's formula used already in LDLt.
378  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
379  }
380 
388  inline Index nonzeroPivots() const
389  {
390  eigen_assert(m_isInitialized && "LU is not initialized.");
391  return m_nonzero_pivots;
392  }
393 
397  RealScalar maxPivot() const { return m_maxpivot; }
398 
399  #ifndef EIGEN_PARSED_BY_DOXYGEN
400  template<typename RhsType, typename DstType>
401  EIGEN_DEVICE_FUNC
402  void _solve_impl(const RhsType &rhs, DstType &dst) const;
403  #endif
404 
405  protected:
406 
407  static void check_template_parameters()
408  {
409  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
410  }
411 
412  void computeInPlace();
413 
414  MatrixType m_qr;
415  HCoeffsType m_hCoeffs;
416  IntDiagSizeVectorType m_rows_transpositions;
417  IntDiagSizeVectorType m_cols_transpositions;
418  PermutationType m_cols_permutation;
419  RowVectorType m_temp;
420  bool m_isInitialized, m_usePrescribedThreshold;
421  RealScalar m_prescribedThreshold, m_maxpivot;
422  Index m_nonzero_pivots;
423  RealScalar m_precision;
424  Index m_det_pq;
425 };
426 
427 template<typename MatrixType>
428 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
429 {
430  using std::abs;
431  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
432  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
433  return abs(m_qr.diagonal().prod());
434 }
435 
436 template<typename MatrixType>
437 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
438 {
439  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
440  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
441  return m_qr.diagonal().cwiseAbs().array().log().sum();
442 }
443 
450 template<typename MatrixType>
451 template<typename InputType>
453 {
454  m_qr = matrix.derived();
455  computeInPlace();
456  return *this;
457 }
458 
459 template<typename MatrixType>
461 {
462  check_template_parameters();
463 
464  using std::abs;
465  Index rows = m_qr.rows();
466  Index cols = m_qr.cols();
467  Index size = (std::min)(rows,cols);
468 
469 
470  m_hCoeffs.resize(size);
471 
472  m_temp.resize(cols);
473 
474  m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
475 
476  m_rows_transpositions.resize(size);
477  m_cols_transpositions.resize(size);
478  Index number_of_transpositions = 0;
479 
480  RealScalar biggest(0);
481 
482  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
483  m_maxpivot = RealScalar(0);
484 
485  for (Index k = 0; k < size; ++k)
486  {
487  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
488  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
489  typedef typename Scoring::result_type Score;
490 
491  Score score = m_qr.bottomRightCorner(rows-k, cols-k)
492  .unaryExpr(Scoring())
493  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
494  row_of_biggest_in_corner += k;
495  col_of_biggest_in_corner += k;
496  RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
497  if(k==0) biggest = biggest_in_corner;
498 
499  // if the corner is negligible, then we have less than full rank, and we can finish early
500  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
501  {
502  m_nonzero_pivots = k;
503  for(Index i = k; i < size; i++)
504  {
505  m_rows_transpositions.coeffRef(i) = i;
506  m_cols_transpositions.coeffRef(i) = i;
507  m_hCoeffs.coeffRef(i) = Scalar(0);
508  }
509  break;
510  }
511 
512  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
513  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
514  if(k != row_of_biggest_in_corner) {
515  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
516  ++number_of_transpositions;
517  }
518  if(k != col_of_biggest_in_corner) {
519  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
520  ++number_of_transpositions;
521  }
522 
523  RealScalar beta;
524  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
525  m_qr.coeffRef(k,k) = beta;
526 
527  // remember the maximum absolute value of diagonal coefficients
528  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
529 
530  m_qr.bottomRightCorner(rows-k, cols-k-1)
531  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
532  }
533 
534  m_cols_permutation.setIdentity(cols);
535  for(Index k = 0; k < size; ++k)
536  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
537 
538  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
539  m_isInitialized = true;
540 }
541 
542 #ifndef EIGEN_PARSED_BY_DOXYGEN
543 template<typename _MatrixType>
544 template<typename RhsType, typename DstType>
545 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
546 {
547  eigen_assert(rhs.rows() == rows());
548  const Index l_rank = rank();
549 
550  // FIXME introduce nonzeroPivots() and use it here. and more generally,
551  // make the same improvements in this dec as in FullPivLU.
552  if(l_rank==0)
553  {
554  dst.setZero();
555  return;
556  }
557 
558  typename RhsType::PlainObject c(rhs);
559 
561  for (Index k = 0; k < l_rank; ++k)
562  {
563  Index remainingSize = rows()-k;
564  c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
565  c.bottomRightCorner(remainingSize, rhs.cols())
566  .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
567  m_hCoeffs.coeff(k), &temp.coeffRef(0));
568  }
569 
570  m_qr.topLeftCorner(l_rank, l_rank)
571  .template triangularView<Upper>()
572  .solveInPlace(c.topRows(l_rank));
573 
574  for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
575  for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
576 }
577 #endif
578 
579 namespace internal {
580 
581 template<typename DstXprType, typename MatrixType, typename Scalar>
582 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar,Scalar>, Dense2Dense>
583 {
584  typedef FullPivHouseholderQR<MatrixType> QrType;
585  typedef Inverse<QrType> SrcXprType;
586  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
587  {
588  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
589  }
590 };
591 
598 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
599  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
600 {
601 public:
603  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
604  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
605  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
606 
607  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
608  const HCoeffsType& hCoeffs,
609  const IntDiagSizeVectorType& rowsTranspositions)
610  : m_qr(qr),
611  m_hCoeffs(hCoeffs),
612  m_rowsTranspositions(rowsTranspositions)
613  {}
614 
615  template <typename ResultType>
616  void evalTo(ResultType& result) const
617  {
618  const Index rows = m_qr.rows();
619  WorkVectorType workspace(rows);
620  evalTo(result, workspace);
621  }
622 
623  template <typename ResultType>
624  void evalTo(ResultType& result, WorkVectorType& workspace) const
625  {
626  using numext::conj;
627  // compute the product H'_0 H'_1 ... H'_n-1,
628  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
629  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
630  const Index rows = m_qr.rows();
631  const Index cols = m_qr.cols();
632  const Index size = (std::min)(rows, cols);
633  workspace.resize(rows);
634  result.setIdentity(rows, rows);
635  for (Index k = size-1; k >= 0; k--)
636  {
637  result.block(k, k, rows-k, rows-k)
638  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
639  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
640  }
641  }
642 
643  Index rows() const { return m_qr.rows(); }
644  Index cols() const { return m_qr.rows(); }
645 
646 protected:
647  typename MatrixType::Nested m_qr;
648  typename HCoeffsType::Nested m_hCoeffs;
649  typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
650 };
651 
652 // template<typename MatrixType>
653 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
654 // : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
655 // {};
656 
657 } // end namespace internal
658 
659 template<typename MatrixType>
660 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
661 {
662  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
663  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
664 }
665 
666 #ifndef __CUDACC__
667 
671 template<typename Derived>
674 {
675  return FullPivHouseholderQR<PlainObject>(eval());
676 }
677 #endif // __CUDACC__
678 
679 } // end namespace Eigen
680 
681 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:362
bool isInvertible() const
Definition: FullPivHouseholderQR.h:304
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:208
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: ForwardDeclarations.h:256
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivHouseholderQR.h:179
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:191
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:437
Namespace containing all symbols from the Eigen library.
Definition: Core:271
bool isSurjective() const
Definition: FullPivHouseholderQR.h:292
Definition: Half.h:502
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:167
RowXpr row(Index i)
Definition: DenseBase.h:802
Index rank() const
Definition: FullPivHouseholderQR.h:249
Derived & derived()
Definition: EigenBase.h:44
Definition: EigenBase.h:28
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:103
Expression of the inverse of another expression.
Definition: Inverse.h:43
bool isInjective() const
Definition: FullPivHouseholderQR.h:279
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:660
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:87
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
FullPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:146
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:347
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:201
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:126
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:315
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:372
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:397
Definition: Constants.h:322
Pseudo expression representing a solving operation.
Definition: Solve.h:62
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:428
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:328
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:673
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:388
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:266