Main MRPT website > C++ reference
MRPT logo

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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011