diff options
| -rw-r--r-- | Eigen/src/Core/util/Constants.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Householder/HouseholderSequence.h | 132 | ||||
| -rw-r--r-- | doc/snippets/HouseholderSequence_HouseholderSequence.cpp | 31 |
3 files changed, 159 insertions, 10 deletions
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 03cadc107..7ba91d6db 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -210,9 +210,11 @@ enum { DontAlign = 0x2 }; +/** \brief Enum for specifying whether to apply or solve on the left or right. + */ enum { - OnTheLeft = 1, - OnTheRight = 2 + OnTheLeft = 1, /**< \brief Apply transformation on the left. */ + OnTheRight = 2 /**< \brief Apply transformation on the right. */ }; /* the following could as well be written: diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index d3616ed70..4b4def02e 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -29,12 +29,28 @@ /** \ingroup Householder_Module * \householder_module * \class HouseholderSequence - * \brief Represents a sequence of householder reflections with decreasing size + * \brief Sequence of Householder reflections acting on subspaces with decreasing size + * \tparam VectorsType type of matrix containing the Householder vectors + * \tparam CoeffsType type of vector containing the Householder coefficients + * \tparam Side either OnTheLeft (the default) or OnTheRight * - * This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$ - * where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$, - * \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$ - * and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i]. + * This class represents a product sequence of Householder reflections where the first Householder reflection + * acts on the whole space, the second Householder reflection leaves the one-dimensional subspace spanned by + * the first unit vector invariant, the third Householder reflection leaves the two-dimensional subspace + * spanned by the first two unit vectors invariant, and so on up to the last reflection which leaves all but + * one dimensions invariant and acts only on the last dimension. Such sequences of Householder reflections + * are used in several algorithms to zero out certain parts of a matrix. Indeed, the methods + * HessenbergDecomposition::matrixQ(), Tridiagonalization::matrixQ(), HouseholderQR::householderQ(), + * and ColPivHouseholderQR::householderQ() all return a %HouseholderSequence. + * + * More precisely, the class %HouseholderSequence represents an \f$ n \times n \f$ matrix \f$ H \f$ of the + * 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 + * v_i^* \f$. The i-th Householder coefficient \f$ h_i \f$ is a scalar and the i-th Householder vector \f$ + * v_i \f$ is a vector of the form + * \f[ + * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. + * \f] + * The last \f$ n-i \f$ entries of \f$ v_i \f$ are called the essential part of the Householder vector. * * Typical usages are listed below, where H is a HouseholderSequence: * \code @@ -46,6 +62,8 @@ * \endcode * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate operators. * + * See the documentation for HouseholderSequence(const VectorsType&, const CoeffsType&) for an example. + * * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ @@ -129,12 +147,30 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS Side > ConjugateReturnType; + /** \brief Constructor. + * \param[in] v %Matrix containing the essential parts of the Householder vectors + * \param[in] h Vector containing the Householder coefficients + * + * Constructs the Householder sequence with coefficients given by \p h and vectors given by \p v. The + * i-th Householder coefficient \f$ h_i \f$ is given by \p h(i) and the essential part of the i-th + * Householder vector \f$ v_i \f$ is given by \p v(k,i) with \p k > \p i (the subdiagonal part of the + * i-th column). If \p v has fewer columns than rows, then the Householder sequence contains as many + * Householder reflections as there are columns. + * + * \note The %HouseholderSequence object stores \p v and \p h by reference. + * + * Example: \include HouseholderSequence_HouseholderSequence.cpp + * Output: \verbinclude HouseholderSequence_HouseholderSequence.out + * + * \sa setLength(), setShift() + */ HouseholderSequence(const VectorsType& v, const CoeffsType& h) : m_vectors(v), m_coeffs(h), m_trans(false), m_length(v.diagonalSize()), m_shift(0) { } + /** \brief Copy constructor. */ HouseholderSequence(const HouseholderSequence& other) : m_vectors(other.m_vectors), m_coeffs(other.m_coeffs), @@ -144,20 +180,45 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS { } + /** \brief Number of rows of transformation viewed as a matrix. + * \returns Number of rows + * \details This equals the dimension of the space that the transformation acts on. + */ Index rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } + + /** \brief Number of columns of transformation viewed as a matrix. + * \returns Number of columns + * \details This equals the dimension of the space that the transformation acts on. + */ Index cols() const { return rows(); } + /** \brief Essential part of a Householder vector. + * \param[in] k Index of Householder reflection + * \returns Vector containing non-trivial entries of k-th Householder vector + * + * This function returns the essential part of the Householder vector \f$ v_i \f$. This is a vector of + * length \f$ n-i \f$ containing the last \f$ n-i \f$ entries of the vector + * \f[ + * v_i = [\underbrace{0, \ldots, 0}_{i-1\mbox{ zeros}}, 1, \underbrace{*, \ldots,*}_{n-i\mbox{ arbitrary entries}} ]. + * \f] + * The index \f$ i \f$ equals \p k + shift(), corresponding to the k-th column of the matrix \p v + * passed to the constructor. + * + * \sa setShift(), shift() + */ const EssentialVectorType essentialVector(Index k) const { eigen_assert(k >= 0 && k < m_length); return internal::hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k); } + /** \brief %Transpose of the Householder sequence. */ HouseholderSequence transpose() const { return HouseholderSequence(*this).setTrans(!m_trans); } + /** \brief Complex conjugate of the Householder sequence. */ ConjugateReturnType conjugate() const { return ConjugateReturnType(m_vectors, m_coeffs.conjugate()) @@ -166,11 +227,13 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS .setShift(m_shift); } + /** \brief Adjoint (conjugate transpose) of the Householder sequence. */ ConjugateReturnType adjoint() const { return conjugate().setTrans(!m_trans); } + /** \brief Inverse of the Householder sequence (equals the adjoint). */ ConjugateReturnType inverse() const { return adjoint(); } /** \internal */ @@ -243,6 +306,13 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS } } + /** \brief Computes the product of a Householder sequence with a matrix. + * \param[in] other %Matrix being multiplied. + * \returns Expression object representing the product. + * + * This function computes \f$ HM \f$ where \f$ H \f$ is the Householder sequence represented by \p *this + * and \f$ M \f$ is the matrix \p other. + */ template<typename OtherDerived> typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other) const { @@ -252,6 +322,14 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS return res; } + /** \brief Computes the product of a matrix with a Householder sequence. + * \param[in] other %Matrix being multiplied. + * \param[in] h %HouseholderSequence being multiplied. + * \returns Expression object representing the product. + * + * This function computes \f$ MH \f$ where \f$ M \f$ is the matrix \p other and \f$ H \f$ is the + * Householder sequence represented by \p h. + */ template<typename OtherDerived> friend typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h) { @@ -263,27 +341,55 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct internal::hseq_side_dependent_impl; + /** \brief Sets the transpose flag. + * \param [in] trans New value of the transpose flag. + * + * By default, the transpose flag is not set. If the transpose flag is set, then this object represents + * \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$. + * + * \sa trans() + */ HouseholderSequence& setTrans(bool trans) { m_trans = trans; return *this; } + /** \brief Sets the length of the Householder sequence. + * \param [in] length New value for the length. + * + * By default, the length \f$ n \f$ of the Householder sequence \f$ H = H_0 H_1 \ldots H_{n-1} \f$ is set + * to the number of columns of the matrix \p v passed to the constructor, or the number of rows if that + * is smaller. After this function is called, the length equals \p length. + * + * \sa length() + */ HouseholderSequence& setLength(Index length) { m_length = length; return *this; } + /** \brief Sets the shift of the Householder sequence. + * \param [in] shift New value for the shift. + * + * By default, a %HouseholderSequence object represents \f$ H = H_0 H_1 \ldots H_{n-1} \f$ and the i-th + * column of the matrix \p v passed to the constructor corresponds to the i-th Householder + * reflection. After this function is called, the object represents \f$ H = H_{\mathrm{shift}} + * H_{\mathrm{shift}+1} \ldots H_{n-1} \f$ and the i-th column of \p v corresponds to the (shift+i)-th + * Householder reflection. + * + * \sa shift() + */ HouseholderSequence& setShift(Index shift) { m_shift = shift; return *this; } - bool trans() const { return m_trans; } - Index length() const { return m_length; } - Index shift() const { return m_shift; } + bool trans() const { return m_trans; } /**< \brief Returns the transpose flag. */ + Index length() const { return m_length; } /**< \brief Returns the length of the Householder sequence. */ + Index shift() const { return m_shift; } /**< \brief Returns the shift of the Householder sequence. */ protected: typename VectorsType::Nested m_vectors; @@ -293,12 +399,22 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS Index m_shift; }; +/** \ingroup Householder_Module \householder_module + * \brief Convenience function for constructing a Householder sequence. + * \returns A HouseholderSequence constructed from the specified arguments. + */ template<typename VectorsType, typename CoeffsType> HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h) { return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h); } +/** \ingroup Householder_Module \householder_module + * \brief Convenience function for constructing a Householder sequence. + * \returns A HouseholderSequence constructed from the specified arguments. + * \details This function differs from householderSequence() in that the template argument \p OnTheSide of + * the constructed HouseholderSequence is set to OnTheRight, instead of the default OnTheLeft. + */ template<typename VectorsType, typename CoeffsType> HouseholderSequence<VectorsType,CoeffsType,OnTheRight> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h) { diff --git a/doc/snippets/HouseholderSequence_HouseholderSequence.cpp b/doc/snippets/HouseholderSequence_HouseholderSequence.cpp new file mode 100644 index 000000000..2632b83b9 --- /dev/null +++ b/doc/snippets/HouseholderSequence_HouseholderSequence.cpp @@ -0,0 +1,31 @@ +Matrix3d v = Matrix3d::Random(); +cout << "The matrix v is:" << endl; +cout << v << endl; + +Vector3d v0(1, v(1,0), v(2,0)); +cout << "The first Householder vector is: v_0 = " << v0.transpose() << endl; +Vector3d v1(0, 1, v(2,1)); +cout << "The second Householder vector is: v_1 = " << v1.transpose() << endl; +Vector3d v2(0, 0, 1); +cout << "The third Householder vector is: v_2 = " << v2.transpose() << endl; + +Vector3d h = Vector3d::Random(); +cout << "The Householder coefficients are: h = " << h.transpose() << endl; + +Matrix3d H0 = Matrix3d::Identity() - h(0) * v0 * v0.adjoint(); +cout << "The first Householder reflection is represented by H_0 = " << endl; +cout << H0 << endl; +Matrix3d H1 = Matrix3d::Identity() - h(1) * v1 * v1.adjoint(); +cout << "The second Householder reflection is represented by H_1 = " << endl; +cout << H1 << endl; +Matrix3d H2 = Matrix3d::Identity() - h(2) * v2 * v2.adjoint(); +cout << "The third Householder reflection is represented by H_2 = " << endl; +cout << H2 << endl; +cout << "Their product is H_0 H_1 H_2 = " << endl; +cout << H0 * H1 * H2 << endl; + +HouseholderSequence<Matrix3d, Vector3d> hhSeq(v, h); +Matrix3d hhSeqAsMatrix(hhSeq); +cout << "If we construct a HouseholderSequence from v and h" << endl; +cout << "and convert it to a matrix, we get:" << endl; +cout << hhSeqAsMatrix << endl; |