Eigen  3.2.92
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  template<typename InputType>
123  : m_qr(matrix.rows(), matrix.cols()),
124  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
125  m_colsPermutation(PermIndexType(matrix.cols())),
126  m_colsTranspositions(matrix.cols()),
127  m_temp(matrix.cols()),
128  m_colSqNorms(matrix.cols()),
129  m_isInitialized(false),
130  m_usePrescribedThreshold(false)
131  {
132  compute(matrix.derived());
133  }
134 
152  template<typename Rhs>
154  solve(const MatrixBase<Rhs>& b) const
155  {
156  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
157  return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
158  }
159 
160  HouseholderSequenceType householderQ() const;
161  HouseholderSequenceType matrixQ() const
162  {
163  return householderQ();
164  }
165 
168  const MatrixType& matrixQR() const
169  {
170  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
171  return m_qr;
172  }
173 
183  const MatrixType& matrixR() const
184  {
185  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
186  return m_qr;
187  }
188 
189  template<typename InputType>
190  ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
191 
193  const PermutationType& colsPermutation() const
194  {
195  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
196  return m_colsPermutation;
197  }
198 
212  typename MatrixType::RealScalar absDeterminant() const;
213 
226  typename MatrixType::RealScalar logAbsDeterminant() const;
227 
234  inline Index rank() const
235  {
236  using std::abs;
237  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
238  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
239  Index result = 0;
240  for(Index i = 0; i < m_nonzero_pivots; ++i)
241  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
242  return result;
243  }
244 
251  inline Index dimensionOfKernel() const
252  {
253  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
254  return cols() - rank();
255  }
256 
264  inline bool isInjective() const
265  {
266  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
267  return rank() == cols();
268  }
269 
277  inline bool isSurjective() const
278  {
279  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
280  return rank() == rows();
281  }
282 
289  inline bool isInvertible() const
290  {
291  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
292  return isInjective() && isSurjective();
293  }
294 
301  {
302  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
303  return Inverse<ColPivHouseholderQR>(*this);
304  }
305 
306  inline Index rows() const { return m_qr.rows(); }
307  inline Index cols() const { return m_qr.cols(); }
308 
313  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
314 
333  {
334  m_usePrescribedThreshold = true;
335  m_prescribedThreshold = threshold;
336  return *this;
337  }
338 
348  {
349  m_usePrescribedThreshold = false;
350  return *this;
351  }
352 
357  RealScalar threshold() const
358  {
359  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
360  return m_usePrescribedThreshold ? m_prescribedThreshold
361  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
362  // and turns out to be identical to Higham's formula used already in LDLt.
363  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
364  }
365 
373  inline Index nonzeroPivots() const
374  {
375  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
376  return m_nonzero_pivots;
377  }
378 
382  RealScalar maxPivot() const { return m_maxpivot; }
383 
391  {
392  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
393  return Success;
394  }
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 
404  static void check_template_parameters()
405  {
406  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
407  }
408 
409  void computeInPlace();
410 
411  MatrixType m_qr;
412  HCoeffsType m_hCoeffs;
413  PermutationType m_colsPermutation;
414  IntRowVectorType m_colsTranspositions;
415  RowVectorType m_temp;
416  RealRowVectorType m_colSqNorms;
417  bool m_isInitialized, m_usePrescribedThreshold;
418  RealScalar m_prescribedThreshold, m_maxpivot;
419  Index m_nonzero_pivots;
420  Index m_det_pq;
421 };
422 
423 template<typename MatrixType>
424 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
425 {
426  using std::abs;
427  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
428  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
429  return abs(m_qr.diagonal().prod());
430 }
431 
432 template<typename MatrixType>
433 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
434 {
435  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
436  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
437  return m_qr.diagonal().cwiseAbs().array().log().sum();
438 }
439 
446 template<typename MatrixType>
447 template<typename InputType>
449 {
450  check_template_parameters();
451 
452  // the column permutation is stored as int indices, so just to be sure:
453  eigen_assert(matrix.cols()<=NumTraits<int>::highest());
454 
455  m_qr = matrix;
456 
457  computeInPlace();
458 
459  return *this;
460 }
461 
462 template<typename MatrixType>
464 {
465  using std::abs;
466  Index rows = m_qr.rows();
467  Index cols = m_qr.cols();
468  Index size = m_qr.diagonalSize();
469 
470  m_hCoeffs.resize(size);
471 
472  m_temp.resize(cols);
473 
474  m_colsTranspositions.resize(m_qr.cols());
475  Index number_of_transpositions = 0;
476 
477  m_colSqNorms.resize(cols);
478  for(Index k = 0; k < cols; ++k)
479  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
480 
481  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
482 
483  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
484  m_maxpivot = RealScalar(0);
485 
486  for(Index k = 0; k < size; ++k)
487  {
488  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
489  Index biggest_col_index;
490  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
491  biggest_col_index += k;
492 
493  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
494  // the actual squared norm of the selected column.
495  // Note that not doing so does result in solve() sometimes returning inf/nan values
496  // when running the unit test with 1000 repetitions.
497  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
498 
499  // we store that back into our table: it can't hurt to correct our table.
500  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
501 
502  // Track the number of meaningful pivots but do not stop the decomposition to make
503  // sure that the initial matrix is properly reproduced. See bug 941.
504  if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
505  m_nonzero_pivots = k;
506 
507  // apply the transposition to the columns
508  m_colsTranspositions.coeffRef(k) = biggest_col_index;
509  if(k != biggest_col_index) {
510  m_qr.col(k).swap(m_qr.col(biggest_col_index));
511  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
512  ++number_of_transpositions;
513  }
514 
515  // generate the householder vector, store it below the diagonal
516  RealScalar beta;
517  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
518 
519  // apply the householder transformation to the diagonal coefficient
520  m_qr.coeffRef(k,k) = beta;
521 
522  // remember the maximum absolute value of diagonal coefficients
523  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
524 
525  // apply the householder transformation
526  m_qr.bottomRightCorner(rows-k, cols-k-1)
527  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
528 
529  // update our table of squared norms of the columns
530  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
531  }
532 
533  m_colsPermutation.setIdentity(PermIndexType(cols));
534  for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
535  m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
536 
537  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
538  m_isInitialized = true;
539 }
540 
541 #ifndef EIGEN_PARSED_BY_DOXYGEN
542 template<typename _MatrixType>
543 template<typename RhsType, typename DstType>
544 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
545 {
546  eigen_assert(rhs.rows() == rows());
547 
548  const Index nonzero_pivots = nonzeroPivots();
549 
550  if(nonzero_pivots == 0)
551  {
552  dst.setZero();
553  return;
554  }
555 
556  typename RhsType::PlainObject c(rhs);
557 
558  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
559  c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
560  .setLength(nonzero_pivots)
561  .transpose()
562  );
563 
564  m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
565  .template triangularView<Upper>()
566  .solveInPlace(c.topRows(nonzero_pivots));
567 
568  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
569  for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
570 }
571 #endif
572 
573 namespace internal {
574 
575 template<typename DstXprType, typename MatrixType, typename Scalar>
576 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
577 {
578  typedef ColPivHouseholderQR<MatrixType> QrType;
579  typedef Inverse<QrType> SrcXprType;
580  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
581  {
582  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
583  }
584 };
585 
586 } // end namespace internal
587 
591 template<typename MatrixType>
592 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
594 {
595  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
596  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
597 }
598 
599 #ifndef __CUDACC__
600 
604 template<typename Derived>
607 {
608  return ColPivHouseholderQR<PlainObject>(eval());
609 }
610 #endif // __CUDACC__
611 
612 } // end namespace Eigen
613 
614 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:193
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:373
const Inverse< ColPivHouseholderQR > inverse() const
Definition: ColPivHouseholderQR.h:300
bool isInjective() const
Definition: ColPivHouseholderQR.h:264
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:154
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:433
Index rank() const
Definition: ColPivHouseholderQR.h:234
bool isInvertible() const
Definition: ColPivHouseholderQR.h:289
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
Derived & derived()
Definition: EigenBase.h:44
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:390
Definition: EigenBase.h:28
ColPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:122
Expression of the inverse of another expression.
Definition: Inverse.h:43
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:347
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:357
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:183
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ForwardDeclarations.h:253
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:168
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:593
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:313
const ColPivHouseholderQR< PlainObject > colPivHouseholderQr() const
Definition: ColPivHouseholderQR.h:606
Definition: Constants.h:432
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:332
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:382
Definition: Eigen_Colamd.h:54
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:424
bool isSurjective() const
Definition: ColPivHouseholderQR.h:277
Index cols() const
Definition: EigenBase.h:61
Pseudo expression representing a solving operation.
Definition: Solve.h:62
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:251