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 
37 template<typename _MatrixType> class ColPivHouseholderQR
38 {
39  public:
40 
41  typedef _MatrixType MatrixType;
42  enum {
43  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
44  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
45  Options = MatrixType::Options,
46  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
47  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
48  };
49  typedef typename MatrixType::Scalar Scalar;
50  typedef typename MatrixType::RealScalar RealScalar;
51  typedef typename MatrixType::Index Index;
52  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
53  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
54  typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
55  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
56  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
57  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
58  typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
59 
60  private:
61 
62  typedef typename PermutationType::Index PermIndexType;
63 
64  public:
65 
73  : m_qr(),
74  m_hCoeffs(),
75  m_colsPermutation(),
76  m_colsTranspositions(),
77  m_temp(),
78  m_colSqNorms(),
79  m_isInitialized(false) {}
80 
87  ColPivHouseholderQR(Index rows, Index cols)
88  : m_qr(rows, cols),
89  m_hCoeffs((std::min)(rows,cols)),
90  m_colsPermutation(PermIndexType(cols)),
91  m_colsTranspositions(cols),
92  m_temp(cols),
93  m_colSqNorms(cols),
94  m_isInitialized(false),
95  m_usePrescribedThreshold(false) {}
96 
97  ColPivHouseholderQR(const MatrixType& matrix)
98  : m_qr(matrix.rows(), matrix.cols()),
99  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
100  m_colsPermutation(PermIndexType(matrix.cols())),
101  m_colsTranspositions(matrix.cols()),
102  m_temp(matrix.cols()),
103  m_colSqNorms(matrix.cols()),
104  m_isInitialized(false),
105  m_usePrescribedThreshold(false)
106  {
107  compute(matrix);
108  }
109 
127  template<typename Rhs>
128  inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
129  solve(const MatrixBase<Rhs>& b) const
130  {
131  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
132  return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
133  }
134 
135  HouseholderSequenceType householderQ(void) const;
136 
139  const MatrixType& matrixQR() const
140  {
141  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
142  return m_qr;
143  }
144 
145  ColPivHouseholderQR& compute(const MatrixType& matrix);
146 
147  const PermutationType& colsPermutation() const
148  {
149  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
150  return m_colsPermutation;
151  }
152 
166  typename MatrixType::RealScalar absDeterminant() const;
167 
180  typename MatrixType::RealScalar logAbsDeterminant() const;
181 
188  inline Index rank() const
189  {
190  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
191  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
192  Index result = 0;
193  for(Index i = 0; i < m_nonzero_pivots; ++i)
194  result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
195  return result;
196  }
197 
204  inline Index dimensionOfKernel() const
205  {
206  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
207  return cols() - rank();
208  }
209 
217  inline bool isInjective() const
218  {
219  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
220  return rank() == cols();
221  }
222 
230  inline bool isSurjective() const
231  {
232  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
233  return rank() == rows();
234  }
235 
242  inline bool isInvertible() const
243  {
244  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
245  return isInjective() && isSurjective();
246  }
247 
253  inline const
254  internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
255  inverse() const
256  {
257  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
258  return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
259  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
260  }
261 
262  inline Index rows() const { return m_qr.rows(); }
263  inline Index cols() const { return m_qr.cols(); }
264  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
265 
284  {
285  m_usePrescribedThreshold = true;
286  m_prescribedThreshold = threshold;
287  return *this;
288  }
289 
299  {
300  m_usePrescribedThreshold = false;
301  return *this;
302  }
303 
308  RealScalar threshold() const
309  {
310  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
311  return m_usePrescribedThreshold ? m_prescribedThreshold
312  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
313  // and turns out to be identical to Higham's formula used already in LDLt.
314  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
315  }
316 
324  inline Index nonzeroPivots() const
325  {
326  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
327  return m_nonzero_pivots;
328  }
329 
333  RealScalar maxPivot() const { return m_maxpivot; }
334 
335  protected:
336  MatrixType m_qr;
337  HCoeffsType m_hCoeffs;
338  PermutationType m_colsPermutation;
339  IntRowVectorType m_colsTranspositions;
340  RowVectorType m_temp;
341  RealRowVectorType m_colSqNorms;
342  bool m_isInitialized, m_usePrescribedThreshold;
343  RealScalar m_prescribedThreshold, m_maxpivot;
344  Index m_nonzero_pivots;
345  Index m_det_pq;
346 };
347 
348 template<typename MatrixType>
349 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
350 {
351  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
352  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
353  return internal::abs(m_qr.diagonal().prod());
354 }
355 
356 template<typename MatrixType>
357 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
358 {
359  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
360  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
361  return m_qr.diagonal().cwiseAbs().array().log().sum();
362 }
363 
364 template<typename MatrixType>
366 {
367  Index rows = matrix.rows();
368  Index cols = matrix.cols();
369  Index size = matrix.diagonalSize();
370 
371  m_qr = matrix;
372  m_hCoeffs.resize(size);
373 
374  m_temp.resize(cols);
375 
376  m_colsTranspositions.resize(matrix.cols());
377  Index number_of_transpositions = 0;
378 
379  m_colSqNorms.resize(cols);
380  for(Index k = 0; k < cols; ++k)
381  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
382 
383  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
384 
385  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
386  m_maxpivot = RealScalar(0);
387 
388  for(Index k = 0; k < size; ++k)
389  {
390  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
391  Index biggest_col_index;
392  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
393  biggest_col_index += k;
394 
395  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
396  // the actual squared norm of the selected column.
397  // Note that not doing so does result in solve() sometimes returning inf/nan values
398  // when running the unit test with 1000 repetitions.
399  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
400 
401  // we store that back into our table: it can't hurt to correct our table.
402  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
403 
404  // if the current biggest column is smaller than epsilon times the initial biggest column,
405  // terminate to avoid generating nan/inf values.
406  // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
407  // repetitions of the unit test, with the result of solve() filled with large values of the order
408  // of 1/(size*epsilon).
409  if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
410  {
411  m_nonzero_pivots = k;
412  m_hCoeffs.tail(size-k).setZero();
413  m_qr.bottomRightCorner(rows-k,cols-k)
414  .template triangularView<StrictlyLower>()
415  .setZero();
416  break;
417  }
418 
419  // apply the transposition to the columns
420  m_colsTranspositions.coeffRef(k) = biggest_col_index;
421  if(k != biggest_col_index) {
422  m_qr.col(k).swap(m_qr.col(biggest_col_index));
423  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
424  ++number_of_transpositions;
425  }
426 
427  // generate the householder vector, store it below the diagonal
428  RealScalar beta;
429  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
430 
431  // apply the householder transformation to the diagonal coefficient
432  m_qr.coeffRef(k,k) = beta;
433 
434  // remember the maximum absolute value of diagonal coefficients
435  if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
436 
437  // apply the householder transformation
438  m_qr.bottomRightCorner(rows-k, cols-k-1)
439  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
440 
441  // update our table of squared norms of the columns
442  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
443  }
444 
445  m_colsPermutation.setIdentity(PermIndexType(cols));
446  for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)
447  m_colsPermutation.applyTranspositionOnTheRight(PermIndexType(k), PermIndexType(m_colsTranspositions.coeff(k)));
448 
449  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
450  m_isInitialized = true;
451 
452  return *this;
453 }
454 
455 namespace internal {
456 
457 template<typename _MatrixType, typename Rhs>
458 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
459  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
460 {
461  EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
462 
463  template<typename Dest> void evalTo(Dest& dst) const
464  {
465  eigen_assert(rhs().rows() == dec().rows());
466 
467  const int cols = dec().cols(),
468  nonzero_pivots = dec().nonzeroPivots();
469 
470  if(nonzero_pivots == 0)
471  {
472  dst.setZero();
473  return;
474  }
475 
476  typename Rhs::PlainObject c(rhs());
477 
478  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
479  c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
480  .setLength(dec().nonzeroPivots())
481  .transpose()
482  );
483 
484  dec().matrixQR()
485  .topLeftCorner(nonzero_pivots, nonzero_pivots)
486  .template triangularView<Upper>()
487  .solveInPlace(c.topRows(nonzero_pivots));
488 
489 
490  typename Rhs::PlainObject d(c);
491  d.topRows(nonzero_pivots)
492  = dec().matrixQR()
493  .topLeftCorner(nonzero_pivots, nonzero_pivots)
494  .template triangularView<Upper>()
495  * c.topRows(nonzero_pivots);
496 
497  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
498  for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
499  }
500 };
501 
502 } // end namespace internal
503 
505 template<typename MatrixType>
506 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
508 {
509  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
510  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
511 }
512 
517 template<typename Derived>
520 {
521  return ColPivHouseholderQR<PlainObject>(eval());
522 }
523 
524 } // end namespace Eigen
525 
526 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H