FullPivHouseholderQR.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_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>
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;
70  // FIXME should be int
71  typedef typename MatrixType::StorageIndex StorageIndex;
74  typedef Matrix<StorageIndex, 1,
75  EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
76  EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
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 
174  template<typename Rhs>
176  solve(const MatrixBase<Rhs>& b) const
177  {
178  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
179  return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
180  }
181 
184  MatrixQReturnType matrixQ(void) const;
185 
188  const MatrixType& matrixQR() const
189  {
190  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
191  return m_qr;
192  }
193 
194  template<typename InputType>
196 
198  const PermutationType& colsPermutation() const
199  {
200  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
201  return m_cols_permutation;
202  }
203 
206  {
207  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
208  return m_rows_transpositions;
209  }
210 
224  typename MatrixType::RealScalar absDeterminant() const;
225 
238  typename MatrixType::RealScalar logAbsDeterminant() const;
239 
246  inline Index rank() const
247  {
248  using std::abs;
249  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
250  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
251  Index result = 0;
252  for(Index i = 0; i < m_nonzero_pivots; ++i)
253  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
254  return result;
255  }
256 
263  inline Index dimensionOfKernel() const
264  {
265  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
266  return cols() - rank();
267  }
268 
276  inline bool isInjective() const
277  {
278  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
279  return rank() == cols();
280  }
281 
289  inline bool isSurjective() const
290  {
291  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
292  return rank() == rows();
293  }
294 
301  inline bool isInvertible() const
302  {
303  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
304  return isInjective() && isSurjective();
305  }
306 
313  {
314  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
315  return Inverse<FullPivHouseholderQR>(*this);
316  }
317 
318  inline Index rows() const { return m_qr.rows(); }
319  inline Index cols() const { return m_qr.cols(); }
320 
325  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
326 
344  FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
345  {
346  m_usePrescribedThreshold = true;
347  m_prescribedThreshold = threshold;
348  return *this;
349  }
350 
360  {
361  m_usePrescribedThreshold = false;
362  return *this;
363  }
364 
369  RealScalar threshold() const
370  {
371  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
372  return m_usePrescribedThreshold ? m_prescribedThreshold
373  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
374  // and turns out to be identical to Higham's formula used already in LDLt.
375  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
376  }
377 
385  inline Index nonzeroPivots() const
386  {
387  eigen_assert(m_isInitialized && "LU is not initialized.");
388  return m_nonzero_pivots;
389  }
390 
394  RealScalar maxPivot() const { return m_maxpivot; }
395 
396  #ifndef EIGEN_PARSED_BY_DOXYGEN
397  template<typename RhsType, typename DstType>
398  EIGEN_DEVICE_FUNC
399  void _solve_impl(const RhsType &rhs, DstType &dst) const;
400  #endif
401 
402  protected:
403 
405  {
407  }
408 
409  void computeInPlace();
410 
411  MatrixType m_qr;
412  HCoeffsType m_hCoeffs;
415  PermutationType m_cols_permutation;
416  RowVectorType m_temp;
417  bool m_isInitialized, m_usePrescribedThreshold;
418  RealScalar m_prescribedThreshold, m_maxpivot;
420  RealScalar m_precision;
422 };
423 
424 template<typename MatrixType>
426 {
427  using std::abs;
428  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
429  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
430  return abs(m_qr.diagonal().prod());
431 }
432 
433 template<typename MatrixType>
435 {
436  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
437  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
438  return m_qr.diagonal().cwiseAbs().array().log().sum();
439 }
440 
447 template<typename MatrixType>
448 template<typename InputType>
450 {
451  m_qr = matrix.derived();
452  computeInPlace();
453  return *this;
454 }
455 
456 template<typename MatrixType>
458 {
459  check_template_parameters();
460 
461  using std::abs;
462  Index rows = m_qr.rows();
463  Index cols = m_qr.cols();
464  Index size = (std::min)(rows,cols);
465 
466 
467  m_hCoeffs.resize(size);
468 
469  m_temp.resize(cols);
470 
471  m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
472 
473  m_rows_transpositions.resize(size);
474  m_cols_transpositions.resize(size);
475  Index number_of_transpositions = 0;
476 
477  RealScalar biggest(0);
478 
479  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
480  m_maxpivot = RealScalar(0);
481 
482  for (Index k = 0; k < size; ++k)
483  {
484  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
486  typedef typename Scoring::result_type Score;
487 
488  Score score = m_qr.bottomRightCorner(rows-k, cols-k)
489  .unaryExpr(Scoring())
490  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
491  row_of_biggest_in_corner += k;
492  col_of_biggest_in_corner += k;
493  RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
494  if(k==0) biggest = biggest_in_corner;
495 
496  // if the corner is negligible, then we have less than full rank, and we can finish early
497  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
498  {
499  m_nonzero_pivots = k;
500  for(Index i = k; i < size; i++)
501  {
502  m_rows_transpositions.coeffRef(i) = i;
503  m_cols_transpositions.coeffRef(i) = i;
504  m_hCoeffs.coeffRef(i) = Scalar(0);
505  }
506  break;
507  }
508 
509  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
510  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
511  if(k != row_of_biggest_in_corner) {
512  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
513  ++number_of_transpositions;
514  }
515  if(k != col_of_biggest_in_corner) {
516  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
517  ++number_of_transpositions;
518  }
519 
520  RealScalar beta;
521  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
522  m_qr.coeffRef(k,k) = beta;
523 
524  // remember the maximum absolute value of diagonal coefficients
525  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
526 
527  m_qr.bottomRightCorner(rows-k, cols-k-1)
528  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
529  }
530 
531  m_cols_permutation.setIdentity(cols);
532  for(Index k = 0; k < size; ++k)
533  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
534 
535  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
536  m_isInitialized = true;
537 }
538 
539 #ifndef EIGEN_PARSED_BY_DOXYGEN
540 template<typename _MatrixType>
541 template<typename RhsType, typename DstType>
542 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
543 {
544  eigen_assert(rhs.rows() == rows());
545  const Index l_rank = rank();
546 
547  // FIXME introduce nonzeroPivots() and use it here. and more generally,
548  // make the same improvements in this dec as in FullPivLU.
549  if(l_rank==0)
550  {
551  dst.setZero();
552  return;
553  }
554 
555  typename RhsType::PlainObject c(rhs);
556 
558  for (Index k = 0; k < l_rank; ++k)
559  {
560  Index remainingSize = rows()-k;
561  c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
562  c.bottomRightCorner(remainingSize, rhs.cols())
563  .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
564  m_hCoeffs.coeff(k), &temp.coeffRef(0));
565  }
566 
567  m_qr.topLeftCorner(l_rank, l_rank)
568  .template triangularView<Upper>()
569  .solveInPlace(c.topRows(l_rank));
570 
571  for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
572  for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
573 }
574 #endif
575 
576 namespace internal {
577 
578 template<typename DstXprType, typename MatrixType>
579 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
580 {
583  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
584  {
585  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
586  }
587 };
588 
595 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
596  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
597 {
598 public:
601  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
602  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
603 
605  const HCoeffsType& hCoeffs,
606  const IntDiagSizeVectorType& rowsTranspositions)
607  : m_qr(qr),
608  m_hCoeffs(hCoeffs),
609  m_rowsTranspositions(rowsTranspositions)
610  {}
611 
612  template <typename ResultType>
613  void evalTo(ResultType& result) const
614  {
615  const Index rows = m_qr.rows();
616  WorkVectorType workspace(rows);
617  evalTo(result, workspace);
618  }
619 
620  template <typename ResultType>
621  void evalTo(ResultType& result, WorkVectorType& workspace) const
622  {
623  using numext::conj;
624  // compute the product H'_0 H'_1 ... H'_n-1,
625  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
626  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
627  const Index rows = m_qr.rows();
628  const Index cols = m_qr.cols();
629  const Index size = (std::min)(rows, cols);
630  workspace.resize(rows);
631  result.setIdentity(rows, rows);
632  for (Index k = size-1; k >= 0; k--)
633  {
634  result.block(k, k, rows-k, rows-k)
635  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
636  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
637  }
638  }
639 
640  Index rows() const { return m_qr.rows(); }
641  Index cols() const { return m_qr.rows(); }
642 
643 protected:
644  typename MatrixType::Nested m_qr;
645  typename HCoeffsType::Nested m_hCoeffs;
646  typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
647 };
648 
649 // template<typename MatrixType>
650 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
651 // : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
652 // {};
653 
654 } // end namespace internal
655 
656 template<typename MatrixType>
658 {
659  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
660  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
661 }
662 
667 template<typename Derived>
670 {
672 }
673 
674 } // end namespace Eigen
675 
676 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
SCALAR Scalar
Definition: bench_gemm.cpp:33
FullPivHouseholderQR & setThreshold(Default_t)
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
EIGEN_DEVICE_FUNC void _solve_impl(const RhsType &rhs, DstType &dst) const
EIGEN_DEVICE_FUNC Index rows() const
Definition: Inverse.h:58
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Scalar * b
Definition: benchVecAdd.cpp:17
MatrixType::StorageIndex StorageIndex
#define min(a, b)
Definition: datatypes.h:19
MatrixQReturnType matrixQ(void) const
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: Half.h:150
MatrixXf MatrixType
Default_t
Definition: Constants.h:352
MatrixType::RealScalar RealScalar
FullPivHouseholderQRMatrixQReturnType(const MatrixType &qr, const HCoeffsType &hCoeffs, const IntDiagSizeVectorType &rowsTranspositions)
#define EIGEN_SIZE_MIN_PREFER_FIXED(a, b)
Definition: Macros.h:889
void evalTo(ResultType &result, WorkVectorType &workspace) const
HouseholderQR< MatrixXf > qr(A)
static double epsilon
Definition: testRot3.cpp:39
FullPivHouseholderQR & compute(const EigenBase< InputType > &matrix)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
internal::plain_diag_type< MatrixType >::type HCoeffsType
Expression of the inverse of another expression.
Definition: Inverse.h:43
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Matrix< typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, MatrixType::MaxRowsAtCompileTime > WorkVectorType
const HCoeffsType & hCoeffs() const
FullPivHouseholderQR()
Default Constructor.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
internal::plain_row_type< MatrixType >::type RowVectorType
const PermutationType & colsPermutation() const
Values result
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
IntDiagSizeVectorType m_rows_transpositions
EIGEN_DEVICE_FUNC Index cols() const
Definition: Inverse.h:59
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:184
internal::FullPivHouseholderQRMatrixQReturnType< MatrixType > MatrixQReturnType
Expression type for return value of FullPivHouseholderQR::matrixQ()
FullPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
internal::plain_col_type< MatrixType >::type ColVectorType
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
PermutationMatrix< ColsAtCompileTime, MaxColsAtCompileTime > PermutationType
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
MatrixType::RealScalar logAbsDeterminant() const
FullPivHouseholderQR< MatrixType >::IntDiagSizeVectorType IntDiagSizeVectorType
Matrix< StorageIndex, 1, EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime, RowsAtCompileTime), RowMajor, 1, EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime, MaxRowsAtCompileTime)> IntDiagSizeVectorType
MatrixType::RealScalar absDeterminant() const
const MatrixType & matrixQR() const
MatrixType::PlainObject PlainObject
IntDiagSizeVectorType m_cols_transpositions
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Pseudo expression representing a solving operation.
Definition: Solve.h:62
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
The matrix class, also used for vectors and row-vectors.
EIGEN_DEVICE_FUNC const XprTypeNestedCleaned & nestedExpression() const
Definition: Inverse.h:61
#define EIGEN_SIZE_MIN_PREFER_DYNAMIC(a, b)
Definition: Macros.h:881
#define abs(x)
Definition: datatypes.h:17
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
internal::plain_diag_type< MatrixType >::type HCoeffsType
ScalarWithExceptions conj(const ScalarWithExceptions &x)
Definition: exceptions.cpp:74
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
const Inverse< FullPivHouseholderQR > inverse() const
v setZero(3)


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:04