Main MRPT website > C++ reference
MRPT logo

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 /** \ingroup QR_Module
00030   *
00031   * \class ColPivHouseholderQR
00032   *
00033   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
00034   *
00035   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
00036   *
00037   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
00038   * such that 
00039   * \f[
00040   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
00041   * \f]
00042   * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 
00043   * upper triangular matrix.
00044   *
00045   * This decomposition performs column pivoting in order to be rank-revealing and improve
00046   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
00047   *
00048   * \sa MatrixBase::colPivHouseholderQr()
00049   */
00050 template<typename _MatrixType> class ColPivHouseholderQR
00051 {
00052   public:
00053 
00054     typedef _MatrixType MatrixType;
00055     enum {
00056       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00057       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00058       Options = MatrixType::Options,
00059       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00060       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00061     };
00062     typedef typename MatrixType::Scalar Scalar;
00063     typedef typename MatrixType::RealScalar RealScalar;
00064     typedef typename MatrixType::Index Index;
00065     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00066     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00067     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00068     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00069     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00070     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00071     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00072 
00073     /**
00074     * \brief Default Constructor.
00075     *
00076     * The default constructor is useful in cases in which the user intends to
00077     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
00078     */
00079     ColPivHouseholderQR()
00080       : m_qr(),
00081         m_hCoeffs(),
00082         m_colsPermutation(),
00083         m_colsTranspositions(),
00084         m_temp(),
00085         m_colSqNorms(),
00086         m_isInitialized(false) {}
00087 
00088     /** \brief Default Constructor with memory preallocation
00089       *
00090       * Like the default constructor but with preallocation of the internal data
00091       * according to the specified problem \a size.
00092       * \sa ColPivHouseholderQR()
00093       */
00094     ColPivHouseholderQR(Index rows, Index cols)
00095       : m_qr(rows, cols),
00096         m_hCoeffs(std::min(rows,cols)),
00097         m_colsPermutation(cols),
00098         m_colsTranspositions(cols),
00099         m_temp(cols),
00100         m_colSqNorms(cols),
00101         m_isInitialized(false),
00102         m_usePrescribedThreshold(false) {}
00103 
00104     ColPivHouseholderQR(const MatrixType& matrix)
00105       : m_qr(matrix.rows(), matrix.cols()),
00106         m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
00107         m_colsPermutation(matrix.cols()),
00108         m_colsTranspositions(matrix.cols()),
00109         m_temp(matrix.cols()),
00110         m_colSqNorms(matrix.cols()),
00111         m_isInitialized(false),
00112         m_usePrescribedThreshold(false)
00113     {
00114       compute(matrix);
00115     }
00116 
00117     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
00118       * *this is the QR decomposition, if any exists.
00119       *
00120       * \param b the right-hand-side of the equation to solve.
00121       *
00122       * \returns a solution.
00123       *
00124       * \note The case where b is a matrix is not yet implemented. Also, this
00125       *       code is space inefficient.
00126       *
00127       * \note_about_checking_solutions
00128       *
00129       * \note_about_arbitrary_choice_of_solution
00130       *
00131       * Example: \include ColPivHouseholderQR_solve.cpp
00132       * Output: \verbinclude ColPivHouseholderQR_solve.out
00133       */
00134     template<typename Rhs>
00135     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00136     solve(const MatrixBase<Rhs>& b) const
00137     {
00138       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00139       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00140     }
00141 
00142     HouseholderSequenceType householderQ(void) const;
00143 
00144     /** \returns a reference to the matrix where the Householder QR decomposition is stored
00145       */
00146     const MatrixType& matrixQR() const
00147     {
00148       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00149       return m_qr;
00150     }
00151 
00152     ColPivHouseholderQR& compute(const MatrixType& matrix);
00153 
00154     const PermutationType& colsPermutation() const
00155     {
00156       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00157       return m_colsPermutation;
00158     }
00159 
00160     /** \returns the absolute value of the determinant of the matrix of which
00161       * *this is the QR decomposition. It has only linear complexity
00162       * (that is, O(n) where n is the dimension of the square matrix)
00163       * as the QR decomposition has already been computed.
00164       *
00165       * \note This is only for square matrices.
00166       *
00167       * \warning a determinant can be very big or small, so for matrices
00168       * of large enough dimension, there is a risk of overflow/underflow.
00169       * One way to work around that is to use logAbsDeterminant() instead.
00170       *
00171       * \sa logAbsDeterminant(), MatrixBase::determinant()
00172       */
00173     typename MatrixType::RealScalar absDeterminant() const;
00174 
00175     /** \returns the natural log of the absolute value of the determinant of the matrix of which
00176       * *this is the QR decomposition. It has only linear complexity
00177       * (that is, O(n) where n is the dimension of the square matrix)
00178       * as the QR decomposition has already been computed.
00179       *
00180       * \note This is only for square matrices.
00181       *
00182       * \note This method is useful to work around the risk of overflow/underflow that's inherent
00183       * to determinant computation.
00184       *
00185       * \sa absDeterminant(), MatrixBase::determinant()
00186       */
00187     typename MatrixType::RealScalar logAbsDeterminant() const;
00188 
00189     /** \returns the rank of the matrix of which *this is the QR decomposition.
00190       *
00191       * \note This method has to determine which pivots should be considered nonzero.
00192       *       For that, it uses the threshold value that you can control by calling
00193       *       setThreshold(const RealScalar&).
00194       */
00195     inline Index rank() const
00196     {
00197       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00198       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00199       Index result = 0;
00200       for(Index i = 0; i < m_nonzero_pivots; ++i)
00201         result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00202       return result;
00203     }
00204 
00205     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
00206       *
00207       * \note This method has to determine which pivots should be considered nonzero.
00208       *       For that, it uses the threshold value that you can control by calling
00209       *       setThreshold(const RealScalar&).
00210       */
00211     inline Index dimensionOfKernel() const
00212     {
00213       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00214       return cols() - rank();
00215     }
00216 
00217     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
00218       *          linear map, i.e. has trivial kernel; false otherwise.
00219       *
00220       * \note This method has to determine which pivots should be considered nonzero.
00221       *       For that, it uses the threshold value that you can control by calling
00222       *       setThreshold(const RealScalar&).
00223       */
00224     inline bool isInjective() const
00225     {
00226       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00227       return rank() == cols();
00228     }
00229 
00230     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
00231       *          linear map; false otherwise.
00232       *
00233       * \note This method has to determine which pivots should be considered nonzero.
00234       *       For that, it uses the threshold value that you can control by calling
00235       *       setThreshold(const RealScalar&).
00236       */
00237     inline bool isSurjective() const
00238     {
00239       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00240       return rank() == rows();
00241     }
00242 
00243     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
00244       *
00245       * \note This method has to determine which pivots should be considered nonzero.
00246       *       For that, it uses the threshold value that you can control by calling
00247       *       setThreshold(const RealScalar&).
00248       */
00249     inline bool isInvertible() const
00250     {
00251       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00252       return isInjective() && isSurjective();
00253     }
00254 
00255     /** \returns the inverse of the matrix of which *this is the QR decomposition.
00256       *
00257       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
00258       *       Use isInvertible() to first determine whether this matrix is invertible.
00259       */
00260     inline const
00261     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00262     inverse() const
00263     {
00264       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00265       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00266                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00267     }
00268 
00269     inline Index rows() const { return m_qr.rows(); }
00270     inline Index cols() const { return m_qr.cols(); }
00271     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00272 
00273     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
00274       * who need to determine when pivots are to be considered nonzero. This is not used for the
00275       * QR decomposition itself.
00276       *
00277       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
00278       * uses a formula to automatically determine a reasonable threshold.
00279       * Once you have called the present method setThreshold(const RealScalar&),
00280       * your value is used instead.
00281       *
00282       * \param threshold The new value to use as the threshold.
00283       *
00284       * A pivot will be considered nonzero if its absolute value is strictly greater than
00285       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
00286       * where maxpivot is the biggest pivot.
00287       *
00288       * If you want to come back to the default behavior, call setThreshold(Default_t)
00289       */
00290     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00291     {
00292       m_usePrescribedThreshold = true;
00293       m_prescribedThreshold = threshold;
00294       return *this;
00295     }
00296 
00297     /** Allows to come back to the default behavior, letting Eigen use its default formula for
00298       * determining the threshold.
00299       *
00300       * You should pass the special object Eigen::Default as parameter here.
00301       * \code qr.setThreshold(Eigen::Default); \endcode
00302       *
00303       * See the documentation of setThreshold(const RealScalar&).
00304       */
00305     ColPivHouseholderQR& setThreshold(Default_t)
00306     {
00307       m_usePrescribedThreshold = false;
00308       return *this;
00309     }
00310 
00311     /** Returns the threshold that will be used by certain methods such as rank().
00312       *
00313       * See the documentation of setThreshold(const RealScalar&).
00314       */
00315     RealScalar threshold() const
00316     {
00317       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00318       return m_usePrescribedThreshold ? m_prescribedThreshold
00319       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00320       // and turns out to be identical to Higham's formula used already in LDLt.
00321                                       : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00322     }
00323 
00324     /** \returns the number of nonzero pivots in the QR decomposition.
00325       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
00326       * So that notion isn't really intrinsically interesting, but it is
00327       * still useful when implementing algorithms.
00328       *
00329       * \sa rank()
00330       */
00331     inline Index nonzeroPivots() const
00332     {
00333       eigen_assert(m_isInitialized && "LU is not initialized.");
00334       return m_nonzero_pivots;
00335     }
00336 
00337     /** \returns the absolute value of the biggest pivot, i.e. the biggest
00338       *          diagonal coefficient of U.
00339       */
00340     RealScalar maxPivot() const { return m_maxpivot; }
00341 
00342   protected:
00343     MatrixType m_qr;
00344     HCoeffsType m_hCoeffs;
00345     PermutationType m_colsPermutation;
00346     IntRowVectorType m_colsTranspositions;
00347     RowVectorType m_temp;
00348     RealRowVectorType m_colSqNorms;
00349     bool m_isInitialized, m_usePrescribedThreshold;
00350     RealScalar m_prescribedThreshold, m_maxpivot;
00351     Index m_nonzero_pivots;
00352     Index m_det_pq;
00353 };
00354 
00355 template<typename MatrixType>
00356 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
00357 {
00358   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00359   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00360   return internal::abs(m_qr.diagonal().prod());
00361 }
00362 
00363 template<typename MatrixType>
00364 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00365 {
00366   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00367   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00368   return m_qr.diagonal().cwiseAbs().array().log().sum();
00369 }
00370 
00371 template<typename MatrixType>
00372 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00373 {
00374   Index rows = matrix.rows();
00375   Index cols = matrix.cols();
00376   Index size = matrix.diagonalSize();
00377 
00378   m_qr = matrix;
00379   m_hCoeffs.resize(size);
00380 
00381   m_temp.resize(cols);
00382 
00383   m_colsTranspositions.resize(matrix.cols());
00384   Index number_of_transpositions = 0;
00385 
00386   m_colSqNorms.resize(cols);
00387   for(Index k = 0; k < cols; ++k)
00388     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00389 
00390   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / rows;
00391 
00392   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00393   m_maxpivot = RealScalar(0);
00394 
00395   for(Index k = 0; k < size; ++k)
00396   {
00397     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
00398     Index biggest_col_index;
00399     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00400     biggest_col_index += k;
00401 
00402     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
00403     // the actual squared norm of the selected column.
00404     // Note that not doing so does result in solve() sometimes returning inf/nan values
00405     // when running the unit test with 1000 repetitions.
00406     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00407 
00408     // we store that back into our table: it can't hurt to correct our table.
00409     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00410 
00411     // if the current biggest column is smaller than epsilon times the initial biggest column,
00412     // terminate to avoid generating nan/inf values.
00413     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
00414     // repetitions of the unit test, with the result of solve() filled with large values of the order
00415     // of 1/(size*epsilon).
00416     if(biggest_col_sq_norm < threshold_helper * (rows-k))
00417     {
00418       m_nonzero_pivots = k;
00419       m_hCoeffs.tail(size-k).setZero();
00420       m_qr.bottomRightCorner(rows-k,cols-k)
00421           .template triangularView<StrictlyLower>()
00422           .setZero();
00423       break;
00424     }
00425 
00426     // apply the transposition to the columns
00427     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00428     if(k != biggest_col_index) {
00429       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00430       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00431       ++number_of_transpositions;
00432     }
00433 
00434     // generate the householder vector, store it below the diagonal
00435     RealScalar beta;
00436     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00437 
00438     // apply the householder transformation to the diagonal coefficient
00439     m_qr.coeffRef(k,k) = beta;
00440 
00441     // remember the maximum absolute value of diagonal coefficients
00442     if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00443 
00444     // apply the householder transformation
00445     m_qr.bottomRightCorner(rows-k, cols-k-1)
00446         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00447 
00448     // update our table of squared norms of the columns
00449     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00450   }
00451 
00452   m_colsPermutation.setIdentity(cols);
00453   for(Index k = 0; k < m_nonzero_pivots; ++k)
00454     m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00455 
00456   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00457   m_isInitialized = true;
00458 
00459   return *this;
00460 }
00461 
00462 namespace internal {
00463 
00464 template<typename _MatrixType, typename Rhs>
00465 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00466   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00467 {
00468   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00469 
00470   template<typename Dest> void evalTo(Dest& dst) const
00471   {
00472     eigen_assert(rhs().rows() == dec().rows());
00473 
00474     const int cols = dec().cols(),
00475     nonzero_pivots = dec().nonzeroPivots();
00476 
00477     if(nonzero_pivots == 0)
00478     {
00479       dst.setZero();
00480       return;
00481     }
00482 
00483     typename Rhs::PlainObject c(rhs());
00484 
00485     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00486     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00487                      .setTrans(true)
00488                      .setLength(dec().nonzeroPivots())
00489       );
00490 
00491     dec().matrixQR()
00492        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00493        .template triangularView<Upper>()
00494        .solveInPlace(c.topRows(nonzero_pivots));
00495 
00496 
00497     typename Rhs::PlainObject d(c);
00498     d.topRows(nonzero_pivots)
00499       = dec().matrixQR()
00500        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00501        .template triangularView<Upper>()
00502        * c.topRows(nonzero_pivots);
00503 
00504     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00505     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00506   }
00507 };
00508 
00509 } // end namespace internal
00510 
00511 /** \returns the matrix Q as a sequence of householder transformations */
00512 template<typename MatrixType>
00513 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00514   ::householderQ() const
00515 {
00516   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00517   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()).setLength(m_nonzero_pivots);
00518 }
00519 
00520 /** \return the column-pivoting Householder QR decomposition of \c *this.
00521   *
00522   * \sa class ColPivHouseholderQR
00523   */
00524 template<typename Derived>
00525 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00526 MatrixBase<Derived>::colPivHouseholderQr() const
00527 {
00528   return ColPivHouseholderQR<PlainObject>(eval());
00529 }
00530 
00531 
00532 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011