Main MRPT website > C++ reference
MRPT logo

HouseholderQR.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-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 // Copyright (C) 2010 Vincent Lejeune
00007 //
00008 // Eigen is free software; you can redistribute it and/or
00009 // modify it under the terms of the GNU Lesser General Public
00010 // License as published by the Free Software Foundation; either
00011 // version 3 of the License, or (at your option) any later version.
00012 //
00013 // Alternatively, you can redistribute it and/or
00014 // modify it under the terms of the GNU General Public License as
00015 // published by the Free Software Foundation; either version 2 of
00016 // the License, or (at your option) any later version.
00017 //
00018 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00019 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00020 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00021 // GNU General Public License for more details.
00022 //
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License and a copy of the GNU General Public License along with
00025 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 #ifndef EIGEN_QR_H
00028 #define EIGEN_QR_H
00029 
00030 /** \ingroup QR_Module
00031   *
00032   *
00033   * \class HouseholderQR
00034   *
00035   * \brief Householder QR decomposition of a matrix
00036   *
00037   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
00038   *
00039   * This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R
00040   * such that 
00041   * \f[
00042   *  \mathbf{A} = \mathbf{Q} \, \mathbf{R}
00043   * \f]
00044   * by using Householder transformations. Here, \b Q a unitary matrix and \b R an upper triangular matrix.
00045   * The result is stored in a compact way compatible with LAPACK.
00046   *
00047   * Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
00048   * If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead.
00049   *
00050   * This Householder QR decomposition is faster, but less numerically stable and less feature-full than
00051   * FullPivHouseholderQR or ColPivHouseholderQR.
00052   *
00053   * \sa MatrixBase::householderQr()
00054   */
00055 template<typename _MatrixType> class HouseholderQR
00056 {
00057   public:
00058 
00059     typedef _MatrixType MatrixType;
00060     enum {
00061       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00062       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00063       Options = MatrixType::Options,
00064       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00065       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00066     };
00067     typedef typename MatrixType::Scalar Scalar;
00068     typedef typename MatrixType::RealScalar RealScalar;
00069     typedef typename MatrixType::Index Index;
00070     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00071     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00072     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00073     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00074 
00075     /**
00076     * \brief Default Constructor.
00077     *
00078     * The default constructor is useful in cases in which the user intends to
00079     * perform decompositions via HouseholderQR::compute(const MatrixType&).
00080     */
00081     HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
00082 
00083     /** \brief Default Constructor with memory preallocation
00084       *
00085       * Like the default constructor but with preallocation of the internal data
00086       * according to the specified problem \a size.
00087       * \sa HouseholderQR()
00088       */
00089     HouseholderQR(Index rows, Index cols)
00090       : m_qr(rows, cols),
00091         m_hCoeffs(std::min(rows,cols)),
00092         m_temp(cols),
00093         m_isInitialized(false) {}
00094 
00095     HouseholderQR(const MatrixType& matrix)
00096       : m_qr(matrix.rows(), matrix.cols()),
00097         m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
00098         m_temp(matrix.cols()),
00099         m_isInitialized(false)
00100     {
00101       compute(matrix);
00102     }
00103 
00104     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
00105       * *this is the QR decomposition, if any exists.
00106       *
00107       * \param b the right-hand-side of the equation to solve.
00108       *
00109       * \returns a solution.
00110       *
00111       * \note The case where b is a matrix is not yet implemented. Also, this
00112       *       code is space inefficient.
00113       *
00114       * \note_about_checking_solutions
00115       *
00116       * \note_about_arbitrary_choice_of_solution
00117       *
00118       * Example: \include HouseholderQR_solve.cpp
00119       * Output: \verbinclude HouseholderQR_solve.out
00120       */
00121     template<typename Rhs>
00122     inline const internal::solve_retval<HouseholderQR, Rhs>
00123     solve(const MatrixBase<Rhs>& b) const
00124     {
00125       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00126       return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
00127     }
00128 
00129     HouseholderSequenceType householderQ() const
00130     {
00131       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00132       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00133     }
00134 
00135     /** \returns a reference to the matrix where the Householder QR decomposition is stored
00136       * in a LAPACK-compatible way.
00137       */
00138     const MatrixType& matrixQR() const
00139     {
00140         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00141         return m_qr;
00142     }
00143 
00144     HouseholderQR& compute(const MatrixType& matrix);
00145 
00146     /** \returns the absolute value of the determinant of the matrix of which
00147       * *this is the QR decomposition. It has only linear complexity
00148       * (that is, O(n) where n is the dimension of the square matrix)
00149       * as the QR decomposition has already been computed.
00150       *
00151       * \note This is only for square matrices.
00152       *
00153       * \warning a determinant can be very big or small, so for matrices
00154       * of large enough dimension, there is a risk of overflow/underflow.
00155       * One way to work around that is to use logAbsDeterminant() instead.
00156       *
00157       * \sa logAbsDeterminant(), MatrixBase::determinant()
00158       */
00159     typename MatrixType::RealScalar absDeterminant() const;
00160 
00161     /** \returns the natural log of the absolute value of the determinant of the matrix of which
00162       * *this is the QR decomposition. It has only linear complexity
00163       * (that is, O(n) where n is the dimension of the square matrix)
00164       * as the QR decomposition has already been computed.
00165       *
00166       * \note This is only for square matrices.
00167       *
00168       * \note This method is useful to work around the risk of overflow/underflow that's inherent
00169       * to determinant computation.
00170       *
00171       * \sa absDeterminant(), MatrixBase::determinant()
00172       */
00173     typename MatrixType::RealScalar logAbsDeterminant() const;
00174 
00175     inline Index rows() const { return m_qr.rows(); }
00176     inline Index cols() const { return m_qr.cols(); }
00177     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00178 
00179   protected:
00180     MatrixType m_qr;
00181     HCoeffsType m_hCoeffs;
00182     RowVectorType m_temp;
00183     bool m_isInitialized;
00184 };
00185 
00186 template<typename MatrixType>
00187 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00188 {
00189   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00190   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00191   return internal::abs(m_qr.diagonal().prod());
00192 }
00193 
00194 template<typename MatrixType>
00195 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00196 {
00197   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00198   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00199   return m_qr.diagonal().cwiseAbs().array().log().sum();
00200 }
00201 
00202 namespace internal {
00203 
00204 /** \internal */
00205 template<typename MatrixQR, typename HCoeffs>
00206 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00207 {
00208   typedef typename MatrixQR::Index Index;
00209   typedef typename MatrixQR::Scalar Scalar;
00210   typedef typename MatrixQR::RealScalar RealScalar;
00211   Index rows = mat.rows();
00212   Index cols = mat.cols();
00213   Index size = std::min(rows,cols);
00214 
00215   eigen_assert(hCoeffs.size() == size);
00216 
00217   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00218   TempType tempVector;
00219   if(tempData==0)
00220   {
00221     tempVector.resize(cols);
00222     tempData = tempVector.data();
00223   }
00224 
00225   for(Index k = 0; k < size; ++k)
00226   {
00227     Index remainingRows = rows - k;
00228     Index remainingCols = cols - k - 1;
00229 
00230     RealScalar beta;
00231     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00232     mat.coeffRef(k,k) = beta;
00233 
00234     // apply H to remaining part of m_qr from the left
00235     mat.bottomRightCorner(remainingRows, remainingCols)
00236         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00237   }
00238 }
00239 
00240 /** \internal */
00241 template<typename MatrixQR, typename HCoeffs>
00242 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
00243                                        typename MatrixQR::Index maxBlockSize=32,
00244                                        typename MatrixQR::Scalar* tempData = 0)
00245 {
00246   typedef typename MatrixQR::Index Index;
00247   typedef typename MatrixQR::Scalar Scalar;
00248   typedef typename MatrixQR::RealScalar RealScalar;
00249   typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00250 
00251   Index rows = mat.rows();
00252   Index cols = mat.cols();
00253   Index size = std::min(rows, cols);
00254 
00255   typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00256   TempType tempVector;
00257   if(tempData==0)
00258   {
00259     tempVector.resize(cols);
00260     tempData = tempVector.data();
00261   }
00262 
00263   Index blockSize = std::min(maxBlockSize,size);
00264 
00265   int k = 0;
00266   for (k = 0; k < size; k += blockSize)
00267   {
00268     Index bs = std::min(size-k,blockSize);  // actual size of the block
00269     Index tcols = cols - k - bs;            // trailing columns
00270     Index brows = rows-k;                   // rows of the block
00271 
00272     // partition the matrix:
00273     //        A00 | A01 | A02
00274     // mat  = A10 | A11 | A12
00275     //        A20 | A21 | A22
00276     // and performs the qr dec of [A11^T A12^T]^T
00277     // and update [A21^T A22^T]^T using level 3 operations.
00278     // Finally, the algorithm continue on A22
00279 
00280     BlockType A11_21 = mat.block(k,k,brows,bs);
00281     Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00282 
00283     householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00284 
00285     if(tcols)
00286     {
00287       BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00288       apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
00289     }
00290   }
00291 }
00292 
00293 template<typename _MatrixType, typename Rhs>
00294 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
00295   : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
00296 {
00297   EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
00298 
00299   template<typename Dest> void evalTo(Dest& dst) const
00300   {
00301     const Index rows = dec().rows(), cols = dec().cols();
00302     const Index rank = std::min(rows, cols);
00303     eigen_assert(rhs().rows() == rows);
00304 
00305     typename Rhs::PlainObject c(rhs());
00306 
00307     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00308     c.applyOnTheLeft(householderSequence(
00309       dec().matrixQR().leftCols(rank),
00310       dec().hCoeffs().head(rank)).transpose()
00311     );
00312 
00313     dec().matrixQR()
00314        .topLeftCorner(rank, rank)
00315        .template triangularView<Upper>()
00316        .solveInPlace(c.topRows(rank));
00317 
00318     dst.topRows(rank) = c.topRows(rank);
00319     dst.bottomRows(cols-rank).setZero();
00320   }
00321 };
00322 
00323 } // end namespace internal
00324 
00325 template<typename MatrixType>
00326 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00327 {
00328   Index rows = matrix.rows();
00329   Index cols = matrix.cols();
00330   Index size = std::min(rows,cols);
00331 
00332   m_qr = matrix;
00333   m_hCoeffs.resize(size);
00334 
00335   m_temp.resize(cols);
00336 
00337   internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
00338 
00339   m_isInitialized = true;
00340   return *this;
00341 }
00342 
00343 /** \return the Householder QR decomposition of \c *this.
00344   *
00345   * \sa class HouseholderQR
00346   */
00347 template<typename Derived>
00348 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00349 MatrixBase<Derived>::householderQr() const
00350 {
00351   return HouseholderQR<PlainObject>(eval());
00352 }
00353 
00354 
00355 #endif // EIGEN_QR_H



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