ColPivHouseholderQR.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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
27 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
28 
29 namespace Eigen {
30 
52 template<typename _MatrixType> class ColPivHouseholderQR
53 {
54  public:
55 
56  typedef _MatrixType MatrixType;
57  enum {
62  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
63  };
64  typedef typename MatrixType::Scalar Scalar;
65  typedef typename MatrixType::RealScalar RealScalar;
66  typedef typename MatrixType::Index Index;
68  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
70  typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
71  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
72  typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
74 
82  : m_qr(),
83  m_hCoeffs(),
86  m_temp(),
87  m_colSqNorms(),
88  m_isInitialized(false) {}
89 
97  : m_qr(rows, cols),
98  m_hCoeffs((std::min)(rows,cols)),
99  m_colsPermutation(cols),
100  m_colsTranspositions(cols),
101  m_temp(cols),
102  m_colSqNorms(cols),
103  m_isInitialized(false),
104  m_usePrescribedThreshold(false) {}
105 
107  : m_qr(matrix.rows(), matrix.cols()),
108  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
109  m_colsPermutation(matrix.cols()),
110  m_colsTranspositions(matrix.cols()),
111  m_temp(matrix.cols()),
112  m_colSqNorms(matrix.cols()),
113  m_isInitialized(false),
115  {
116  compute(matrix);
117  }
118 
136  template<typename Rhs>
137  inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
138  solve(const MatrixBase<Rhs>& b) const
139  {
140  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
141  return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
142  }
143 
145 
148  const MatrixType& matrixQR() const
149  {
150  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
151  return m_qr;
152  }
153 
154  ColPivHouseholderQR& compute(const MatrixType& matrix);
155 
157  {
158  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
159  return m_colsPermutation;
160  }
161 
175  typename MatrixType::RealScalar absDeterminant() const;
176 
189  typename MatrixType::RealScalar logAbsDeterminant() const;
190 
197  inline Index rank() const
198  {
199  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
200  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
201  Index result = 0;
202  for(Index i = 0; i < m_nonzero_pivots; ++i)
203  result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
204  return result;
205  }
206 
213  inline Index dimensionOfKernel() const
214  {
215  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
216  return cols() - rank();
217  }
218 
226  inline bool isInjective() const
227  {
228  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
229  return rank() == cols();
230  }
231 
239  inline bool isSurjective() const
240  {
241  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
242  return rank() == rows();
243  }
244 
251  inline bool isInvertible() const
252  {
253  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
254  return isInjective() && isSurjective();
255  }
256 
262  inline const
263  internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
264  inverse() const
265  {
266  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
267  return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
268  (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
269  }
270 
271  inline Index rows() const { return m_qr.rows(); }
272  inline Index cols() const { return m_qr.cols(); }
273  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
274 
293  {
296  return *this;
297  }
298 
308  {
309  m_usePrescribedThreshold = false;
310  return *this;
311  }
312 
318  {
321  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
322  // and turns out to be identical to Higham's formula used already in LDLt.
323  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
324  }
325 
333  inline Index nonzeroPivots() const
334  {
335  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
336  return m_nonzero_pivots;
337  }
338 
342  RealScalar maxPivot() const { return m_maxpivot; }
343 
344  protected:
355 };
356 
357 template<typename MatrixType>
358 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
359 {
360  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
361  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
362  return internal::abs(m_qr.diagonal().prod());
363 }
364 
365 template<typename MatrixType>
366 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
367 {
368  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
369  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
370  return m_qr.diagonal().cwiseAbs().array().log().sum();
371 }
372 
373 template<typename MatrixType>
375 {
376  Index rows = matrix.rows();
377  Index cols = matrix.cols();
378  Index size = matrix.diagonalSize();
379 
380  m_qr = matrix;
381  m_hCoeffs.resize(size);
382 
383  m_temp.resize(cols);
384 
385  m_colsTranspositions.resize(matrix.cols());
386  Index number_of_transpositions = 0;
387 
388  m_colSqNorms.resize(cols);
389  for(Index k = 0; k < cols; ++k)
390  m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
391 
392  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
393 
394  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
395  m_maxpivot = RealScalar(0);
396 
397  for(Index k = 0; k < size; ++k)
398  {
399  // first, we look up in our table m_colSqNorms which column has the biggest squared norm
400  Index biggest_col_index;
401  RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
402  biggest_col_index += k;
403 
404  // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
405  // the actual squared norm of the selected column.
406  // Note that not doing so does result in solve() sometimes returning inf/nan values
407  // when running the unit test with 1000 repetitions.
408  biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
409 
410  // we store that back into our table: it can't hurt to correct our table.
411  m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
412 
413  // if the current biggest column is smaller than epsilon times the initial biggest column,
414  // terminate to avoid generating nan/inf values.
415  // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
416  // repetitions of the unit test, with the result of solve() filled with large values of the order
417  // of 1/(size*epsilon).
418  if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
419  {
420  m_nonzero_pivots = k;
421  m_hCoeffs.tail(size-k).setZero();
422  m_qr.bottomRightCorner(rows-k,cols-k)
423  .template triangularView<StrictlyLower>()
424  .setZero();
425  break;
426  }
427 
428  // apply the transposition to the columns
429  m_colsTranspositions.coeffRef(k) = biggest_col_index;
430  if(k != biggest_col_index) {
431  m_qr.col(k).swap(m_qr.col(biggest_col_index));
432  std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
433  ++number_of_transpositions;
434  }
435 
436  // generate the householder vector, store it below the diagonal
437  RealScalar beta;
438  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
439 
440  // apply the householder transformation to the diagonal coefficient
441  m_qr.coeffRef(k,k) = beta;
442 
443  // remember the maximum absolute value of diagonal coefficients
444  if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
445 
446  // apply the householder transformation
447  m_qr.bottomRightCorner(rows-k, cols-k-1)
448  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
449 
450  // update our table of squared norms of the columns
451  m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
452  }
453 
454  m_colsPermutation.setIdentity(cols);
455  for(Index k = 0; k < m_nonzero_pivots; ++k)
456  m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
457 
458  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
459  m_isInitialized = true;
460 
461  return *this;
462 }
463 
464 namespace internal {
465 
466 template<typename _MatrixType, typename Rhs>
467 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
468  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
469 {
471 
472  template<typename Dest> void evalTo(Dest& dst) const
473  {
474  eigen_assert(rhs().rows() == dec().rows());
475 
476  const int cols = dec().cols(),
477  nonzero_pivots = dec().nonzeroPivots();
478 
479  if(nonzero_pivots == 0)
480  {
481  dst.setZero();
482  return;
483  }
484 
485  typename Rhs::PlainObject c(rhs());
486 
487  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
488  c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
489  .setLength(dec().nonzeroPivots())
490  .transpose()
491  );
492 
493  dec().matrixQR()
494  .topLeftCorner(nonzero_pivots, nonzero_pivots)
495  .template triangularView<Upper>()
496  .solveInPlace(c.topRows(nonzero_pivots));
497 
498 
499  typename Rhs::PlainObject d(c);
500  d.topRows(nonzero_pivots)
501  = dec().matrixQR()
502  .topLeftCorner(nonzero_pivots, nonzero_pivots)
503  .template triangularView<Upper>()
504  * c.topRows(nonzero_pivots);
505 
506  for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
507  for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
508  }
509 };
510 
511 } // end namespace internal
512 
514 template<typename MatrixType>
515 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
517 {
518  eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
519  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
520 }
521 
526 template<typename Derived>
529 {
530  return ColPivHouseholderQR<PlainObject>(eval());
531 }
532 
533 } // end namespace Eigen
534 
535 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H