// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012-2016 Gael Guennebaud // Copyright (C) 2010,2012 Jitse Niesen // Copyright (C) 2016 Tobias Wood // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H #define EIGEN_GENERALIZEDEIGENSOLVER_H #include "./RealQZ.h" // IWYU pragma: private #include "./InternalHeaderCheck.h" namespace Eigen { /** \eigenvalues_module \ingroup Eigenvalues_Module * * * \class GeneralizedEigenSolver * * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices * * \tparam MatrixType_ the type of the matrices of which we are computing the * eigen-decomposition; this is expected to be an instantiation of the Matrix * class template. Currently, only real matrices are supported. * * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. * * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: * \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 * called the left eigenvector. * * Call the function compute() to compute the generalized eigenvalues and eigenvectors of * a given matrix pair. Alternatively, you can use the * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the * eigenvalues and eigenvectors at construction time. Once the eigenvalue and * eigenvectors are computed, they can be retrieved with the eigenvalues() and * eigenvectors() functions. * * Here is an usage example of this class: * Example: \include GeneralizedEigenSolver.cpp * Output: \verbinclude GeneralizedEigenSolver.out * * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver */ template class GeneralizedEigenSolver { public: /** \brief Synonym for the template parameter \p MatrixType_. */ typedef MatrixType_ MatrixType; enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = internal::traits::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 /** \brief Complex scalar type for #MatrixType. * * This is \c std::complex if #Scalar is real (e.g., * \c float or \c double) and just \c Scalar if #Scalar is * complex. */ typedef std::complex ComplexScalar; /** \brief Type for vector of real scalar values eigenvalues as returned by betas(). * * This is a column vector with entries of type #Scalar. * The length of the vector is the size of #MatrixType. */ typedef Matrix VectorType; /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). * * This is a column vector with entries of type #ComplexScalar. * The length of the vector is the size of #MatrixType. */ typedef Matrix ComplexVectorType; /** \brief Expression type for the eigenvalues as returned by eigenvalues(). */ typedef CwiseBinaryOp, ComplexVectorType, VectorType> EigenvalueType; /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). * * This is a square matrix with entries of type #ComplexScalar. * The size is the same as the size of #MatrixType. */ typedef Matrix EigenvectorsType; /** \brief Default constructor. * * The default constructor is useful in cases in which the user intends to * perform decompositions via EigenSolver::compute(const MatrixType&, bool). * * \sa compute() for an example. */ GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_computeEigenvectors(false), m_isInitialized(false), m_realQZ() {} /** \brief Default constructor with memory preallocation * * Like the default constructor but with preallocation of the internal data * according to the specified problem \a size. * \sa GeneralizedEigenSolver() */ explicit GeneralizedEigenSolver(Index size) : m_eivec(size, size), m_alphas(size), m_betas(size), m_computeEigenvectors(false), m_isInitialized(false), m_realQZ(size), m_tmp(size) {} /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. * * \param[in] A Square matrix whose eigendecomposition is to be computed. * \param[in] B Square matrix whose eigendecomposition is to be computed. * \param[in] computeEigenvectors If true, both the eigenvectors and the * eigenvalues are computed; if false, only the eigenvalues are computed. * * This constructor calls compute() to compute the generalized eigenvalues * and eigenvectors. * * \sa compute() */ GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) : m_eivec(A.rows(), A.cols()), m_alphas(A.cols()), m_betas(A.cols()), m_computeEigenvectors(false), m_isInitialized(false), m_realQZ(A.cols()), m_tmp(A.cols()) { compute(A, B, computeEigenvectors); } /** \brief Returns the computed generalized eigenvectors. * * \returns %Matrix whose columns are the (possibly complex) right eigenvectors. * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. * * \pre Either the constructor * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function * compute(const MatrixType&, const MatrixType& bool) has been called before, and * \p computeEigenvectors was set to true (the default). * * \sa eigenvalues() */ EigenvectorsType eigenvectors() const { eigen_assert(info() == Success && "GeneralizedEigenSolver failed to compute eigenvectors"); eigen_assert(m_computeEigenvectors && "Eigenvectors for GeneralizedEigenSolver were not calculated"); return m_eivec; } /** \brief Returns an expression of the computed generalized eigenvalues. * * \returns An expression of the column vector containing the eigenvalues. * * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode * Not that betas might contain zeros. It is therefore not recommended to use this function, * but rather directly deal with the alphas and betas vectors. * * \pre Either the constructor * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function * compute(const MatrixType&,const MatrixType&,bool) has been called before. * * The eigenvalues are repeated according to their algebraic multiplicity, * so there are as many eigenvalues as rows in the matrix. The eigenvalues * are not sorted in any particular order. * * \sa alphas(), betas(), eigenvectors() */ EigenvalueType eigenvalues() const { eigen_assert(info() == Success && "GeneralizedEigenSolver failed to compute eigenvalues."); return EigenvalueType(m_alphas, m_betas); } /** \returns A const reference to the vectors containing the alpha values * * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). * * \sa betas(), eigenvalues() */ const ComplexVectorType& alphas() const { eigen_assert(info() == Success && "GeneralizedEigenSolver failed to compute alphas."); return m_alphas; } /** \returns A const reference to the vectors containing the beta values * * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). * * \sa alphas(), eigenvalues() */ const VectorType& betas() const { eigen_assert(info() == Success && "GeneralizedEigenSolver failed to compute betas."); return m_betas; } /** \brief Computes generalized eigendecomposition of given matrix. * * \param[in] A Square matrix whose eigendecomposition is to be computed. * \param[in] B Square matrix whose eigendecomposition is to be computed. * \param[in] computeEigenvectors If true, both the eigenvectors and the * eigenvalues are computed; if false, only the eigenvalues are * computed. * \returns Reference to \c *this * * This function computes the eigenvalues of the real matrix \p matrix. * The eigenvalues() function can be used to retrieve them. If * \p computeEigenvectors is true, then the eigenvectors are also computed * and can be retrieved by calling eigenvectors(). * * The matrix is first reduced to real generalized Schur form using the RealQZ * class. The generalized Schur decomposition is then used to compute the eigenvalues * and eigenvectors. * * The cost of the computation is dominated by the cost of the * generalized Schur decomposition. * * This method reuses of the allocated data in the GeneralizedEigenSolver object. */ GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); ComputationInfo info() const { eigen_assert(m_isInitialized && "EigenSolver is not initialized."); return m_realQZ.info(); } /** Sets the maximal number of iterations allowed. */ GeneralizedEigenSolver& setMaxIterations(Index maxIters) { m_realQZ.setMaxIterations(maxIters); return *this; } protected: EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) EIGEN_STATIC_ASSERT(!NumTraits::IsComplex, NUMERIC_TYPE_MUST_BE_REAL) EigenvectorsType m_eivec; ComplexVectorType m_alphas; VectorType m_betas; bool m_computeEigenvectors; bool m_isInitialized; RealQZ m_realQZ; ComplexVectorType m_tmp; }; template GeneralizedEigenSolver& GeneralizedEigenSolver::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) { using std::abs; using std::sqrt; eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); Index size = A.cols(); // Reduce to generalized real Schur form: // A = Q S Z and B = Q T Z m_realQZ.compute(A, B, computeEigenvectors); if (m_realQZ.info() == Success) { // Resize storage m_alphas.resize(size); m_betas.resize(size); if (computeEigenvectors) { m_eivec.resize(size, size); m_tmp.resize(size); } // Aliases: Map v(reinterpret_cast(m_tmp.data()), size); ComplexVectorType& cv = m_tmp; const MatrixType& mS = m_realQZ.matrixS(); const MatrixType& mT = m_realQZ.matrixT(); Index i = 0; while (i < size) { if (i == size - 1 || mS.coeff(i + 1, i) == Scalar(0)) { // Real eigenvalue m_alphas.coeffRef(i) = mS.diagonal().coeff(i); m_betas.coeffRef(i) = mT.diagonal().coeff(i); if (computeEigenvectors) { v.setConstant(Scalar(0.0)); v.coeffRef(i) = Scalar(1.0); // For singular eigenvalues do nothing more if (abs(m_betas.coeffRef(i)) >= (std::numeric_limits::min)()) { // Non-singular eigenvalue const Scalar alpha = real(m_alphas.coeffRef(i)); const Scalar beta = m_betas.coeffRef(i); for (Index j = i - 1; j >= 0; j--) { const Index st = j + 1; const Index sz = i - j; if (j > 0 && mS.coeff(j, j - 1) != Scalar(0)) { // 2x2 block Matrix 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)); Matrix lhs = beta * mS.template block<2, 2>(j - 1, j - 1) - alpha * mT.template block<2, 2>(j - 1, j - 1); v.template segment<2>(j - 1) = lhs.partialPivLu().solve(rhs); j--; } else { 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)); } } } m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v; m_eivec.col(i).real().normalize(); m_eivec.col(i).imag().setConstant(0); } ++i; } else { // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal // 2x2 block T 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): // T = [a 0] // [0 b] RealScalar a = mT.diagonal().coeff(i), b = mT.diagonal().coeff(i + 1); const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i + 1) = a * b; // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. Matrix S2 = mS.template block<2, 2>(i, i) * Matrix(b, a).asDiagonal(); Scalar p = Scalar(0.5) * (S2.coeff(0, 0) - S2.coeff(1, 1)); Scalar z = sqrt(abs(p * p + S2.coeff(1, 0) * S2.coeff(0, 1))); const ComplexScalar alpha = ComplexScalar(S2.coeff(1, 1) + p, (beta > 0) ? z : -z); m_alphas.coeffRef(i) = conj(alpha); m_alphas.coeffRef(i + 1) = alpha; if (computeEigenvectors) { // Compute eigenvector in position (i+1) and then position (i) is just the conjugate cv.setZero(); cv.coeffRef(i + 1) = Scalar(1.0); // here, the "static_cast" workaround expression template issues. cv.coeffRef(i) = -(static_cast(beta * mS.coeffRef(i, i + 1)) - alpha * mT.coeffRef(i, i + 1)) / (static_cast(beta * mS.coeffRef(i, i)) - alpha * mT.coeffRef(i, i)); for (Index j = i - 1; j >= 0; j--) { const Index st = j + 1; const Index sz = i + 1 - j; if (j > 0 && mS.coeff(j, j - 1) != Scalar(0)) { // 2x2 block Matrix 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)); Matrix lhs = beta * mS.template block<2, 2>(j - 1, j - 1) - alpha * mT.template block<2, 2>(j - 1, j - 1); cv.template segment<2>(j - 1) = lhs.partialPivLu().solve(rhs); j--; } else { cv.coeffRef(j) = cv.segment(st, sz) .transpose() .cwiseProduct(beta * mS.block(j, st, 1, sz) - alpha * mT.block(j, st, 1, sz)) .sum() / (alpha * mT.coeffRef(j, j) - static_cast(beta * mS.coeffRef(j, j))); } } m_eivec.col(i + 1).noalias() = (m_realQZ.matrixZ().transpose() * cv); m_eivec.col(i + 1).normalize(); m_eivec.col(i) = m_eivec.col(i + 1).conjugate(); } i += 2; } } } m_computeEigenvectors = computeEigenvectors; m_isInitialized = true; return *this; } } // end namespace Eigen #endif // EIGEN_GENERALIZEDEIGENSOLVER_H