Eigen  3.2.91
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 
46 template<typename _MatrixType> class ColPivHouseholderQR
47 {
48  public:
49 
50  typedef _MatrixType MatrixType;
51  enum {
52  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
53  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
54  Options = MatrixType::Options,
55  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
56  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
57  };
58  typedef typename MatrixType::Scalar Scalar;
59  typedef typename MatrixType::RealScalar RealScalar;
60  // FIXME should be int
61  typedef typename MatrixType::StorageIndex StorageIndex;
62  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
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_colSqNorms(),
90  m_isInitialized(false),
91  m_usePrescribedThreshold(false) {}
92 
99  ColPivHouseholderQR(Index rows, Index cols)
100  : m_qr(rows, cols),
101  m_hCoeffs((std::min)(rows,cols)),
102  m_colsPermutation(PermIndexType(cols)),
103  m_colsTranspositions(cols),
104  m_temp(cols),
105  m_colSqNorms(cols),
106  m_isInitialized(false),
107  m_usePrescribedThreshold(false) {}
108 
121  explicit ColPivHouseholderQR(const MatrixType& matrix)
122  : m_qr(matrix.rows(), matrix.cols()),
123  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
124  m_colsPermutation(PermIndexType(matrix.cols())),
125  m_colsTranspositions(matrix.cols()),
126  m_temp(matrix.cols()),
127  m_colSqNorms(matrix.cols()),
128  m_isInitialized(false),
129  m_usePrescribedThreshold(false)
130  {
131  compute(matrix);
132  }
133 
151  template<typename Rhs>
153  solve(const MatrixBase<Rhs>& b) const
154  {
155  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
156  return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
157  }
158 
159  HouseholderSequenceType householderQ() const;
160  HouseholderSequenceType matrixQ() const
161  {
162  return householderQ();
163  }
164 
167  const MatrixType& matrixQR() const
168  {
169  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
170  return m_qr;
171  }
172 
182  const MatrixType& matrixR() const
183  {
184  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
185  return m_qr;
186  }
187 
188  ColPivHouseholderQR& compute(const MatrixType& matrix);
189 
191  const PermutationType& colsPermutation() const
192  {
193  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
194  return m_colsPermutation;
195  }
196 
210  typename MatrixType::RealScalar absDeterminant() const;
211 
224  typename MatrixType::RealScalar logAbsDeterminant() const;
225 
232  inline Index rank() const
233  {
234  using std::abs;
235  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
236  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
237  Index result = 0;
238  for(Index i = 0; i < m_nonzero_pivots; ++i)
239  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
240  return result;
241  }
242 
249  inline Index dimensionOfKernel() const
250  {
251  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
252  return cols() - rank();
253  }
254 
262  inline bool isInjective() const
263  {
264  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
265  return rank() == cols();
266  }
267 
275  inline bool isSurjective() const
276  {
277  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
278  return rank() == rows();
279  }
280 
287  inline bool isInvertible() const
288  {
289  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
290  return isInjective() && isSurjective();
291  }
292 
299  {
300  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
301  return Inverse<ColPivHouseholderQR>(*this);
302  }
303 
304  inline Index rows() const { return m_qr.rows(); }
305  inline Index cols() const { return m_qr.cols(); }
306 
311  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
312 
331  {
332  m_usePrescribedThreshold = true;
333  m_prescribedThreshold = threshold;
334  return *this;
335  }
336 
346  {
347  m_usePrescribedThreshold = false;
348  return *this;
349  }
350 
355  RealScalar threshold() const
356  {
357  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
358  return m_usePrescribedThreshold ? m_prescribedThreshold
359  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
360  // and turns out to be identical to Higham's formula used already in LDLt.
361  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
362  }
363 
371  inline Index nonzeroPivots() const
372  {
373  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
374  return m_nonzero_pivots;
375  }
376 
380  RealScalar maxPivot() const { return m_maxpivot; }
381 
389  {
390  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
391  return Success;
392  }
393 
394  #ifndef EIGEN_PARSED_BY_DOXYGEN
395  template<typename RhsType, typename DstType>
396  EIGEN_DEVICE_FUNC
397  void _solve_impl(const RhsType &rhs, DstType &dst) const;
398  #endif
399 
400  protected:
401 
402  static void check_template_parameters()
403  {
404  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
405  }
406 
407  MatrixType m_qr;
408  HCoeffsType m_hCoeffs;
409  PermutationType m_colsPermutation;
410  IntRowVectorType m_colsTranspositions;
411  RowVectorType m_temp;
412  RealRowVectorType m_colSqNorms;
413  bool m_isInitialized, m_usePrescribedThreshold;
414  RealScalar m_prescribedThreshold, m_maxpivot;
415  Index m_nonzero_pivots;
416  Index m_det_pq;
417 };
418 
419 template<typename MatrixType>
420 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
421 {
422  using std::abs;
423  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
424  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
425  return abs(m_qr.diagonal().prod());
426 }
427 
428 template<typename MatrixType>
429 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
430 {
431  eigen_assert(m_isInitialized && "ColPivHouseholderQR 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 m_qr.diagonal().cwiseAbs().array().log().sum();
434 }
435 
442 template<typename MatrixType>
444 {
445  check_template_parameters();
446 
447  using std::abs;
448  Index rows = matrix.rows();
449  Index cols = matrix.cols();
450  Index size = matrix.diagonalSize();
451 
452  // the column permutation is stored as int indices, so just to be sure:
453  eigen_assert(cols<=NumTraits<int>::highest());
454 
455  m_qr = matrix;
456  m_hCoeffs.resize(size);
457 
458  m_temp.resize(cols);
459 
460  m_colsTranspositions.resize(matrix.cols());
461  Index number_of_transpositions = 0;
462 
463  m_colSqNorms.resize(cols);
464  for(Index k = 0; k < cols; ++k)
465  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
466 
467  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
468 
469  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
470  m_maxpivot = RealScalar(0);
471 
472  for(Index k = 0; k < size; ++k)
473  {
474  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
475  Index biggest_col_index;
476  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
477  biggest_col_index += k;
478 
479  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
480  // the actual squared norm of the selected column.
481  // Note that not doing so does result in solve() sometimes returning inf/nan values
482  // when running the unit test with 1000 repetitions.
483  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
484 
485  // we store that back into our table: it can't hurt to correct our table.
486  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
487 
488  // Track the number of meaningful pivots but do not stop the decomposition to make
489  // sure that the initial matrix is properly reproduced. See bug 941.
490  if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
491  m_nonzero_pivots = k;
492 
493  // apply the transposition to the columns
494  m_colsTranspositions.coeffRef(k) = biggest_col_index;
495  if(k != biggest_col_index) {
496  m_qr.col(k).swap(m_qr.col(biggest_col_index));
497  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
498  ++number_of_transpositions;
499  }
500 
501  // generate the householder vector, store it below the diagonal
502  RealScalar beta;
503  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
504 
505  // apply the householder transformation to the diagonal coefficient
506  m_qr.coeffRef(k,k) = beta;
507 
508  // remember the maximum absolute value of diagonal coefficients
509  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
510 
511  // apply the householder transformation
512  m_qr.bottomRightCorner(rows-k, cols-k-1)
513  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
514 
515  // update our table of squared norms of the columns
516  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
517  }
518 
519  m_colsPermutation.setIdentity(PermIndexType(cols));
520  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
521  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
522 
523  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
524  m_isInitialized = true;
525 
526  return *this;
527 }
528 
529 #ifndef EIGEN_PARSED_BY_DOXYGEN
530 template<typename _MatrixType>
531 template<typename RhsType, typename DstType>
532 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
533 {
534  eigen_assert(rhs.rows() == rows());
535 
536  const Index nonzero_pivots = nonzeroPivots();
537 
538  if(nonzero_pivots == 0)
539  {
540  dst.setZero();
541  return;
542  }
543 
544  typename RhsType::PlainObject c(rhs);
545 
546  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
547  c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
548  .setLength(nonzero_pivots)
549  .transpose()
550  );
551 
552  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
553  .template triangularView<Upper>()
554  .solveInPlace(c.topRows(nonzero_pivots));
555 
556  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
557  for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
558 }
559 #endif
560 
561 namespace internal {
562 
563 template<typename DstXprType, typename MatrixType, typename Scalar>
564 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
565 {
566  typedef ColPivHouseholderQR<MatrixType> QrType;
567  typedef Inverse<QrType> SrcXprType;
568  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
569  {
570  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
571  }
572 };
573 
574 } // end namespace internal
575 
579 template<typename MatrixType>
580 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
582 {
583  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
584  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
585 }
586 
587 #ifndef __CUDACC__
588 
592 template<typename Derived>
595 {
596  return ColPivHouseholderQR<PlainObject>(eval());
597 }
598 #endif // __CUDACC__
599 
600 } // end namespace Eigen
601 
602 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:191
ColPivHouseholderQR & compute(const MatrixType &matrix)
Definition: ColPivHouseholderQR.h:443
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:371
const Inverse< ColPivHouseholderQR > inverse() const
Definition: ColPivHouseholderQR.h:298
bool isInjective() const
Definition: ColPivHouseholderQR.h:262
ColPivHouseholderQR()
Default Constructor.
Definition: ColPivHouseholderQR.h:83
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: ColPivHouseholderQR.h:99
const Solve< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: ColPivHouseholderQR.h:153
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:429
Index rank() const
Definition: ColPivHouseholderQR.h:232
bool isInvertible() const
Definition: ColPivHouseholderQR.h:287
Definition: LDLT.h:16
Definition: StdDeque.h:58
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:452
ComputationInfo info() const
Reports whether the QR factorization was succesful.
Definition: ColPivHouseholderQR.h:388
Expression of the inverse of another expression.
Definition: Inverse.h:45
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:345
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:355
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:182
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ForwardDeclarations.h:252
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:167
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:581
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:311
const ColPivHouseholderQR< PlainObject > colPivHouseholderQr() const
Definition: ColPivHouseholderQR.h:594
Definition: Constants.h:424
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:330
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:380
Definition: Eigen_Colamd.h:54
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:420
bool isSurjective() const
Definition: ColPivHouseholderQR.h:275
ColPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:121
Pseudo expression representing a solving operation.
Definition: Solve.h:63
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:249