Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
GeneralizedEigenSolver.h
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 #include "./InternalHeaderCheck.h"
18 
19 namespace Eigen {
20 
60 template<typename MatrixType_> class GeneralizedEigenSolver
61 {
62  public:
63 
65  typedef MatrixType_ MatrixType;
66 
67  enum {
68  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
69  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
70  Options = MatrixType::Options,
71  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
72  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
73  };
74 
76  typedef typename MatrixType::Scalar Scalar;
77  typedef typename NumTraits<Scalar>::Real RealScalar;
78  typedef Eigen::Index Index;
79 
86  typedef std::complex<RealScalar> ComplexScalar;
87 
94 
101 
105 
112 
121  : m_eivec(),
122  m_alphas(),
123  m_betas(),
124  m_valuesOkay(false),
125  m_vectorsOkay(false),
126  m_realQZ()
127  {}
128 
136  : m_eivec(size, size),
137  m_alphas(size),
138  m_betas(size),
139  m_valuesOkay(false),
140  m_vectorsOkay(false),
141  m_realQZ(size),
142  m_tmp(size)
143  {}
144 
157  GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
158  : m_eivec(A.rows(), A.cols()),
159  m_alphas(A.cols()),
160  m_betas(A.cols()),
161  m_valuesOkay(false),
162  m_vectorsOkay(false),
163  m_realQZ(A.cols()),
164  m_tmp(A.cols())
165  {
166  compute(A, B, computeEigenvectors);
167  }
168 
169  /* \brief Returns the computed generalized eigenvectors.
170  *
171  * \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
172  * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
173  *
174  * \pre Either the constructor
175  * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
176  * compute(const MatrixType&, const MatrixType& bool) has been called before, and
177  * \p computeEigenvectors was set to true (the default).
178  *
179  * \sa eigenvalues()
180  */
181  EigenvectorsType eigenvectors() const {
182  eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
183  return m_eivec;
184  }
185 
205  {
206  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
207  return EigenvalueType(m_alphas,m_betas);
208  }
209 
215  const ComplexVectorType& alphas() const
216  {
217  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
218  return m_alphas;
219  }
220 
226  const VectorType& betas() const
227  {
228  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
229  return m_betas;
230  }
231 
255  GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
256 
257  ComputationInfo info() const
258  {
259  eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
260  return m_realQZ.info();
261  }
262 
266  {
267  m_realQZ.setMaxIterations(maxIters);
268  return *this;
269  }
270 
271  protected:
272 
273  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
274  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL)
275 
276  EigenvectorsType m_eivec;
277  ComplexVectorType m_alphas;
278  VectorType m_betas;
279  bool m_valuesOkay, m_vectorsOkay;
280  RealQZ<MatrixType> m_realQZ;
281  ComplexVectorType m_tmp;
282 };
283 
284 template<typename MatrixType>
285 GeneralizedEigenSolver<MatrixType>&
286 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
287 {
288  using std::sqrt;
289  using std::abs;
290  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
291  Index size = A.cols();
292  m_valuesOkay = false;
293  m_vectorsOkay = false;
294  // Reduce to generalized real Schur form:
295  // A = Q S Z and B = Q T Z
296  m_realQZ.compute(A, B, computeEigenvectors);
297  if (m_realQZ.info() == Success)
298  {
299  // Resize storage
300  m_alphas.resize(size);
301  m_betas.resize(size);
302  if (computeEigenvectors)
303  {
304  m_eivec.resize(size,size);
305  m_tmp.resize(size);
306  }
307 
308  // Aliases:
309  Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
310  ComplexVectorType &cv = m_tmp;
311  const MatrixType &mS = m_realQZ.matrixS();
312  const MatrixType &mT = m_realQZ.matrixT();
313 
314  Index i = 0;
315  while (i < size)
316  {
317  if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
318  {
319  // Real eigenvalue
320  m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
321  m_betas.coeffRef(i) = mT.diagonal().coeff(i);
322  if (computeEigenvectors)
323  {
324  v.setConstant(Scalar(0.0));
325  v.coeffRef(i) = Scalar(1.0);
326  // For singular eigenvalues do nothing more
327  if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
328  {
329  // Non-singular eigenvalue
330  const Scalar alpha = real(m_alphas.coeffRef(i));
331  const Scalar beta = m_betas.coeffRef(i);
332  for (Index j = i-1; j >= 0; j--)
333  {
334  const Index st = j+1;
335  const Index sz = i-j;
336  if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
337  {
338  // 2x2 block
339  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) );
340  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);
341  v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
342  j--;
343  }
344  else
345  {
346  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));
347  }
348  }
349  }
350  m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
351  m_eivec.col(i).real().normalize();
352  m_eivec.col(i).imag().setConstant(0);
353  }
354  ++i;
355  }
356  else
357  {
358  // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
359  // 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):
360 
361  // T = [a 0]
362  // [0 b]
363  RealScalar a = mT.diagonal().coeff(i),
364  b = mT.diagonal().coeff(i+1);
365  const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
366 
367  // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
368  Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
369 
370  Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
371  Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
372  const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
373  m_alphas.coeffRef(i) = conj(alpha);
374  m_alphas.coeffRef(i+1) = alpha;
375 
376  if (computeEigenvectors) {
377  // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
378  cv.setZero();
379  cv.coeffRef(i+1) = Scalar(1.0);
380  // here, the "static_cast" workaound expression template issues.
381  cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
382  / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
383  for (Index j = i-1; j >= 0; j--)
384  {
385  const Index st = j+1;
386  const Index sz = i+1-j;
387  if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
388  {
389  // 2x2 block
390  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) );
391  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);
392  cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
393  j--;
394  } else {
395  cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
396  / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
397  }
398  }
399  m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
400  m_eivec.col(i+1).normalize();
401  m_eivec.col(i) = m_eivec.col(i+1).conjugate();
402  }
403  i += 2;
404  }
405  }
406 
407  m_valuesOkay = true;
408  m_vectorsOkay = computeEigenvectors;
409  }
410  return *this;
411 }
412 
413 } // end namespace Eigen
414 
415 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:86
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:61
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:286
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:265
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:86
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by alphas().
Definition: GeneralizedEigenSolver.h:100
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:93
const VectorType & betas() const
Definition: GeneralizedEigenSolver.h:226
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:104
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:76
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:120
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:78
const ComplexVectorType & alphas() const
Definition: GeneralizedEigenSolver.h:215
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:157
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:204
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:135
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:111
MatrixType_ MatrixType
Synonym for the template parameter MatrixType_.
Definition: GeneralizedEigenSolver.h:65
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:327
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:599
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:164
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:187
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:564
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:188
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:171
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231