GeneralizedEigenSolver.h 13.6 KB
Newer Older
Don Gagne's avatar
Don Gagne committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
// 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 { 

/** \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<typename _MatrixType> class GeneralizedEigenSolver
{
  public:

    /** \brief Synonym for the template parameter \p _MatrixType. */
    typedef _MatrixType MatrixType;

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

    /** \brief Scalar type for matrices of type #MatrixType. */
    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef typename MatrixType::Index Index;

    /** \brief Complex scalar type for #MatrixType. 
      *
      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
      * \c float or \c double) and just \c Scalar if #Scalar is
      * complex.
      */
    typedef std::complex<RealScalar> 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<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;

    /** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
      *
      * This is a column vector with entries of type #ComplexScalar.
      * The length of the vector is the size of #MatrixType.
      */
    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;

    /** \brief Expression type for the eigenvalues as returned by eigenvalues().
      */
    typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,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<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> 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_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}

    /** \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()
      */
    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)
    {}

    /** \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_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;

    /** \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(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      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() */
    ComplexVectorType alphas() const
    {
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      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() */
    VectorType betas() const
    {
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
      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:
266 267 268 269 270 271 272
    
    static void check_template_parameters()
    {
      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
    }
    
Don Gagne's avatar
Don Gagne committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
    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)
{
300 301
  check_template_parameters();
  
Don Gagne's avatar
Don Gagne committed
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
  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
      {
        Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
        Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
        m_alphas.coeffRef(i)   = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
        m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);

        m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
        m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
        i += 2;
      }
    }
  }

  m_isInitialized = true;
  m_eigenvectorsOk = false;//computeEigenvectors;

  return *this;
}

} // end namespace Eigen

#endif // EIGEN_GENERALIZEDEIGENSOLVER_H