ColPivHouseholderQR.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00028 
00029 namespace Eigen { 
00030 
00052 template<typename _MatrixType> class ColPivHouseholderQR
00053 {
00054   public:
00055 
00056     typedef _MatrixType MatrixType;
00057     enum {
00058       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00059       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00060       Options = MatrixType::Options,
00061       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00062       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00063     };
00064     typedef typename MatrixType::Scalar Scalar;
00065     typedef typename MatrixType::RealScalar RealScalar;
00066     typedef typename MatrixType::Index Index;
00067     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00068     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00069     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00070     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00071     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00072     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00073     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00074 
00081     ColPivHouseholderQR()
00082       : m_qr(),
00083         m_hCoeffs(),
00084         m_colsPermutation(),
00085         m_colsTranspositions(),
00086         m_temp(),
00087         m_colSqNorms(),
00088         m_isInitialized(false) {}
00089 
00096     ColPivHouseholderQR(Index rows, Index cols)
00097       : m_qr(rows, cols),
00098         m_hCoeffs((std::min)(rows,cols)),
00099         m_colsPermutation(cols),
00100         m_colsTranspositions(cols),
00101         m_temp(cols),
00102         m_colSqNorms(cols),
00103         m_isInitialized(false),
00104         m_usePrescribedThreshold(false) {}
00105 
00106     ColPivHouseholderQR(const MatrixType& matrix)
00107       : m_qr(matrix.rows(), matrix.cols()),
00108         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00109         m_colsPermutation(matrix.cols()),
00110         m_colsTranspositions(matrix.cols()),
00111         m_temp(matrix.cols()),
00112         m_colSqNorms(matrix.cols()),
00113         m_isInitialized(false),
00114         m_usePrescribedThreshold(false)
00115     {
00116       compute(matrix);
00117     }
00118 
00136     template<typename Rhs>
00137     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00138     solve(const MatrixBase<Rhs>& b) const
00139     {
00140       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00141       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00142     }
00143 
00144     HouseholderSequenceType householderQ(void) const;
00145 
00148     const MatrixType& matrixQR() const
00149     {
00150       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00151       return m_qr;
00152     }
00153 
00154     ColPivHouseholderQR& compute(const MatrixType& matrix);
00155 
00156     const PermutationType& colsPermutation() const
00157     {
00158       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00159       return m_colsPermutation;
00160     }
00161 
00175     typename MatrixType::RealScalar absDeterminant() const;
00176 
00189     typename MatrixType::RealScalar logAbsDeterminant() const;
00190 
00197     inline Index rank() const
00198     {
00199       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00200       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00201       Index result = 0;
00202       for(Index i = 0; i < m_nonzero_pivots; ++i)
00203         result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00204       return result;
00205     }
00206 
00213     inline Index dimensionOfKernel() const
00214     {
00215       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00216       return cols() - rank();
00217     }
00218 
00226     inline bool isInjective() const
00227     {
00228       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00229       return rank() == cols();
00230     }
00231 
00239     inline bool isSurjective() const
00240     {
00241       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00242       return rank() == rows();
00243     }
00244 
00251     inline bool isInvertible() const
00252     {
00253       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00254       return isInjective() && isSurjective();
00255     }
00256 
00262     inline const
00263     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00264     inverse() const
00265     {
00266       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00267       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00268                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00269     }
00270 
00271     inline Index rows() const { return m_qr.rows(); }
00272     inline Index cols() const { return m_qr.cols(); }
00273     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00274 
00292     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00293     {
00294       m_usePrescribedThreshold = true;
00295       m_prescribedThreshold = threshold;
00296       return *this;
00297     }
00298 
00307     ColPivHouseholderQR& setThreshold(Default_t)
00308     {
00309       m_usePrescribedThreshold = false;
00310       return *this;
00311     }
00312 
00317     RealScalar threshold() const
00318     {
00319       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00320       return m_usePrescribedThreshold ? m_prescribedThreshold
00321       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00322       // and turns out to be identical to Higham's formula used already in LDLt.
00323                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00324     }
00325 
00333     inline Index nonzeroPivots() const
00334     {
00335       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00336       return m_nonzero_pivots;
00337     }
00338 
00342     RealScalar maxPivot() const { return m_maxpivot; }
00343 
00344   protected:
00345     MatrixType m_qr;
00346     HCoeffsType m_hCoeffs;
00347     PermutationType m_colsPermutation;
00348     IntRowVectorType m_colsTranspositions;
00349     RowVectorType m_temp;
00350     RealRowVectorType m_colSqNorms;
00351     bool m_isInitialized, m_usePrescribedThreshold;
00352     RealScalar m_prescribedThreshold, m_maxpivot;
00353     Index m_nonzero_pivots;
00354     Index m_det_pq;
00355 };
00356 
00357 template<typename MatrixType>
00358 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00359 {
00360   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00361   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00362   return internal::abs(m_qr.diagonal().prod());
00363 }
00364 
00365 template<typename MatrixType>
00366 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00367 {
00368   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00369   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00370   return m_qr.diagonal().cwiseAbs().array().log().sum();
00371 }
00372 
00373 template<typename MatrixType>
00374 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00375 {
00376   Index rows = matrix.rows();
00377   Index cols = matrix.cols();
00378   Index size = matrix.diagonalSize();
00379 
00380   m_qr = matrix;
00381   m_hCoeffs.resize(size);
00382 
00383   m_temp.resize(cols);
00384 
00385   m_colsTranspositions.resize(matrix.cols());
00386   Index number_of_transpositions = 0;
00387 
00388   m_colSqNorms.resize(cols);
00389   for(Index k = 0; k < cols; ++k)
00390     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00391 
00392   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00393 
00394   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00395   m_maxpivot = RealScalar(0);
00396 
00397   for(Index k = 0; k < size; ++k)
00398   {
00399     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
00400     Index biggest_col_index;
00401     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00402     biggest_col_index += k;
00403 
00404     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
00405     // the actual squared norm of the selected column.
00406     // Note that not doing so does result in solve() sometimes returning inf/nan values
00407     // when running the unit test with 1000 repetitions.
00408     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00409 
00410     // we store that back into our table: it can't hurt to correct our table.
00411     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00412 
00413     // if the current biggest column is smaller than epsilon times the initial biggest column,
00414     // terminate to avoid generating nan/inf values.
00415     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
00416     // repetitions of the unit test, with the result of solve() filled with large values of the order
00417     // of 1/(size*epsilon).
00418     if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00419     {
00420       m_nonzero_pivots = k;
00421       m_hCoeffs.tail(size-k).setZero();
00422       m_qr.bottomRightCorner(rows-k,cols-k)
00423           .template triangularView<StrictlyLower>()
00424           .setZero();
00425       break;
00426     }
00427 
00428     // apply the transposition to the columns
00429     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00430     if(k != biggest_col_index) {
00431       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00432       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00433       ++number_of_transpositions;
00434     }
00435 
00436     // generate the householder vector, store it below the diagonal
00437     RealScalar beta;
00438     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00439 
00440     // apply the householder transformation to the diagonal coefficient
00441     m_qr.coeffRef(k,k) = beta;
00442 
00443     // remember the maximum absolute value of diagonal coefficients
00444     if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00445 
00446     // apply the householder transformation
00447     m_qr.bottomRightCorner(rows-k, cols-k-1)
00448         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00449 
00450     // update our table of squared norms of the columns
00451     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00452   }
00453 
00454   m_colsPermutation.setIdentity(cols);
00455   for(Index k = 0; k < m_nonzero_pivots; ++k)
00456     m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00457 
00458   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00459   m_isInitialized = true;
00460 
00461   return *this;
00462 }
00463 
00464 namespace internal {
00465 
00466 template<typename _MatrixType, typename Rhs>
00467 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00468   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00469 {
00470   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00471 
00472   template<typename Dest> void evalTo(Dest& dst) const
00473   {
00474     eigen_assert(rhs().rows() == dec().rows());
00475 
00476     const int cols = dec().cols(),
00477     nonzero_pivots = dec().nonzeroPivots();
00478 
00479     if(nonzero_pivots == 0)
00480     {
00481       dst.setZero();
00482       return;
00483     }
00484 
00485     typename Rhs::PlainObject c(rhs());
00486 
00487     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00488     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00489                      .setLength(dec().nonzeroPivots())
00490                      .transpose()
00491       );
00492 
00493     dec().matrixQR()
00494        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00495        .template triangularView<Upper>()
00496        .solveInPlace(c.topRows(nonzero_pivots));
00497 
00498 
00499     typename Rhs::PlainObject d(c);
00500     d.topRows(nonzero_pivots)
00501       = dec().matrixQR()
00502        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00503        .template triangularView<Upper>()
00504        * c.topRows(nonzero_pivots);
00505 
00506     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00507     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00508   }
00509 };
00510 
00511 } // end namespace internal
00512 
00514 template<typename MatrixType>
00515 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00516   ::householderQ() const
00517 {
00518   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00519   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
00520 }
00521 
00526 template<typename Derived>
00527 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00528 MatrixBase<Derived>::colPivHouseholderQr() const
00529 {
00530   return ColPivHouseholderQR<PlainObject>(eval());
00531 }
00532 
00533 } // end namespace Eigen
00534 
00535 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H