Eigen  3.3.7
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 namespace Eigen {
15 
16 namespace internal {
17 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
18  : traits<_MatrixType>
19 {
20  enum { Flags = 0 };
21 };
22 
23 } // end namespace internal
24 
48 template<typename _MatrixType> class ColPivHouseholderQR
49 {
50  public:
51 
52  typedef _MatrixType MatrixType;
53  enum {
54  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
55  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
56  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
58  };
59  typedef typename MatrixType::Scalar Scalar;
60  typedef typename MatrixType::RealScalar RealScalar;
61  // FIXME should be int
62  typedef typename MatrixType::StorageIndex StorageIndex;
63  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
64  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
65  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
66  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
67  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
68  typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
69  typedef typename MatrixType::PlainObject PlainObject;
70 
71  private:
72 
73  typedef typename PermutationType::StorageIndex PermIndexType;
74 
75  public:
76 
84  : m_qr(),
85  m_hCoeffs(),
86  m_colsPermutation(),
87  m_colsTranspositions(),
88  m_temp(),
89  m_colNormsUpdated(),
90  m_colNormsDirect(),
91  m_isInitialized(false),
92  m_usePrescribedThreshold(false) {}
93 
101  : m_qr(rows, cols),
102  m_hCoeffs((std::min)(rows,cols)),
103  m_colsPermutation(PermIndexType(cols)),
104  m_colsTranspositions(cols),
105  m_temp(cols),
106  m_colNormsUpdated(cols),
107  m_colNormsDirect(cols),
108  m_isInitialized(false),
109  m_usePrescribedThreshold(false) {}
110 
123  template<typename InputType>
125  : m_qr(matrix.rows(), matrix.cols()),
126  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
127  m_colsPermutation(PermIndexType(matrix.cols())),
128  m_colsTranspositions(matrix.cols()),
129  m_temp(matrix.cols()),
130  m_colNormsUpdated(matrix.cols()),
131  m_colNormsDirect(matrix.cols()),
132  m_isInitialized(false),
133  m_usePrescribedThreshold(false)
134  {
135  compute(matrix.derived());
136  }
137 
144  template<typename InputType>
146  : m_qr(matrix.derived()),
147  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
148  m_colsPermutation(PermIndexType(matrix.cols())),
149  m_colsTranspositions(matrix.cols()),
150  m_temp(matrix.cols()),
151  m_colNormsUpdated(matrix.cols()),
152  m_colNormsDirect(matrix.cols()),
153  m_isInitialized(false),
154  m_usePrescribedThreshold(false)
155  {
156  computeInPlace();
157  }
158 
173  template<typename Rhs>
175  solve(const MatrixBase<Rhs>& b) const
176  {
177  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
178  return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
179  }
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>
211  ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
212 
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 
378  RealScalar threshold() const
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  EIGEN_DEVICE_FUNC
420  void _solve_impl(const RhsType &rhs, DstType &dst) const;
421  #endif
422 
423  protected:
424 
425  friend class CompleteOrthogonalDecomposition<MatrixType>;
426 
427  static void check_template_parameters()
428  {
429  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
430  }
431 
432  void computeInPlace();
433 
434  MatrixType m_qr;
435  HCoeffsType m_hCoeffs;
436  PermutationType m_colsPermutation;
437  IntRowVectorType m_colsTranspositions;
438  RowVectorType m_temp;
439  RealRowVectorType m_colNormsUpdated;
440  RealRowVectorType m_colNormsDirect;
441  bool m_isInitialized, m_usePrescribedThreshold;
442  RealScalar m_prescribedThreshold, m_maxpivot;
443  Index m_nonzero_pivots;
444  Index m_det_pq;
445 };
446 
447 template<typename MatrixType>
448 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
449 {
450  using std::abs;
451  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
452  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
453  return abs(m_qr.diagonal().prod());
454 }
455 
456 template<typename MatrixType>
457 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
458 {
459  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
460  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
461  return m_qr.diagonal().cwiseAbs().array().log().sum();
462 }
463 
470 template<typename MatrixType>
471 template<typename InputType>
473 {
474  m_qr = matrix.derived();
475  computeInPlace();
476  return *this;
477 }
478 
479 template<typename MatrixType>
481 {
482  check_template_parameters();
483 
484  // the column permutation is stored as int indices, so just to be sure:
485  eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
486 
487  using std::abs;
488 
489  Index rows = m_qr.rows();
490  Index cols = m_qr.cols();
491  Index size = m_qr.diagonalSize();
492 
493  m_hCoeffs.resize(size);
494 
495  m_temp.resize(cols);
496 
497  m_colsTranspositions.resize(m_qr.cols());
498  Index number_of_transpositions = 0;
499 
500  m_colNormsUpdated.resize(cols);
501  m_colNormsDirect.resize(cols);
502  for (Index k = 0; k < cols; ++k) {
503  // colNormsDirect(k) caches the most recent directly computed norm of
504  // column k.
505  m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
506  m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
507  }
508 
509  RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
510  RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
511 
512  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
513  m_maxpivot = RealScalar(0);
514 
515  for(Index k = 0; k < size; ++k)
516  {
517  // first, we look up in our table m_colNormsUpdated which column has the biggest norm
518  Index biggest_col_index;
519  RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
520  biggest_col_index += k;
521 
522  // Track the number of meaningful pivots but do not stop the decomposition to make
523  // sure that the initial matrix is properly reproduced. See bug 941.
524  if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
525  m_nonzero_pivots = k;
526 
527  // apply the transposition to the columns
528  m_colsTranspositions.coeffRef(k) = biggest_col_index;
529  if(k != biggest_col_index) {
530  m_qr.col(k).swap(m_qr.col(biggest_col_index));
531  std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
532  std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
533  ++number_of_transpositions;
534  }
535 
536  // generate the householder vector, store it below the diagonal
537  RealScalar beta;
538  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
539 
540  // apply the householder transformation to the diagonal coefficient
541  m_qr.coeffRef(k,k) = beta;
542 
543  // remember the maximum absolute value of diagonal coefficients
544  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
545 
546  // apply the householder transformation
547  m_qr.bottomRightCorner(rows-k, cols-k-1)
548  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
549 
550  // update our table of norms of the columns
551  for (Index j = k + 1; j < cols; ++j) {
552  // The following implements the stable norm downgrade step discussed in
553  // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
554  // and used in LAPACK routines xGEQPF and xGEQP3.
555  // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
556  if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
557  RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
558  temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
559  temp = temp < RealScalar(0) ? RealScalar(0) : temp;
560  RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
561  m_colNormsDirect.coeffRef(j));
562  if (temp2 <= norm_downdate_threshold) {
563  // The updated norm has become too inaccurate so re-compute the column
564  // norm directly.
565  m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
566  m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
567  } else {
568  m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
569  }
570  }
571  }
572  }
573 
574  m_colsPermutation.setIdentity(PermIndexType(cols));
575  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
576  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
577 
578  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
579  m_isInitialized = true;
580 }
581 
582 #ifndef EIGEN_PARSED_BY_DOXYGEN
583 template<typename _MatrixType>
584 template<typename RhsType, typename DstType>
585 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
586 {
587  eigen_assert(rhs.rows() == rows());
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  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
600  c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
601  .setLength(nonzero_pivots)
602  .transpose()
603  );
604 
605  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
606  .template triangularView<Upper>()
607  .solveInPlace(c.topRows(nonzero_pivots));
608 
609  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
610  for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
611 }
612 #endif
613 
614 namespace internal {
615 
616 template<typename DstXprType, typename MatrixType>
617 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
618 {
619  typedef ColPivHouseholderQR<MatrixType> QrType;
620  typedef Inverse<QrType> SrcXprType;
621  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
622  {
623  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
624  }
625 };
626 
627 } // end namespace internal
628 
632 template<typename MatrixType>
633 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
634  ::householderQ() const
635 {
636  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
637  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
638 }
639 
644 template<typename Derived>
647 {
648  return ColPivHouseholderQR<PlainObject>(eval());
649 }
650 
651 } // end namespace Eigen
652 
653 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
Index rank() const
Definition: ColPivHouseholderQR.h:255
bool isSurjective() const
Definition: ColPivHouseholderQR.h:298
ColPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:145
ColPivHouseholderQR()
Default Constructor.
Definition: ColPivHouseholderQR.h:83
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: ColPivHouseholderQR.h:100
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:334
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Definition: Half.h:150
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
ComputationInfo info() const
Reports whether the QR factorization was succesful.
Definition: ColPivHouseholderQR.h:411
Derived & derived()
Definition: EigenBase.h:45
Complete orthogonal decomposition (COD) of a matrix.
Definition: ForwardDeclarations.h:253
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:451
Definition: EigenBase.h:29
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: ForwardDeclarations.h:258
ColPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:124
Expression of the inverse of another expression.
Definition: Inverse.h:43
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:634
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:403
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:368
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ForwardDeclarations.h:251
Derived & derived()
Definition: EigenBase.h:45
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:204
bool isInjective() const
Definition: ColPivHouseholderQR.h:285
const Solve< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: ColPivHouseholderQR.h:175
bool isInvertible() const
Definition: ColPivHouseholderQR.h:310
Definition: Constants.h:432
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:448
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:353
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:378
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:457
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:272
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:394
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:189
ComputationInfo
Definition: Constants.h:430
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:214
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Inverse< ColPivHouseholderQR > inverse() const
Definition: ColPivHouseholderQR.h:321