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 // 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_FULLPIVOTINGHOUSEHOLDERQR_H
27 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
28 
29 namespace Eigen {
30 
31 namespace internal {
32 
33 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
34 
35 template<typename MatrixType>
36 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
37 {
38  typedef typename MatrixType::PlainObject ReturnType;
39 };
40 
41 }
42 
64 template<typename _MatrixType> class FullPivHouseholderQR
65 {
66  public:
67 
68  typedef _MatrixType MatrixType;
69  enum {
74  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
75  };
76  typedef typename MatrixType::Scalar Scalar;
77  typedef typename MatrixType::RealScalar RealScalar;
78  typedef typename MatrixType::Index Index;
79  typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
80  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
83  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
84  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
85  typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
86 
93  : m_qr(),
94  m_hCoeffs(),
98  m_temp(),
99  m_isInitialized(false),
100  m_usePrescribedThreshold(false) {}
101 
109  : m_qr(rows, cols),
110  m_hCoeffs((std::min)(rows,cols)),
111  m_rows_transpositions(rows),
112  m_cols_transpositions(cols),
113  m_cols_permutation(cols),
114  m_temp((std::min)(rows,cols)),
115  m_isInitialized(false),
116  m_usePrescribedThreshold(false) {}
117 
119  : m_qr(matrix.rows(), matrix.cols()),
120  m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
121  m_rows_transpositions(matrix.rows()),
122  m_cols_transpositions(matrix.cols()),
123  m_cols_permutation(matrix.cols()),
124  m_temp((std::min)(matrix.rows(), matrix.cols())),
125  m_isInitialized(false),
127  {
128  compute(matrix);
129  }
130 
148  template<typename Rhs>
149  inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
150  solve(const MatrixBase<Rhs>& b) const
151  {
152  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
153  return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
154  }
155 
158  MatrixQReturnType matrixQ(void) const;
159 
162  const MatrixType& matrixQR() const
163  {
164  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
165  return m_qr;
166  }
167 
168  FullPivHouseholderQR& compute(const MatrixType& matrix);
169 
171  {
172  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
173  return m_cols_permutation;
174  }
175 
177  {
178  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
179  return m_rows_transpositions;
180  }
181 
195  typename MatrixType::RealScalar absDeterminant() const;
196 
209  typename MatrixType::RealScalar logAbsDeterminant() const;
210 
217  inline Index rank() const
218  {
219  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
220  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
221  Index result = 0;
222  for(Index i = 0; i < m_nonzero_pivots; ++i)
223  result += (internal::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
282  internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
283  inverse() const
284  {
285  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
286  return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
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  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
293 
312  {
315  return *this;
316  }
317 
327  {
328  m_usePrescribedThreshold = false;
329  return *this;
330  }
331 
337  {
340  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
341  // and turns out to be identical to Higham's formula used already in LDLt.
342  : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
343  }
344 
352  inline Index nonzeroPivots() const
353  {
354  eigen_assert(m_isInitialized && "LU is not initialized.");
355  return m_nonzero_pivots;
356  }
357 
361  RealScalar maxPivot() const { return m_maxpivot; }
362 
363  protected:
375 };
376 
377 template<typename MatrixType>
378 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
379 {
380  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
381  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
382  return internal::abs(m_qr.diagonal().prod());
383 }
384 
385 template<typename MatrixType>
386 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
387 {
388  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
389  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
390  return m_qr.diagonal().cwiseAbs().array().log().sum();
391 }
392 
393 template<typename MatrixType>
395 {
396  Index rows = matrix.rows();
397  Index cols = matrix.cols();
398  Index size = (std::min)(rows,cols);
399 
400  m_qr = matrix;
401  m_hCoeffs.resize(size);
402 
403  m_temp.resize(cols);
404 
405  m_precision = NumTraits<Scalar>::epsilon() * size;
406 
407  m_rows_transpositions.resize(matrix.rows());
408  m_cols_transpositions.resize(matrix.cols());
409  Index number_of_transpositions = 0;
410 
411  RealScalar biggest(0);
412 
413  m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
414  m_maxpivot = RealScalar(0);
415 
416  for (Index k = 0; k < size; ++k)
417  {
418  Index row_of_biggest_in_corner, col_of_biggest_in_corner;
419  RealScalar biggest_in_corner;
420 
421  biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
422  .cwiseAbs()
423  .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
424  row_of_biggest_in_corner += k;
425  col_of_biggest_in_corner += k;
426  if(k==0) biggest = biggest_in_corner;
427 
428  // if the corner is negligible, then we have less than full rank, and we can finish early
429  if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
430  {
431  m_nonzero_pivots = k;
432  for(Index i = k; i < size; i++)
433  {
434  m_rows_transpositions.coeffRef(i) = i;
435  m_cols_transpositions.coeffRef(i) = i;
436  m_hCoeffs.coeffRef(i) = Scalar(0);
437  }
438  break;
439  }
440 
441  m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
442  m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
443  if(k != row_of_biggest_in_corner) {
444  m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
445  ++number_of_transpositions;
446  }
447  if(k != col_of_biggest_in_corner) {
448  m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
449  ++number_of_transpositions;
450  }
451 
452  RealScalar beta;
453  m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
454  m_qr.coeffRef(k,k) = beta;
455 
456  // remember the maximum absolute value of diagonal coefficients
457  if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
458 
459  m_qr.bottomRightCorner(rows-k, cols-k-1)
460  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
461  }
462 
463  m_cols_permutation.setIdentity(cols);
464  for(Index k = 0; k < size; ++k)
465  m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
466 
467  m_det_pq = (number_of_transpositions%2) ? -1 : 1;
468  m_isInitialized = true;
469 
470  return *this;
471 }
472 
473 namespace internal {
474 
475 template<typename _MatrixType, typename Rhs>
476 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
477  : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
478 {
480 
481  template<typename Dest> void evalTo(Dest& dst) const
482  {
483  const Index rows = dec().rows(), cols = dec().cols();
484  eigen_assert(rhs().rows() == rows);
485 
486  // FIXME introduce nonzeroPivots() and use it here. and more generally,
487  // make the same improvements in this dec as in FullPivLU.
488  if(dec().rank()==0)
489  {
490  dst.setZero();
491  return;
492  }
493 
494  typename Rhs::PlainObject c(rhs());
495 
496  Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
497  for (Index k = 0; k < dec().rank(); ++k)
498  {
499  Index remainingSize = rows-k;
500  c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
501  c.bottomRightCorner(remainingSize, rhs().cols())
502  .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
503  dec().hCoeffs().coeff(k), &temp.coeffRef(0));
504  }
505 
506  if(!dec().isSurjective())
507  {
508  // is c is in the image of R ?
509  RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
510  RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
511  // FIXME brain dead
512  const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
513  // this internal:: prefix is needed by at least gcc 3.4 and ICC
514  if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
515  return;
516  }
517  dec().matrixQR()
518  .topLeftCorner(dec().rank(), dec().rank())
519  .template triangularView<Upper>()
520  .solveInPlace(c.topRows(dec().rank()));
521 
522  for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
523  for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
524  }
525 };
526 
533 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
534  : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
535 {
536 public:
537  typedef typename MatrixType::Index Index;
538  typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
539  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
540  typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
541  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
542 
543  FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
544  const HCoeffsType& hCoeffs,
545  const IntColVectorType& rowsTranspositions)
546  : m_qr(qr),
547  m_hCoeffs(hCoeffs),
548  m_rowsTranspositions(rowsTranspositions)
549  {}
550 
551  template <typename ResultType>
552  void evalTo(ResultType& result) const
553  {
554  const Index rows = m_qr.rows();
555  WorkVectorType workspace(rows);
556  evalTo(result, workspace);
557  }
558 
559  template <typename ResultType>
560  void evalTo(ResultType& result, WorkVectorType& workspace) const
561  {
562  // compute the product H'_0 H'_1 ... H'_n-1,
563  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
564  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
565  const Index rows = m_qr.rows();
566  const Index cols = m_qr.cols();
567  const Index size = (std::min)(rows, cols);
568  workspace.resize(rows);
569  result.setIdentity(rows, rows);
570  for (Index k = size-1; k >= 0; k--)
571  {
572  result.block(k, k, rows-k, rows-k)
573  .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
574  result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
575  }
576  }
577 
578  Index rows() const { return m_qr.rows(); }
579  Index cols() const { return m_qr.rows(); }
580 
581 protected:
582  typename MatrixType::Nested m_qr;
583  typename HCoeffsType::Nested m_hCoeffs;
584  typename IntColVectorType::Nested m_rowsTranspositions;
585 };
586 
587 } // end namespace internal
588 
589 template<typename MatrixType>
591 {
592  eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
593  return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
594 }
595 
600 template<typename Derived>
603 {
604  return FullPivHouseholderQR<PlainObject>(eval());
605 }
606 
607 } // end namespace Eigen
608 
609 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H