Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
ColPivHouseholderQR.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_COLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
13 
14 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 template<typename MatrixType_> struct traits<ColPivHouseholderQR<MatrixType_> >
20  : traits<MatrixType_>
21 {
22  typedef MatrixXpr XprKind;
23  typedef SolverStorage StorageKind;
24  typedef int StorageIndex;
25  enum { Flags = 0 };
26 };
27 
28 } // end namespace internal
29 
53 template<typename MatrixType_> class ColPivHouseholderQR
54  : public SolverBase<ColPivHouseholderQR<MatrixType_> >
55 {
56  public:
57 
58  typedef MatrixType_ MatrixType;
60  friend class SolverBase<ColPivHouseholderQR>;
61 
62  EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR)
63  enum {
64  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
66  };
67  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
69  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
70  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
71  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
73  typedef typename MatrixType::PlainObject PlainObject;
74 
75  private:
76 
77  typedef typename PermutationType::StorageIndex PermIndexType;
78 
79  public:
80 
88  : m_qr(),
89  m_hCoeffs(),
90  m_colsPermutation(),
91  m_colsTranspositions(),
92  m_temp(),
93  m_colNormsUpdated(),
94  m_colNormsDirect(),
95  m_isInitialized(false),
96  m_usePrescribedThreshold(false) {}
97 
105  : m_qr(rows, cols),
106  m_hCoeffs((std::min)(rows,cols)),
107  m_colsPermutation(PermIndexType(cols)),
108  m_colsTranspositions(cols),
109  m_temp(cols),
110  m_colNormsUpdated(cols),
111  m_colNormsDirect(cols),
112  m_isInitialized(false),
113  m_usePrescribedThreshold(false) {}
114 
127  template<typename InputType>
129  : m_qr(matrix.rows(), matrix.cols()),
130  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
131  m_colsPermutation(PermIndexType(matrix.cols())),
132  m_colsTranspositions(matrix.cols()),
133  m_temp(matrix.cols()),
134  m_colNormsUpdated(matrix.cols()),
135  m_colNormsDirect(matrix.cols()),
136  m_isInitialized(false),
137  m_usePrescribedThreshold(false)
138  {
139  compute(matrix.derived());
140  }
141 
148  template<typename InputType>
150  : m_qr(matrix.derived()),
151  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
152  m_colsPermutation(PermIndexType(matrix.cols())),
153  m_colsTranspositions(matrix.cols()),
154  m_temp(matrix.cols()),
155  m_colNormsUpdated(matrix.cols()),
156  m_colNormsDirect(matrix.cols()),
157  m_isInitialized(false),
158  m_usePrescribedThreshold(false)
159  {
160  computeInPlace();
161  }
162 
163  #ifdef EIGEN_PARSED_BY_DOXYGEN
178  template<typename Rhs>
180  solve(const MatrixBase<Rhs>& b) const;
181  #endif
182 
184  HouseholderSequenceType matrixQ() const
185  {
186  return householderQ();
187  }
188 
191  const MatrixType& matrixQR() const
192  {
193  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
194  return m_qr;
195  }
196 
206  const MatrixType& matrixR() const
207  {
208  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
209  return m_qr;
210  }
211 
212  template<typename InputType>
213  ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
214 
217  {
218  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
219  return m_colsPermutation;
220  }
221 
235  typename MatrixType::RealScalar absDeterminant() const;
236 
249  typename MatrixType::RealScalar logAbsDeterminant() const;
250 
257  inline Index rank() const
258  {
259  using std::abs;
260  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
261  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
262  Index result = 0;
263  for(Index i = 0; i < m_nonzero_pivots; ++i)
264  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
265  return result;
266  }
267 
274  inline Index dimensionOfKernel() const
275  {
276  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
277  return cols() - rank();
278  }
279 
287  inline bool isInjective() const
288  {
289  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
290  return rank() == cols();
291  }
292 
300  inline bool isSurjective() const
301  {
302  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
303  return rank() == rows();
304  }
305 
312  inline bool isInvertible() const
313  {
314  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
315  return isInjective() && isSurjective();
316  }
317 
324  {
325  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
326  return Inverse<ColPivHouseholderQR>(*this);
327  }
328 
329  inline Index rows() const { return m_qr.rows(); }
330  inline Index cols() const { return m_qr.cols(); }
331 
336  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
337 
356  {
357  m_usePrescribedThreshold = true;
358  m_prescribedThreshold = threshold;
359  return *this;
360  }
361 
371  {
372  m_usePrescribedThreshold = false;
373  return *this;
374  }
375 
380  RealScalar threshold() const
381  {
382  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
383  return m_usePrescribedThreshold ? m_prescribedThreshold
384  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
385  // and turns out to be identical to Higham's formula used already in LDLt.
386  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
387  }
388 
396  inline Index nonzeroPivots() const
397  {
398  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
399  return m_nonzero_pivots;
400  }
401 
405  RealScalar maxPivot() const { return m_maxpivot; }
406 
414  {
415  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
416  return Success;
417  }
418 
419  #ifndef EIGEN_PARSED_BY_DOXYGEN
420  template<typename RhsType, typename DstType>
421  void _solve_impl(const RhsType &rhs, DstType &dst) const;
422 
423  template<bool Conjugate, typename RhsType, typename DstType>
424  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
425  #endif
426 
427  protected:
428 
429  friend class CompleteOrthogonalDecomposition<MatrixType>;
430 
431  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
432 
433  void computeInPlace();
434 
435  MatrixType m_qr;
436  HCoeffsType m_hCoeffs;
437  PermutationType m_colsPermutation;
438  IntRowVectorType m_colsTranspositions;
439  RowVectorType m_temp;
440  RealRowVectorType m_colNormsUpdated;
441  RealRowVectorType m_colNormsDirect;
442  bool m_isInitialized, m_usePrescribedThreshold;
443  RealScalar m_prescribedThreshold, m_maxpivot;
444  Index m_nonzero_pivots;
445  Index m_det_pq;
446 };
447 
448 template<typename MatrixType>
449 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
450 {
451  using std::abs;
452  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
453  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
454  return abs(m_qr.diagonal().prod());
455 }
456 
457 template<typename MatrixType>
458 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
459 {
460  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
461  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
462  return m_qr.diagonal().cwiseAbs().array().log().sum();
463 }
464 
471 template<typename MatrixType>
472 template<typename InputType>
474 {
475  m_qr = matrix.derived();
476  computeInPlace();
477  return *this;
478 }
479 
480 template<typename MatrixType>
482 {
483  // the column permutation is stored as int indices, so just to be sure:
484  eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
485 
486  using std::abs;
487 
488  Index rows = m_qr.rows();
489  Index cols = m_qr.cols();
490  Index size = m_qr.diagonalSize();
491 
492  m_hCoeffs.resize(size);
493 
494  m_temp.resize(cols);
495 
496  m_colsTranspositions.resize(m_qr.cols());
497  Index number_of_transpositions = 0;
498 
499  m_colNormsUpdated.resize(cols);
500  m_colNormsDirect.resize(cols);
501  for (Index k = 0; k < cols; ++k) {
502  // colNormsDirect(k) caches the most recent directly computed norm of
503  // column k.
504  m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
505  m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
506  }
507 
508  RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
509  RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
510 
511  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
512  m_maxpivot = RealScalar(0);
513 
514  for(Index k = 0; k < size; ++k)
515  {
516  // first, we look up in our table m_colNormsUpdated which column has the biggest norm
517  Index biggest_col_index;
518  RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
519  biggest_col_index += k;
520 
521  // Track the number of meaningful pivots but do not stop the decomposition to make
522  // sure that the initial matrix is properly reproduced. See bug 941.
523  if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
524  m_nonzero_pivots = k;
525 
526  // apply the transposition to the columns
527  m_colsTranspositions.coeffRef(k) = biggest_col_index;
528  if(k != biggest_col_index) {
529  m_qr.col(k).swap(m_qr.col(biggest_col_index));
530  std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
531  std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
532  ++number_of_transpositions;
533  }
534 
535  // generate the householder vector, store it below the diagonal
536  RealScalar beta;
537  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
538 
539  // apply the householder transformation to the diagonal coefficient
540  m_qr.coeffRef(k,k) = beta;
541 
542  // remember the maximum absolute value of diagonal coefficients
543  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
544 
545  // apply the householder transformation
546  m_qr.bottomRightCorner(rows-k, cols-k-1)
547  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
548 
549  // update our table of norms of the columns
550  for (Index j = k + 1; j < cols; ++j) {
551  // The following implements the stable norm downgrade step discussed in
552  // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
553  // and used in LAPACK routines xGEQPF and xGEQP3.
554  // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
555  if (!numext::is_exactly_zero(m_colNormsUpdated.coeffRef(j))) {
556  RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
557  temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
558  temp = temp < RealScalar(0) ? RealScalar(0) : temp;
559  RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
560  m_colNormsDirect.coeffRef(j));
561  if (temp2 <= norm_downdate_threshold) {
562  // The updated norm has become too inaccurate so re-compute the column
563  // norm directly.
564  m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
565  m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
566  } else {
567  m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
568  }
569  }
570  }
571  }
572 
573  m_colsPermutation.setIdentity(PermIndexType(cols));
574  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
575  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
576 
577  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
578  m_isInitialized = true;
579 }
580 
581 #ifndef EIGEN_PARSED_BY_DOXYGEN
582 template<typename MatrixType_>
583 template<typename RhsType, typename DstType>
584 void ColPivHouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
585 {
586  const Index nonzero_pivots = nonzeroPivots();
587 
588  if(nonzero_pivots == 0)
589  {
590  dst.setZero();
591  return;
592  }
593 
594  typename RhsType::PlainObject c(rhs);
595 
596  c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
597 
598  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
599  .template triangularView<Upper>()
600  .solveInPlace(c.topRows(nonzero_pivots));
601 
602  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
603  for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
604 }
605 
606 template<typename MatrixType_>
607 template<bool Conjugate, typename RhsType, typename DstType>
608 void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
609 {
610  const Index nonzero_pivots = nonzeroPivots();
611 
612  if(nonzero_pivots == 0)
613  {
614  dst.setZero();
615  return;
616  }
617 
618  typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs);
619 
620  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
621  .template triangularView<Upper>()
622  .transpose().template conjugateIf<Conjugate>()
623  .solveInPlace(c.topRows(nonzero_pivots));
624 
625  dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
626  dst.bottomRows(rows()-nonzero_pivots).setZero();
627 
628  dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).template conjugateIf<!Conjugate>() );
629 }
630 #endif
631 
632 namespace internal {
633 
634 template<typename DstXprType, typename MatrixType>
635 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
636 {
637  typedef ColPivHouseholderQR<MatrixType> QrType;
638  typedef Inverse<QrType> SrcXprType;
639  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
640  {
641  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
642  }
643 };
644 
645 } // end namespace internal
646 
650 template<typename MatrixType>
651 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
652  ::householderQ() const
653 {
654  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
655  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
656 }
657 
662 template<typename Derived>
665 {
666  return ColPivHouseholderQR<PlainObject>(eval());
667 }
668 
669 } // end namespace Eigen
670 
671 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ColPivHouseholderQR.h:55
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:336
ColPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:149
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:396
bool isInvertible() const
Definition: ColPivHouseholderQR.h:312
ColPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:128
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:652
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:191
ComputationInfo info() const
Reports whether the QR factorization was successful.
Definition: ColPivHouseholderQR.h:413
Index rank() const
Definition: ColPivHouseholderQR.h:257
const Inverse< ColPivHouseholderQR > inverse() const
Definition: ColPivHouseholderQR.h:323
ColPivHouseholderQR()
Default Constructor.
Definition: ColPivHouseholderQR.h:87
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: ColPivHouseholderQR.h:104
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:380
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:216
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:405
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:206
bool isSurjective() const
Definition: ColPivHouseholderQR.h:300
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:449
const Solve< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:370
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:274
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:355
bool isInjective() const
Definition: ColPivHouseholderQR.h:287
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:458
Complete orthogonal decomposition (COD) of a matrix.
Definition: CompleteOrthogonalDecomposition.h:54
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
ColPivHouseholderQR< MatrixType_ > & derived()
Definition: EigenBase.h:48
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:107
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231