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 
16 namespace internal {
17 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
18  : traits<_MatrixType>
19 {
20  typedef MatrixXpr XprKind;
22  typedef int StorageIndex;
23  enum { Flags = 0 };
24 };
25 
26 } // end namespace internal
27 
51 template<typename _MatrixType> class ColPivHouseholderQR
52  : public SolverBase<ColPivHouseholderQR<_MatrixType> >
53 {
54  public:
55 
56  typedef _MatrixType MatrixType;
59 
61  enum {
62  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
63  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64  };
71  typedef typename MatrixType::PlainObject PlainObject;
72 
73  private:
74 
76 
77  public:
78 
86  : m_qr(),
87  m_hCoeffs(),
88  m_colsPermutation(),
89  m_colsTranspositions(),
90  m_temp(),
91  m_colNormsUpdated(),
92  m_colNormsDirect(),
93  m_isInitialized(false),
94  m_usePrescribedThreshold(false) {}
95 
103  : m_qr(rows, cols),
104  m_hCoeffs((std::min)(rows,cols)),
105  m_colsPermutation(PermIndexType(cols)),
106  m_colsTranspositions(cols),
107  m_temp(cols),
108  m_colNormsUpdated(cols),
109  m_colNormsDirect(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_colsPermutation(PermIndexType(matrix.cols())),
130  m_colsTranspositions(matrix.cols()),
131  m_temp(matrix.cols()),
132  m_colNormsUpdated(matrix.cols()),
133  m_colNormsDirect(matrix.cols()),
134  m_isInitialized(false),
135  m_usePrescribedThreshold(false)
136  {
137  compute(matrix.derived());
138  }
139 
146  template<typename InputType>
148  : m_qr(matrix.derived()),
149  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
150  m_colsPermutation(PermIndexType(matrix.cols())),
151  m_colsTranspositions(matrix.cols()),
152  m_temp(matrix.cols()),
153  m_colNormsUpdated(matrix.cols()),
154  m_colNormsDirect(matrix.cols()),
155  m_isInitialized(false),
156  m_usePrescribedThreshold(false)
157  {
158  computeInPlace();
159  }
160 
161  #ifdef EIGEN_PARSED_BY_DOXYGEN
162 
176  template<typename Rhs>
178  solve(const MatrixBase<Rhs>& b) const;
179  #endif
180 
181  HouseholderSequenceType householderQ() const;
182  HouseholderSequenceType matrixQ() const
183  {
184  return householderQ();
185  }
186 
189  const MatrixType& matrixQR() const
190  {
191  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
192  return m_qr;
193  }
194 
204  const MatrixType& matrixR() const
205  {
206  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
207  return m_qr;
208  }
209 
210  template<typename InputType>
212 
214  const PermutationType& colsPermutation() const
215  {
216  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
217  return m_colsPermutation;
218  }
219 
233  typename MatrixType::RealScalar absDeterminant() const;
234 
247  typename MatrixType::RealScalar logAbsDeterminant() const;
248 
255  inline Index rank() const
256  {
257  using std::abs;
258  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
259  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
260  Index result = 0;
261  for(Index i = 0; i < m_nonzero_pivots; ++i)
262  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
263  return result;
264  }
265 
272  inline Index dimensionOfKernel() const
273  {
274  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
275  return cols() - rank();
276  }
277 
285  inline bool isInjective() const
286  {
287  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
288  return rank() == cols();
289  }
290 
298  inline bool isSurjective() const
299  {
300  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
301  return rank() == rows();
302  }
303 
310  inline bool isInvertible() const
311  {
312  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
313  return isInjective() && isSurjective();
314  }
315 
322  {
323  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
324  return Inverse<ColPivHouseholderQR>(*this);
325  }
326 
327  inline Index rows() const { return m_qr.rows(); }
328  inline Index cols() const { return m_qr.cols(); }
329 
334  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
335 
354  {
355  m_usePrescribedThreshold = true;
356  m_prescribedThreshold = threshold;
357  return *this;
358  }
359 
369  {
370  m_usePrescribedThreshold = false;
371  return *this;
372  }
373 
379  {
380  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
381  return m_usePrescribedThreshold ? m_prescribedThreshold
382  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
383  // and turns out to be identical to Higham's formula used already in LDLt.
384  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
385  }
386 
394  inline Index nonzeroPivots() const
395  {
396  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
397  return m_nonzero_pivots;
398  }
399 
403  RealScalar maxPivot() const { return m_maxpivot; }
404 
412  {
413  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
414  return Success;
415  }
416 
417  #ifndef EIGEN_PARSED_BY_DOXYGEN
418  template<typename RhsType, typename DstType>
419  void _solve_impl(const RhsType &rhs, DstType &dst) const;
420 
421  template<bool Conjugate, typename RhsType, typename DstType>
422  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
423  #endif
424 
425  protected:
426 
427  friend class CompleteOrthogonalDecomposition<MatrixType>;
428 
430  {
432  }
433 
434  void computeInPlace();
435 
436  MatrixType m_qr;
437  HCoeffsType m_hCoeffs;
438  PermutationType m_colsPermutation;
439  IntRowVectorType m_colsTranspositions;
440  RowVectorType m_temp;
441  RealRowVectorType m_colNormsUpdated;
442  RealRowVectorType m_colNormsDirect;
443  bool m_isInitialized, m_usePrescribedThreshold;
447 };
448 
449 template<typename MatrixType>
451 {
452  using std::abs;
453  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
454  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
455  return abs(m_qr.diagonal().prod());
456 }
457 
458 template<typename MatrixType>
460 {
461  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
462  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
463  return m_qr.diagonal().cwiseAbs().array().log().sum();
464 }
465 
472 template<typename MatrixType>
473 template<typename InputType>
475 {
476  m_qr = matrix.derived();
477  computeInPlace();
478  return *this;
479 }
480 
481 template<typename MatrixType>
483 {
484  check_template_parameters();
485 
486  // the column permutation is stored as int indices, so just to be sure:
487  eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
488 
489  using std::abs;
490 
491  Index rows = m_qr.rows();
492  Index cols = m_qr.cols();
493  Index size = m_qr.diagonalSize();
494 
495  m_hCoeffs.resize(size);
496 
497  m_temp.resize(cols);
498 
499  m_colsTranspositions.resize(m_qr.cols());
500  Index number_of_transpositions = 0;
501 
502  m_colNormsUpdated.resize(cols);
503  m_colNormsDirect.resize(cols);
504  for (Index k = 0; k < cols; ++k) {
505  // colNormsDirect(k) caches the most recent directly computed norm of
506  // column k.
507  m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
508  m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
509  }
510 
511  RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
512  RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
513 
514  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
515  m_maxpivot = RealScalar(0);
516 
517  for(Index k = 0; k < size; ++k)
518  {
519  // first, we look up in our table m_colNormsUpdated which column has the biggest norm
520  Index biggest_col_index;
521  RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
522  biggest_col_index += k;
523 
524  // Track the number of meaningful pivots but do not stop the decomposition to make
525  // sure that the initial matrix is properly reproduced. See bug 941.
526  if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
527  m_nonzero_pivots = k;
528 
529  // apply the transposition to the columns
530  m_colsTranspositions.coeffRef(k) = biggest_col_index;
531  if(k != biggest_col_index) {
532  m_qr.col(k).swap(m_qr.col(biggest_col_index));
533  std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
534  std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
535  ++number_of_transpositions;
536  }
537 
538  // generate the householder vector, store it below the diagonal
539  RealScalar beta;
540  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
541 
542  // apply the householder transformation to the diagonal coefficient
543  m_qr.coeffRef(k,k) = beta;
544 
545  // remember the maximum absolute value of diagonal coefficients
546  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
547 
548  // apply the householder transformation
549  m_qr.bottomRightCorner(rows-k, cols-k-1)
550  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
551 
552  // update our table of norms of the columns
553  for (Index j = k + 1; j < cols; ++j) {
554  // The following implements the stable norm downgrade step discussed in
555  // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
556  // and used in LAPACK routines xGEQPF and xGEQP3.
557  // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
558  if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
559  RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
560  temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
561  temp = temp < RealScalar(0) ? RealScalar(0) : temp;
562  RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
563  m_colNormsDirect.coeffRef(j));
564  if (temp2 <= norm_downdate_threshold) {
565  // The updated norm has become too inaccurate so re-compute the column
566  // norm directly.
567  m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
568  m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
569  } else {
570  m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
571  }
572  }
573  }
574  }
575 
576  m_colsPermutation.setIdentity(PermIndexType(cols));
577  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
578  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
579 
580  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
581  m_isInitialized = true;
582 }
583 
584 #ifndef EIGEN_PARSED_BY_DOXYGEN
585 template<typename _MatrixType>
586 template<typename RhsType, typename DstType>
587 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
588 {
589  const Index nonzero_pivots = nonzeroPivots();
590 
591  if(nonzero_pivots == 0)
592  {
593  dst.setZero();
594  return;
595  }
596 
597  typename RhsType::PlainObject c(rhs);
598 
599  c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
600 
601  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
602  .template triangularView<Upper>()
603  .solveInPlace(c.topRows(nonzero_pivots));
604 
605  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
606  for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
607 }
608 
609 template<typename _MatrixType>
610 template<bool Conjugate, typename RhsType, typename DstType>
611 void ColPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
612 {
613  const Index nonzero_pivots = nonzeroPivots();
614 
615  if(nonzero_pivots == 0)
616  {
617  dst.setZero();
618  return;
619  }
620 
621  typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs);
622 
623  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
624  .template triangularView<Upper>()
625  .transpose().template conjugateIf<Conjugate>()
626  .solveInPlace(c.topRows(nonzero_pivots));
627 
628  dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
629  dst.bottomRows(rows()-nonzero_pivots).setZero();
630 
631  dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).template conjugateIf<!Conjugate>() );
632 }
633 #endif
634 
635 namespace internal {
636 
637 template<typename DstXprType, typename MatrixType>
638 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
639 {
642  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
643  {
644  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
645  }
646 };
647 
648 } // end namespace internal
649 
653 template<typename MatrixType>
656 {
657  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
658  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
659 }
660 
665 template<typename Derived>
668 {
670 }
671 
672 } // end namespace Eigen
673 
674 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
PermutationMatrix< ColsAtCompileTime, MaxColsAtCompileTime > PermutationType
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: Inverse.h:57
SolverBase< ColPivHouseholderQR > Base
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:1264
ColPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Scalar * b
Definition: benchVecAdd.cpp:17
ColPivHouseholderQR()
Default Constructor.
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
internal::traits< ColPivHouseholderQR< TransposeTypeWithSameStorageOrder > >::Scalar Scalar
Definition: SolverBase.h:73
const HCoeffsType & hCoeffs() const
PermutationType::StorageIndex PermIndexType
const ColPivHouseholderQR< PlainObject > colPivHouseholderQr() const
internal::plain_diag_type< MatrixType >::type HCoeffsType
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Definition: BFloat16.h:88
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Default_t
Definition: Constants.h:362
ComputationInfo info() const
Reports whether the QR factorization was successful.
Complete orthogonal decomposition (COD) of a matrix.
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
RealRowVectorType m_colNormsUpdated
static double epsilon
Definition: testRot3.cpp:37
Sequence of Householder reflections acting on subspaces with decreasing size.
ColPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Expression of the inverse of another expression.
Definition: Inverse.h:43
HouseholderSequenceType householderQ() const
ColPivHouseholderQR & compute(const EigenBase< InputType > &matrix)
ColPivHouseholderQR & setThreshold(Default_t)
internal::plain_row_type< MatrixType >::type RowVectorType
Values result
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
const MatrixType & matrixR() const
#define eigen_assert(x)
Definition: Macros.h:1037
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
void _solve_impl(const RhsType &rhs, DstType &dst) const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
MatrixType::PlainObject PlainObject
internal::plain_row_type< MatrixType, RealScalar >::type RealRowVectorType
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
MatrixType::RealScalar absDeterminant() const
RealRowVectorType m_colNormsDirect
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
HouseholderSequenceType matrixQ() const
MatrixType::RealScalar logAbsDeterminant() const
internal::plain_row_type< MatrixType, Index >::type IntRowVectorType
internal::nested_eval< T, 1 >::type eval(const T &xpr)
void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
Pseudo expression representing a solving operation.
Definition: Solve.h:62
IntRowVectorType m_colsTranspositions
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
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)
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sqrt(const float &x)
const MatrixType & matrixQR() const
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:68
const PermutationType & colsPermutation() const
EIGEN_DEVICE_FUNC bool abs2(bool x)
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Inverse< ColPivHouseholderQR > inverse() const
std::ptrdiff_t j
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: Inverse.h:58
EIGEN_DEVICE_FUNC const XprTypeNestedCleaned & nestedExpression() const
Definition: Inverse.h:60
v setZero(3)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:01