Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
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
Generated on Tue Jul 12 2022 17:46:55 by 1.7.2