GeneralizedEigenSolver.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// 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"
 
namespace Eigen { 

template<typename _MatrixType> class GeneralizedEigenSolver
{
  public:

    typedef _MatrixType MatrixType;
 
    enum {
      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
      Options = MatrixType::Options,
      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
    };

    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef typename MatrixType::Index Index;

    typedef std::complex<RealScalar> ComplexScalar;

    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;

    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;

    typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType;

    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;

    GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}

    GeneralizedEigenSolver(Index size)
      : m_eivec(size, size),
        m_alphas(size),
        m_betas(size),
        m_isInitialized(false),
        m_eigenvectorsOk(false),
        m_realQZ(size),
        m_matS(size, size),
        m_tmp(size)
    {}

    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_isInitialized(false),
        m_eigenvectorsOk(false),
        m_realQZ(A.cols()),
        m_matS(A.rows(), A.cols()),
        m_tmp(A.cols())
    {
      compute(A, B, computeEigenvectors);
    }
 
    /* \brief Returns the computed generalized eigenvectors.
      *
      * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
      *
      * \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).
      *
      * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
      * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
      * eigenvectors are normalized to have (Euclidean) norm equal to one. The
      * matrix returned by this function is the matrix \f$ V \f$ in the
      * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
      *
      * \sa eigenvalues()
      */
//    EigenvectorsType eigenvectors() const;

    EigenvalueType eigenvalues() const
    {
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      return EigenvalueType(m_alphas,m_betas);
    }

    ComplexVectorType alphas() const
    {
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      return m_alphas;
    }

    VectorType betas() const
    {
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      return m_betas;
    }

    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();
    }

    GeneralizedEigenSolver& setMaxIterations(Index maxIters)
    {
      m_realQZ.setMaxIterations(maxIters);
      return *this;
    }
 
  protected:
    
    static void check_template_parameters()
    {
      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
    }
    
    MatrixType m_eivec;
    ComplexVectorType m_alphas;
    VectorType m_betas;
    bool m_isInitialized;
    bool m_eigenvectorsOk;
    RealQZ<MatrixType> m_realQZ;
    MatrixType m_matS;
 
    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
    ColumnVectorType m_tmp;
};
 
//template<typename MatrixType>
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
//{
//  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
//  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
//  Index n = m_eivec.cols();
//  EigenvectorsType matV(n,n);
//  // TODO
//  return matV;
//}
 
template<typename MatrixType>
GeneralizedEigenSolver<MatrixType>&
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
{
  check_template_parameters();
  
  using std::sqrt;
  using std::abs;
  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
 
  // 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)
  {
    m_matS = m_realQZ.matrixS();
    if (computeEigenvectors)
      m_eivec = m_realQZ.matrixZ().transpose();
  
    // Compute eigenvalues from matS
    m_alphas.resize(A.cols());
    m_betas.resize(A.cols());
    Index i = 0;
    while (i < A.cols())
    {
      if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
      {
        m_alphas.coeffRef(i) = m_matS.coeff(i, i);
        m_betas.coeffRef(i)  = m_realQZ.matrixT().coeff(i,i);
        ++i;
      }
      else
      {
        // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a triangular 2x2 block T
        // From the eigen decomposition of T = U * E * U^-1,
        // we can extract the eigenvalues of (U^-1 * S * U) / E
        // Here, we can take advantage that E = diag(T), and U = [ 1 T_01 ; 0 T_11-T_00], and U^-1 = [1 -T_11/(T_11-T_00) ; 0 1/(T_11-T_00)].
        // Then taking beta=T_00*T_11*(T_11-T_00), we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00) * (T_11-T_00):
 
        // T =  [a b ; 0 c]
        // S =  [e f ; g h]
        RealScalar a = m_realQZ.matrixT().coeff(i, i), b = m_realQZ.matrixT().coeff(i, i+1), c = m_realQZ.matrixT().coeff(i+1, i+1);
        RealScalar e = m_matS.coeff(i, i), f = m_matS.coeff(i, i+1), g = m_matS.coeff(i+1, i), h = m_matS.coeff(i+1, i+1);
        RealScalar d = c-a;
        RealScalar gb = g*b;
        Matrix<RealScalar,2,2> A;
        A <<  (e*d-gb)*c,  ((e*b+f*d-h*b)*d-gb*b)*a,
              g*c       ,  (gb+h*d)*a;
 
        // NOTE, we could also compute the SVD of T's block during the QZ factorization so that the respective T block is guaranteed to be diagonal,
        // and then we could directly apply the formula below (while taking care of scaling S columns by T11,T00):
 
        Scalar p = Scalar(0.5) * (A.coeff(i, i) - A.coeff(i+1, i+1));
        Scalar z = sqrt(abs(p * p + A.coeff(i+1, i) * A.coeff(i, i+1)));
        m_alphas.coeffRef(i)   = ComplexScalar(A.coeff(i+1, i+1) + p, z);
        m_alphas.coeffRef(i+1) = ComplexScalar(A.coeff(i+1, i+1) + p, -z);
 
        m_betas.coeffRef(i)   =
        m_betas.coeffRef(i+1) = a*c*d;
 
        i += 2;
      }
    }
  }
 
  m_isInitialized = true;
  m_eigenvectorsOk = false;//computeEigenvectors;
 
  return *this;
}
 
} // end namespace Eigen
 
#endif // EIGEN_GENERALIZEDEIGENSOLVER_H