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 namespace Eigen { 
00031 
00057 template<typename _MatrixType> class HouseholderQR
00058 {
00059   public:
00060 
00061     typedef _MatrixType MatrixType;
00062     enum {
00063       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065       Options = MatrixType::Options,
00066       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00067       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00068     };
00069     typedef typename MatrixType::Scalar Scalar;
00070     typedef typename MatrixType::RealScalar RealScalar;
00071     typedef typename MatrixType::Index Index;
00072     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00073     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00074     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00075     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00076 
00083     HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
00084 
00091     HouseholderQR(Index rows, Index cols)
00092       : m_qr(rows, cols),
00093         m_hCoeffs((std::min)(rows,cols)),
00094         m_temp(cols),
00095         m_isInitialized(false) {}
00096 
00097     HouseholderQR(const MatrixType& matrix)
00098       : m_qr(matrix.rows(), matrix.cols()),
00099         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00100         m_temp(matrix.cols()),
00101         m_isInitialized(false)
00102     {
00103       compute(matrix);
00104     }
00105 
00123     template<typename Rhs>
00124     inline const internal::solve_retval<HouseholderQR, Rhs>
00125     solve(const MatrixBase<Rhs>& b) const
00126     {
00127       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00128       return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
00129     }
00130 
00131     HouseholderSequenceType householderQ() const
00132     {
00133       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00134       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00135     }
00136 
00140     const MatrixType& matrixQR() const
00141     {
00142         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00143         return m_qr;
00144     }
00145 
00146     HouseholderQR& compute(const MatrixType& matrix);
00147 
00161     typename MatrixType::RealScalar absDeterminant() const;
00162 
00175     typename MatrixType::RealScalar logAbsDeterminant() const;
00176 
00177     inline Index rows() const { return m_qr.rows(); }
00178     inline Index cols() const { return m_qr.cols(); }
00179     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00180 
00181   protected:
00182     MatrixType m_qr;
00183     HCoeffsType m_hCoeffs;
00184     RowVectorType m_temp;
00185     bool m_isInitialized;
00186 };
00187 
00188 template<typename MatrixType>
00189 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00190 {
00191   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00192   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00193   return internal::abs(m_qr.diagonal().prod());
00194 }
00195 
00196 template<typename MatrixType>
00197 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00198 {
00199   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00200   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00201   return m_qr.diagonal().cwiseAbs().array().log().sum();
00202 }
00203 
00204 namespace internal {
00205 
00207 template<typename MatrixQR, typename HCoeffs>
00208 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00209 {
00210   typedef typename MatrixQR::Index Index;
00211   typedef typename MatrixQR::Scalar Scalar;
00212   typedef typename MatrixQR::RealScalar RealScalar;
00213   Index rows = mat.rows();
00214   Index cols = mat.cols();
00215   Index size = (std::min)(rows,cols);
00216 
00217   eigen_assert(hCoeffs.size() == size);
00218 
00219   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00220   TempType tempVector;
00221   if(tempData==0)
00222   {
00223     tempVector.resize(cols);
00224     tempData = tempVector.data();
00225   }
00226 
00227   for(Index k = 0; k < size; ++k)
00228   {
00229     Index remainingRows = rows - k;
00230     Index remainingCols = cols - k - 1;
00231 
00232     RealScalar beta;
00233     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00234     mat.coeffRef(k,k) = beta;
00235 
00236     // apply H to remaining part of m_qr from the left
00237     mat.bottomRightCorner(remainingRows, remainingCols)
00238         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00239   }
00240 }
00241 
00243 template<typename MatrixQR, typename HCoeffs>
00244 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
00245                                        typename MatrixQR::Index maxBlockSize=32,
00246                                        typename MatrixQR::Scalar* tempData = 0)
00247 {
00248   typedef typename MatrixQR::Index Index;
00249   typedef typename MatrixQR::Scalar Scalar;
00250   typedef typename MatrixQR::RealScalar RealScalar;
00251   typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00252 
00253   Index rows = mat.rows();
00254   Index cols = mat.cols();
00255   Index size = (std::min)(rows, cols);
00256 
00257   typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00258   TempType tempVector;
00259   if(tempData==0)
00260   {
00261     tempVector.resize(cols);
00262     tempData = tempVector.data();
00263   }
00264 
00265   Index blockSize = (std::min)(maxBlockSize,size);
00266 
00267   Index k = 0;
00268   for (k = 0; k < size; k += blockSize)
00269   {
00270     Index bs = (std::min)(size-k,blockSize);  // actual size of the block
00271     Index tcols = cols - k - bs;            // trailing columns
00272     Index brows = rows-k;                   // rows of the block
00273 
00274     // partition the matrix:
00275     //        A00 | A01 | A02
00276     // mat  = A10 | A11 | A12
00277     //        A20 | A21 | A22
00278     // and performs the qr dec of [A11^T A12^T]^T
00279     // and update [A21^T A22^T]^T using level 3 operations.
00280     // Finally, the algorithm continue on A22
00281 
00282     BlockType A11_21 = mat.block(k,k,brows,bs);
00283     Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00284 
00285     householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00286 
00287     if(tcols)
00288     {
00289       BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00290       apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
00291     }
00292   }
00293 }
00294 
00295 template<typename _MatrixType, typename Rhs>
00296 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
00297   : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
00298 {
00299   EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
00300 
00301   template<typename Dest> void evalTo(Dest& dst) const
00302   {
00303     const Index rows = dec().rows(), cols = dec().cols();
00304     const Index rank = (std::min)(rows, cols);
00305     eigen_assert(rhs().rows() == rows);
00306 
00307     typename Rhs::PlainObject c(rhs());
00308 
00309     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00310     c.applyOnTheLeft(householderSequence(
00311       dec().matrixQR().leftCols(rank),
00312       dec().hCoeffs().head(rank)).transpose()
00313     );
00314 
00315     dec().matrixQR()
00316        .topLeftCorner(rank, rank)
00317        .template triangularView<Upper>()
00318        .solveInPlace(c.topRows(rank));
00319 
00320     dst.topRows(rank) = c.topRows(rank);
00321     dst.bottomRows(cols-rank).setZero();
00322   }
00323 };
00324 
00325 } // end namespace internal
00326 
00327 template<typename MatrixType>
00328 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00329 {
00330   Index rows = matrix.rows();
00331   Index cols = matrix.cols();
00332   Index size = (std::min)(rows,cols);
00333 
00334   m_qr = matrix;
00335   m_hCoeffs.resize(size);
00336 
00337   m_temp.resize(cols);
00338 
00339   internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
00340 
00341   m_isInitialized = true;
00342   return *this;
00343 }
00344 
00349 template<typename Derived>
00350 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00351 MatrixBase<Derived>::householderQr() const
00352 {
00353   return HouseholderQR<PlainObject>(eval());
00354 }
00355 
00356 } // end namespace Eigen
00357 
00358 #endif // EIGEN_QR_H