Main MRPT website > C++ reference
MRPT logo

HessenbergDecomposition.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) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
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_HESSENBERGDECOMPOSITION_H
00027 #define EIGEN_HESSENBERGDECOMPOSITION_H
00028 
00029 namespace internal {
00030   
00031 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
00032 template<typename MatrixType>
00033 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00034 {
00035   typedef MatrixType ReturnType;
00036 };
00037 
00038 }
00039 
00040 /** \eigenvalues_module \ingroup Eigenvalues_Module
00041   *
00042   *
00043   * \class HessenbergDecomposition
00044   *
00045   * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
00046   *
00047   * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
00048   *
00049   * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
00050   * the real case, the Hessenberg decomposition consists of an orthogonal
00051   * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
00052   * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
00053   * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
00054   * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
00055   * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
00056   * \f$ Q^{-1} = Q^* \f$).
00057   *
00058   * Call the function compute() to compute the Hessenberg decomposition of a
00059   * given matrix. Alternatively, you can use the
00060   * HessenbergDecomposition(const MatrixType&) constructor which computes the
00061   * Hessenberg decomposition at construction time. Once the decomposition is
00062   * computed, you can use the matrixH() and matrixQ() functions to construct
00063   * the matrices H and Q in the decomposition.
00064   *
00065   * The documentation for matrixH() contains an example of the typical use of
00066   * this class.
00067   *
00068   * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
00069   */
00070 template<typename _MatrixType> class HessenbergDecomposition
00071 {
00072   public:
00073 
00074     /** \brief Synonym for the template parameter \p _MatrixType. */
00075     typedef _MatrixType MatrixType;
00076 
00077     enum {
00078       Size = MatrixType::RowsAtCompileTime,
00079       SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
00080       Options = MatrixType::Options,
00081       MaxSize = MatrixType::MaxRowsAtCompileTime,
00082       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
00083     };
00084 
00085     /** \brief Scalar type for matrices of type #MatrixType. */
00086     typedef typename MatrixType::Scalar Scalar;
00087     typedef typename MatrixType::Index Index;
00088 
00089     /** \brief Type for vector of Householder coefficients.
00090       *
00091       * This is column vector with entries of type #Scalar. The length of the
00092       * vector is one less than the size of #MatrixType, if it is a fixed-side
00093       * type.
00094       */
00095     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00096 
00097     /** \brief Return type of matrixQ() */
00098     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00099     
00100     typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
00101 
00102     /** \brief Default constructor; the decomposition will be computed later.
00103       *
00104       * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
00105       *
00106       * The default constructor is useful in cases in which the user intends to
00107       * perform decompositions via compute().  The \p size parameter is only
00108       * used as a hint. It is not an error to give a wrong \p size, but it may
00109       * impair performance.
00110       *
00111       * \sa compute() for an example.
00112       */
00113     HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
00114       : m_matrix(size,size),
00115         m_temp(size),
00116         m_isInitialized(false)
00117     {
00118       if(size>1)
00119         m_hCoeffs.resize(size-1);
00120     }
00121 
00122     /** \brief Constructor; computes Hessenberg decomposition of given matrix.
00123       *
00124       * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
00125       *
00126       * This constructor calls compute() to compute the Hessenberg
00127       * decomposition.
00128       *
00129       * \sa matrixH() for an example.
00130       */
00131     HessenbergDecomposition(const MatrixType& matrix)
00132       : m_matrix(matrix),
00133         m_temp(matrix.rows()),
00134         m_isInitialized(false)
00135     {
00136       if(matrix.rows()<2)
00137       {
00138         m_isInitialized = true;
00139         return;
00140       }
00141       m_hCoeffs.resize(matrix.rows()-1,1);
00142       _compute(m_matrix, m_hCoeffs, m_temp);
00143       m_isInitialized = true;
00144     }
00145 
00146     /** \brief Computes Hessenberg decomposition of given matrix.
00147       *
00148       * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
00149       * \returns    Reference to \c *this
00150       *
00151       * The Hessenberg decomposition is computed by bringing the columns of the
00152       * matrix successively in the required form using Householder reflections
00153       * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
00154       * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
00155       * denotes the size of the given matrix.
00156       *
00157       * This method reuses of the allocated data in the HessenbergDecomposition
00158       * object.
00159       *
00160       * Example: \include HessenbergDecomposition_compute.cpp
00161       * Output: \verbinclude HessenbergDecomposition_compute.out
00162       */
00163     HessenbergDecomposition& compute(const MatrixType& matrix)
00164     {
00165       m_matrix = matrix;
00166       if(matrix.rows()<2)
00167       {
00168         m_isInitialized = true;
00169         return *this;
00170       }
00171       m_hCoeffs.resize(matrix.rows()-1,1);
00172       _compute(m_matrix, m_hCoeffs, m_temp);
00173       m_isInitialized = true;
00174       return *this;
00175     }
00176 
00177     /** \brief Returns the Householder coefficients.
00178       *
00179       * \returns a const reference to the vector of Householder coefficients
00180       *
00181       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
00182       * or the member function compute(const MatrixType&) has been called
00183       * before to compute the Hessenberg decomposition of a matrix.
00184       *
00185       * The Householder coefficients allow the reconstruction of the matrix
00186       * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
00187       *
00188       * \sa packedMatrix(), \ref Householder_Module "Householder module"
00189       */
00190     const CoeffVectorType& householderCoefficients() const
00191     {
00192       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00193       return m_hCoeffs;
00194     }
00195 
00196     /** \brief Returns the internal representation of the decomposition
00197       *
00198       * \returns a const reference to a matrix with the internal representation
00199       *          of the decomposition.
00200       *
00201       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
00202       * or the member function compute(const MatrixType&) has been called
00203       * before to compute the Hessenberg decomposition of a matrix.
00204       *
00205       * The returned matrix contains the following information:
00206       *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
00207       *  - the rest of the lower part contains the Householder vectors that, combined with
00208       *    Householder coefficients returned by householderCoefficients(),
00209       *    allows to reconstruct the matrix Q as
00210       *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
00211       *    Here, the matrices \f$ H_i \f$ are the Householder transformations
00212       *       \f$ H_i = (I - h_i v_i v_i^T) \f$
00213       *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
00214       *    \f$ v_i \f$ is the Householder vector defined by
00215       *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
00216       *    with M the matrix returned by this function.
00217       *
00218       * See LAPACK for further details on this packed storage.
00219       *
00220       * Example: \include HessenbergDecomposition_packedMatrix.cpp
00221       * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
00222       *
00223       * \sa householderCoefficients()
00224       */
00225     const MatrixType& packedMatrix() const
00226     {
00227       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00228       return m_matrix;
00229     }
00230 
00231     /** \brief Reconstructs the orthogonal matrix Q in the decomposition
00232       *
00233       * \returns object representing the matrix Q
00234       *
00235       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
00236       * or the member function compute(const MatrixType&) has been called
00237       * before to compute the Hessenberg decomposition of a matrix.
00238       *
00239       * This function returns a light-weight object of template class
00240       * HouseholderSequence. You can either apply it directly to a matrix or
00241       * you can convert it to a matrix of type #MatrixType.
00242       *
00243       * \sa matrixH() for an example, class HouseholderSequence
00244       */
00245     HouseholderSequenceType matrixQ() const
00246     {
00247       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00248       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00249              .setLength(m_matrix.rows() - 1)
00250              .setShift(1);
00251     }
00252 
00253     /** \brief Constructs the Hessenberg matrix H in the decomposition
00254       *
00255       * \returns expression object representing the matrix H
00256       *
00257       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
00258       * or the member function compute(const MatrixType&) has been called
00259       * before to compute the Hessenberg decomposition of a matrix.
00260       *
00261       * The object returned by this function constructs the Hessenberg matrix H
00262       * when it is assigned to a matrix or otherwise evaluated. The matrix H is
00263       * constructed from the packed matrix as returned by packedMatrix(): The
00264       * upper part (including the subdiagonal) of the packed matrix contains
00265       * the matrix H. It may sometimes be better to directly use the packed
00266       * matrix instead of constructing the matrix H.
00267       *
00268       * Example: \include HessenbergDecomposition_matrixH.cpp
00269       * Output: \verbinclude HessenbergDecomposition_matrixH.out
00270       *
00271       * \sa matrixQ(), packedMatrix()
00272       */
00273     MatrixHReturnType matrixH() const
00274     {
00275       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00276       return MatrixHReturnType(*this);
00277     }
00278 
00279   private:
00280 
00281     typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
00282     typedef typename NumTraits<Scalar>::Real RealScalar;
00283     static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
00284 
00285   protected:
00286     MatrixType m_matrix;
00287     CoeffVectorType m_hCoeffs;
00288     VectorType m_temp;
00289     bool m_isInitialized;
00290 };
00291 
00292 /** \internal
00293   * Performs a tridiagonal decomposition of \a matA in place.
00294   *
00295   * \param matA the input selfadjoint matrix
00296   * \param hCoeffs returned Householder coefficients
00297   *
00298   * The result is written in the lower triangular part of \a matA.
00299   *
00300   * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
00301   *
00302   * \sa packedMatrix()
00303   */
00304 template<typename MatrixType>
00305 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
00306 {
00307   assert(matA.rows()==matA.cols());
00308   Index n = matA.rows();
00309   temp.resize(n);
00310   for (Index i = 0; i<n-1; ++i)
00311   {
00312     // let's consider the vector v = i-th column starting at position i+1
00313     Index remainingSize = n-i-1;
00314     RealScalar beta;
00315     Scalar h;
00316     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00317     matA.col(i).coeffRef(i+1) = beta;
00318     hCoeffs.coeffRef(i) = h;
00319 
00320     // Apply similarity transformation to remaining columns,
00321     // i.e., compute A = H A H'
00322 
00323     // A = H A
00324     matA.bottomRightCorner(remainingSize, remainingSize)
00325         .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
00326 
00327     // A = A H'
00328     matA.rightCols(remainingSize)
00329         .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
00330   }
00331 }
00332 
00333 namespace internal {
00334 
00335 /** \eigenvalues_module \ingroup Eigenvalues_Module
00336   *
00337   *
00338   * \brief Expression type for return value of HessenbergDecomposition::matrixH()
00339   *
00340   * \tparam MatrixType type of matrix in the Hessenberg decomposition
00341   *
00342   * Objects of this type represent the Hessenberg matrix in the Hessenberg
00343   * decomposition of some matrix. The object holds a reference to the
00344   * HessenbergDecomposition class until the it is assigned or evaluated for
00345   * some other reason (the reference should remain valid during the life time
00346   * of this object). This class is the return type of
00347   * HessenbergDecomposition::matrixH(); there is probably no other use for this
00348   * class.
00349   */
00350 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
00351 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00352 {
00353     typedef typename MatrixType::Index Index;
00354   public:
00355     /** \brief Constructor.
00356       *
00357       * \param[in] hess  Hessenberg decomposition
00358       */
00359     HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
00360 
00361     /** \brief Hessenberg matrix in decomposition.
00362       *
00363       * \param[out] result  Hessenberg matrix in decomposition \p hess which
00364       *                     was passed to the constructor
00365       */
00366     template <typename ResultType>
00367     inline void evalTo(ResultType& result) const
00368     {
00369       result = m_hess.packedMatrix();
00370       Index n = result.rows();
00371       if (n>2)
00372         result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
00373     }
00374 
00375     Index rows() const { return m_hess.packedMatrix().rows(); }
00376     Index cols() const { return m_hess.packedMatrix().cols(); }
00377 
00378   protected:
00379     const HessenbergDecomposition<MatrixType>& m_hess;
00380 };
00381 
00382 }
00383 
00384 #endif // EIGEN_HESSENBERGDECOMPOSITION_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011