Home | History | Annotate | Download | only in Householder
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1 (at) gmail.com>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
     12 #define EIGEN_HOUSEHOLDER_SEQUENCE_H
     13 
     14 namespace Eigen {
     15 
     16 /** \ingroup Householder_Module
     17   * \householder_module
     18   * \class HouseholderSequence
     19   * \brief Sequence of Householder reflections acting on subspaces with decreasing size
     20   * \tparam VectorsType type of matrix containing the Householder vectors
     21   * \tparam CoeffsType  type of vector containing the Householder coefficients
     22   * \tparam Side        either OnTheLeft (the default) or OnTheRight
     23   *
     24   * This class represents a product sequence of Householder reflections where the first Householder reflection
     25   * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by
     26   * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace
     27   * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but
     28   * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections
     29   * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods
     30   * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(),
     31   * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence.
     32   *
     33   * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the
     34   * 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
     35   * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$
     36   * v_i \f$ is a vector of the form
     37   * \f[
     38   * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ].
     39   * \f]
     40   * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector.
     41   *
     42   * Typical usages are listed below, where H is a HouseholderSequence:
     43   * \code
     44   * A.applyOnTheRight(H);             // A = A * H
     45   * A.applyOnTheLeft(H);              // A = H * A
     46   * A.applyOnTheRight(H.adjoint());   // A = A * H^*
     47   * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
     48   * MatrixXd Q = H;                   // conversion to a dense matrix
     49   * \endcode
     50   * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators.
     51   *
     52   * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example.
     53   *
     54   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
     55   */
     56 
     57 namespace internal {
     58 
     59 template<typename VectorsType, typename CoeffsType, int Side>
     60 struct traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
     61 {
     62   typedef typename VectorsType::Scalar Scalar;
     63   typedef typename VectorsType::Index Index;
     64   typedef typename VectorsType::StorageKind StorageKind;
     65   enum {
     66     RowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::RowsAtCompileTime
     67                                         : traits<VectorsType>::ColsAtCompileTime,
     68     ColsAtCompileTime = RowsAtCompileTime,
     69     MaxRowsAtCompileTime = Side==OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime
     70                                            : traits<VectorsType>::MaxColsAtCompileTime,
     71     MaxColsAtCompileTime = MaxRowsAtCompileTime,
     72     Flags = 0
     73   };
     74 };
     75 
     76 template<typename VectorsType, typename CoeffsType, int Side>
     77 struct hseq_side_dependent_impl
     78 {
     79   typedef Block<const VectorsType, Dynamic, 1> EssentialVectorType;
     80   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
     81   typedef typename VectorsType::Index Index;
     82   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
     83   {
     84     Index start = k+1+h.m_shift;
     85     return Block<const VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
     86   }
     87 };
     88 
     89 template<typename VectorsType, typename CoeffsType>
     90 struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
     91 {
     92   typedef Transpose<Block<const VectorsType, 1, Dynamic> > EssentialVectorType;
     93   typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
     94   typedef typename VectorsType::Index Index;
     95   static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k)
     96   {
     97     Index start = k+1+h.m_shift;
     98     return Block<const VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
     99   }
    100 };
    101 
    102 template<typename OtherScalarType, typename MatrixType> struct matrix_type_times_scalar_type
    103 {
    104   typedef typename scalar_product_traits<OtherScalarType, typename MatrixType::Scalar>::ReturnType
    105     ResultScalar;
    106   typedef Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
    107                  0, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime> Type;
    108 };
    109 
    110 } // end namespace internal
    111 
    112 template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
    113   : public EigenBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
    114 {
    115     typedef typename internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType EssentialVectorType;
    116 
    117   public:
    118     enum {
    119       RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
    120       ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
    121       MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
    122       MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
    123     };
    124     typedef typename internal::traits<HouseholderSequence>::Scalar Scalar;
    125     typedef typename VectorsType::Index Index;
    126 
    127     typedef HouseholderSequence<
    128       typename internal::conditional<NumTraits<Scalar>::IsComplex,
    129         typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
    130         VectorsType>::type,
    131       typename internal::conditional<NumTraits<Scalar>::IsComplex,
    132         typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
    133         CoeffsType>::type,
    134       Side
    135     > ConjugateReturnType;
    136 
    137     /** \brief Constructor.
    138       * \param[in]  v      %Matrix containing the essential parts of the Householder vectors
    139       * \param[in]  h      Vector containing the Householder coefficients
    140       *
    141       * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The
    142       * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th
    143       * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the
    144       * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many
    145       * Householder reflections as there are columns.
    146       *
    147       * \note The %HouseholderSequence object stores \p v and \p h by reference.
    148       *
    149       * Example: \include HouseholderSequence_HouseholderSequence.cpp
    150       * Output: \verbinclude HouseholderSequence_HouseholderSequence.out
    151       *
    152       * \sa setLength(), setShift()
    153       */
    154     HouseholderSequence(const VectorsType& v, const CoeffsType& h)
    155       : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()),
    156         m_shift(0)
    157     {
    158     }
    159 
    160     /** \brief Copy constructor. */
    161     HouseholderSequence(const HouseholderSequence& other)
    162       : m_vectors(other.m_vectors),
    163         m_coeffs(other.m_coeffs),
    164         m_trans(other.m_trans),
    165         m_length(other.m_length),
    166         m_shift(other.m_shift)
    167     {
    168     }
    169 
    170     /** \brief Number of rows of transformation viewed as a matrix.
    171       * \returns Number of rows
    172       * \details This equals the dimension of the space that the transformation acts on.
    173       */
    174     Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
    175 
    176     /** \brief Number of columns of transformation viewed as a matrix.
    177       * \returns Number of columns
    178       * \details This equals the dimension of the space that the transformation acts on.
    179       */
    180     Index cols() const { return rows(); }
    181 
    182     /** \brief Essential part of a Householder vector.
    183       * \param[in]  k  Index of Householder reflection
    184       * \returns    Vector containing non-trivial entries of k-th Householder vector
    185       *
    186       * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of
    187       * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector
    188       * \f[
    189       * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ].
    190       * \f]
    191       * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v
    192       * passed to the constructor.
    193       *
    194       * \sa setShift(), shift()
    195       */
    196     const EssentialVectorType essentialVector(Index k) const
    197     {
    198       eigen_assert(k >= 0 && k < m_length);
    199       return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
    200     }
    201 
    202     /** \brief %Transpose of the Householder sequence. */
    203     HouseholderSequence transpose() const
    204     {
    205       return HouseholderSequence(*this).setTrans(!m_trans);
    206     }
    207 
    208     /** \brief Complex conjugate of the Householder sequence. */
    209     ConjugateReturnType conjugate() const
    210     {
    211       return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
    212              .setTrans(m_trans)
    213              .setLength(m_length)
    214              .setShift(m_shift);
    215     }
    216 
    217     /** \brief Adjoint (conjugate transpose) of the Householder sequence. */
    218     ConjugateReturnType adjoint() const
    219     {
    220       return conjugate().setTrans(!m_trans);
    221     }
    222 
    223     /** \brief Inverse of the Householder sequence (equals the adjoint). */
    224     ConjugateReturnType inverse() const { return adjoint(); }
    225 
    226     /** \internal */
    227     template<typename DestType> inline void evalTo(DestType& dst) const
    228     {
    229       Matrix<Scalar, DestType::RowsAtCompileTime, 1,
    230              AutoAlign|ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(rows());
    231       evalTo(dst, workspace);
    232     }
    233 
    234     /** \internal */
    235     template<typename Dest, typename Workspace>
    236     void evalTo(Dest& dst, Workspace& workspace) const
    237     {
    238       workspace.resize(rows());
    239       Index vecs = m_length;
    240       if(    internal::is_same<typename internal::remove_all<VectorsType>::type,Dest>::value
    241           && internal::extract_data(dst) == internal::extract_data(m_vectors))
    242       {
    243         // in-place
    244         dst.diagonal().setOnes();
    245         dst.template triangularView<StrictlyUpper>().setZero();
    246         for(Index k = vecs-1; k >= 0; --k)
    247         {
    248           Index cornerSize = rows() - k - m_shift;
    249           if(m_trans)
    250             dst.bottomRightCorner(cornerSize, cornerSize)
    251                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
    252           else
    253             dst.bottomRightCorner(cornerSize, cornerSize)
    254                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
    255 
    256           // clear the off diagonal vector
    257           dst.col(k).tail(rows()-k-1).setZero();
    258         }
    259         // clear the remaining columns if needed
    260         for(Index k = 0; k<cols()-vecs ; ++k)
    261           dst.col(k).tail(rows()-k-1).setZero();
    262       }
    263       else
    264       {
    265         dst.setIdentity(rows(), rows());
    266         for(Index k = vecs-1; k >= 0; --k)
    267         {
    268           Index cornerSize = rows() - k - m_shift;
    269           if(m_trans)
    270             dst.bottomRightCorner(cornerSize, cornerSize)
    271                .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
    272           else
    273             dst.bottomRightCorner(cornerSize, cornerSize)
    274                .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
    275         }
    276       }
    277     }
    278 
    279     /** \internal */
    280     template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
    281     {
    282       Matrix<Scalar,1,Dest::RowsAtCompileTime,RowMajor,1,Dest::MaxRowsAtCompileTime> workspace(dst.rows());
    283       applyThisOnTheRight(dst, workspace);
    284     }
    285 
    286     /** \internal */
    287     template<typename Dest, typename Workspace>
    288     inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const
    289     {
    290       workspace.resize(dst.rows());
    291       for(Index k = 0; k < m_length; ++k)
    292       {
    293         Index actual_k = m_trans ? m_length-k-1 : k;
    294         dst.rightCols(rows()-m_shift-actual_k)
    295            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
    296       }
    297     }
    298 
    299     /** \internal */
    300     template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
    301     {
    302       Matrix<Scalar,1,Dest::ColsAtCompileTime,RowMajor,1,Dest::MaxColsAtCompileTime> workspace(dst.cols());
    303       applyThisOnTheLeft(dst, workspace);
    304     }
    305 
    306     /** \internal */
    307     template<typename Dest, typename Workspace>
    308     inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace) const
    309     {
    310       workspace.resize(dst.cols());
    311       for(Index k = 0; k < m_length; ++k)
    312       {
    313         Index actual_k = m_trans ? k : m_length-k-1;
    314         dst.bottomRows(rows()-m_shift-actual_k)
    315            .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
    316       }
    317     }
    318 
    319     /** \brief Computes the product of a Householder sequence with a matrix.
    320       * \param[in]  other  %Matrix being multiplied.
    321       * \returns    Expression object representing the product.
    322       *
    323       * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this
    324       * and \f$ M \f$ is the matrix \p other.
    325       */
    326     template<typename OtherDerived>
    327     typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const
    328     {
    329       typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type
    330         res(other.template cast<typename internal::matrix_type_times_scalar_type<Scalar,OtherDerived>::ResultScalar>());
    331       applyThisOnTheLeft(res);
    332       return res;
    333     }
    334 
    335     template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl;
    336 
    337     /** \brief Sets the length of the Householder sequence.
    338       * \param [in]  length  New value for the length.
    339       *
    340       * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set
    341       * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that
    342       * is smaller. After this function is called, the length equals \p length.
    343       *
    344       * \sa length()
    345       */
    346     HouseholderSequence& setLength(Index length)
    347     {
    348       m_length = length;
    349       return *this;
    350     }
    351 
    352     /** \brief Sets the shift of the Householder sequence.
    353       * \param [in]  shift  New value for the shift.
    354       *
    355       * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th
    356       * column of the matrix \p v passed to the constructor corresponds to the i-th Householder
    357       * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}}
    358       * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th
    359       * Householder reflection.
    360       *
    361       * \sa shift()
    362       */
    363     HouseholderSequence& setShift(Index shift)
    364     {
    365       m_shift = shift;
    366       return *this;
    367     }
    368 
    369     Index length() const { return m_length; }  /**< \brief Returns the length of the Householder sequence. */
    370     Index shift() const { return m_shift; }    /**< \brief Returns the shift of the Householder sequence. */
    371 
    372     /* Necessary for .adjoint() and .conjugate() */
    373     template <typename VectorsType2, typename CoeffsType2, int Side2> friend class HouseholderSequence;
    374 
    375   protected:
    376 
    377     /** \brief Sets the transpose flag.
    378       * \param [in]  trans  New value of the transpose flag.
    379       *
    380       * By default, the transpose flag is not set. If the transpose flag is set, then this object represents
    381       * \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$.
    382       *
    383       * \sa trans()
    384       */
    385     HouseholderSequence& setTrans(bool trans)
    386     {
    387       m_trans = trans;
    388       return *this;
    389     }
    390 
    391     bool trans() const { return m_trans; }     /**< \brief Returns the transpose flag. */
    392 
    393     typename VectorsType::Nested m_vectors;
    394     typename CoeffsType::Nested m_coeffs;
    395     bool m_trans;
    396     Index m_length;
    397     Index m_shift;
    398 };
    399 
    400 /** \brief Computes the product of a matrix with a Householder sequence.
    401   * \param[in]  other  %Matrix being multiplied.
    402   * \param[in]  h      %HouseholderSequence being multiplied.
    403   * \returns    Expression object representing the product.
    404   *
    405   * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the
    406   * Householder sequence represented by \p h.
    407   */
    408 template<typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
    409 typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType,CoeffsType,Side>& h)
    410 {
    411   typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::Type
    412     res(other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,OtherDerived>::ResultScalar>());
    413   h.applyThisOnTheRight(res);
    414   return res;
    415 }
    416 
    417 /** \ingroup Householder_Module \householder_module
    418   * \brief Convenience function for constructing a Householder sequence.
    419   * \returns A HouseholderSequence constructed from the specified arguments.
    420   */
    421 template<typename VectorsType, typename CoeffsType>
    422 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h)
    423 {
    424   return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h);
    425 }
    426 
    427 /** \ingroup Householder_Module \householder_module
    428   * \brief Convenience function for constructing a Householder sequence.
    429   * \returns A HouseholderSequence constructed from the specified arguments.
    430   * \details This function differs from householderSequence() in that the template argument \p OnTheSide of
    431   * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft.
    432   */
    433 template<typename VectorsType, typename CoeffsType>
    434 HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h)
    435 {
    436   return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h);
    437 }
    438 
    439 } // end namespace Eigen
    440 
    441 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
    442