HouseholderSequence.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
00012 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
00013 
00014 namespace Eigen { 
00015 
00016 /** \ingroup Householder_Module
00017   * \householder_module
00018   * \class HouseholderSequence
00019   * \brief Sequence of Householder reflections acting on subspaces with decreasing size
00020   * \tparam VectorsType type of matrix containing the Householder vectors
00021   * \tparam CoeffsType  type of vector containing the Householder coefficients
00022   * \tparam Side        either OnTheLeft (the default) or OnTheRight
00023   *
00024   * This class represents a product sequence of Householder reflections where the first Householder reflection
00025   * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by
00026   * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace
00027   * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but
00028   * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections
00029   * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods
00030   * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(),
00031   * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence.
00032   *
00033   * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the
00034   * form \f$ H = \prod_{i=0}^{n-1} H_i \f$ where the i-th Householder reflection is \f$ H_i = I - h_i v_i
00035   * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$
00036   * v_i \f$ is a vector of the form
00037   * \f[ 
00038   * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
00039   * \f]
00040   * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector.
00041   *
00042   * Typical usages are listed below, where H is a HouseholderSequence:
00043   * \code
00044   * A.applyOnTheRight(H);             // A = A * H
00045   * A.applyOnTheLeft(H);              // A = H * A
00046   * A.applyOnTheRight(H.adjoint());   // A = A * H^*
00047   * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
00048   * MatrixXd Q = H;                   // conversion to a dense matrix
00049   * \endcode
00050   * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators.
00051   *
00052   * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example.
00053   *
00054   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
00055   */
00056 
00057 namespace internal {
00058 
00059 template<typename VectorsType, typename CoeffsType, int Side>
00060 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
00061 {
00062   typedef typename VectorsType::Scalar Scalar;
00063   typedef typename VectorsType::Index Index;
00064   typedef typename VectorsType::StorageKind StorageKind;
00065   enum {
00066     RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
00067                                         : traits<VectorsType>::ColsAtCompileTime,
00068     ColsAtCompileTime = RowsAtCompileTime,
00069     MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
00070                                            : traits<VectorsType>::MaxColsAtCompileTime,
00071     MaxColsAtCompileTime = MaxRowsAtCompileTime,
00072     Flags = 0
00073   };
00074 };
00075 
00076 template<typename VectorsType, typename CoeffsType, int Side>
00077 struct hseq_side_dependent_impl
00078 {
00079   typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
00080   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
00081   typedef typename VectorsType::Index Index;
00082   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00083   {
00084     Index start = k+1+h.m_shift;
00085     return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
00086   }
00087 };
00088 
00089 template<typename VectorsType, typename CoeffsType>
00090 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
00091 {
00092   typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
00093   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
00094   typedef typename VectorsType::Index Index;
00095   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
00096   {
00097     Index start = k+1+h.m_shift;
00098     return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
00099   }
00100 };
00101 
00102 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
00103 {
00104   typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
00105     ResultScalar;
00106   typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
00107                  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
00108 };
00109 
00110 } // end namespace internal
00111 
00112 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence 
00113   : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
00114 {
00115     typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType;
00116   
00117   public:
00118     enum {
00119       RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
00120       ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
00121       MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
00122       MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
00123     };
00124     typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
00125     typedef typename VectorsType::Index Index;
00126 
00127     typedef HouseholderSequence <
00128       typename internal::conditional<NumTraits<Scalar>::IsComplex,
00129         typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
00130         VectorsType>::type,
00131       typename internal::conditional<NumTraits<Scalar>::IsComplex,
00132         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
00133         CoeffsType>::type,
00134       Side
00135     > ConjugateReturnType ;
00136 
00137     /** \brief Constructor.
00138       * \param[in]  v      %Matrix containing the essential parts of the Householder vectors
00139       * \param[in]  h      Vector containing the Householder coefficients
00140       *
00141       * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The
00142       * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th
00143       * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the
00144       * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many
00145       * Householder reflections as there are columns.
00146       *
00147       * \note The %HouseholderSequence object stores \p v and \p h by reference.
00148       *
00149       * Example: \include HouseholderSequence_HouseholderSequence.cpp
00150       * Output: \verbinclude HouseholderSequence_HouseholderSequence.out
00151       *
00152       * \sa setLength(), setShift()
00153       */
00154     HouseholderSequence (const VectorsType& v, const CoeffsType& h)
00155       : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
00156         m_shift(0)
00157     {
00158     }
00159 
00160     /** \brief Copy constructor. */
00161     HouseholderSequence (const HouseholderSequence & other)
00162       : m_vectors(other.m_vectors),
00163         m_coeffs(other.m_coeffs),
00164         m_trans(other.m_trans),
00165         m_length(other.m_length),
00166         m_shift(other.m_shift)
00167     {
00168     }
00169 
00170     /** \brief Number of rows of transformation viewed as a matrix.
00171       * \returns Number of rows 
00172       * \details This equals the dimension of the space that the transformation acts on.
00173       */
00174     Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
00175 
00176     /** \brief Number of columns of transformation viewed as a matrix.
00177       * \returns Number of columns
00178       * \details This equals the dimension of the space that the transformation acts on.
00179       */
00180     Index cols() const { return rows(); }
00181 
00182     /** \brief Essential part of a Householder vector.
00183       * \param[in]  k  Index of Householder reflection
00184       * \returns    Vector containing non-trivial entries of k-th Householder vector
00185       *
00186       * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of
00187       * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector
00188       * \f[ 
00189       * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. 
00190       * \f]
00191       * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v
00192       * passed to the constructor.
00193       *
00194       * \sa setShift(), shift()
00195       */
00196     const EssentialVectorType essentialVector(Index k) const
00197     {
00198       eigen_assert(k >= 0 && k < m_length);
00199       return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
00200     }
00201 
00202     /** \brief %Transpose of the Householder sequence. */
00203     HouseholderSequence  transpose() const
00204     {
00205       return HouseholderSequence (*this).setTrans(!m_trans);
00206     }
00207 
00208     /** \brief Complex conjugate of the Householder sequence. */
00209     ConjugateReturnType  conjugate() const
00210     {
00211       return ConjugateReturnType (m_vectors.conjugate(), m_coeffs.conjugate())
00212              .setTrans(m_trans)
00213              .setLength(m_length)
00214              .setShift(m_shift);
00215     }
00216 
00217     /** \brief Adjoint (conjugate transpose) of the Householder sequence. */
00218     ConjugateReturnType  adjoint() const
00219     {
00220       return conjugate().setTrans(!m_trans);
00221     }
00222 
00223     /** \brief Inverse of the Householder sequence (equals the adjoint). */
00224     ConjugateReturnType  inverse() const { return adjoint(); }
00225 
00226     /** \internal */
00227     template<typename DestType> inline void evalTo(DestType& dst) const
00228     {
00229       Matrix<Scalar, DestType::RowsAtCompileTime, 1,
00230              AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
00231       evalTo(dst, workspace);
00232     }
00233 
00234     /** \internal */
00235     template<typename Dest, typename Workspace>
00236     void evalTo(Dest& dst, Workspace& workspace) const
00237     {
00238       workspace.resize(rows());
00239       Index vecs = m_length;
00240       if(    internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
00241           && internal::extract_data(dst) == internal::extract_data(m_vectors))
00242       {
00243         // in-place
00244         dst.diagonal().setOnes();
00245         dst.template triangularView<StrictlyUpper>().setZero();
00246         for(Index k = vecs-1; k >= 0; --k)
00247         {
00248           Index cornerSize = rows() - k - m_shift;
00249           if(m_trans)
00250             dst.bottomRightCorner(cornerSize, cornerSize)
00251                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00252           else
00253             dst.bottomRightCorner(cornerSize, cornerSize)
00254                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
00255 
00256           // clear the off diagonal vector
00257           dst.col(k).tail(rows()-k-1).setZero();
00258         }
00259         // clear the remaining columns if needed
00260         for(Index k = 0; k<cols()-vecs ; ++k)
00261           dst.col(k).tail(rows()-k-1).setZero();
00262       }
00263       else
00264       {
00265         dst.setIdentity(rows(), rows());
00266         for(Index k = vecs-1; k >= 0; --k)
00267         {
00268           Index cornerSize = rows() - k - m_shift;
00269           if(m_trans)
00270             dst.bottomRightCorner(cornerSize, cornerSize)
00271                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00272           else
00273             dst.bottomRightCorner(cornerSize, cornerSize)
00274                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
00275         }
00276       }
00277     }
00278 
00279     /** \internal */
00280     template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
00281     {
00282       Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
00283       applyThisOnTheRight(dst, workspace);
00284     }
00285 
00286     /** \internal */
00287     template<typename Dest, typename Workspace>
00288     inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
00289     {
00290       workspace.resize(dst.rows());
00291       for(Index k = 0; k < m_length; ++k)
00292       {
00293         Index actual_k = m_trans ? m_length-k-1 : k;
00294         dst.rightCols(rows()-m_shift-actual_k)
00295            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00296       }
00297     }
00298 
00299     /** \internal */
00300     template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
00301     {
00302       Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
00303       applyThisOnTheLeft(dst, workspace);
00304     }
00305 
00306     /** \internal */
00307     template<typename Dest, typename Workspace>
00308     inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
00309     {
00310       workspace.resize(dst.cols());
00311       for(Index k = 0; k < m_length; ++k)
00312       {
00313         Index actual_k = m_trans ? k : m_length-k-1;
00314         dst.bottomRows(rows()-m_shift-actual_k)
00315            .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
00316       }
00317     }
00318 
00319     /** \brief Computes the product of a Householder sequence with a matrix.
00320       * \param[in]  other  %Matrix being multiplied.
00321       * \returns    Expression object representing the product.
00322       *
00323       * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
00324       * and \f$ M \f$ is the matrix \p other.
00325       */
00326     template<typename OtherDerived>
00327     typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
00328     {
00329       typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
00330         res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
00331       applyThisOnTheLeft(res);
00332       return res;
00333     }
00334 
00335     template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
00336 
00337     /** \brief Sets the length of the Householder sequence.
00338       * \param [in]  length  New value for the length.
00339       *
00340       * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set
00341       * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that
00342       * is smaller. After this function is called, the length equals \p length.
00343       *
00344       * \sa length()
00345       */
00346     HouseholderSequence & setLength(Index length)
00347     {
00348       m_length = length;
00349       return *this;
00350     }
00351 
00352     /** \brief Sets the shift of the Householder sequence.
00353       * \param [in]  shift  New value for the shift.
00354       *
00355       * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th
00356       * column of the matrix \p v passed to the constructor corresponds to the i-th Householder
00357       * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}}
00358       * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th
00359       * Householder reflection.
00360       *
00361       * \sa shift()
00362       */
00363     HouseholderSequence & setShift(Index shift)
00364     {
00365       m_shift = shift;
00366       return *this;
00367     }
00368 
00369     Index length() const { return m_length; }  /**< \brief Returns the length of the Householder sequence. */
00370     Index shift() const { return m_shift; }    /**< \brief Returns the shift of the Householder sequence. */
00371 
00372     /* Necessary for .adjoint() and .conjugate() */
00373     template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence ;
00374 
00375   protected:
00376 
00377     /** \brief Sets the transpose flag.
00378       * \param [in]  trans  New value of the transpose flag.
00379       *
00380       * By default, the transpose flag is not set. If the transpose flag is set, then this object represents 
00381       * \f$ H^T = H_{n-1}^T \ldots H_1^T H_0^T \f$ instead of \f$ H = H_0 H_1 \ldots H_{n-1} \f$.
00382       *
00383       * \sa trans()
00384       */
00385     HouseholderSequence & setTrans(bool trans)
00386     {
00387       m_trans = trans;
00388       return *this;
00389     }
00390 
00391     bool trans() const { return m_trans; }     /**< \brief Returns the transpose flag. */
00392 
00393     typename VectorsType::Nested m_vectors;
00394     typename CoeffsType::Nested m_coeffs;
00395     bool m_trans;
00396     Index m_length;
00397     Index m_shift;
00398 };
00399 
00400 /** \brief Computes the product of a matrix with a Householder sequence.
00401   * \param[in]  other  %Matrix being multiplied.
00402   * \param[in]  h      %HouseholderSequence being multiplied.
00403   * \returns    Expression object representing the product.
00404   *
00405   * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the
00406   * Householder sequence represented by \p h.
00407   */
00408 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
00409 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator* (const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side> & h)
00410 {
00411   typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
00412     res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
00413   h.applyThisOnTheRight(res);
00414   return res;
00415 }
00416 
00417 /** \ingroup Householder_Module \householder_module
00418   * \brief Convenience function for constructing a Householder sequence. 
00419   * \returns A HouseholderSequence constructed from the specified arguments.
00420   */
00421 template<typename VectorsType, typename CoeffsType>
00422 HouseholderSequence<VectorsType,CoeffsType>  householderSequence(const VectorsType& v, const CoeffsType& h)
00423 {
00424   return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft> (v, h);
00425 }
00426 
00427 /** \ingroup Householder_Module \householder_module
00428   * \brief Convenience function for constructing a Householder sequence. 
00429   * \returns A HouseholderSequence constructed from the specified arguments.
00430   * \details This function differs from householderSequence() in that the template argument \p OnTheSide of
00431   * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft.
00432   */
00433 template<typename VectorsType, typename CoeffsType>
00434 HouseholderSequence<VectorsType,CoeffsType,OnTheRight>  rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
00435 {
00436   return HouseholderSequence<VectorsType,CoeffsType,OnTheRight> (v, h);
00437 }
00438 
00439 } // end namespace Eigen
00440 
00441 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H