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 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00012 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00013 
00014 namespace Eigen { 
00015 
00037 template<typename _MatrixType> class ColPivHouseholderQR
00038 {
00039   public:
00040 
00041     typedef _MatrixType MatrixType;
00042     enum {
00043       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00044       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00045       Options = MatrixType::Options,
00046       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00047       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00048     };
00049     typedef typename MatrixType::Scalar Scalar;
00050     typedef typename MatrixType::RealScalar RealScalar;
00051     typedef typename MatrixType::Index Index;
00052     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00053     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00054     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00055     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00056     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00057     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00058     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00059 
00066     ColPivHouseholderQR()
00067       : m_qr(),
00068         m_hCoeffs(),
00069         m_colsPermutation(),
00070         m_colsTranspositions(),
00071         m_temp(),
00072         m_colSqNorms(),
00073         m_isInitialized(false) {}
00074 
00081     ColPivHouseholderQR(Index rows, Index cols)
00082       : m_qr(rows, cols),
00083         m_hCoeffs((std::min)(rows,cols)),
00084         m_colsPermutation(cols),
00085         m_colsTranspositions(cols),
00086         m_temp(cols),
00087         m_colSqNorms(cols),
00088         m_isInitialized(false),
00089         m_usePrescribedThreshold(false) {}
00090 
00091     ColPivHouseholderQR(const MatrixType& matrix)
00092       : m_qr(matrix.rows(), matrix.cols()),
00093         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00094         m_colsPermutation(matrix.cols()),
00095         m_colsTranspositions(matrix.cols()),
00096         m_temp(matrix.cols()),
00097         m_colSqNorms(matrix.cols()),
00098         m_isInitialized(false),
00099         m_usePrescribedThreshold(false)
00100     {
00101       compute(matrix);
00102     }
00103 
00121     template<typename Rhs>
00122     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00123     solve(const MatrixBase<Rhs>& b) const
00124     {
00125       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00126       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00127     }
00128 
00129     HouseholderSequenceType householderQ(void) const;
00130 
00133     const MatrixType& matrixQR() const
00134     {
00135       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00136       return m_qr;
00137     }
00138 
00139     ColPivHouseholderQR& compute(const MatrixType& matrix);
00140 
00141     const PermutationType& colsPermutation() const
00142     {
00143       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00144       return m_colsPermutation;
00145     }
00146 
00160     typename MatrixType::RealScalar absDeterminant() const;
00161 
00174     typename MatrixType::RealScalar logAbsDeterminant() const;
00175 
00182     inline Index rank() const
00183     {
00184       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00185       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00186       Index result = 0;
00187       for(Index i = 0; i < m_nonzero_pivots; ++i)
00188         result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00189       return result;
00190     }
00191 
00198     inline Index dimensionOfKernel() const
00199     {
00200       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00201       return cols() - rank();
00202     }
00203 
00211     inline bool isInjective() const
00212     {
00213       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00214       return rank() == cols();
00215     }
00216 
00224     inline bool isSurjective() const
00225     {
00226       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00227       return rank() == rows();
00228     }
00229 
00236     inline bool isInvertible() const
00237     {
00238       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00239       return isInjective() && isSurjective();
00240     }
00241 
00247     inline const
00248     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00249     inverse() const
00250     {
00251       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00252       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00253                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00254     }
00255 
00256     inline Index rows() const { return m_qr.rows(); }
00257     inline Index cols() const { return m_qr.cols(); }
00258     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00259 
00277     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00278     {
00279       m_usePrescribedThreshold = true;
00280       m_prescribedThreshold = threshold;
00281       return *this;
00282     }
00283 
00292     ColPivHouseholderQR& setThreshold(Default_t)
00293     {
00294       m_usePrescribedThreshold = false;
00295       return *this;
00296     }
00297 
00302     RealScalar threshold() const
00303     {
00304       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00305       return m_usePrescribedThreshold ? m_prescribedThreshold
00306       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00307       // and turns out to be identical to Higham's formula used already in LDLt.
00308                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00309     }
00310 
00318     inline Index nonzeroPivots() const
00319     {
00320       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00321       return m_nonzero_pivots;
00322     }
00323 
00327     RealScalar maxPivot() const { return m_maxpivot; }
00328 
00329   protected:
00330     MatrixType m_qr;
00331     HCoeffsType m_hCoeffs;
00332     PermutationType m_colsPermutation;
00333     IntRowVectorType m_colsTranspositions;
00334     RowVectorType m_temp;
00335     RealRowVectorType m_colSqNorms;
00336     bool m_isInitialized, m_usePrescribedThreshold;
00337     RealScalar m_prescribedThreshold, m_maxpivot;
00338     Index m_nonzero_pivots;
00339     Index m_det_pq;
00340 };
00341 
00342 template<typename MatrixType>
00343 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00344 {
00345   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00346   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00347   return internal::abs(m_qr.diagonal().prod());
00348 }
00349 
00350 template<typename MatrixType>
00351 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00352 {
00353   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00354   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00355   return m_qr.diagonal().cwiseAbs().array().log().sum();
00356 }
00357 
00358 template<typename MatrixType>
00359 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00360 {
00361   Index rows = matrix.rows();
00362   Index cols = matrix.cols();
00363   Index size = matrix.diagonalSize();
00364 
00365   m_qr = matrix;
00366   m_hCoeffs.resize(size);
00367 
00368   m_temp.resize(cols);
00369 
00370   m_colsTranspositions.resize(matrix.cols());
00371   Index number_of_transpositions = 0;
00372 
00373   m_colSqNorms.resize(cols);
00374   for(Index k = 0; k < cols; ++k)
00375     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00376 
00377   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00378 
00379   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00380   m_maxpivot = RealScalar(0);
00381 
00382   for(Index k = 0; k < size; ++k)
00383   {
00384     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
00385     Index biggest_col_index;
00386     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00387     biggest_col_index += k;
00388 
00389     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
00390     // the actual squared norm of the selected column.
00391     // Note that not doing so does result in solve() sometimes returning inf/nan values
00392     // when running the unit test with 1000 repetitions.
00393     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00394 
00395     // we store that back into our table: it can't hurt to correct our table.
00396     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00397 
00398     // if the current biggest column is smaller than epsilon times the initial biggest column,
00399     // terminate to avoid generating nan/inf values.
00400     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
00401     // repetitions of the unit test, with the result of solve() filled with large values of the order
00402     // of 1/(size*epsilon).
00403     if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00404     {
00405       m_nonzero_pivots = k;
00406       m_hCoeffs.tail(size-k).setZero();
00407       m_qr.bottomRightCorner(rows-k,cols-k)
00408           .template triangularView<StrictlyLower>()
00409           .setZero();
00410       break;
00411     }
00412 
00413     // apply the transposition to the columns
00414     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00415     if(k != biggest_col_index) {
00416       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00417       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00418       ++number_of_transpositions;
00419     }
00420 
00421     // generate the householder vector, store it below the diagonal
00422     RealScalar beta;
00423     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00424 
00425     // apply the householder transformation to the diagonal coefficient
00426     m_qr.coeffRef(k,k) = beta;
00427 
00428     // remember the maximum absolute value of diagonal coefficients
00429     if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00430 
00431     // apply the householder transformation
00432     m_qr.bottomRightCorner(rows-k, cols-k-1)
00433         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00434 
00435     // update our table of squared norms of the columns
00436     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00437   }
00438 
00439   m_colsPermutation.setIdentity(cols);
00440   for(Index k = 0; k < m_nonzero_pivots; ++k)
00441     m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00442 
00443   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00444   m_isInitialized = true;
00445 
00446   return *this;
00447 }
00448 
00449 namespace internal {
00450 
00451 template<typename _MatrixType, typename Rhs>
00452 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00453   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00454 {
00455   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00456 
00457   template<typename Dest> void evalTo(Dest& dst) const
00458   {
00459     eigen_assert(rhs().rows() == dec().rows());
00460 
00461     const int cols = dec().cols(),
00462     nonzero_pivots = dec().nonzeroPivots();
00463 
00464     if(nonzero_pivots == 0)
00465     {
00466       dst.setZero();
00467       return;
00468     }
00469 
00470     typename Rhs::PlainObject c(rhs());
00471 
00472     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00473     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00474                      .setLength(dec().nonzeroPivots())
00475                      .transpose()
00476       );
00477 
00478     dec().matrixQR()
00479        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00480        .template triangularView<Upper>()
00481        .solveInPlace(c.topRows(nonzero_pivots));
00482 
00483 
00484     typename Rhs::PlainObject d(c);
00485     d.topRows(nonzero_pivots)
00486       = dec().matrixQR()
00487        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00488        .template triangularView<Upper>()
00489        * c.topRows(nonzero_pivots);
00490 
00491     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00492     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00493   }
00494 };
00495 
00496 } // end namespace internal
00497 
00499 template<typename MatrixType>
00500 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00501   ::householderQ() const
00502 {
00503   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00504   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
00505 }
00506 
00511 template<typename Derived>
00512 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00513 MatrixBase<Derived>::colPivHouseholderQr() const
00514 {
00515   return ColPivHouseholderQR<PlainObject>(eval());
00516 }
00517 
00518 } // end namespace Eigen
00519 
00520 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:17