• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@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 Eigen::Index Index; ///< \deprecated since Eigen 3.3
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 HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> 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     explicit 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     template<typename InputType>
119     explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
120       : m_matrix(matrix.derived()),
121         m_temp(matrix.rows()),
122         m_isInitialized(false)
123     {
124       if(matrix.rows()<2)
125       {
126         m_isInitialized = true;
127         return;
128       }
129       m_hCoeffs.resize(matrix.rows()-1,1);
130       _compute(m_matrix, m_hCoeffs, m_temp);
131       m_isInitialized = true;
132     }
133 
134     /** \brief Computes Hessenberg decomposition of given matrix.
135       *
136       * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
137       * \returns    Reference to \c *this
138       *
139       * The Hessenberg decomposition is computed by bringing the columns of the
140       * matrix successively in the required form using Householder reflections
141       * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
142       * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
143       * denotes the size of the given matrix.
144       *
145       * This method reuses of the allocated data in the HessenbergDecomposition
146       * object.
147       *
148       * Example: \include HessenbergDecomposition_compute.cpp
149       * Output: \verbinclude HessenbergDecomposition_compute.out
150       */
151     template<typename InputType>
152     HessenbergDecomposition& compute(const EigenBase<InputType>& matrix)
153     {
154       m_matrix = matrix.derived();
155       if(matrix.rows()<2)
156       {
157         m_isInitialized = true;
158         return *this;
159       }
160       m_hCoeffs.resize(matrix.rows()-1,1);
161       _compute(m_matrix, m_hCoeffs, m_temp);
162       m_isInitialized = true;
163       return *this;
164     }
165 
166     /** \brief Returns the Householder coefficients.
167       *
168       * \returns a const reference to the vector of Householder coefficients
169       *
170       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
171       * or the member function compute(const MatrixType&) has been called
172       * before to compute the Hessenberg decomposition of a matrix.
173       *
174       * The Householder coefficients allow the reconstruction of the matrix
175       * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
176       *
177       * \sa packedMatrix(), \ref Householder_Module "Householder module"
178       */
179     const CoeffVectorType& householderCoefficients() const
180     {
181       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
182       return m_hCoeffs;
183     }
184 
185     /** \brief Returns the internal representation of the decomposition
186       *
187       *	\returns a const reference to a matrix with the internal representation
188       *	         of the decomposition.
189       *
190       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
191       * or the member function compute(const MatrixType&) has been called
192       * before to compute the Hessenberg decomposition of a matrix.
193       *
194       * The returned matrix contains the following information:
195       *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
196       *  - the rest of the lower part contains the Householder vectors that, combined with
197       *    Householder coefficients returned by householderCoefficients(),
198       *    allows to reconstruct the matrix Q as
199       *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
200       *    Here, the matrices \f$ H_i \f$ are the Householder transformations
201       *       \f$ H_i = (I - h_i v_i v_i^T) \f$
202       *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
203       *    \f$ v_i \f$ is the Householder vector defined by
204       *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
205       *    with M the matrix returned by this function.
206       *
207       * See LAPACK for further details on this packed storage.
208       *
209       * Example: \include HessenbergDecomposition_packedMatrix.cpp
210       * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
211       *
212       * \sa householderCoefficients()
213       */
214     const MatrixType& packedMatrix() const
215     {
216       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
217       return m_matrix;
218     }
219 
220     /** \brief Reconstructs the orthogonal matrix Q in the decomposition
221       *
222       * \returns object representing the matrix Q
223       *
224       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
225       * or the member function compute(const MatrixType&) has been called
226       * before to compute the Hessenberg decomposition of a matrix.
227       *
228       * This function returns a light-weight object of template class
229       * HouseholderSequence. You can either apply it directly to a matrix or
230       * you can convert it to a matrix of type #MatrixType.
231       *
232       * \sa matrixH() for an example, class HouseholderSequence
233       */
234     HouseholderSequenceType matrixQ() const
235     {
236       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
237       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
238              .setLength(m_matrix.rows() - 1)
239              .setShift(1);
240     }
241 
242     /** \brief Constructs the Hessenberg matrix H in the decomposition
243       *
244       * \returns expression object representing the matrix H
245       *
246       * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
247       * or the member function compute(const MatrixType&) has been called
248       * before to compute the Hessenberg decomposition of a matrix.
249       *
250       * The object returned by this function constructs the Hessenberg matrix H
251       * when it is assigned to a matrix or otherwise evaluated. The matrix H is
252       * constructed from the packed matrix as returned by packedMatrix(): The
253       * upper part (including the subdiagonal) of the packed matrix contains
254       * the matrix H. It may sometimes be better to directly use the packed
255       * matrix instead of constructing the matrix H.
256       *
257       * Example: \include HessenbergDecomposition_matrixH.cpp
258       * Output: \verbinclude HessenbergDecomposition_matrixH.out
259       *
260       * \sa matrixQ(), packedMatrix()
261       */
262     MatrixHReturnType matrixH() const
263     {
264       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
265       return MatrixHReturnType(*this);
266     }
267 
268   private:
269 
270     typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
271     typedef typename NumTraits<Scalar>::Real RealScalar;
272     static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
273 
274   protected:
275     MatrixType m_matrix;
276     CoeffVectorType m_hCoeffs;
277     VectorType m_temp;
278     bool m_isInitialized;
279 };
280 
281 /** \internal
282   * Performs a tridiagonal decomposition of \a matA in place.
283   *
284   * \param matA the input selfadjoint matrix
285   * \param hCoeffs returned Householder coefficients
286   *
287   * The result is written in the lower triangular part of \a matA.
288   *
289   * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
290   *
291   * \sa packedMatrix()
292   */
293 template<typename MatrixType>
294 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
295 {
296   eigen_assert(matA.rows()==matA.cols());
297   Index n = matA.rows();
298   temp.resize(n);
299   for (Index i = 0; i<n-1; ++i)
300   {
301     // let's consider the vector v = i-th column starting at position i+1
302     Index remainingSize = n-i-1;
303     RealScalar beta;
304     Scalar h;
305     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
306     matA.col(i).coeffRef(i+1) = beta;
307     hCoeffs.coeffRef(i) = h;
308 
309     // Apply similarity transformation to remaining columns,
310     // i.e., compute A = H A H'
311 
312     // A = H A
313     matA.bottomRightCorner(remainingSize, remainingSize)
314         .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
315 
316     // A = A H'
317     matA.rightCols(remainingSize)
318         .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
319   }
320 }
321 
322 namespace internal {
323 
324 /** \eigenvalues_module \ingroup Eigenvalues_Module
325   *
326   *
327   * \brief Expression type for return value of HessenbergDecomposition::matrixH()
328   *
329   * \tparam MatrixType type of matrix in the Hessenberg decomposition
330   *
331   * Objects of this type represent the Hessenberg matrix in the Hessenberg
332   * decomposition of some matrix. The object holds a reference to the
333   * HessenbergDecomposition class until the it is assigned or evaluated for
334   * some other reason (the reference should remain valid during the life time
335   * of this object). This class is the return type of
336   * HessenbergDecomposition::matrixH(); there is probably no other use for this
337   * class.
338   */
339 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
340 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
341 {
342   public:
343     /** \brief Constructor.
344       *
345       * \param[in] hess  Hessenberg decomposition
346       */
347     HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
348 
349     /** \brief Hessenberg matrix in decomposition.
350       *
351       * \param[out] result  Hessenberg matrix in decomposition \p hess which
352       *                     was passed to the constructor
353       */
354     template <typename ResultType>
355     inline void evalTo(ResultType& result) const
356     {
357       result = m_hess.packedMatrix();
358       Index n = result.rows();
359       if (n>2)
360         result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
361     }
362 
363     Index rows() const { return m_hess.packedMatrix().rows(); }
364     Index cols() const { return m_hess.packedMatrix().cols(); }
365 
366   protected:
367     const HessenbergDecomposition<MatrixType>& m_hess;
368 };
369 
370 } // end namespace internal
371 
372 } // end namespace Eigen
373 
374 #endif // EIGEN_HESSENBERGDECOMPOSITION_H
375