src/Eigenvalues/GeneralizedEigenSolver.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
0006 // Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
0013 #define EIGEN_GENERALIZEDEIGENSOLVER_H
0014
0015 #include "./RealQZ.h"
0016
0017 namespace Eigen {
0018
0019 /** \eigenvalues_module \ingroup Eigenvalues_Module
0020   *
0021   *
0022   * \class GeneralizedEigenSolver
0023   *
0024   * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices
0025   *
0026   * \tparam _MatrixType the type of the matrices of which we are computing the
0027   * eigen-decomposition; this is expected to be an instantiation of the Matrix
0028   * class template. Currently, only real matrices are supported.
0029   *
0030   * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars
0031   * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$.  If
0032   * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
0033   * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
0034   * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
0035   * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition.
0036   *
0037   * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the
0038   * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is
0039   * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$
0040   * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero,
0041   * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that:
0042   * \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
0043   * called the left eigenvector.
0044   *
0045   * Call the function compute() to compute the generalized eigenvalues and eigenvectors of
0046   * a given matrix pair. Alternatively, you can use the
0047   * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the
0048   * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
0049   * eigenvectors are computed, they can be retrieved with the eigenvalues() and
0050   * eigenvectors() functions.
0051   *
0052   * Here is an usage example of this class:
0053   * Example: \include GeneralizedEigenSolver.cpp
0054   * Output: \verbinclude GeneralizedEigenSolver.out
0055   *
0056   * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
0057   */
0058 template<typename _MatrixType> class GeneralizedEigenSolver
0059 {
0060   public:
0061
0062     /** \brief Synonym for the template parameter \p _MatrixType. */
0063     typedef _MatrixType MatrixType;
0064
0065     enum {
0066       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0067       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0068       Options = MatrixType::Options,
0069       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0070       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0071     };
0072
0073     /** \brief Scalar type for matrices of type #MatrixType. */
0074     typedef typename MatrixType::Scalar Scalar;
0075     typedef typename NumTraits<Scalar>::Real RealScalar;
0076     typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
0077
0078     /** \brief Complex scalar type for #MatrixType.
0079       *
0080       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
0081       * \c float or \c double) and just \c Scalar if #Scalar is
0082       * complex.
0083       */
0084     typedef std::complex<RealScalar> ComplexScalar;
0085
0086     /** \brief Type for vector of real scalar values eigenvalues as returned by betas().
0087       *
0088       * This is a column vector with entries of type #Scalar.
0089       * The length of the vector is the size of #MatrixType.
0090       */
0091     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
0092
0093     /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
0094       *
0095       * This is a column vector with entries of type #ComplexScalar.
0096       * The length of the vector is the size of #MatrixType.
0097       */
0098     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;
0099
0100     /** \brief Expression type for the eigenvalues as returned by eigenvalues().
0101       */
0102     typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType;
0103
0104     /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
0105       *
0106       * This is a square matrix with entries of type #ComplexScalar.
0107       * The size is the same as the size of #MatrixType.
0108       */
0109     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
0110
0111     /** \brief Default constructor.
0112       *
0113       * The default constructor is useful in cases in which the user intends to
0114       * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
0115       *
0116       * \sa compute() for an example.
0117       */
0118     GeneralizedEigenSolver()
0119       : m_eivec(),
0120         m_alphas(),
0121         m_betas(),
0122         m_valuesOkay(false),
0123         m_vectorsOkay(false),
0124         m_realQZ()
0125     {}
0126
0127     /** \brief Default constructor with memory preallocation
0128       *
0129       * Like the default constructor but with preallocation of the internal data
0130       * according to the specified problem \a size.
0131       * \sa GeneralizedEigenSolver()
0132       */
0133     explicit GeneralizedEigenSolver(Index size)
0134       : m_eivec(size, size),
0135         m_alphas(size),
0136         m_betas(size),
0137         m_valuesOkay(false),
0138         m_vectorsOkay(false),
0139         m_realQZ(size),
0140         m_tmp(size)
0141     {}
0142
0143     /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
0144       *
0145       * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
0146       * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
0147       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
0148       *    eigenvalues are computed; if false, only the eigenvalues are computed.
0149       *
0150       * This constructor calls compute() to compute the generalized eigenvalues
0151       * and eigenvectors.
0152       *
0153       * \sa compute()
0154       */
0155     GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
0156       : m_eivec(A.rows(), A.cols()),
0157         m_alphas(A.cols()),
0158         m_betas(A.cols()),
0159         m_valuesOkay(false),
0160         m_vectorsOkay(false),
0161         m_realQZ(A.cols()),
0162         m_tmp(A.cols())
0163     {
0164       compute(A, B, computeEigenvectors);
0165     }
0166
0167     /* \brief Returns the computed generalized eigenvectors.
0168       *
0169       * \returns  %Matrix whose columns are the (possibly complex) right eigenvectors.
0170       * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
0171       *
0172       * \pre Either the constructor
0173       * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
0174       * compute(const MatrixType&, const MatrixType& bool) has been called before, and
0175       * \p computeEigenvectors was set to true (the default).
0176       *
0177       * \sa eigenvalues()
0178       */
0179     EigenvectorsType eigenvectors() const {
0180       eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
0181       return m_eivec;
0182     }
0183
0184     /** \brief Returns an expression of the computed generalized eigenvalues.
0185       *
0186       * \returns An expression of the column vector containing the eigenvalues.
0187       *
0188       * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode
0189       * Not that betas might contain zeros. It is therefore not recommended to use this function,
0190       * but rather directly deal with the alphas and betas vectors.
0191       *
0192       * \pre Either the constructor
0193       * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
0194       * compute(const MatrixType&,const MatrixType&,bool) has been called before.
0195       *
0196       * The eigenvalues are repeated according to their algebraic multiplicity,
0197       * so there are as many eigenvalues as rows in the matrix. The eigenvalues
0198       * are not sorted in any particular order.
0199       *
0200       * \sa alphas(), betas(), eigenvectors()
0201       */
0202     EigenvalueType eigenvalues() const
0203     {
0204       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
0205       return EigenvalueType(m_alphas,m_betas);
0206     }
0207
0208     /** \returns A const reference to the vectors containing the alpha values
0209       *
0210       * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
0211       *
0212       * \sa betas(), eigenvalues() */
0213     ComplexVectorType alphas() const
0214     {
0215       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
0216       return m_alphas;
0217     }
0218
0219     /** \returns A const reference to the vectors containing the beta values
0220       *
0221       * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
0222       *
0223       * \sa alphas(), eigenvalues() */
0224     VectorType betas() const
0225     {
0226       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
0227       return m_betas;
0228     }
0229
0230     /** \brief Computes generalized eigendecomposition of given matrix.
0231       *
0232       * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
0233       * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
0234       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
0235       *    eigenvalues are computed; if false, only the eigenvalues are
0236       *    computed.
0237       * \returns    Reference to \c *this
0238       *
0239       * This function computes the eigenvalues of the real matrix \p matrix.
0240       * The eigenvalues() function can be used to retrieve them.  If
0241       * \p computeEigenvectors is true, then the eigenvectors are also computed
0242       * and can be retrieved by calling eigenvectors().
0243       *
0244       * The matrix is first reduced to real generalized Schur form using the RealQZ
0245       * class. The generalized Schur decomposition is then used to compute the eigenvalues
0246       * and eigenvectors.
0247       *
0248       * The cost of the computation is dominated by the cost of the
0249       * generalized Schur decomposition.
0250       *
0251       * This method reuses of the allocated data in the GeneralizedEigenSolver object.
0252       */
0253     GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
0254
0255     ComputationInfo info() const
0256     {
0257       eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
0258       return m_realQZ.info();
0259     }
0260
0261     /** Sets the maximal number of iterations allowed.
0262     */
0263     GeneralizedEigenSolver& setMaxIterations(Index maxIters)
0264     {
0265       m_realQZ.setMaxIterations(maxIters);
0266       return *this;
0267     }
0268
0269   protected:
0270
0271     static void check_template_parameters()
0272     {
0273       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0274       EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
0275     }
0276
0277     EigenvectorsType m_eivec;
0278     ComplexVectorType m_alphas;
0279     VectorType m_betas;
0280     bool m_valuesOkay, m_vectorsOkay;
0281     RealQZ<MatrixType> m_realQZ;
0282     ComplexVectorType m_tmp;
0283 };
0284
0285 template<typename MatrixType>
0286 GeneralizedEigenSolver<MatrixType>&
0287 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
0288 {
0289   check_template_parameters();
0290
0291   using std::sqrt;
0292   using std::abs;
0293   eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
0294   Index size = A.cols();
0295   m_valuesOkay = false;
0296   m_vectorsOkay = false;
0297   // Reduce to generalized real Schur form:
0298   // A = Q S Z and B = Q T Z
0299   m_realQZ.compute(A, B, computeEigenvectors);
0300   if (m_realQZ.info() == Success)
0301   {
0302     // Resize storage
0303     m_alphas.resize(size);
0304     m_betas.resize(size);
0305     if (computeEigenvectors)
0306     {
0307       m_eivec.resize(size,size);
0308       m_tmp.resize(size);
0309     }
0310
0311     // Aliases:
0312     Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
0313     ComplexVectorType &cv = m_tmp;
0314     const MatrixType &mS = m_realQZ.matrixS();
0315     const MatrixType &mT = m_realQZ.matrixT();
0316
0317     Index i = 0;
0318     while (i < size)
0319     {
0320       if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
0321       {
0322         // Real eigenvalue
0323         m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
0324         m_betas.coeffRef(i)  = mT.diagonal().coeff(i);
0325         if (computeEigenvectors)
0326         {
0327           v.setConstant(Scalar(0.0));
0328           v.coeffRef(i) = Scalar(1.0);
0329           // For singular eigenvalues do nothing more
0330           if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
0331           {
0332             // Non-singular eigenvalue
0333             const Scalar alpha = real(m_alphas.coeffRef(i));
0334             const Scalar beta = m_betas.coeffRef(i);
0335             for (Index j = i-1; j >= 0; j--)
0336             {
0337               const Index st = j+1;
0338               const Index sz = i-j;
0339               if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
0340               {
0341                 // 2x2 block
0342                 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) );
0343                 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);
0344                 v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
0345                 j--;
0346               }
0347               else
0348               {
0349                 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));
0350               }
0351             }
0352           }
0353           m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
0354           m_eivec.col(i).real().normalize();
0355           m_eivec.col(i).imag().setConstant(0);
0356         }
0357         ++i;
0358       }
0359       else
0360       {
0361         // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
0362         // 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):
0363
0364         // T =  [a 0]
0365         //      [0 b]
0366         RealScalar a = mT.diagonal().coeff(i),
0367                    b = mT.diagonal().coeff(i+1);
0368         const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
0369
0370         // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
0371         Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
0372
0373         Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
0374         Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
0375         const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
0376         m_alphas.coeffRef(i)   = conj(alpha);
0377         m_alphas.coeffRef(i+1) = alpha;
0378
0379         if (computeEigenvectors) {
0380           // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
0381           cv.setZero();
0382           cv.coeffRef(i+1) = Scalar(1.0);
0383           // here, the "static_cast" workaound expression template issues.
0384           cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
0385                           / (static_cast<Scalar>(beta*mS.coeffRef(i,i))   - alpha*mT.coeffRef(i,i));
0386           for (Index j = i-1; j >= 0; j--)
0387           {
0388             const Index st = j+1;
0389             const Index sz = i+1-j;
0390             if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
0391             {
0392               // 2x2 block
0393               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) );
0394               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);
0395               cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
0396               j--;
0397             } else {
0398               cv.coeffRef(j) =  cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
0399                               / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
0400             }
0401           }
0402           m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
0403           m_eivec.col(i+1).normalize();
0404           m_eivec.col(i) = m_eivec.col(i+1).conjugate();
0405         }
0406         i += 2;
0407       }
0408     }
0409
0410     m_valuesOkay = true;
0411     m_vectorsOkay = computeEigenvectors;
0412   }
0413   return *this;
0414 }
0415
0416 } // end namespace Eigen
0417
0418 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H