Home | History | Annotate | Download | only in Eigenvalues
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 // Copyright (C) 2010 Jitse Niesen <jitse (at) maths.leeds.ac.uk>
      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_HESSENBERGDECOMPOSITION_H
     12 #define EIGEN_HESSENBERGDECOMPOSITION_H
     13 
     14 namespace Eigen {
     15 
     16 namespace internal {
     17 
     18 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
     19 template<typename MatrixType>
     20 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
     21 {
     22   typedef MatrixType ReturnType;
     23 };
     24 
     25 }
     26 
     27 /** \eigenvalues_module \ingroup Eigenvalues_Module
     28   *
     29   *
     30   * \class HessenbergDecomposition
     31   *
     32   * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
     33   *
     34   * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
     35   *
     36   * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
     37   * the real case, the Hessenberg decomposition consists of an orthogonal
     38   * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
     39   * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
     40   * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
     41   * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
     42   * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
     43   * \f$ Q^{-1} = Q^* \f$).
     44   *
     45   * Call the function compute() to compute the Hessenberg decomposition of a
     46   * given matrix. Alternatively, you can use the
     47   * HessenbergDecomposition(const MatrixType&) constructor which computes the
     48   * Hessenberg decomposition at construction time. Once the decomposition is
     49   * computed, you can use the matrixH() and matrixQ() functions to construct
     50   * the matrices H and Q in the decomposition.
     51   *
     52   * The documentation for matrixH() contains an example of the typical use of
     53   * this class.
     54   *
     55   * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
     56   */
     57 template<typename _MatrixType> class HessenbergDecomposition
     58 {
     59   public:
     60 
     61     /** \brief Synonym for the template parameter \p _MatrixType. */
     62     typedef _MatrixType MatrixType;
     63 
     64     enum {
     65       Size = MatrixType::RowsAtCompileTime,
     66       SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
     67       Options = MatrixType::Options,
     68       MaxSize = MatrixType::MaxRowsAtCompileTime,
     69       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
     70     };
     71 
     72     /** \brief Scalar type for matrices of type #MatrixType. */
     73     typedef typename MatrixType::Scalar Scalar;
     74     typedef typename MatrixType::Index Index;
     75 
     76     /** \brief Type for vector of Householder coefficients.
     77       *
     78       * This is column vector with entries of type #Scalar. The length of the
     79       * vector is one less than the size of #MatrixType, if it is a fixed-side
     80       * type.
     81       */
     82     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
     83 
     84     /** \brief Return type of matrixQ() */
     85     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
     86 
     87     typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
     88 
     89     /** \brief Default constructor; the decomposition will be computed later.
     90       *
     91       * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
     92       *
     93       * The default constructor is useful in cases in which the user intends to
     94       * perform decompositions via compute().  The \p size parameter is only
     95       * used as a hint. It is not an error to give a wrong \p size, but it may
     96       * impair performance.
     97       *
     98       * \sa compute() for an example.
     99       */
    100     HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
    101       : m_matrix(size,size),
    102         m_temp(size),
    103         m_isInitialized(false)
    104     {
    105       if(size>1)
    106         m_hCoeffs.resize(size-1);
    107     }
    108 
    109     /** \brief Constructor; computes Hessenberg decomposition of given matrix.
    110       *
    111       * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
    112       *
    113       * This constructor calls compute() to compute the Hessenberg
    114       * decomposition.
    115       *
    116       * \sa matrixH() for an example.
    117       */
    118     HessenbergDecomposition(const MatrixType& matrix)
    119       : m_matrix(matrix),
    120         m_temp(matrix.rows()),
    121         m_isInitialized(false)
    122     {
    123       if(matrix.rows()<2)
    124       {
    125         m_isInitialized = true;
    126         return;
    127       }
    128       m_hCoeffs.resize(matrix.rows()-1,1);
    129       _compute(m_matrix, m_hCoeffs, m_temp);
    130       m_isInitialized = true;
    131     }
    132 
    133     /** \brief Computes Hessenberg decomposition of given matrix.
    134       *
    135       * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
    136       * \returns    Reference to \c *this
    137       *
    138       * The Hessenberg decomposition is computed by bringing the columns of the
    139       * matrix successively in the required form using Householder reflections
    140       * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
    141       * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
    142       * denotes the size of the given matrix.
    143       *
    144       * This method reuses of the allocated data in the HessenbergDecomposition
    145       * object.
    146       *
    147       * Example: \include HessenbergDecomposition_compute.cpp
    148       * Output: \verbinclude HessenbergDecomposition_compute.out
    149       */
    150     HessenbergDecomposition& compute(const MatrixType& matrix)
    151     {
    152       m_matrix = matrix;
    153       if(matrix.rows()<2)
    154       {
    155         m_isInitialized = true;
    156         return *this;
    157       }
    158       m_hCoeffs.resize(matrix.rows()-1,1);
    159       _compute(m_matrix, m_hCoeffs, m_temp);
    160       m_isInitialized = true;
    161       return *this;
    162     }
    163 
    164     /** \brief Returns the Householder coefficients.
    165       *
    166       * \returns a const reference to the vector of Householder coefficients
    167       *
    168       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
    169       * or the member function compute(const MatrixType&) has been called
    170       * before to compute the Hessenberg decomposition of a matrix.
    171       *
    172       * The Householder coefficients allow the reconstruction of the matrix
    173       * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
    174       *
    175       * \sa packedMatrix(), \ref Householder_Module "Householder module"
    176       */
    177     const CoeffVectorType& householderCoefficients() const
    178     {
    179       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    180       return m_hCoeffs;
    181     }
    182 
    183     /** \brief Returns the internal representation of the decomposition
    184       *
    185       *	\returns a const reference to a matrix with the internal representation
    186       *	         of the decomposition.
    187       *
    188       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
    189       * or the member function compute(const MatrixType&) has been called
    190       * before to compute the Hessenberg decomposition of a matrix.
    191       *
    192       * The returned matrix contains the following information:
    193       *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
    194       *  - the rest of the lower part contains the Householder vectors that, combined with
    195       *    Householder coefficients returned by householderCoefficients(),
    196       *    allows to reconstruct the matrix Q as
    197       *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
    198       *    Here, the matrices \f$ H_i \f$ are the Householder transformations
    199       *       \f$ H_i = (I - h_i v_i v_i^T) \f$
    200       *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
    201       *    \f$ v_i \f$ is the Householder vector defined by
    202       *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
    203       *    with M the matrix returned by this function.
    204       *
    205       * See LAPACK for further details on this packed storage.
    206       *
    207       * Example: \include HessenbergDecomposition_packedMatrix.cpp
    208       * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
    209       *
    210       * \sa householderCoefficients()
    211       */
    212     const MatrixType& packedMatrix() const
    213     {
    214       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    215       return m_matrix;
    216     }
    217 
    218     /** \brief Reconstructs the orthogonal matrix Q in the decomposition
    219       *
    220       * \returns object representing the matrix Q
    221       *
    222       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
    223       * or the member function compute(const MatrixType&) has been called
    224       * before to compute the Hessenberg decomposition of a matrix.
    225       *
    226       * This function returns a light-weight object of template class
    227       * HouseholderSequence. You can either apply it directly to a matrix or
    228       * you can convert it to a matrix of type #MatrixType.
    229       *
    230       * \sa matrixH() for an example, class HouseholderSequence
    231       */
    232     HouseholderSequenceType matrixQ() const
    233     {
    234       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    235       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
    236              .setLength(m_matrix.rows() - 1)
    237              .setShift(1);
    238     }
    239 
    240     /** \brief Constructs the Hessenberg matrix H in the decomposition
    241       *
    242       * \returns expression object representing the matrix H
    243       *
    244       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
    245       * or the member function compute(const MatrixType&) has been called
    246       * before to compute the Hessenberg decomposition of a matrix.
    247       *
    248       * The object returned by this function constructs the Hessenberg matrix H
    249       * when it is assigned to a matrix or otherwise evaluated. The matrix H is
    250       * constructed from the packed matrix as returned by packedMatrix(): The
    251       * upper part (including the subdiagonal) of the packed matrix contains
    252       * the matrix H. It may sometimes be better to directly use the packed
    253       * matrix instead of constructing the matrix H.
    254       *
    255       * Example: \include HessenbergDecomposition_matrixH.cpp
    256       * Output: \verbinclude HessenbergDecomposition_matrixH.out
    257       *
    258       * \sa matrixQ(), packedMatrix()
    259       */
    260     MatrixHReturnType matrixH() const
    261     {
    262       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    263       return MatrixHReturnType(*this);
    264     }
    265 
    266   private:
    267 
    268     typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
    269     typedef typename NumTraits<Scalar>::Real RealScalar;
    270     static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
    271 
    272   protected:
    273     MatrixType m_matrix;
    274     CoeffVectorType m_hCoeffs;
    275     VectorType m_temp;
    276     bool m_isInitialized;
    277 };
    278 
    279 /** \internal
    280   * Performs a tridiagonal decomposition of \a matA in place.
    281   *
    282   * \param matA the input selfadjoint matrix
    283   * \param hCoeffs returned Householder coefficients
    284   *
    285   * The result is written in the lower triangular part of \a matA.
    286   *
    287   * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
    288   *
    289   * \sa packedMatrix()
    290   */
    291 template<typename MatrixType>
    292 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
    293 {
    294   assert(matA.rows()==matA.cols());
    295   Index n = matA.rows();
    296   temp.resize(n);
    297   for (Index i = 0; i<n-1; ++i)
    298   {
    299     // let's consider the vector v = i-th column starting at position i+1
    300     Index remainingSize = n-i-1;
    301     RealScalar beta;
    302     Scalar h;
    303     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
    304     matA.col(i).coeffRef(i+1) = beta;
    305     hCoeffs.coeffRef(i) = h;
    306 
    307     // Apply similarity transformation to remaining columns,
    308     // i.e., compute A = H A H'
    309 
    310     // A = H A
    311     matA.bottomRightCorner(remainingSize, remainingSize)
    312         .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
    313 
    314     // A = A H'
    315     matA.rightCols(remainingSize)
    316         .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
    317   }
    318 }
    319 
    320 namespace internal {
    321 
    322 /** \eigenvalues_module \ingroup Eigenvalues_Module
    323   *
    324   *
    325   * \brief Expression type for return value of HessenbergDecomposition::matrixH()
    326   *
    327   * \tparam MatrixType type of matrix in the Hessenberg decomposition
    328   *
    329   * Objects of this type represent the Hessenberg matrix in the Hessenberg
    330   * decomposition of some matrix. The object holds a reference to the
    331   * HessenbergDecomposition class until the it is assigned or evaluated for
    332   * some other reason (the reference should remain valid during the life time
    333   * of this object). This class is the return type of
    334   * HessenbergDecomposition::matrixH(); there is probably no other use for this
    335   * class.
    336   */
    337 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
    338 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
    339 {
    340     typedef typename MatrixType::Index Index;
    341   public:
    342     /** \brief Constructor.
    343       *
    344       * \param[in] hess  Hessenberg decomposition
    345       */
    346     HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
    347 
    348     /** \brief Hessenberg matrix in decomposition.
    349       *
    350       * \param[out] result  Hessenberg matrix in decomposition \p hess which
    351       *                     was passed to the constructor
    352       */
    353     template <typename ResultType>
    354     inline void evalTo(ResultType& result) const
    355     {
    356       result = m_hess.packedMatrix();
    357       Index n = result.rows();
    358       if (n>2)
    359         result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
    360     }
    361 
    362     Index rows() const { return m_hess.packedMatrix().rows(); }
    363     Index cols() const { return m_hess.packedMatrix().cols(); }
    364 
    365   protected:
    366     const HessenbergDecomposition<MatrixType>& m_hess;
    367 };
    368 
    369 } // end namespace internal
    370 
    371 } // end namespace Eigen
    372 
    373 #endif // EIGEN_HESSENBERGDECOMPOSITION_H
    374