FullPivHouseholderQR.h
Go to the documentation of this file.
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 FullPivHouseholderQRMatrixQReturnType;
19 
20 template<typename MatrixType>
22 {
23  typedef typename MatrixType::PlainObject ReturnType;
24 };
25 
26 }
27 
49 template<typename _MatrixType> class FullPivHouseholderQR
50 {
51  public:
52 
53  typedef _MatrixType MatrixType;
54  enum {
55  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
57  Options = MatrixType::Options,
58  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
59  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
60  };
61  typedef typename MatrixType::Scalar Scalar;
62  typedef typename MatrixType::RealScalar RealScalar;
63  typedef typename MatrixType::Index Index;
71 
78  : m_qr(),
79  m_hCoeffs(),
80  m_rows_transpositions(),
81  m_cols_transpositions(),
82  m_cols_permutation(),
83  m_temp(),
84  m_isInitialized(false),
85  m_usePrescribedThreshold(false) {}
86 
93  FullPivHouseholderQR(Index rows, Index cols)
94  : m_qr(rows, cols),
95  m_hCoeffs((std::min)(rows,cols)),
96  m_rows_transpositions(rows),
97  m_cols_transpositions(cols),
98  m_cols_permutation(cols),
99  m_temp((std::min)(rows,cols)),
100  m_isInitialized(false),
101  m_usePrescribedThreshold(false) {}
102 
115  FullPivHouseholderQR(const MatrixType& matrix)
116  : m_qr(matrix.rows(), matrix.cols()),
117  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
118  m_rows_transpositions(matrix.rows()),
119  m_cols_transpositions(matrix.cols()),
120  m_cols_permutation(matrix.cols()),
121  m_temp((std::min)(matrix.rows(), matrix.cols())),
122  m_isInitialized(false),
123  m_usePrescribedThreshold(false)
124  {
125  compute(matrix);
126  }
127 
145  template<typename Rhs>
147  solve(const MatrixBase<Rhs>& b) const
148  {
149  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
150  return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
151  }
152 
155  MatrixQReturnType matrixQ(void) const;
156 
159  const MatrixType& matrixQR() const
160  {
161  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
162  return m_qr;
163  }
164 
165  FullPivHouseholderQR& compute(const MatrixType& matrix);
166 
168  const PermutationType& colsPermutation() const
169  {
170  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
171  return m_cols_permutation;
172  }
173 
175  const IntColVectorType& rowsTranspositions() const
176  {
177  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
178  return m_rows_transpositions;
179  }
180 
194  typename MatrixType::RealScalar absDeterminant() const;
195 
208  typename MatrixType::RealScalar logAbsDeterminant() const;
209 
216  inline Index rank() const
217  {
218  using std::abs;
219  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
220  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
221  Index result = 0;
222  for(Index i = 0; i < m_nonzero_pivots; ++i)
223  result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
224  return result;
225  }
226 
233  inline Index dimensionOfKernel() const
234  {
235  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
236  return cols() - rank();
237  }
238 
246  inline bool isInjective() const
247  {
248  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
249  return rank() == cols();
250  }
251 
259  inline bool isSurjective() const
260  {
261  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
262  return rank() == rows();
263  }
264 
271  inline bool isInvertible() const
272  {
273  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
274  return isInjective() && isSurjective();
275  }
276  inline const
283  inverse() const
284  {
285  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
287  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
288  }
289 
290  inline Index rows() const { return m_qr.rows(); }
291  inline Index cols() const { return m_qr.cols(); }
292 
297  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
298 
316  FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
317  {
318  m_usePrescribedThreshold = true;
319  m_prescribedThreshold = threshold;
320  return *this;
321  }
322 
332  {
333  m_usePrescribedThreshold = false;
334  return *this;
335  }
336 
341  RealScalar threshold() const
342  {
343  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
344  return m_usePrescribedThreshold ? m_prescribedThreshold
345  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
346  // and turns out to be identical to Higham's formula used already in LDLt.
347  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
348  }
349 
357  inline Index nonzeroPivots() const
358  {
359  eigen_assert(m_isInitialized && "LU is not initialized.");
360  return m_nonzero_pivots;
361  }
362 
366  RealScalar maxPivot() const { return m_maxpivot; }
367 
368  protected:
369  MatrixType m_qr;
370  HCoeffsType m_hCoeffs;
371  IntColVectorType m_rows_transpositions;
372  IntRowVectorType m_cols_transpositions;
373  PermutationType m_cols_permutation;
374  RowVectorType m_temp;
375  bool m_isInitialized, m_usePrescribedThreshold;
376  RealScalar m_prescribedThreshold, m_maxpivot;
378  RealScalar m_precision;
379  Index m_det_pq;
380 };
381 
382 template<typename MatrixType>
383 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
384 {
385  using std::abs;
386  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
387  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
388  return abs(m_qr.diagonal().prod());
389 }
390 
391 template<typename MatrixType>
392 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
393 {
394  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
395  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
396  return m_qr.diagonal().cwiseAbs().array().log().sum();
397 }
398 
405 template<typename MatrixType>
407 {
408  using std::abs;
409  Index rows = matrix.rows();
410  Index cols = matrix.cols();
411  Index size = (std::min)(rows,cols);
412 
413  m_qr = matrix;
414  m_hCoeffs.resize(size);
415 
416  m_temp.resize(cols);
417 
418  m_precision = NumTraits<Scalar>::epsilon() * size;
419 
420  m_rows_transpositions.resize(matrix.rows());
421  m_cols_transpositions.resize(matrix.cols());
422  Index number_of_transpositions = 0;
423 
424  RealScalar biggest(0);
425 
426  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
427  m_maxpivot = RealScalar(0);
428 
429  for (Index k = 0; k < size; ++k)
430  {
431  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
432  RealScalar biggest_in_corner;
433 
434  biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
435  .cwiseAbs()
436  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
437  row_of_biggest_in_corner += k;
438  col_of_biggest_in_corner += k;
439  if(k==0) biggest = biggest_in_corner;
440 
441  // if the corner is negligible, then we have less than full rank, and we can finish early
442  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
443  {
444  m_nonzero_pivots = k;
445  for(Index i = k; i < size; i++)
446  {
447  m_rows_transpositions.coeffRef(i) = i;
448  m_cols_transpositions.coeffRef(i) = i;
449  m_hCoeffs.coeffRef(i) = Scalar(0);
450  }
451  break;
452  }
453 
454  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
455  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
456  if(k != row_of_biggest_in_corner) {
457  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
458  ++number_of_transpositions;
459  }
460  if(k != col_of_biggest_in_corner) {
461  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
462  ++number_of_transpositions;
463  }
464 
465  RealScalar beta;
466  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
467  m_qr.coeffRef(k,k) = beta;
468 
469  // remember the maximum absolute value of diagonal coefficients
470  if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
471 
472  m_qr.bottomRightCorner(rows-k, cols-k-1)
473  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
474  }
475 
476  m_cols_permutation.setIdentity(cols);
477  for(Index k = 0; k < size; ++k)
478  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
479 
480  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
481  m_isInitialized = true;
482 
483  return *this;
484 }
485 
486 namespace internal {
487 
488 template<typename _MatrixType, typename Rhs>
489 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
490  : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
491 {
493 
494  template<typename Dest> void evalTo(Dest& dst) const
495  {
496  const Index rows = dec().rows(), cols = dec().cols();
497  eigen_assert(rhs().rows() == rows);
498 
499  // FIXME introduce nonzeroPivots() and use it here. and more generally,
500  // make the same improvements in this dec as in FullPivLU.
501  if(dec().rank()==0)
502  {
503  dst.setZero();
504  return;
505  }
506 
507  typename Rhs::PlainObject c(rhs());
508 
510  for (Index k = 0; k < dec().rank(); ++k)
511  {
512  Index remainingSize = rows-k;
513  c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
514  c.bottomRightCorner(remainingSize, rhs().cols())
515  .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
516  dec().hCoeffs().coeff(k), &temp.coeffRef(0));
517  }
518 
519  if(!dec().isSurjective())
520  {
521  // is c is in the image of R ?
522  RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
523  RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
524  // FIXME brain dead
525  const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
526  // this internal:: prefix is needed by at least gcc 3.4 and ICC
527  if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
528  return;
529  }
530  dec().matrixQR()
531  .topLeftCorner(dec().rank(), dec().rank())
532  .template triangularView<Upper>()
533  .solveInPlace(c.topRows(dec().rank()));
534 
535  for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
536  for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
537  }
538 };
539 
546 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
547  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
548 {
549 public:
550  typedef typename MatrixType::Index Index;
553  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
554  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
555 
557  const HCoeffsType& hCoeffs,
558  const IntColVectorType& rowsTranspositions)
559  : m_qr(qr),
560  m_hCoeffs(hCoeffs),
561  m_rowsTranspositions(rowsTranspositions)
562  {}
563 
564  template <typename ResultType>
565  void evalTo(ResultType& result) const
566  {
567  const Index rows = m_qr.rows();
568  WorkVectorType workspace(rows);
569  evalTo(result, workspace);
570  }
571 
572  template <typename ResultType>
573  void evalTo(ResultType& result, WorkVectorType& workspace) const
574  {
575  using numext::conj;
576  // compute the product H'_0 H'_1 ... H'_n-1,
577  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
578  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
579  const Index rows = m_qr.rows();
580  const Index cols = m_qr.cols();
581  const Index size = (std::min)(rows, cols);
582  workspace.resize(rows);
583  result.setIdentity(rows, rows);
584  for (Index k = size-1; k >= 0; k--)
585  {
586  result.block(k, k, rows-k, rows-k)
587  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
588  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
589  }
590  }
591 
592  Index rows() const { return m_qr.rows(); }
593  Index cols() const { return m_qr.rows(); }
594 
595 protected:
596  typename MatrixType::Nested m_qr;
597  typename HCoeffsType::Nested m_hCoeffs;
598  typename IntColVectorType::Nested m_rowsTranspositions;
599 };
600 
601 } // end namespace internal
602 
603 template<typename MatrixType>
605 {
606  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
607  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
608 }
609 
614 template<typename Derived>
617 {
618  return FullPivHouseholderQR<PlainObject>(eval());
619 }
620 
621 } // end namespace Eigen
622 
623 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
FullPivHouseholderQR & setThreshold(Default_t)
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Householder rank-revealing QR decomposition of a matrix with full pivoting.
internal::plain_col_type< MatrixType, Index >::type IntColVectorType
const internal::solve_retval< FullPivHouseholderQR, typename MatrixType::IdentityReturnType > inverse() const
const IntColVectorType & rowsTranspositions() const
MatrixQReturnType matrixQ(void) const
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
const internal::solve_retval< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
MatrixType::RealScalar RealScalar
bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, typename NumTraits< Scalar >::Real precision=NumTraits< Scalar >::dummy_precision())
void evalTo(ResultType &result, WorkVectorType &workspace) const
FullPivHouseholderQRMatrixQReturnType(const MatrixType &qr, const HCoeffsType &hCoeffs, const IntColVectorType &rowsTranspositions)
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
internal::plain_diag_type< MatrixType >::type HCoeffsType
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
Matrix< typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, MatrixType::MaxRowsAtCompileTime > WorkVectorType
const HCoeffsType & hCoeffs() const
FullPivHouseholderQR()
Default Constructor.
internal::plain_row_type< MatrixType >::type RowVectorType
const PermutationType & colsPermutation() const
internal::plain_col_type< MatrixType, Index >::type IntColVectorType
FullPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
SegmentReturnType tail(Index vecSize)
Definition: BlockMethods.h:810
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
internal::FullPivHouseholderQRMatrixQReturnType< MatrixType > MatrixQReturnType
Expression type for return value of FullPivHouseholderQR::matrixQ()
internal::plain_col_type< MatrixType >::type ColVectorType
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
PermutationMatrix< ColsAtCompileTime, MaxColsAtCompileTime > PermutationType
void rhs(const real_t *x, real_t *f)
MatrixType::RealScalar logAbsDeterminant() const
MatrixType::RealScalar absDeterminant() const
const MatrixType & matrixQR() const
Matrix< Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime > IntRowVectorType
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
ColXpr col(Index i)
Definition: BlockMethods.h:708
#define eigen_assert(x)
FullPivHouseholderQR & compute(const MatrixType &matrix)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
internal::plain_diag_type< MatrixType >::type HCoeffsType
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:35