1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
5 | // Copyright (C) 2010,2012 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_GENERALIZEDEIGENSOLVER_H
|
---|
12 | #define EIGEN_GENERALIZEDEIGENSOLVER_H
|
---|
13 |
|
---|
14 | #include "./RealQZ.h"
|
---|
15 |
|
---|
16 | namespace Eigen {
|
---|
17 |
|
---|
18 | /** \eigenvalues_module \ingroup Eigenvalues_Module
|
---|
19 | *
|
---|
20 | *
|
---|
21 | * \class GeneralizedEigenSolver
|
---|
22 | *
|
---|
23 | * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices
|
---|
24 | *
|
---|
25 | * \tparam _MatrixType the type of the matrices of which we are computing the
|
---|
26 | * eigen-decomposition; this is expected to be an instantiation of the Matrix
|
---|
27 | * class template. Currently, only real matrices are supported.
|
---|
28 | *
|
---|
29 | * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars
|
---|
30 | * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If
|
---|
31 | * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
|
---|
32 | * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
|
---|
33 | * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
|
---|
34 | * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition.
|
---|
35 | *
|
---|
36 | * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the
|
---|
37 | * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is
|
---|
38 | * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$
|
---|
39 | * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero,
|
---|
40 | * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that:
|
---|
41 | * \f$ \mu_i A v_i = B v_i \f$, or even \f$ \mu_i u_i^T A = u_i^T B \f$ where \f$ u_i \f$ is
|
---|
42 | * called the left eigenvector.
|
---|
43 | *
|
---|
44 | * Call the function compute() to compute the generalized eigenvalues and eigenvectors of
|
---|
45 | * a given matrix pair. Alternatively, you can use the
|
---|
46 | * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the
|
---|
47 | * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
|
---|
48 | * eigenvectors are computed, they can be retrieved with the eigenvalues() and
|
---|
49 | * eigenvectors() functions.
|
---|
50 | *
|
---|
51 | * Here is an usage example of this class:
|
---|
52 | * Example: \include GeneralizedEigenSolver.cpp
|
---|
53 | * Output: \verbinclude GeneralizedEigenSolver.out
|
---|
54 | *
|
---|
55 | * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
|
---|
56 | */
|
---|
57 | template<typename _MatrixType> class GeneralizedEigenSolver
|
---|
58 | {
|
---|
59 | public:
|
---|
60 |
|
---|
61 | /** \brief Synonym for the template parameter \p _MatrixType. */
|
---|
62 | typedef _MatrixType MatrixType;
|
---|
63 |
|
---|
64 | enum {
|
---|
65 | RowsAtCompileTime = MatrixType::RowsAtCompileTime,
|
---|
66 | ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
---|
67 | Options = MatrixType::Options,
|
---|
68 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
|
---|
69 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
|
---|
70 | };
|
---|
71 |
|
---|
72 | /** \brief Scalar type for matrices of type #MatrixType. */
|
---|
73 | typedef typename MatrixType::Scalar Scalar;
|
---|
74 | typedef typename NumTraits<Scalar>::Real RealScalar;
|
---|
75 | typedef typename MatrixType::Index Index;
|
---|
76 |
|
---|
77 | /** \brief Complex scalar type for #MatrixType.
|
---|
78 | *
|
---|
79 | * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
|
---|
80 | * \c float or \c double) and just \c Scalar if #Scalar is
|
---|
81 | * complex.
|
---|
82 | */
|
---|
83 | typedef std::complex<RealScalar> ComplexScalar;
|
---|
84 |
|
---|
85 | /** \brief Type for vector of real scalar values eigenvalues as returned by betas().
|
---|
86 | *
|
---|
87 | * This is a column vector with entries of type #Scalar.
|
---|
88 | * The length of the vector is the size of #MatrixType.
|
---|
89 | */
|
---|
90 | typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
|
---|
91 |
|
---|
92 | /** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
|
---|
93 | *
|
---|
94 | * This is a column vector with entries of type #ComplexScalar.
|
---|
95 | * The length of the vector is the size of #MatrixType.
|
---|
96 | */
|
---|
97 | typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;
|
---|
98 |
|
---|
99 | /** \brief Expression type for the eigenvalues as returned by eigenvalues().
|
---|
100 | */
|
---|
101 | typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType;
|
---|
102 |
|
---|
103 | /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
|
---|
104 | *
|
---|
105 | * This is a square matrix with entries of type #ComplexScalar.
|
---|
106 | * The size is the same as the size of #MatrixType.
|
---|
107 | */
|
---|
108 | typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
|
---|
109 |
|
---|
110 | /** \brief Default constructor.
|
---|
111 | *
|
---|
112 | * The default constructor is useful in cases in which the user intends to
|
---|
113 | * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
|
---|
114 | *
|
---|
115 | * \sa compute() for an example.
|
---|
116 | */
|
---|
117 | GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
|
---|
118 |
|
---|
119 | /** \brief Default constructor with memory preallocation
|
---|
120 | *
|
---|
121 | * Like the default constructor but with preallocation of the internal data
|
---|
122 | * according to the specified problem \a size.
|
---|
123 | * \sa GeneralizedEigenSolver()
|
---|
124 | */
|
---|
125 | GeneralizedEigenSolver(Index size)
|
---|
126 | : m_eivec(size, size),
|
---|
127 | m_alphas(size),
|
---|
128 | m_betas(size),
|
---|
129 | m_isInitialized(false),
|
---|
130 | m_eigenvectorsOk(false),
|
---|
131 | m_realQZ(size),
|
---|
132 | m_matS(size, size),
|
---|
133 | m_tmp(size)
|
---|
134 | {}
|
---|
135 |
|
---|
136 | /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
|
---|
137 | *
|
---|
138 | * \param[in] A Square matrix whose eigendecomposition is to be computed.
|
---|
139 | * \param[in] B Square matrix whose eigendecomposition is to be computed.
|
---|
140 | * \param[in] computeEigenvectors If true, both the eigenvectors and the
|
---|
141 | * eigenvalues are computed; if false, only the eigenvalues are computed.
|
---|
142 | *
|
---|
143 | * This constructor calls compute() to compute the generalized eigenvalues
|
---|
144 | * and eigenvectors.
|
---|
145 | *
|
---|
146 | * \sa compute()
|
---|
147 | */
|
---|
148 | GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
|
---|
149 | : m_eivec(A.rows(), A.cols()),
|
---|
150 | m_alphas(A.cols()),
|
---|
151 | m_betas(A.cols()),
|
---|
152 | m_isInitialized(false),
|
---|
153 | m_eigenvectorsOk(false),
|
---|
154 | m_realQZ(A.cols()),
|
---|
155 | m_matS(A.rows(), A.cols()),
|
---|
156 | m_tmp(A.cols())
|
---|
157 | {
|
---|
158 | compute(A, B, computeEigenvectors);
|
---|
159 | }
|
---|
160 |
|
---|
161 | /* \brief Returns the computed generalized eigenvectors.
|
---|
162 | *
|
---|
163 | * \returns %Matrix whose columns are the (possibly complex) eigenvectors.
|
---|
164 | *
|
---|
165 | * \pre Either the constructor
|
---|
166 | * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
|
---|
167 | * compute(const MatrixType&, const MatrixType& bool) has been called before, and
|
---|
168 | * \p computeEigenvectors was set to true (the default).
|
---|
169 | *
|
---|
170 | * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
|
---|
171 | * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
|
---|
172 | * eigenvectors are normalized to have (Euclidean) norm equal to one. The
|
---|
173 | * matrix returned by this function is the matrix \f$ V \f$ in the
|
---|
174 | * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
|
---|
175 | *
|
---|
176 | * \sa eigenvalues()
|
---|
177 | */
|
---|
178 | // EigenvectorsType eigenvectors() const;
|
---|
179 |
|
---|
180 | /** \brief Returns an expression of the computed generalized eigenvalues.
|
---|
181 | *
|
---|
182 | * \returns An expression of the column vector containing the eigenvalues.
|
---|
183 | *
|
---|
184 | * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode
|
---|
185 | * Not that betas might contain zeros. It is therefore not recommended to use this function,
|
---|
186 | * but rather directly deal with the alphas and betas vectors.
|
---|
187 | *
|
---|
188 | * \pre Either the constructor
|
---|
189 | * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
|
---|
190 | * compute(const MatrixType&,const MatrixType&,bool) has been called before.
|
---|
191 | *
|
---|
192 | * The eigenvalues are repeated according to their algebraic multiplicity,
|
---|
193 | * so there are as many eigenvalues as rows in the matrix. The eigenvalues
|
---|
194 | * are not sorted in any particular order.
|
---|
195 | *
|
---|
196 | * \sa alphas(), betas(), eigenvectors()
|
---|
197 | */
|
---|
198 | EigenvalueType eigenvalues() const
|
---|
199 | {
|
---|
200 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
|
---|
201 | return EigenvalueType(m_alphas,m_betas);
|
---|
202 | }
|
---|
203 |
|
---|
204 | /** \returns A const reference to the vectors containing the alpha values
|
---|
205 | *
|
---|
206 | * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
|
---|
207 | *
|
---|
208 | * \sa betas(), eigenvalues() */
|
---|
209 | ComplexVectorType alphas() const
|
---|
210 | {
|
---|
211 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
|
---|
212 | return m_alphas;
|
---|
213 | }
|
---|
214 |
|
---|
215 | /** \returns A const reference to the vectors containing the beta values
|
---|
216 | *
|
---|
217 | * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
|
---|
218 | *
|
---|
219 | * \sa alphas(), eigenvalues() */
|
---|
220 | VectorType betas() const
|
---|
221 | {
|
---|
222 | eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
|
---|
223 | return m_betas;
|
---|
224 | }
|
---|
225 |
|
---|
226 | /** \brief Computes generalized eigendecomposition of given matrix.
|
---|
227 | *
|
---|
228 | * \param[in] A Square matrix whose eigendecomposition is to be computed.
|
---|
229 | * \param[in] B Square matrix whose eigendecomposition is to be computed.
|
---|
230 | * \param[in] computeEigenvectors If true, both the eigenvectors and the
|
---|
231 | * eigenvalues are computed; if false, only the eigenvalues are
|
---|
232 | * computed.
|
---|
233 | * \returns Reference to \c *this
|
---|
234 | *
|
---|
235 | * This function computes the eigenvalues of the real matrix \p matrix.
|
---|
236 | * The eigenvalues() function can be used to retrieve them. If
|
---|
237 | * \p computeEigenvectors is true, then the eigenvectors are also computed
|
---|
238 | * and can be retrieved by calling eigenvectors().
|
---|
239 | *
|
---|
240 | * The matrix is first reduced to real generalized Schur form using the RealQZ
|
---|
241 | * class. The generalized Schur decomposition is then used to compute the eigenvalues
|
---|
242 | * and eigenvectors.
|
---|
243 | *
|
---|
244 | * The cost of the computation is dominated by the cost of the
|
---|
245 | * generalized Schur decomposition.
|
---|
246 | *
|
---|
247 | * This method reuses of the allocated data in the GeneralizedEigenSolver object.
|
---|
248 | */
|
---|
249 | GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
|
---|
250 |
|
---|
251 | ComputationInfo info() const
|
---|
252 | {
|
---|
253 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
254 | return m_realQZ.info();
|
---|
255 | }
|
---|
256 |
|
---|
257 | /** Sets the maximal number of iterations allowed.
|
---|
258 | */
|
---|
259 | GeneralizedEigenSolver& setMaxIterations(Index maxIters)
|
---|
260 | {
|
---|
261 | m_realQZ.setMaxIterations(maxIters);
|
---|
262 | return *this;
|
---|
263 | }
|
---|
264 |
|
---|
265 | protected:
|
---|
266 |
|
---|
267 | static void check_template_parameters()
|
---|
268 | {
|
---|
269 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
|
---|
270 | EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
|
---|
271 | }
|
---|
272 |
|
---|
273 | MatrixType m_eivec;
|
---|
274 | ComplexVectorType m_alphas;
|
---|
275 | VectorType m_betas;
|
---|
276 | bool m_isInitialized;
|
---|
277 | bool m_eigenvectorsOk;
|
---|
278 | RealQZ<MatrixType> m_realQZ;
|
---|
279 | MatrixType m_matS;
|
---|
280 |
|
---|
281 | typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
|
---|
282 | ColumnVectorType m_tmp;
|
---|
283 | };
|
---|
284 |
|
---|
285 | //template<typename MatrixType>
|
---|
286 | //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
|
---|
287 | //{
|
---|
288 | // eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
289 | // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
---|
290 | // Index n = m_eivec.cols();
|
---|
291 | // EigenvectorsType matV(n,n);
|
---|
292 | // // TODO
|
---|
293 | // return matV;
|
---|
294 | //}
|
---|
295 |
|
---|
296 | template<typename MatrixType>
|
---|
297 | GeneralizedEigenSolver<MatrixType>&
|
---|
298 | GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
|
---|
299 | {
|
---|
300 | check_template_parameters();
|
---|
301 |
|
---|
302 | using std::sqrt;
|
---|
303 | using std::abs;
|
---|
304 | eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
|
---|
305 |
|
---|
306 | // Reduce to generalized real Schur form:
|
---|
307 | // A = Q S Z and B = Q T Z
|
---|
308 | m_realQZ.compute(A, B, computeEigenvectors);
|
---|
309 |
|
---|
310 | if (m_realQZ.info() == Success)
|
---|
311 | {
|
---|
312 | m_matS = m_realQZ.matrixS();
|
---|
313 | if (computeEigenvectors)
|
---|
314 | m_eivec = m_realQZ.matrixZ().transpose();
|
---|
315 |
|
---|
316 | // Compute eigenvalues from matS
|
---|
317 | m_alphas.resize(A.cols());
|
---|
318 | m_betas.resize(A.cols());
|
---|
319 | Index i = 0;
|
---|
320 | while (i < A.cols())
|
---|
321 | {
|
---|
322 | if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
|
---|
323 | {
|
---|
324 | m_alphas.coeffRef(i) = m_matS.coeff(i, i);
|
---|
325 | m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
|
---|
326 | ++i;
|
---|
327 | }
|
---|
328 | else
|
---|
329 | {
|
---|
330 | // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a triangular 2x2 block T
|
---|
331 | // From the eigen decomposition of T = U * E * U^-1,
|
---|
332 | // we can extract the eigenvalues of (U^-1 * S * U) / E
|
---|
333 | // Here, we can take advantage that E = diag(T), and U = [ 1 T_01 ; 0 T_11-T_00], and U^-1 = [1 -T_11/(T_11-T_00) ; 0 1/(T_11-T_00)].
|
---|
334 | // Then taking beta=T_00*T_11*(T_11-T_00), we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00) * (T_11-T_00):
|
---|
335 |
|
---|
336 | // T = [a b ; 0 c]
|
---|
337 | // S = [e f ; g h]
|
---|
338 | RealScalar a = m_realQZ.matrixT().coeff(i, i), b = m_realQZ.matrixT().coeff(i, i+1), c = m_realQZ.matrixT().coeff(i+1, i+1);
|
---|
339 | RealScalar e = m_matS.coeff(i, i), f = m_matS.coeff(i, i+1), g = m_matS.coeff(i+1, i), h = m_matS.coeff(i+1, i+1);
|
---|
340 | RealScalar d = c-a;
|
---|
341 | RealScalar gb = g*b;
|
---|
342 | Matrix<RealScalar,2,2> A;
|
---|
343 | A << (e*d-gb)*c, ((e*b+f*d-h*b)*d-gb*b)*a,
|
---|
344 | g*c , (gb+h*d)*a;
|
---|
345 |
|
---|
346 | // NOTE, we could also compute the SVD of T's block during the QZ factorization so that the respective T block is guaranteed to be diagonal,
|
---|
347 | // and then we could directly apply the formula below (while taking care of scaling S columns by T11,T00):
|
---|
348 |
|
---|
349 | Scalar p = Scalar(0.5) * (A.coeff(i, i) - A.coeff(i+1, i+1));
|
---|
350 | Scalar z = sqrt(abs(p * p + A.coeff(i+1, i) * A.coeff(i, i+1)));
|
---|
351 | m_alphas.coeffRef(i) = ComplexScalar(A.coeff(i+1, i+1) + p, z);
|
---|
352 | m_alphas.coeffRef(i+1) = ComplexScalar(A.coeff(i+1, i+1) + p, -z);
|
---|
353 |
|
---|
354 | m_betas.coeffRef(i) =
|
---|
355 | m_betas.coeffRef(i+1) = a*c*d;
|
---|
356 |
|
---|
357 | i += 2;
|
---|
358 | }
|
---|
359 | }
|
---|
360 | }
|
---|
361 |
|
---|
362 | m_isInitialized = true;
|
---|
363 | m_eigenvectorsOk = false;//computeEigenvectors;
|
---|
364 |
|
---|
365 | return *this;
|
---|
366 | }
|
---|
367 |
|
---|
368 | } // end namespace Eigen
|
---|
369 |
|
---|
370 | #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
|
---|