GeneralizedEigenSolver.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
00012 #define EIGEN_GENERALIZEDEIGENSOLVER_H
00013 
00014 #include "./RealQZ.h"
00015 
00016 namespace Eigen { 
00017 
00018 /** \eigenvalues_module \ingroup Eigenvalues_Module
00019   *
00020   *
00021   * \class GeneralizedEigenSolver
00022   *
00023   * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices
00024   *
00025   * \tparam _MatrixType the type of the matrices of which we are computing the
00026   * eigen-decomposition; this is expected to be an instantiation of the Matrix
00027   * class template. Currently, only real matrices are supported.
00028   *
00029   * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars
00030   * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$.  If
00031   * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
00032   * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
00033   * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
00034   * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition.
00035   *
00036   * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the
00037   * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is
00038   * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$
00039   * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero,
00040   * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that:
00041   * \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
00042   * called the left eigenvector.
00043   *
00044   * Call the function compute() to compute the generalized eigenvalues and eigenvectors of
00045   * a given matrix pair. Alternatively, you can use the
00046   * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the
00047   * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
00048   * eigenvectors are computed, they can be retrieved with the eigenvalues() and
00049   * eigenvectors() functions.
00050   *
00051   * Here is an usage example of this class:
00052   * Example: \include GeneralizedEigenSolver.cpp
00053   * Output: \verbinclude GeneralizedEigenSolver.out
00054   *
00055   * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
00056   */
00057 template<typename _MatrixType> class GeneralizedEigenSolver 
00058 {
00059   public:
00060 
00061     /** \brief Synonym for the template parameter \p _MatrixType. */
00062     typedef _MatrixType MatrixType;
00063 
00064     enum {
00065       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00066       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00067       Options = MatrixType::Options,
00068       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00069       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00070     };
00071 
00072     /** \brief Scalar type for matrices of type #MatrixType. */
00073     typedef typename MatrixType::Scalar Scalar;
00074     typedef typename NumTraits<Scalar>::Real RealScalar;
00075     typedef typename MatrixType::Index Index;
00076 
00077     /** \brief Complex scalar type for #MatrixType. 
00078       *
00079       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
00080       * \c float or \c double) and just \c Scalar if #Scalar is
00081       * complex.
00082       */
00083     typedef std::complex<RealScalar> ComplexScalar;
00084 
00085     /** \brief Type for vector of real scalar values eigenvalues as returned by betas().
00086       *
00087       * This is a column vector with entries of type #Scalar.
00088       * The length of the vector is the size of #MatrixType.
00089       */
00090     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1>  VectorType;
00091 
00092     /** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
00093       *
00094       * This is a column vector with entries of type #ComplexScalar.
00095       * The length of the vector is the size of #MatrixType.
00096       */
00097     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1>  ComplexVectorType;
00098 
00099     /** \brief Expression type for the eigenvalues as returned by eigenvalues().
00100       */
00101     typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType > EigenvalueType;
00102 
00103     /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 
00104       *
00105       * This is a square matrix with entries of type #ComplexScalar. 
00106       * The size is the same as the size of #MatrixType.
00107       */
00108     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime>  EigenvectorsType;
00109 
00110     /** \brief Default constructor.
00111       *
00112       * The default constructor is useful in cases in which the user intends to
00113       * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
00114       *
00115       * \sa compute() for an example.
00116       */
00117     GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
00118 
00119     /** \brief Default constructor with memory preallocation
00120       *
00121       * Like the default constructor but with preallocation of the internal data
00122       * according to the specified problem \a size.
00123       * \sa GeneralizedEigenSolver()
00124       */
00125     GeneralizedEigenSolver(Index size)
00126       : m_eivec(size, size),
00127         m_alphas(size),
00128         m_betas(size),
00129         m_isInitialized(false),
00130         m_eigenvectorsOk(false),
00131         m_realQZ(size),
00132         m_matS(size, size),
00133         m_tmp(size)
00134     {}
00135 
00136     /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
00137       * 
00138       * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
00139       * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
00140       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00141       *    eigenvalues are computed; if false, only the eigenvalues are computed.
00142       *
00143       * This constructor calls compute() to compute the generalized eigenvalues
00144       * and eigenvectors.
00145       *
00146       * \sa compute()
00147       */
00148     GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
00149       : m_eivec(A.rows(), A.cols()),
00150         m_alphas(A.cols()),
00151         m_betas(A.cols()),
00152         m_isInitialized(false),
00153         m_eigenvectorsOk(false),
00154         m_realQZ(A.cols()),
00155         m_matS(A.rows(), A.cols()),
00156         m_tmp(A.cols())
00157     {
00158       compute(A, B, computeEigenvectors);
00159     }
00160 
00161     /* \brief Returns the computed generalized eigenvectors.
00162       *
00163       * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
00164       *
00165       * \pre Either the constructor 
00166       * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
00167       * compute(const MatrixType&, const MatrixType& bool) has been called before, and
00168       * \p computeEigenvectors was set to true (the default).
00169       *
00170       * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
00171       * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
00172       * eigenvectors are normalized to have (Euclidean) norm equal to one. The
00173       * matrix returned by this function is the matrix \f$ V \f$ in the
00174       * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
00175       *
00176       * \sa eigenvalues()
00177       */
00178 //    EigenvectorsType eigenvectors() const;
00179 
00180     /** \brief Returns an expression of the computed generalized eigenvalues.
00181       *
00182       * \returns An expression of the column vector containing the eigenvalues.
00183       *
00184       * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode
00185       * Not that betas might contain zeros. It is therefore not recommended to use this function,
00186       * but rather directly deal with the alphas and betas vectors.
00187       *
00188       * \pre Either the constructor 
00189       * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
00190       * compute(const MatrixType&,const MatrixType&,bool) has been called before.
00191       *
00192       * The eigenvalues are repeated according to their algebraic multiplicity,
00193       * so there are as many eigenvalues as rows in the matrix. The eigenvalues 
00194       * are not sorted in any particular order.
00195       *
00196       * \sa alphas(), betas(), eigenvectors()
00197       */
00198     EigenvalueType eigenvalues() const
00199     {
00200       eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
00201       return EigenvalueType(m_alphas,m_betas);
00202     }
00203 
00204     /** \returns A const reference to the vectors containing the alpha values
00205       *
00206       * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
00207       *
00208       * \sa betas(), eigenvalues() */
00209     ComplexVectorType  alphas () const
00210     {
00211       eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
00212       return m_alphas;
00213     }
00214 
00215     /** \returns A const reference to the vectors containing the beta values
00216       *
00217       * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
00218       *
00219       * \sa alphas(), eigenvalues() */
00220     VectorType  betas () const
00221     {
00222       eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
00223       return m_betas;
00224     }
00225 
00226     /** \brief Computes generalized eigendecomposition of given matrix.
00227       * 
00228       * \param[in]  A  Square matrix whose eigendecomposition is to be computed.
00229       * \param[in]  B  Square matrix whose eigendecomposition is to be computed.
00230       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00231       *    eigenvalues are computed; if false, only the eigenvalues are
00232       *    computed. 
00233       * \returns    Reference to \c *this
00234       *
00235       * This function computes the eigenvalues of the real matrix \p matrix.
00236       * The eigenvalues() function can be used to retrieve them.  If 
00237       * \p computeEigenvectors is true, then the eigenvectors are also computed
00238       * and can be retrieved by calling eigenvectors().
00239       *
00240       * The matrix is first reduced to real generalized Schur form using the RealQZ
00241       * class. The generalized Schur decomposition is then used to compute the eigenvalues
00242       * and eigenvectors.
00243       *
00244       * The cost of the computation is dominated by the cost of the
00245       * generalized Schur decomposition.
00246       *
00247       * This method reuses of the allocated data in the GeneralizedEigenSolver object.
00248       */
00249     GeneralizedEigenSolver & compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
00250 
00251     ComputationInfo info() const
00252     {
00253       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00254       return m_realQZ.info();
00255     }
00256 
00257     /** Sets the maximal number of iterations allowed.
00258     */
00259     GeneralizedEigenSolver & setMaxIterations(Index maxIters)
00260     {
00261       m_realQZ.setMaxIterations(maxIters);
00262       return *this;
00263     }
00264 
00265   protected:
00266     
00267     static void check_template_parameters()
00268     {
00269       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00270       EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
00271     }
00272     
00273     MatrixType m_eivec;
00274     ComplexVectorType m_alphas;
00275     VectorType m_betas;
00276     bool m_isInitialized;
00277     bool m_eigenvectorsOk;
00278     RealQZ<MatrixType> m_realQZ;
00279     MatrixType m_matS;
00280 
00281     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00282     ColumnVectorType m_tmp;
00283 };
00284 
00285 //template<typename MatrixType>
00286 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
00287 //{
00288 //  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00289 //  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00290 //  Index n = m_eivec.cols();
00291 //  EigenvectorsType matV(n,n);
00292 //  // TODO
00293 //  return matV;
00294 //}
00295 
00296 template<typename MatrixType>
00297 GeneralizedEigenSolver<MatrixType>&
00298 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
00299 {
00300   check_template_parameters();
00301   
00302   using std::sqrt;
00303   using std::abs;
00304   eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
00305 
00306   // Reduce to generalized real Schur form:
00307   // A = Q S Z and B = Q T Z
00308   m_realQZ.compute(A, B, computeEigenvectors);
00309 
00310   if (m_realQZ.info() == Success)
00311   {
00312     m_matS = m_realQZ.matrixS();
00313     if (computeEigenvectors)
00314       m_eivec = m_realQZ.matrixZ().transpose();
00315   
00316     // Compute eigenvalues from matS
00317     m_alphas.resize(A.cols());
00318     m_betas.resize(A.cols());
00319     Index i = 0;
00320     while (i < A.cols())
00321     {
00322       if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
00323       {
00324         m_alphas.coeffRef(i) = m_matS.coeff(i, i);
00325         m_betas.coeffRef(i)  = m_realQZ.matrixT().coeff(i,i);
00326         ++i;
00327       }
00328       else
00329       {
00330         Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
00331         Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
00332         m_alphas.coeffRef(i)   = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
00333         m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
00334 
00335         m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
00336         m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
00337         i += 2;
00338       }
00339     }
00340   }
00341 
00342   m_isInitialized = true;
00343   m_eigenvectorsOk = false;//computeEigenvectors;
00344 
00345   return *this;
00346 }
00347 
00348 } // end namespace Eigen
00349 
00350 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H