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 */
GeneralizedEigenSolver()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 */
GeneralizedEigenSolver(Index size)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 */
eigenvalues()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() */
alphas()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() */
betas()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
info()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 */
setMaxIterations(Index maxIters)259 GeneralizedEigenSolver& setMaxIterations(Index maxIters)
260 {
261 m_realQZ.setMaxIterations(maxIters);
262 return *this;
263 }
264
265 protected:
266 MatrixType m_eivec;
267 ComplexVectorType m_alphas;
268 VectorType m_betas;
269 bool m_isInitialized;
270 bool m_eigenvectorsOk;
271 RealQZ<MatrixType> m_realQZ;
272 MatrixType m_matS;
273
274 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
275 ColumnVectorType m_tmp;
276 };
277
278 //template<typename MatrixType>
279 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
280 //{
281 // eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
282 // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
283 // Index n = m_eivec.cols();
284 // EigenvectorsType matV(n,n);
285 // // TODO
286 // return matV;
287 //}
288
289 template<typename MatrixType>
290 GeneralizedEigenSolver<MatrixType>&
compute(const MatrixType & A,const MatrixType & B,bool computeEigenvectors)291 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
292 {
293 using std::sqrt;
294 using std::abs;
295 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
296
297 // Reduce to generalized real Schur form:
298 // A = Q S Z and B = Q T Z
299 m_realQZ.compute(A, B, computeEigenvectors);
300
301 if (m_realQZ.info() == Success)
302 {
303 m_matS = m_realQZ.matrixS();
304 if (computeEigenvectors)
305 m_eivec = m_realQZ.matrixZ().transpose();
306
307 // Compute eigenvalues from matS
308 m_alphas.resize(A.cols());
309 m_betas.resize(A.cols());
310 Index i = 0;
311 while (i < A.cols())
312 {
313 if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
314 {
315 m_alphas.coeffRef(i) = m_matS.coeff(i, i);
316 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
317 ++i;
318 }
319 else
320 {
321 Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
322 Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
323 m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
324 m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
325
326 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
327 m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
328 i += 2;
329 }
330 }
331 }
332
333 m_isInitialized = true;
334 m_eigenvectorsOk = false;//computeEigenvectors;
335
336 return *this;
337 }
338
339 } // end namespace Eigen
340
341 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
342