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 
check_template_parameters()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>&
compute(const MatrixType & A,const MatrixType & B,bool computeEigenvectors)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         Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
331         Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
332         m_alphas.coeffRef(i)   = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
333         m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
334 
335         m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
336         m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
337         i += 2;
338       }
339     }
340   }
341 
342   m_isInitialized = true;
343   m_eigenvectorsOk = false;//computeEigenvectors;
344 
345   return *this;
346 }
347 
348 } // end namespace Eigen
349 
350 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
351