WPILibC++ 2023.4.3-108-ge5452e3
HessenbergDecomposition.h
Go to the documentation of this file.
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
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
19template<typename 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 */
57template<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,
67 Options = MatrixType::Options,
68 MaxSize = MatrixType::MaxRowsAtCompileTime,
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 */
83
84 /** \brief Return type of matrixQ() */
86
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 */
101 : m_matrix(size,size),
102 m_temp(size),
103 m_isInitialized(false)
104 {
105 if(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>
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>
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 */
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 */
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 */
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 */
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, int(Options) | int(RowMajor), 1, MaxSize> VectorType;
271 typedef typename NumTraits<Scalar>::Real RealScalar;
272 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
273
274 protected:
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 */
293template<typename MatrixType>
294void 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), numext::conj(h), &temp.coeffRef(0));
319 }
320}
321
322namespace 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 */
339template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
340: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
341{
342 public:
343 /** \brief Constructor.
344 *
345 * \param[in] hess Hessenberg decomposition
346 */
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:
368};
369
370} // end namespace internal
371
372} // end namespace Eigen
373
374#endif // EIGEN_HESSENBERGDECOMPOSITION_H
#define eigen_assert(x)
Definition: Macros.h:1047
constexpr common_return_t< T1, T2 > beta(const T1 a, const T2 b) noexcept
Compile-time beta function.
Definition: beta.hpp:36
\eigenvalues_module
Definition: HessenbergDecomposition.h:58
HessenbergDecomposition & compute(const EigenBase< InputType > &matrix)
Computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:152
HouseholderSequenceType matrixQ() const
Reconstructs the orthogonal matrix Q in the decomposition.
Definition: HessenbergDecomposition.h:234
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: HessenbergDecomposition.h:214
bool m_isInitialized
Definition: HessenbergDecomposition.h:278
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
Type for vector of Householder coefficients.
Definition: HessenbergDecomposition.h:82
VectorType m_temp
Definition: HessenbergDecomposition.h:277
CoeffVectorType m_hCoeffs
Definition: HessenbergDecomposition.h:276
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: HessenbergDecomposition.h:85
Eigen::Index Index
Definition: HessenbergDecomposition.h:74
MatrixHReturnType matrixH() const
Constructs the Hessenberg matrix H in the decomposition.
Definition: HessenbergDecomposition.h:262
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: HessenbergDecomposition.h:62
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: HessenbergDecomposition.h:73
MatrixType m_matrix
Definition: HessenbergDecomposition.h:275
internal::HessenbergDecompositionMatrixHReturnType< MatrixType > MatrixHReturnType
Definition: HessenbergDecomposition.h:87
HessenbergDecomposition(const EigenBase< InputType > &matrix)
Constructor; computes Hessenberg decomposition of given matrix.
Definition: HessenbergDecomposition.h:119
@ MaxSizeMinusOne
Definition: HessenbergDecomposition.h:69
@ Size
Definition: HessenbergDecomposition.h:65
@ Options
Definition: HessenbergDecomposition.h:67
@ SizeMinusOne
Definition: HessenbergDecomposition.h:66
@ MaxSize
Definition: HessenbergDecomposition.h:68
HessenbergDecomposition(Index size=Size==Dynamic ? 2 :Size)
Default constructor; the decomposition will be computed later.
Definition: HessenbergDecomposition.h:100
const CoeffVectorType & householderCoefficients() const
Returns the Householder coefficients.
Definition: HessenbergDecomposition.h:179
\householder_module
Definition: HouseholderSequence.h:121
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Resizes *this to a rows x cols matrix.
Definition: PlainObjectBase.h:271
Definition: ReturnByValue.h:52
@ RowMajor
Storage order is row major (see TopicStorageOrders).
Definition: Constants.h:321
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
This value means that a positive quantity (e.g., a size) is not known at compile-time,...
Definition: Constants.h:22
result
Definition: format.h:2564
Definition: Eigen_Colamd.h:50
static constexpr const unit_t< compound_unit< energy::joule, time::seconds > > h(6.626070040e-34)
Planck constant.
Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor Matrix...
Definition: EigenBase.h:30
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Holds information about the various numeric (i.e.
Definition: NumTraits.h:233
\eigenvalues_module
Definition: HessenbergDecomposition.h:341
Index cols() const
Definition: HessenbergDecomposition.h:364
void evalTo(ResultType &result) const
Hessenberg matrix in decomposition.
Definition: HessenbergDecomposition.h:355
HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition< MatrixType > &hess)
Constructor.
Definition: HessenbergDecomposition.h:347
const HessenbergDecomposition< MatrixType > & m_hess
Definition: HessenbergDecomposition.h:367
Index rows() const
Definition: HessenbergDecomposition.h:363
Definition: ForwardDeclarations.h:17