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;
53  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
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;
59 
67  : m_qr(),
68  m_hCoeffs(),
69  m_colsPermutation(),
70  m_colsTranspositions(),
71  m_temp(),
72  m_colSqNorms(),
73  m_isInitialized(false) {}
74 
81  ColPivHouseholderQR(Index rows, Index cols)
82  : m_qr(rows, cols),
83  m_hCoeffs((std::min)(rows,cols)),
84  m_colsPermutation(cols),
85  m_colsTranspositions(cols),
86  m_temp(cols),
87  m_colSqNorms(cols),
88  m_isInitialized(false),
89  m_usePrescribedThreshold(false) {}
90 
91  ColPivHouseholderQR(const MatrixType& matrix)
92  : m_qr(matrix.rows(), matrix.cols()),
93  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
94  m_colsPermutation(matrix.cols()),
95  m_colsTranspositions(matrix.cols()),
96  m_temp(matrix.cols()),
97  m_colSqNorms(matrix.cols()),
98  m_isInitialized(false),
99  m_usePrescribedThreshold(false)
100  {
101  compute(matrix);
102  }
103 
121  template<typename Rhs>
122  inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
123  solve(const MatrixBase<Rhs>& b) const
124  {
125  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
126  return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
127  }
128 
129  HouseholderSequenceType householderQ(void) const;
130 
133  const MatrixType& matrixQR() const
134  {
135  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
136  return m_qr;
137  }
138 
139  ColPivHouseholderQR& compute(const MatrixType& matrix);
140 
141  const PermutationType& colsPermutation() const
142  {
143  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
144  return m_colsPermutation;
145  }
146 
160  typename MatrixType::RealScalar absDeterminant() const;
161 
174  typename MatrixType::RealScalar logAbsDeterminant() const;
175 
182  inline Index rank() const
183  {
184  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
185  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
186  Index result = 0;
187  for(Index i = 0; i < m_nonzero_pivots; ++i)
188  result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
189  return result;
190  }
191 
198  inline Index dimensionOfKernel() const
199  {
200  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
201  return cols() - rank();
202  }
203 
211  inline bool isInjective() const
212  {
213  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
214  return rank() == cols();
215  }
216 
224  inline bool isSurjective() const
225  {
226  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
227  return rank() == rows();
228  }
229 
236  inline bool isInvertible() const
237  {
238  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
239  return isInjective() && isSurjective();
240  }
241 
247  inline const
248  internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
249  inverse() const
250  {
251  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
252  return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
253  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
254  }
255 
256  inline Index rows() const { return m_qr.rows(); }
257  inline Index cols() const { return m_qr.cols(); }
258  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
259 
278  {
279  m_usePrescribedThreshold = true;
280  m_prescribedThreshold = threshold;
281  return *this;
282  }
283 
293  {
294  m_usePrescribedThreshold = false;
295  return *this;
296  }
297 
302  RealScalar threshold() const
303  {
304  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
305  return m_usePrescribedThreshold ? m_prescribedThreshold
306  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
307  // and turns out to be identical to Higham's formula used already in LDLt.
308  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
309  }
310 
318  inline Index nonzeroPivots() const
319  {
320  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
321  return m_nonzero_pivots;
322  }
323 
327  RealScalar maxPivot() const { return m_maxpivot; }
328 
329  protected:
330  MatrixType m_qr;
331  HCoeffsType m_hCoeffs;
332  PermutationType m_colsPermutation;
333  IntRowVectorType m_colsTranspositions;
334  RowVectorType m_temp;
335  RealRowVectorType m_colSqNorms;
336  bool m_isInitialized, m_usePrescribedThreshold;
337  RealScalar m_prescribedThreshold, m_maxpivot;
338  Index m_nonzero_pivots;
339  Index m_det_pq;
340 };
341 
342 template<typename MatrixType>
343 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
344 {
345  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
346  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
347  return internal::abs(m_qr.diagonal().prod());
348 }
349 
350 template<typename MatrixType>
351 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
352 {
353  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
354  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
355  return m_qr.diagonal().cwiseAbs().array().log().sum();
356 }
357 
358 template<typename MatrixType>
360 {
361  Index rows = matrix.rows();
362  Index cols = matrix.cols();
363  Index size = matrix.diagonalSize();
364 
365  m_qr = matrix;
366  m_hCoeffs.resize(size);
367 
368  m_temp.resize(cols);
369 
370  m_colsTranspositions.resize(matrix.cols());
371  Index number_of_transpositions = 0;
372 
373  m_colSqNorms.resize(cols);
374  for(Index k = 0; k < cols; ++k)
375  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
376 
377  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
378 
379  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
380  m_maxpivot = RealScalar(0);
381 
382  for(Index k = 0; k < size; ++k)
383  {
384  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
385  Index biggest_col_index;
386  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
387  biggest_col_index += k;
388 
389  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
390  // the actual squared norm of the selected column.
391  // Note that not doing so does result in solve() sometimes returning inf/nan values
392  // when running the unit test with 1000 repetitions.
393  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
394 
395  // we store that back into our table: it can't hurt to correct our table.
396  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
397 
398  // if the current biggest column is smaller than epsilon times the initial biggest column,
399  // terminate to avoid generating nan/inf values.
400  // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
401  // repetitions of the unit test, with the result of solve() filled with large values of the order
402  // of 1/(size*epsilon).
403  if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
404  {
405  m_nonzero_pivots = k;
406  m_hCoeffs.tail(size-k).setZero();
407  m_qr.bottomRightCorner(rows-k,cols-k)
408  .template triangularView<StrictlyLower>()
409  .setZero();
410  break;
411  }
412 
413  // apply the transposition to the columns
414  m_colsTranspositions.coeffRef(k) = biggest_col_index;
415  if(k != biggest_col_index) {
416  m_qr.col(k).swap(m_qr.col(biggest_col_index));
417  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
418  ++number_of_transpositions;
419  }
420 
421  // generate the householder vector, store it below the diagonal
422  RealScalar beta;
423  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
424 
425  // apply the householder transformation to the diagonal coefficient
426  m_qr.coeffRef(k,k) = beta;
427 
428  // remember the maximum absolute value of diagonal coefficients
429  if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
430 
431  // apply the householder transformation
432  m_qr.bottomRightCorner(rows-k, cols-k-1)
433  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
434 
435  // update our table of squared norms of the columns
436  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
437  }
438 
439  m_colsPermutation.setIdentity(cols);
440  for(Index k = 0; k < m_nonzero_pivots; ++k)
441  m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
442 
443  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
444  m_isInitialized = true;
445 
446  return *this;
447 }
448 
449 namespace internal {
450 
451 template<typename _MatrixType, typename Rhs>
452 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
453  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
454 {
455  EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
456 
457  template<typename Dest> void evalTo(Dest& dst) const
458  {
459  eigen_assert(rhs().rows() == dec().rows());
460 
461  const int cols = dec().cols(),
462  nonzero_pivots = dec().nonzeroPivots();
463 
464  if(nonzero_pivots == 0)
465  {
466  dst.setZero();
467  return;
468  }
469 
470  typename Rhs::PlainObject c(rhs());
471 
472  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
473  c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
474  .setLength(dec().nonzeroPivots())
475  .transpose()
476  );
477 
478  dec().matrixQR()
479  .topLeftCorner(nonzero_pivots, nonzero_pivots)
480  .template triangularView<Upper>()
481  .solveInPlace(c.topRows(nonzero_pivots));
482 
483 
484  typename Rhs::PlainObject d(c);
485  d.topRows(nonzero_pivots)
486  = dec().matrixQR()
487  .topLeftCorner(nonzero_pivots, nonzero_pivots)
488  .template triangularView<Upper>()
489  * c.topRows(nonzero_pivots);
490 
491  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
492  for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
493  }
494 };
495 
496 } // end namespace internal
497 
499 template<typename MatrixType>
500 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
502 {
503  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
504  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
505 }
506 
511 template<typename Derived>
514 {
515  return ColPivHouseholderQR<PlainObject>(eval());
516 }
517 
518 } // end namespace Eigen
519 
520 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H