1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 // Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
13 #define EIGEN_GENERALIZEDEIGENSOLVER_H
14
15 #include "./RealQZ.h"
16
17 namespace Eigen {
18
19 /** \eigenvalues_module \ingroup Eigenvalues_Module
20 *
21 *
22 * \class GeneralizedEigenSolver
23 *
24 * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices
25 *
26 * \tparam _MatrixType the type of the matrices of which we are computing the
27 * eigen-decomposition; this is expected to be an instantiation of the Matrix
28 * class template. Currently, only real matrices are supported.
29 *
30 * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars
31 * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If
32 * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
33 * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
34 * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
35 * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition.
36 *
37 * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the
38 * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is
39 * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$
40 * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero,
41 * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that:
42 * \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
43 * called the left eigenvector.
44 *
45 * Call the function compute() to compute the generalized eigenvalues and eigenvectors of
46 * a given matrix pair. Alternatively, you can use the
47 * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the
48 * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
49 * eigenvectors are computed, they can be retrieved with the eigenvalues() and
50 * eigenvectors() functions.
51 *
52 * Here is an usage example of this class:
53 * Example: \include GeneralizedEigenSolver.cpp
54 * Output: \verbinclude GeneralizedEigenSolver.out
55 *
56 * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
57 */
58 template<typename _MatrixType> class GeneralizedEigenSolver
59 {
60 public:
61
62 /** \brief Synonym for the template parameter \p _MatrixType. */
63 typedef _MatrixType MatrixType;
64
65 enum {
66 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68 Options = MatrixType::Options,
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71 };
72
73 /** \brief Scalar type for matrices of type #MatrixType. */
74 typedef typename MatrixType::Scalar Scalar;
75 typedef typename NumTraits<Scalar>::Real RealScalar;
76 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
77
78 /** \brief Complex scalar type for #MatrixType.
79 *
80 * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
81 * \c float or \c double) and just \c Scalar if #Scalar is
82 * complex.
83 */
84 typedef std::complex<RealScalar> ComplexScalar;
85
86 /** \brief Type for vector of real scalar values eigenvalues as returned by betas().
87 *
88 * This is a column vector with entries of type #Scalar.
89 * The length of the vector is the size of #MatrixType.
90 */
91 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
92
93 /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
94 *
95 * This is a column vector with entries of type #ComplexScalar.
96 * The length of the vector is the size of #MatrixType.
97 */
98 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;
99
100 /** \brief Expression type for the eigenvalues as returned by eigenvalues().
101 */
102 typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType;
103
104 /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
105 *
106 * This is a square matrix with entries of type #ComplexScalar.
107 * The size is the same as the size of #MatrixType.
108 */
109 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
110
111 /** \brief Default constructor.
112 *
113 * The default constructor is useful in cases in which the user intends to
114 * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
115 *
116 * \sa compute() for an example.
117 */
GeneralizedEigenSolver()118 GeneralizedEigenSolver()
119 : m_eivec(),
120 m_alphas(),
121 m_betas(),
122 m_valuesOkay(false),
123 m_vectorsOkay(false),
124 m_realQZ()
125 {}
126
127 /** \brief Default constructor with memory preallocation
128 *
129 * Like the default constructor but with preallocation of the internal data
130 * according to the specified problem \a size.
131 * \sa GeneralizedEigenSolver()
132 */
GeneralizedEigenSolver(Index size)133 explicit GeneralizedEigenSolver(Index size)
134 : m_eivec(size, size),
135 m_alphas(size),
136 m_betas(size),
137 m_valuesOkay(false),
138 m_vectorsOkay(false),
139 m_realQZ(size),
140 m_tmp(size)
141 {}
142
143 /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
144 *
145 * \param[in] A Square matrix whose eigendecomposition is to be computed.
146 * \param[in] B Square matrix whose eigendecomposition is to be computed.
147 * \param[in] computeEigenvectors If true, both the eigenvectors and the
148 * eigenvalues are computed; if false, only the eigenvalues are computed.
149 *
150 * This constructor calls compute() to compute the generalized eigenvalues
151 * and eigenvectors.
152 *
153 * \sa compute()
154 */
155 GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
156 : m_eivec(A.rows(), A.cols()),
157 m_alphas(A.cols()),
158 m_betas(A.cols()),
159 m_valuesOkay(false),
160 m_vectorsOkay(false),
161 m_realQZ(A.cols()),
162 m_tmp(A.cols())
163 {
164 compute(A, B, computeEigenvectors);
165 }
166
167 /* \brief Returns the computed generalized eigenvectors.
168 *
169 * \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
170 * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
171 *
172 * \pre Either the constructor
173 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
174 * compute(const MatrixType&, const MatrixType& bool) has been called before, and
175 * \p computeEigenvectors was set to true (the default).
176 *
177 * \sa eigenvalues()
178 */
eigenvectors()179 EigenvectorsType eigenvectors() const {
180 eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
181 return m_eivec;
182 }
183
184 /** \brief Returns an expression of the computed generalized eigenvalues.
185 *
186 * \returns An expression of the column vector containing the eigenvalues.
187 *
188 * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode
189 * Not that betas might contain zeros. It is therefore not recommended to use this function,
190 * but rather directly deal with the alphas and betas vectors.
191 *
192 * \pre Either the constructor
193 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
194 * compute(const MatrixType&,const MatrixType&,bool) has been called before.
195 *
196 * The eigenvalues are repeated according to their algebraic multiplicity,
197 * so there are as many eigenvalues as rows in the matrix. The eigenvalues
198 * are not sorted in any particular order.
199 *
200 * \sa alphas(), betas(), eigenvectors()
201 */
eigenvalues()202 EigenvalueType eigenvalues() const
203 {
204 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
205 return EigenvalueType(m_alphas,m_betas);
206 }
207
208 /** \returns A const reference to the vectors containing the alpha values
209 *
210 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
211 *
212 * \sa betas(), eigenvalues() */
alphas()213 ComplexVectorType alphas() const
214 {
215 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
216 return m_alphas;
217 }
218
219 /** \returns A const reference to the vectors containing the beta values
220 *
221 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
222 *
223 * \sa alphas(), eigenvalues() */
betas()224 VectorType betas() const
225 {
226 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
227 return m_betas;
228 }
229
230 /** \brief Computes generalized eigendecomposition of given matrix.
231 *
232 * \param[in] A Square matrix whose eigendecomposition is to be computed.
233 * \param[in] B Square matrix whose eigendecomposition is to be computed.
234 * \param[in] computeEigenvectors If true, both the eigenvectors and the
235 * eigenvalues are computed; if false, only the eigenvalues are
236 * computed.
237 * \returns Reference to \c *this
238 *
239 * This function computes the eigenvalues of the real matrix \p matrix.
240 * The eigenvalues() function can be used to retrieve them. If
241 * \p computeEigenvectors is true, then the eigenvectors are also computed
242 * and can be retrieved by calling eigenvectors().
243 *
244 * The matrix is first reduced to real generalized Schur form using the RealQZ
245 * class. The generalized Schur decomposition is then used to compute the eigenvalues
246 * and eigenvectors.
247 *
248 * The cost of the computation is dominated by the cost of the
249 * generalized Schur decomposition.
250 *
251 * This method reuses of the allocated data in the GeneralizedEigenSolver object.
252 */
253 GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
254
info()255 ComputationInfo info() const
256 {
257 eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
258 return m_realQZ.info();
259 }
260
261 /** Sets the maximal number of iterations allowed.
262 */
setMaxIterations(Index maxIters)263 GeneralizedEigenSolver& setMaxIterations(Index maxIters)
264 {
265 m_realQZ.setMaxIterations(maxIters);
266 return *this;
267 }
268
269 protected:
270
check_template_parameters()271 static void check_template_parameters()
272 {
273 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
274 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
275 }
276
277 EigenvectorsType m_eivec;
278 ComplexVectorType m_alphas;
279 VectorType m_betas;
280 bool m_valuesOkay, m_vectorsOkay;
281 RealQZ<MatrixType> m_realQZ;
282 ComplexVectorType m_tmp;
283 };
284
285 template<typename MatrixType>
286 GeneralizedEigenSolver<MatrixType>&
compute(const MatrixType & A,const MatrixType & B,bool computeEigenvectors)287 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
288 {
289 check_template_parameters();
290
291 using std::sqrt;
292 using std::abs;
293 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
294 Index size = A.cols();
295 m_valuesOkay = false;
296 m_vectorsOkay = false;
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 if (m_realQZ.info() == Success)
301 {
302 // Resize storage
303 m_alphas.resize(size);
304 m_betas.resize(size);
305 if (computeEigenvectors)
306 {
307 m_eivec.resize(size,size);
308 m_tmp.resize(size);
309 }
310
311 // Aliases:
312 Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
313 ComplexVectorType &cv = m_tmp;
314 const MatrixType &mZ = m_realQZ.matrixZ();
315 const MatrixType &mS = m_realQZ.matrixS();
316 const MatrixType &mT = m_realQZ.matrixT();
317
318 Index i = 0;
319 while (i < size)
320 {
321 if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
322 {
323 // Real eigenvalue
324 m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
325 m_betas.coeffRef(i) = mT.diagonal().coeff(i);
326 if (computeEigenvectors)
327 {
328 v.setConstant(Scalar(0.0));
329 v.coeffRef(i) = Scalar(1.0);
330 // For singular eigenvalues do nothing more
331 if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
332 {
333 // Non-singular eigenvalue
334 const Scalar alpha = real(m_alphas.coeffRef(i));
335 const Scalar beta = m_betas.coeffRef(i);
336 for (Index j = i-1; j >= 0; j--)
337 {
338 const Index st = j+1;
339 const Index sz = i-j;
340 if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
341 {
342 // 2x2 block
343 Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) );
344 Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
345 v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
346 j--;
347 }
348 else
349 {
350 v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
351 }
352 }
353 }
354 m_eivec.col(i).real().noalias() = mZ.transpose() * v;
355 m_eivec.col(i).real().normalize();
356 m_eivec.col(i).imag().setConstant(0);
357 }
358 ++i;
359 }
360 else
361 {
362 // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
363 // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
364
365 // T = [a 0]
366 // [0 b]
367 RealScalar a = mT.diagonal().coeff(i),
368 b = mT.diagonal().coeff(i+1);
369 const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
370
371 // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
372 Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
373
374 Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
375 Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
376 const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
377 m_alphas.coeffRef(i) = conj(alpha);
378 m_alphas.coeffRef(i+1) = alpha;
379
380 if (computeEigenvectors) {
381 // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
382 cv.setZero();
383 cv.coeffRef(i+1) = Scalar(1.0);
384 // here, the "static_cast" workaound expression template issues.
385 cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
386 / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
387 for (Index j = i-1; j >= 0; j--)
388 {
389 const Index st = j+1;
390 const Index sz = i+1-j;
391 if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
392 {
393 // 2x2 block
394 Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) );
395 Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
396 cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
397 j--;
398 } else {
399 cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
400 / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
401 }
402 }
403 m_eivec.col(i+1).noalias() = (mZ.transpose() * cv);
404 m_eivec.col(i+1).normalize();
405 m_eivec.col(i) = m_eivec.col(i+1).conjugate();
406 }
407 i += 2;
408 }
409 }
410
411 m_valuesOkay = true;
412 m_vectorsOkay = computeEigenvectors;
413 }
414 return *this;
415 }
416
417 } // end namespace Eigen
418
419 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
420