src/QR/HouseholderQR.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0006 // Copyright (C) 2010 Vincent Lejeune
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012 #ifndef EIGEN_QR_H
0013 #define EIGEN_QR_H
0014
0015 namespace Eigen {
0016
0017 namespace internal {
0018 template<typename _MatrixType> struct traits<HouseholderQR<_MatrixType> >
0019  : traits<_MatrixType>
0020 {
0021   typedef MatrixXpr XprKind;
0022   typedef SolverStorage StorageKind;
0023   typedef int StorageIndex;
0024   enum { Flags = 0 };
0025 };
0026
0027 } // end namespace internal
0028
0029 /** \ingroup QR_Module
0030   *
0031   *
0032   * \class HouseholderQR
0033   *
0034   * \brief Householder QR decomposition of a matrix
0035   *
0036   * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
0037   *
0038   * This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R
0039   * such that
0040   * \f[
0041   *  \mathbf{A} = \mathbf{Q} \, \mathbf{R}
0042   * \f]
0043   * by using Householder transformations. Here, \b Q a unitary matrix and \b R an upper triangular matrix.
0044   * The result is stored in a compact way compatible with LAPACK.
0045   *
0046   * Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
0047   * If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead.
0048   *
0049   * This Householder QR decomposition is faster, but less numerically stable and less feature-full than
0050   * FullPivHouseholderQR or ColPivHouseholderQR.
0051   *
0052   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0053   *
0054   * \sa MatrixBase::householderQr()
0055   */
0056 template<typename _MatrixType> class HouseholderQR
0057         : public SolverBase<HouseholderQR<_MatrixType> >
0058 {
0059   public:
0060
0061     typedef _MatrixType MatrixType;
0062     typedef SolverBase<HouseholderQR> Base;
0063     friend class SolverBase<HouseholderQR>;
0064
0065     EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR)
0066     enum {
0067       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0068       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0069     };
0070     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
0071     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
0072     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
0073     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
0074
0075     /**
0076       * \brief Default Constructor.
0077       *
0078       * The default constructor is useful in cases in which the user intends to
0079       * perform decompositions via HouseholderQR::compute(const MatrixType&).
0080       */
0081     HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
0082
0083     /** \brief Default Constructor with memory preallocation
0084       *
0085       * Like the default constructor but with preallocation of the internal data
0086       * according to the specified problem \a size.
0087       * \sa HouseholderQR()
0088       */
0089     HouseholderQR(Index rows, Index cols)
0090       : m_qr(rows, cols),
0091         m_hCoeffs((std::min)(rows,cols)),
0092         m_temp(cols),
0093         m_isInitialized(false) {}
0094
0095     /** \brief Constructs a QR factorization from a given matrix
0096       *
0097       * This constructor computes the QR factorization of the matrix \a matrix by calling
0098       * the method compute(). It is a short cut for:
0099       *
0100       * \code
0101       * HouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
0102       * qr.compute(matrix);
0103       * \endcode
0104       *
0105       * \sa compute()
0106       */
0107     template<typename InputType>
0108     explicit HouseholderQR(const EigenBase<InputType>& matrix)
0109       : m_qr(matrix.rows(), matrix.cols()),
0110         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
0111         m_temp(matrix.cols()),
0112         m_isInitialized(false)
0113     {
0114       compute(matrix.derived());
0115     }
0116
0117
0118     /** \brief Constructs a QR factorization from a given matrix
0119       *
0120       * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
0121       * \c MatrixType is a Eigen::Ref.
0122       *
0123       * \sa HouseholderQR(const EigenBase&)
0124       */
0125     template<typename InputType>
0126     explicit HouseholderQR(EigenBase<InputType>& matrix)
0127       : m_qr(matrix.derived()),
0128         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
0129         m_temp(matrix.cols()),
0130         m_isInitialized(false)
0131     {
0132       computeInPlace();
0133     }
0134
0135     #ifdef EIGEN_PARSED_BY_DOXYGEN
0136     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
0137       * *this is the QR decomposition, if any exists.
0138       *
0139       * \param b the right-hand-side of the equation to solve.
0140       *
0141       * \returns a solution.
0142       *
0143       * \note_about_checking_solutions
0144       *
0145       * \note_about_arbitrary_choice_of_solution
0146       *
0147       * Example: \include HouseholderQR_solve.cpp
0148       * Output: \verbinclude HouseholderQR_solve.out
0149       */
0150     template<typename Rhs>
0151     inline const Solve<HouseholderQR, Rhs>
0152     solve(const MatrixBase<Rhs>& b) const;
0153     #endif
0154
0155     /** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations.
0156       *
0157       * The returned expression can directly be used to perform matrix products. It can also be assigned to a dense Matrix object.
0158       * Here is an example showing how to recover the full or thin matrix Q, as well as how to perform matrix products using operator*:
0159       *
0160       * Example: \include HouseholderQR_householderQ.cpp
0161       * Output: \verbinclude HouseholderQR_householderQ.out
0162       */
0163     HouseholderSequenceType householderQ() const
0164     {
0165       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
0166       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
0167     }
0168
0169     /** \returns a reference to the matrix where the Householder QR decomposition is stored
0170       * in a LAPACK-compatible way.
0171       */
0172     const MatrixType& matrixQR() const
0173     {
0174         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
0175         return m_qr;
0176     }
0177
0178     template<typename InputType>
0179     HouseholderQR& compute(const EigenBase<InputType>& matrix) {
0180       m_qr = matrix.derived();
0181       computeInPlace();
0182       return *this;
0183     }
0184
0185     /** \returns the absolute value of the determinant of the matrix of which
0186       * *this is the QR decomposition. It has only linear complexity
0187       * (that is, O(n) where n is the dimension of the square matrix)
0188       * as the QR decomposition has already been computed.
0189       *
0190       * \note This is only for square matrices.
0191       *
0192       * \warning a determinant can be very big or small, so for matrices
0193       * of large enough dimension, there is a risk of overflow/underflow.
0194       * One way to work around that is to use logAbsDeterminant() instead.
0195       *
0196       * \sa logAbsDeterminant(), MatrixBase::determinant()
0197       */
0198     typename MatrixType::RealScalar absDeterminant() const;
0199
0200     /** \returns the natural log of the absolute value of the determinant of the matrix of which
0201       * *this is the QR decomposition. It has only linear complexity
0202       * (that is, O(n) where n is the dimension of the square matrix)
0203       * as the QR decomposition has already been computed.
0204       *
0205       * \note This is only for square matrices.
0206       *
0207       * \note This method is useful to work around the risk of overflow/underflow that's inherent
0208       * to determinant computation.
0209       *
0210       * \sa absDeterminant(), MatrixBase::determinant()
0211       */
0212     typename MatrixType::RealScalar logAbsDeterminant() const;
0213
0214     inline Index rows() const { return m_qr.rows(); }
0215     inline Index cols() const { return m_qr.cols(); }
0216
0217     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
0218       *
0219       * For advanced uses only.
0220       */
0221     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
0222
0223     #ifndef EIGEN_PARSED_BY_DOXYGEN
0224     template<typename RhsType, typename DstType>
0225     void _solve_impl(const RhsType &rhs, DstType &dst) const;
0226
0227     template<bool Conjugate, typename RhsType, typename DstType>
0228     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
0229     #endif
0230
0231   protected:
0232
0233     static void check_template_parameters()
0234     {
0235       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0236     }
0237
0238     void computeInPlace();
0239
0240     MatrixType m_qr;
0241     HCoeffsType m_hCoeffs;
0242     RowVectorType m_temp;
0243     bool m_isInitialized;
0244 };
0245
0246 template<typename MatrixType>
0247 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
0248 {
0249   using std::abs;
0250   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
0251   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0252   return abs(m_qr.diagonal().prod());
0253 }
0254
0255 template<typename MatrixType>
0256 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
0257 {
0258   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
0259   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0260   return m_qr.diagonal().cwiseAbs().array().log().sum();
0261 }
0262
0263 namespace internal {
0264
0265 /** \internal */
0266 template<typename MatrixQR, typename HCoeffs>
0267 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
0268 {
0269   typedef typename MatrixQR::Scalar Scalar;
0270   typedef typename MatrixQR::RealScalar RealScalar;
0271   Index rows = mat.rows();
0272   Index cols = mat.cols();
0273   Index size = (std::min)(rows,cols);
0274
0275   eigen_assert(hCoeffs.size() == size);
0276
0277   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
0278   TempType tempVector;
0279   if(tempData==0)
0280   {
0281     tempVector.resize(cols);
0282     tempData = tempVector.data();
0283   }
0284
0285   for(Index k = 0; k < size; ++k)
0286   {
0287     Index remainingRows = rows - k;
0288     Index remainingCols = cols - k - 1;
0289
0290     RealScalar beta;
0291     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
0292     mat.coeffRef(k,k) = beta;
0293
0294     // apply H to remaining part of m_qr from the left
0295     mat.bottomRightCorner(remainingRows, remainingCols)
0296         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
0297   }
0298 }
0299
0300 /** \internal */
0301 template<typename MatrixQR, typename HCoeffs,
0302   typename MatrixQRScalar = typename MatrixQR::Scalar,
0303   bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
0304 struct householder_qr_inplace_blocked
0305 {
0306   // This is specialized for LAPACK-supported Scalar types in HouseholderQR_LAPACKE.h
0307   static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
0308       typename MatrixQR::Scalar* tempData = 0)
0309   {
0310     typedef typename MatrixQR::Scalar Scalar;
0311     typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
0312
0313     Index rows = mat.rows();
0314     Index cols = mat.cols();
0315     Index size = (std::min)(rows, cols);
0316
0317     typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
0318     TempType tempVector;
0319     if(tempData==0)
0320     {
0321       tempVector.resize(cols);
0322       tempData = tempVector.data();
0323     }
0324
0325     Index blockSize = (std::min)(maxBlockSize,size);
0326
0327     Index k = 0;
0328     for (k = 0; k < size; k += blockSize)
0329     {
0330       Index bs = (std::min)(size-k,blockSize);  // actual size of the block
0331       Index tcols = cols - k - bs;              // trailing columns
0332       Index brows = rows-k;                     // rows of the block
0333
0334       // partition the matrix:
0335       //        A00 | A01 | A02
0336       // mat  = A10 | A11 | A12
0337       //        A20 | A21 | A22
0338       // and performs the qr dec of [A11^T A12^T]^T
0339       // and update [A21^T A22^T]^T using level 3 operations.
0340       // Finally, the algorithm continue on A22
0341
0342       BlockType A11_21 = mat.block(k,k,brows,bs);
0343       Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
0344
0345       householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
0346
0347       if(tcols)
0348       {
0349         BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
0350         apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
0351       }
0352     }
0353   }
0354 };
0355
0356 } // end namespace internal
0357
0358 #ifndef EIGEN_PARSED_BY_DOXYGEN
0359 template<typename _MatrixType>
0360 template<typename RhsType, typename DstType>
0361 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
0362 {
0363   const Index rank = (std::min)(rows(), cols());
0364
0365   typename RhsType::PlainObject c(rhs);
0366
0367   c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
0368
0369   m_qr.topLeftCorner(rank, rank)
0370       .template triangularView<Upper>()
0371       .solveInPlace(c.topRows(rank));
0372
0373   dst.topRows(rank) = c.topRows(rank);
0374   dst.bottomRows(cols()-rank).setZero();
0375 }
0376
0377 template<typename _MatrixType>
0378 template<bool Conjugate, typename RhsType, typename DstType>
0379 void HouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
0380 {
0381   const Index rank = (std::min)(rows(), cols());
0382
0383   typename RhsType::PlainObject c(rhs);
0384
0385   m_qr.topLeftCorner(rank, rank)
0386       .template triangularView<Upper>()
0387       .transpose().template conjugateIf<Conjugate>()
0388       .solveInPlace(c.topRows(rank));
0389
0390   dst.topRows(rank) = c.topRows(rank);
0391   dst.bottomRows(rows()-rank).setZero();
0392
0393   dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
0394 }
0395 #endif
0396
0397 /** Performs the QR factorization of the given matrix \a matrix. The result of
0398   * the factorization is stored into \c *this, and a reference to \c *this
0399   * is returned.
0400   *
0401   * \sa class HouseholderQR, HouseholderQR(const MatrixType&)
0402   */
0403 template<typename MatrixType>
0404 void HouseholderQR<MatrixType>::computeInPlace()
0405 {
0406   check_template_parameters();
0407
0408   Index rows = m_qr.rows();
0409   Index cols = m_qr.cols();
0410   Index size = (std::min)(rows,cols);
0411
0412   m_hCoeffs.resize(size);
0413
0414   m_temp.resize(cols);
0415
0416   internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
0417
0418   m_isInitialized = true;
0419 }
0420
0421 /** \return the Householder QR decomposition of \c *this.
0422   *
0423   * \sa class HouseholderQR
0424   */
0425 template<typename Derived>
0426 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
0427 MatrixBase<Derived>::householderQr() const
0428 {
0429   return HouseholderQR<PlainObject>(eval());
0430 }
0431
0432 } // end namespace Eigen
0433
0434 #endif // EIGEN_QR_H