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 FullPivHouseholderQRMatrixQReturnType;
19 
20 template<typename MatrixType>
21 struct traits<FullPivHouseholderQRMatrixQReturnType<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;
64  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
65  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
66  typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
67  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
68  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
69  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
70  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
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 
103  FullPivHouseholderQR(const MatrixType& matrix)
104  : m_qr(matrix.rows(), matrix.cols()),
105  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
106  m_rows_transpositions(matrix.rows()),
107  m_cols_transpositions(matrix.cols()),
108  m_cols_permutation(matrix.cols()),
109  m_temp((std::min)(matrix.rows(), matrix.cols())),
110  m_isInitialized(false),
111  m_usePrescribedThreshold(false)
112  {
113  compute(matrix);
114  }
115 
133  template<typename Rhs>
134  inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
135  solve(const MatrixBase<Rhs>& b) const
136  {
137  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
138  return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
139  }
140 
143  MatrixQReturnType matrixQ(void) const;
144 
147  const MatrixType& matrixQR() const
148  {
149  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
150  return m_qr;
151  }
152 
153  FullPivHouseholderQR& compute(const MatrixType& matrix);
154 
155  const PermutationType& colsPermutation() const
156  {
157  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
158  return m_cols_permutation;
159  }
160 
161  const IntColVectorType& rowsTranspositions() const
162  {
163  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
164  return m_rows_transpositions;
165  }
166 
180  typename MatrixType::RealScalar absDeterminant() const;
181 
194  typename MatrixType::RealScalar logAbsDeterminant() const;
195 
202  inline Index rank() const
203  {
204  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
205  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
206  Index result = 0;
207  for(Index i = 0; i < m_nonzero_pivots; ++i)
208  result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
209  return result;
210  }
211 
218  inline Index dimensionOfKernel() const
219  {
220  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
221  return cols() - rank();
222  }
223 
231  inline bool isInjective() const
232  {
233  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
234  return rank() == cols();
235  }
236 
244  inline bool isSurjective() const
245  {
246  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
247  return rank() == rows();
248  }
249 
256  inline bool isInvertible() const
257  {
258  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
259  return isInjective() && isSurjective();
260  }
261  inline const
267  internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
268  inverse() const
269  {
270  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
271  return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
272  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
273  }
274 
275  inline Index rows() const { return m_qr.rows(); }
276  inline Index cols() const { return m_qr.cols(); }
277  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
278 
297  {
298  m_usePrescribedThreshold = true;
299  m_prescribedThreshold = threshold;
300  return *this;
301  }
302 
312  {
313  m_usePrescribedThreshold = false;
314  return *this;
315  }
316 
321  RealScalar threshold() const
322  {
323  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
324  return m_usePrescribedThreshold ? m_prescribedThreshold
325  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
326  // and turns out to be identical to Higham's formula used already in LDLt.
327  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
328  }
329 
337  inline Index nonzeroPivots() const
338  {
339  eigen_assert(m_isInitialized && "LU is not initialized.");
340  return m_nonzero_pivots;
341  }
342 
346  RealScalar maxPivot() const { return m_maxpivot; }
347 
348  protected:
349  MatrixType m_qr;
350  HCoeffsType m_hCoeffs;
351  IntColVectorType m_rows_transpositions;
352  IntRowVectorType m_cols_transpositions;
353  PermutationType m_cols_permutation;
354  RowVectorType m_temp;
355  bool m_isInitialized, m_usePrescribedThreshold;
356  RealScalar m_prescribedThreshold, m_maxpivot;
357  Index m_nonzero_pivots;
358  RealScalar m_precision;
359  Index m_det_pq;
360 };
361 
362 template<typename MatrixType>
363 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
364 {
365  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
366  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
367  return internal::abs(m_qr.diagonal().prod());
368 }
369 
370 template<typename MatrixType>
371 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
372 {
373  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
374  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
375  return m_qr.diagonal().cwiseAbs().array().log().sum();
376 }
377 
378 template<typename MatrixType>
380 {
381  Index rows = matrix.rows();
382  Index cols = matrix.cols();
383  Index size = (std::min)(rows,cols);
384 
385  m_qr = matrix;
386  m_hCoeffs.resize(size);
387 
388  m_temp.resize(cols);
389 
390  m_precision = NumTraits<Scalar>::epsilon() * size;
391 
392  m_rows_transpositions.resize(matrix.rows());
393  m_cols_transpositions.resize(matrix.cols());
394  Index number_of_transpositions = 0;
395 
396  RealScalar biggest(0);
397 
398  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
399  m_maxpivot = RealScalar(0);
400 
401  for (Index k = 0; k < size; ++k)
402  {
403  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
404  RealScalar biggest_in_corner;
405 
406  biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
407  .cwiseAbs()
408  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
409  row_of_biggest_in_corner += k;
410  col_of_biggest_in_corner += k;
411  if(k==0) biggest = biggest_in_corner;
412 
413  // if the corner is negligible, then we have less than full rank, and we can finish early
414  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
415  {
416  m_nonzero_pivots = k;
417  for(Index i = k; i < size; i++)
418  {
419  m_rows_transpositions.coeffRef(i) = i;
420  m_cols_transpositions.coeffRef(i) = i;
421  m_hCoeffs.coeffRef(i) = Scalar(0);
422  }
423  break;
424  }
425 
426  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
427  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
428  if(k != row_of_biggest_in_corner) {
429  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
430  ++number_of_transpositions;
431  }
432  if(k != col_of_biggest_in_corner) {
433  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
434  ++number_of_transpositions;
435  }
436 
437  RealScalar beta;
438  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
439  m_qr.coeffRef(k,k) = beta;
440 
441  // remember the maximum absolute value of diagonal coefficients
442  if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
443 
444  m_qr.bottomRightCorner(rows-k, cols-k-1)
445  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
446  }
447 
448  m_cols_permutation.setIdentity(cols);
449  for(Index k = 0; k < size; ++k)
450  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
451 
452  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
453  m_isInitialized = true;
454 
455  return *this;
456 }
457 
458 namespace internal {
459 
460 template<typename _MatrixType, typename Rhs>
461 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
462  : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
463 {
464  EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
465 
466  template<typename Dest> void evalTo(Dest& dst) const
467  {
468  const Index rows = dec().rows(), cols = dec().cols();
469  eigen_assert(rhs().rows() == rows);
470 
471  // FIXME introduce nonzeroPivots() and use it here. and more generally,
472  // make the same improvements in this dec as in FullPivLU.
473  if(dec().rank()==0)
474  {
475  dst.setZero();
476  return;
477  }
478 
479  typename Rhs::PlainObject c(rhs());
480 
481  Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
482  for (Index k = 0; k < dec().rank(); ++k)
483  {
484  Index remainingSize = rows-k;
485  c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
486  c.bottomRightCorner(remainingSize, rhs().cols())
487  .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
488  dec().hCoeffs().coeff(k), &temp.coeffRef(0));
489  }
490 
491  if(!dec().isSurjective())
492  {
493  // is c is in the image of R ?
494  RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
495  RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
496  // FIXME brain dead
497  const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
498  // this internal:: prefix is needed by at least gcc 3.4 and ICC
499  if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
500  return;
501  }
502  dec().matrixQR()
503  .topLeftCorner(dec().rank(), dec().rank())
504  .template triangularView<Upper>()
505  .solveInPlace(c.topRows(dec().rank()));
506 
507  for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
508  for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
509  }
510 };
511 
518 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
519  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
520 {
521 public:
522  typedef typename MatrixType::Index Index;
523  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
524  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
525  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
526  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
527 
528  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
529  const HCoeffsType& hCoeffs,
530  const IntColVectorType& rowsTranspositions)
531  : m_qr(qr),
532  m_hCoeffs(hCoeffs),
533  m_rowsTranspositions(rowsTranspositions)
534  {}
535 
536  template <typename ResultType>
537  void evalTo(ResultType& result) const
538  {
539  const Index rows = m_qr.rows();
540  WorkVectorType workspace(rows);
541  evalTo(result, workspace);
542  }
543 
544  template <typename ResultType>
545  void evalTo(ResultType& result, WorkVectorType& workspace) const
546  {
547  // compute the product H'_0 H'_1 ... H'_n-1,
548  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
549  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
550  const Index rows = m_qr.rows();
551  const Index cols = m_qr.cols();
552  const Index size = (std::min)(rows, cols);
553  workspace.resize(rows);
554  result.setIdentity(rows, rows);
555  for (Index k = size-1; k >= 0; k--)
556  {
557  result.block(k, k, rows-k, rows-k)
558  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
559  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
560  }
561  }
562 
563  Index rows() const { return m_qr.rows(); }
564  Index cols() const { return m_qr.rows(); }
565 
566 protected:
567  typename MatrixType::Nested m_qr;
568  typename HCoeffsType::Nested m_hCoeffs;
569  typename IntColVectorType::Nested m_rowsTranspositions;
570 };
571 
572 } // end namespace internal
573 
574 template<typename MatrixType>
575 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
576 {
577  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
578  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
579 }
580 
585 template<typename Derived>
588 {
589  return FullPivHouseholderQR<PlainObject>(eval());
590 }
591 
592 } // end namespace Eigen
593 
594 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H