Eigen  3.2.91
FullPivHouseholderQR.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_FULLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
18 template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
19  : traits<_MatrixType>
20 {
21  enum { Flags = 0 };
22 };
23 
24 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
25 
26 template<typename MatrixType>
27 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
28 {
29  typedef typename MatrixType::PlainObject ReturnType;
30 };
31 
32 } // end namespace internal
33 
55 template<typename _MatrixType> class FullPivHouseholderQR
56 {
57  public:
58 
59  typedef _MatrixType MatrixType;
60  enum {
61  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
62  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
63  Options = MatrixType::Options,
64  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
66  };
67  typedef typename MatrixType::Scalar Scalar;
68  typedef typename MatrixType::RealScalar RealScalar;
69  // FIXME should be int
70  typedef typename MatrixType::StorageIndex StorageIndex;
71  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
72  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
73  typedef Matrix<StorageIndex, 1,
74  EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
75  EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
76  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
77  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
78  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
79  typedef typename MatrixType::PlainObject PlainObject;
80 
87  : m_qr(),
88  m_hCoeffs(),
89  m_rows_transpositions(),
90  m_cols_transpositions(),
91  m_cols_permutation(),
92  m_temp(),
93  m_isInitialized(false),
94  m_usePrescribedThreshold(false) {}
95 
102  FullPivHouseholderQR(Index rows, Index cols)
103  : m_qr(rows, cols),
104  m_hCoeffs((std::min)(rows,cols)),
105  m_rows_transpositions((std::min)(rows,cols)),
106  m_cols_transpositions((std::min)(rows,cols)),
107  m_cols_permutation(cols),
108  m_temp(cols),
109  m_isInitialized(false),
110  m_usePrescribedThreshold(false) {}
111 
124  explicit FullPivHouseholderQR(const MatrixType& matrix)
125  : m_qr(matrix.rows(), matrix.cols()),
126  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
127  m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
128  m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
129  m_cols_permutation(matrix.cols()),
130  m_temp(matrix.cols()),
131  m_isInitialized(false),
132  m_usePrescribedThreshold(false)
133  {
134  compute(matrix);
135  }
136 
155  template<typename Rhs>
157  solve(const MatrixBase<Rhs>& b) const
158  {
159  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
160  return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
161  }
162 
165  MatrixQReturnType matrixQ(void) const;
166 
169  const MatrixType& matrixQR() const
170  {
171  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
172  return m_qr;
173  }
174 
175  FullPivHouseholderQR& compute(const MatrixType& matrix);
176 
178  const PermutationType& colsPermutation() const
179  {
180  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
181  return m_cols_permutation;
182  }
183 
185  const IntDiagSizeVectorType& rowsTranspositions() const
186  {
187  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
188  return m_rows_transpositions;
189  }
190 
204  typename MatrixType::RealScalar absDeterminant() const;
205 
218  typename MatrixType::RealScalar logAbsDeterminant() const;
219 
226  inline Index rank() const
227  {
228  using std::abs;
229  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
230  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
231  Index result = 0;
232  for(Index i = 0; i < m_nonzero_pivots; ++i)
233  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
234  return result;
235  }
236 
243  inline Index dimensionOfKernel() const
244  {
245  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
246  return cols() - rank();
247  }
248 
256  inline bool isInjective() const
257  {
258  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
259  return rank() == cols();
260  }
261 
269  inline bool isSurjective() const
270  {
271  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
272  return rank() == rows();
273  }
274 
281  inline bool isInvertible() const
282  {
283  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
284  return isInjective() && isSurjective();
285  }
286 
293  {
294  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
295  return Inverse<FullPivHouseholderQR>(*this);
296  }
297 
298  inline Index rows() const { return m_qr.rows(); }
299  inline Index cols() const { return m_qr.cols(); }
300 
305  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
306 
325  {
326  m_usePrescribedThreshold = true;
327  m_prescribedThreshold = threshold;
328  return *this;
329  }
330 
340  {
341  m_usePrescribedThreshold = false;
342  return *this;
343  }
344 
349  RealScalar threshold() const
350  {
351  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
352  return m_usePrescribedThreshold ? m_prescribedThreshold
353  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
354  // and turns out to be identical to Higham's formula used already in LDLt.
355  : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
356  }
357 
365  inline Index nonzeroPivots() const
366  {
367  eigen_assert(m_isInitialized && "LU is not initialized.");
368  return m_nonzero_pivots;
369  }
370 
374  RealScalar maxPivot() const { return m_maxpivot; }
375 
376  #ifndef EIGEN_PARSED_BY_DOXYGEN
377  template<typename RhsType, typename DstType>
378  EIGEN_DEVICE_FUNC
379  void _solve_impl(const RhsType &rhs, DstType &dst) const;
380  #endif
381 
382  protected:
383 
384  static void check_template_parameters()
385  {
386  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
387  }
388 
389  MatrixType m_qr;
390  HCoeffsType m_hCoeffs;
391  IntDiagSizeVectorType m_rows_transpositions;
392  IntDiagSizeVectorType m_cols_transpositions;
393  PermutationType m_cols_permutation;
394  RowVectorType m_temp;
395  bool m_isInitialized, m_usePrescribedThreshold;
396  RealScalar m_prescribedThreshold, m_maxpivot;
397  Index m_nonzero_pivots;
398  RealScalar m_precision;
399  Index m_det_pq;
400 };
401 
402 template<typename MatrixType>
403 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
404 {
405  using std::abs;
406  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
407  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
408  return abs(m_qr.diagonal().prod());
409 }
410 
411 template<typename MatrixType>
412 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
413 {
414  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
415  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
416  return m_qr.diagonal().cwiseAbs().array().log().sum();
417 }
418 
425 template<typename MatrixType>
427 {
428  check_template_parameters();
429 
430  using std::abs;
431  Index rows = matrix.rows();
432  Index cols = matrix.cols();
433  Index size = (std::min)(rows,cols);
434 
435  m_qr = matrix;
436  m_hCoeffs.resize(size);
437 
438  m_temp.resize(cols);
439 
440  m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
441 
442  m_rows_transpositions.resize(size);
443  m_cols_transpositions.resize(size);
444  Index number_of_transpositions = 0;
445 
446  RealScalar biggest(0);
447 
448  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
449  m_maxpivot = RealScalar(0);
450 
451  for (Index k = 0; k < size; ++k)
452  {
453  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
454  typedef internal::scalar_score_coeff_op<Scalar> Scoring;
455  typedef typename Scoring::result_type Score;
456 
457  Score score = m_qr.bottomRightCorner(rows-k, cols-k)
458  .unaryExpr(Scoring())
459  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
460  row_of_biggest_in_corner += k;
461  col_of_biggest_in_corner += k;
462  RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
463  if(k==0) biggest = biggest_in_corner;
464 
465  // if the corner is negligible, then we have less than full rank, and we can finish early
466  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
467  {
468  m_nonzero_pivots = k;
469  for(Index i = k; i < size; i++)
470  {
471  m_rows_transpositions.coeffRef(i) = i;
472  m_cols_transpositions.coeffRef(i) = i;
473  m_hCoeffs.coeffRef(i) = Scalar(0);
474  }
475  break;
476  }
477 
478  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
479  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
480  if(k != row_of_biggest_in_corner) {
481  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
482  ++number_of_transpositions;
483  }
484  if(k != col_of_biggest_in_corner) {
485  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
486  ++number_of_transpositions;
487  }
488 
489  RealScalar beta;
490  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
491  m_qr.coeffRef(k,k) = beta;
492 
493  // remember the maximum absolute value of diagonal coefficients
494  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
495 
496  m_qr.bottomRightCorner(rows-k, cols-k-1)
497  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
498  }
499 
500  m_cols_permutation.setIdentity(cols);
501  for(Index k = 0; k < size; ++k)
502  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
503 
504  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
505  m_isInitialized = true;
506 
507  return *this;
508 }
509 
510 #ifndef EIGEN_PARSED_BY_DOXYGEN
511 template<typename _MatrixType>
512 template<typename RhsType, typename DstType>
513 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
514 {
515  eigen_assert(rhs.rows() == rows());
516  const Index l_rank = rank();
517 
518  // FIXME introduce nonzeroPivots() and use it here. and more generally,
519  // make the same improvements in this dec as in FullPivLU.
520  if(l_rank==0)
521  {
522  dst.setZero();
523  return;
524  }
525 
526  typename RhsType::PlainObject c(rhs);
527 
528  Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
529  for (Index k = 0; k < l_rank; ++k)
530  {
531  Index remainingSize = rows()-k;
532  c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
533  c.bottomRightCorner(remainingSize, rhs.cols())
534  .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
535  m_hCoeffs.coeff(k), &temp.coeffRef(0));
536  }
537 
538  m_qr.topLeftCorner(l_rank, l_rank)
539  .template triangularView<Upper>()
540  .solveInPlace(c.topRows(l_rank));
541 
542  for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
543  for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
544 }
545 #endif
546 
547 namespace internal {
548 
549 template<typename DstXprType, typename MatrixType, typename Scalar>
550 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
551 {
552  typedef FullPivHouseholderQR<MatrixType> QrType;
553  typedef Inverse<QrType> SrcXprType;
554  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
555  {
556  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
557  }
558 };
559 
566 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
567  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
568 {
569 public:
570  typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
571  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
572  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
573  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
574 
575  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
576  const HCoeffsType& hCoeffs,
577  const IntDiagSizeVectorType& rowsTranspositions)
578  : m_qr(qr),
579  m_hCoeffs(hCoeffs),
580  m_rowsTranspositions(rowsTranspositions)
581  {}
582 
583  template <typename ResultType>
584  void evalTo(ResultType& result) const
585  {
586  const Index rows = m_qr.rows();
587  WorkVectorType workspace(rows);
588  evalTo(result, workspace);
589  }
590 
591  template <typename ResultType>
592  void evalTo(ResultType& result, WorkVectorType& workspace) const
593  {
594  using numext::conj;
595  // compute the product H'_0 H'_1 ... H'_n-1,
596  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
597  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
598  const Index rows = m_qr.rows();
599  const Index cols = m_qr.cols();
600  const Index size = (std::min)(rows, cols);
601  workspace.resize(rows);
602  result.setIdentity(rows, rows);
603  for (Index k = size-1; k >= 0; k--)
604  {
605  result.block(k, k, rows-k, rows-k)
606  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
607  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
608  }
609  }
610 
611  Index rows() const { return m_qr.rows(); }
612  Index cols() const { return m_qr.rows(); }
613 
614 protected:
615  typename MatrixType::Nested m_qr;
616  typename HCoeffsType::Nested m_hCoeffs;
617  typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
618 };
619 
620 // template<typename MatrixType>
621 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
622 // : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
623 // {};
624 
625 } // end namespace internal
626 
627 template<typename MatrixType>
628 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
629 {
630  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
631  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
632 }
633 
634 #ifndef __CUDACC__
635 
639 template<typename Derived>
642 {
643  return FullPivHouseholderQR<PlainObject>(eval());
644 }
645 #endif // __CUDACC__
646 
647 } // end namespace Eigen
648 
649 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:243
bool isInvertible() const
Definition: FullPivHouseholderQR.h:281
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:339
Definition: Constants.h:314
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: ForwardDeclarations.h:253
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:628
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:374
Definition: LDLT.h:16
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:185
bool isInjective() const
Definition: FullPivHouseholderQR.h:256
Definition: StdDeque.h:58
bool isSurjective() const
Definition: FullPivHouseholderQR.h:269
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:102
Expression of the inverse of another expression.
Definition: Inverse.h:45
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivHouseholderQR.h:157
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:305
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:86
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:178
FullPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:124
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:324
Definition: Eigen_Colamd.h:54
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:412
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:349
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:403
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:169
Index rank() const
Definition: FullPivHouseholderQR.h:226
Pseudo expression representing a solving operation.
Definition: Solve.h:63
FullPivHouseholderQR & compute(const MatrixType &matrix)
Definition: FullPivHouseholderQR.h:426
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:641
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:292
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:365