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 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 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 | 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 | eigen_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(), numext::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
|
---|