Main MRPT website > C++ reference
MRPT logo

FullPivHouseholderQR.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_FULLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00028 
00029 /** \ingroup QR_Module
00030   *
00031   * \class FullPivHouseholderQR
00032   *
00033   * \brief Householder rank-revealing QR decomposition of a matrix with full 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 a very prudent full pivoting in order to be rank-revealing and achieve optimal
00046   * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR.
00047   *
00048   * \sa MatrixBase::fullPivHouseholderQr()
00049   */
00050 template<typename _MatrixType> class FullPivHouseholderQR
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 Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00068     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00069     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00070     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00071     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00072 
00073     /** \brief Default Constructor.
00074       *
00075       * The default constructor is useful in cases in which the user intends to
00076       * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
00077       */
00078     FullPivHouseholderQR()
00079       : m_qr(),
00080         m_hCoeffs(),
00081         m_rows_transpositions(),
00082         m_cols_transpositions(),
00083         m_cols_permutation(),
00084         m_temp(),
00085         m_isInitialized(false) {}
00086 
00087     /** \brief Default Constructor with memory preallocation
00088       *
00089       * Like the default constructor but with preallocation of the internal data
00090       * according to the specified problem \a size.
00091       * \sa FullPivHouseholderQR()
00092       */
00093     FullPivHouseholderQR(Index rows, Index cols)
00094       : m_qr(rows, cols),
00095         m_hCoeffs(std::min(rows,cols)),
00096         m_rows_transpositions(rows),
00097         m_cols_transpositions(cols),
00098         m_cols_permutation(cols),
00099         m_temp(std::min(rows,cols)),
00100         m_isInitialized(false) {}
00101 
00102     FullPivHouseholderQR(const MatrixType& matrix)
00103       : m_qr(matrix.rows(), matrix.cols()),
00104         m_hCoeffs(std::min(matrix.rows(), matrix.cols())),
00105         m_rows_transpositions(matrix.rows()),
00106         m_cols_transpositions(matrix.cols()),
00107         m_cols_permutation(matrix.cols()),
00108         m_temp(std::min(matrix.rows(), matrix.cols())),
00109         m_isInitialized(false)
00110     {
00111       compute(matrix);
00112     }
00113 
00114     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
00115       * *this is the QR decomposition, if any exists.
00116       *
00117       * \param b the right-hand-side of the equation to solve.
00118       *
00119       * \returns a solution.
00120       *
00121       * \note The case where b is a matrix is not yet implemented. Also, this
00122       *       code is space inefficient.
00123       *
00124       * \note_about_checking_solutions
00125       *
00126       * \note_about_arbitrary_choice_of_solution
00127       *
00128       * Example: \include FullPivHouseholderQR_solve.cpp
00129       * Output: \verbinclude FullPivHouseholderQR_solve.out
00130       */
00131     template<typename Rhs>
00132     inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00133     solve(const MatrixBase<Rhs>& b) const
00134     {
00135       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00136       return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00137     }
00138 
00139     MatrixQType matrixQ(void) const;
00140 
00141     /** \returns a reference to the matrix where the Householder QR decomposition is stored
00142       */
00143     const MatrixType& matrixQR() const
00144     {
00145       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00146       return m_qr;
00147     }
00148 
00149     FullPivHouseholderQR& compute(const MatrixType& matrix);
00150 
00151     const PermutationType& colsPermutation() const
00152     {
00153       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00154       return m_cols_permutation;
00155     }
00156 
00157     const IntColVectorType& rowsTranspositions() const
00158     {
00159       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00160       return m_rows_transpositions;
00161     }
00162 
00163     /** \returns the absolute value of the determinant of the matrix of which
00164       * *this is the QR decomposition. It has only linear complexity
00165       * (that is, O(n) where n is the dimension of the square matrix)
00166       * as the QR decomposition has already been computed.
00167       *
00168       * \note This is only for square matrices.
00169       *
00170       * \warning a determinant can be very big or small, so for matrices
00171       * of large enough dimension, there is a risk of overflow/underflow.
00172       * One way to work around that is to use logAbsDeterminant() instead.
00173       *
00174       * \sa logAbsDeterminant(), MatrixBase::determinant()
00175       */
00176     typename MatrixType::RealScalar absDeterminant() const;
00177 
00178     /** \returns the natural log of the absolute value of the determinant of the matrix of which
00179       * *this is the QR decomposition. It has only linear complexity
00180       * (that is, O(n) where n is the dimension of the square matrix)
00181       * as the QR decomposition has already been computed.
00182       *
00183       * \note This is only for square matrices.
00184       *
00185       * \note This method is useful to work around the risk of overflow/underflow that's inherent
00186       * to determinant computation.
00187       *
00188       * \sa absDeterminant(), MatrixBase::determinant()
00189       */
00190     typename MatrixType::RealScalar logAbsDeterminant() const;
00191 
00192     /** \returns the rank of the matrix of which *this is the QR decomposition.
00193       *
00194       * \note This is computed at the time of the construction of the QR decomposition. This
00195       *       method does not perform any further computation.
00196       */
00197     inline Index rank() const
00198     {
00199       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00200       return m_rank;
00201     }
00202 
00203     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
00204       *
00205       * \note Since the rank is computed at the time of the construction of the QR decomposition, this
00206       *       method almost does not perform any further computation.
00207       */
00208     inline Index dimensionOfKernel() const
00209     {
00210       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00211       return m_qr.cols() - m_rank;
00212     }
00213 
00214     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
00215       *          linear map, i.e. has trivial kernel; false otherwise.
00216       *
00217       * \note Since the rank is computed at the time of the construction of the QR decomposition, this
00218       *       method almost does not perform any further computation.
00219       */
00220     inline bool isInjective() const
00221     {
00222       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00223       return m_rank == m_qr.cols();
00224     }
00225 
00226     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
00227       *          linear map; false otherwise.
00228       *
00229       * \note Since the rank is computed at the time of the construction of the QR decomposition, this
00230       *       method almost does not perform any further computation.
00231       */
00232     inline bool isSurjective() const
00233     {
00234       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00235       return m_rank == m_qr.rows();
00236     }
00237 
00238     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
00239       *
00240       * \note Since the rank is computed at the time of the construction of the QR decomposition, this
00241       *       method almost does not perform any further computation.
00242       */
00243     inline bool isInvertible() const
00244     {
00245       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00246       return isInjective() && isSurjective();
00247     }
00248 
00249     /** \returns the inverse of the matrix of which *this is the QR decomposition.
00250       *
00251       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
00252       *       Use isInvertible() to first determine whether this matrix is invertible.
00253       */    inline const
00254     internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
00255     inverse() const
00256     {
00257       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00258       return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
00259                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00260     }
00261 
00262     inline Index rows() const { return m_qr.rows(); }
00263     inline Index cols() const { return m_qr.cols(); }
00264     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00265 
00266   protected:
00267     MatrixType m_qr;
00268     HCoeffsType m_hCoeffs;
00269     IntColVectorType m_rows_transpositions;
00270     IntRowVectorType m_cols_transpositions;
00271     PermutationType m_cols_permutation;
00272     RowVectorType m_temp;
00273     bool m_isInitialized;
00274     RealScalar m_precision;
00275     Index m_rank;
00276     Index m_det_pq;
00277 };
00278 
00279 template<typename MatrixType>
00280 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00281 {
00282   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00283   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00284   return internal::abs(m_qr.diagonal().prod());
00285 }
00286 
00287 template<typename MatrixType>
00288 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00289 {
00290   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00291   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00292   return m_qr.diagonal().cwiseAbs().array().log().sum();
00293 }
00294 
00295 template<typename MatrixType>
00296 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00297 {
00298   Index rows = matrix.rows();
00299   Index cols = matrix.cols();
00300   Index size = std::min(rows,cols);
00301   m_rank = size;
00302 
00303   m_qr = matrix;
00304   m_hCoeffs.resize(size);
00305 
00306   m_temp.resize(cols);
00307 
00308   m_precision = NumTraits<Scalar>::epsilon() * size;
00309 
00310   m_rows_transpositions.resize(matrix.rows());
00311   m_cols_transpositions.resize(matrix.cols());
00312   Index number_of_transpositions = 0;
00313 
00314   RealScalar biggest(0);
00315 
00316   for (Index k = 0; k < size; ++k)
00317   {
00318     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00319     RealScalar biggest_in_corner;
00320 
00321     biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00322                             .cwiseAbs()
00323                             .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00324     row_of_biggest_in_corner += k;
00325     col_of_biggest_in_corner += k;
00326     if(k==0) biggest = biggest_in_corner;
00327 
00328     // if the corner is negligible, then we have less than full rank, and we can finish early
00329     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00330     {
00331       m_rank = k;
00332       for(Index i = k; i < size; i++)
00333       {
00334         m_rows_transpositions.coeffRef(i) = i;
00335         m_cols_transpositions.coeffRef(i) = i;
00336         m_hCoeffs.coeffRef(i) = Scalar(0);
00337       }
00338       break;
00339     }
00340 
00341     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00342     m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00343     if(k != row_of_biggest_in_corner) {
00344       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00345       ++number_of_transpositions;
00346     }
00347     if(k != col_of_biggest_in_corner) {
00348       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00349       ++number_of_transpositions;
00350     }
00351 
00352     RealScalar beta;
00353     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00354     m_qr.coeffRef(k,k) = beta;
00355 
00356     m_qr.bottomRightCorner(rows-k, cols-k-1)
00357         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00358   }
00359 
00360   m_cols_permutation.setIdentity(cols);
00361   for(Index k = 0; k < size; ++k)
00362     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00363 
00364   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00365   m_isInitialized = true;
00366 
00367   return *this;
00368 }
00369 
00370 namespace internal {
00371 
00372 template<typename _MatrixType, typename Rhs>
00373 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00374   : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00375 {
00376   EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00377 
00378   template<typename Dest> void evalTo(Dest& dst) const
00379   {
00380     const Index rows = dec().rows(), cols = dec().cols();
00381     eigen_assert(rhs().rows() == rows);
00382 
00383     // FIXME introduce nonzeroPivots() and use it here. and more generally,
00384     // make the same improvements in this dec as in FullPivLU.
00385     if(dec().rank()==0)
00386     {
00387       dst.setZero();
00388       return;
00389     }
00390 
00391     typename Rhs::PlainObject c(rhs());
00392 
00393     Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00394     for (Index k = 0; k < dec().rank(); ++k)
00395     {
00396       Index remainingSize = rows-k;
00397       c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00398       c.bottomRightCorner(remainingSize, rhs().cols())
00399        .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00400                                   dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00401     }
00402 
00403     if(!dec().isSurjective())
00404     {
00405       // is c is in the image of R ?
00406       RealScalar biggest_in_upper_part_of_c = c.topRows(   dec().rank()     ).cwiseAbs().maxCoeff();
00407       RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00408       // FIXME brain dead
00409       const RealScalar m_precision = NumTraits<Scalar>::epsilon() * std::min(rows,cols);
00410       // this internal:: prefix is needed by at least gcc 3.4 and ICC
00411       if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00412         return;
00413     }
00414     dec().matrixQR()
00415        .topLeftCorner(dec().rank(), dec().rank())
00416        .template triangularView<Upper>()
00417        .solveInPlace(c.topRows(dec().rank()));
00418 
00419     for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00420     for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00421   }
00422 };
00423 
00424 } // end namespace internal
00425 
00426 /** \returns the matrix Q */
00427 template<typename MatrixType>
00428 typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
00429 {
00430   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00431   // compute the product H'_0 H'_1 ... H'_n-1,
00432   // where H_k is the k-th Householder transformation I - h_k v_k v_k'
00433   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
00434   Index rows = m_qr.rows();
00435   Index cols = m_qr.cols();
00436   Index size = std::min(rows,cols);
00437   MatrixQType res = MatrixQType::Identity(rows, rows);
00438   Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
00439   for (Index k = size-1; k >= 0; k--)
00440   {
00441     res.block(k, k, rows-k, rows-k)
00442        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
00443     res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
00444   }
00445   return res;
00446 }
00447 
00448 /** \return the full-pivoting Householder QR decomposition of \c *this.
00449   *
00450   * \sa class FullPivHouseholderQR
00451   */
00452 template<typename Derived>
00453 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00454 MatrixBase<Derived>::fullPivHouseholderQr() const
00455 {
00456   return FullPivHouseholderQR<PlainObject>(eval());
00457 }
00458 
00459 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011