HessenbergDecomposition.h

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