1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2010,2012 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H 12 #define EIGEN_GENERALIZEDEIGENSOLVER_H 13 14 #include "./RealQZ.h" 15 16 namespace Eigen { 17 18 /** \eigenvalues_module \ingroup Eigenvalues_Module 19 * 20 * 21 * \class GeneralizedEigenSolver 22 * 23 * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices 24 * 25 * \tparam _MatrixType the type of the matrices of which we are computing the 26 * eigen-decomposition; this is expected to be an instantiation of the Matrix 27 * class template. Currently, only real matrices are supported. 28 * 29 * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars 30 * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If 31 * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and 32 * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = 33 * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we 34 * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. 35 * 36 * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the 37 * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is 38 * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ 39 * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, 40 * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: 41 * \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 42 * called the left eigenvector. 43 * 44 * Call the function compute() to compute the generalized eigenvalues and eigenvectors of 45 * a given matrix pair. Alternatively, you can use the 46 * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the 47 * eigenvalues and eigenvectors at construction time. Once the eigenvalue and 48 * eigenvectors are computed, they can be retrieved with the eigenvalues() and 49 * eigenvectors() functions. 50 * 51 * Here is an usage example of this class: 52 * Example: \include GeneralizedEigenSolver.cpp 53 * Output: \verbinclude GeneralizedEigenSolver.out 54 * 55 * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver 56 */ 57 template<typename _MatrixType> class GeneralizedEigenSolver 58 { 59 public: 60 61 /** \brief Synonym for the template parameter \p _MatrixType. */ 62 typedef _MatrixType MatrixType; 63 64 enum { 65 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 66 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 67 Options = MatrixType::Options, 68 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 69 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 70 }; 71 72 /** \brief Scalar type for matrices of type #MatrixType. */ 73 typedef typename MatrixType::Scalar Scalar; 74 typedef typename NumTraits<Scalar>::Real RealScalar; 75 typedef typename MatrixType::Index Index; 76 77 /** \brief Complex scalar type for #MatrixType. 78 * 79 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 80 * \c float or \c double) and just \c Scalar if #Scalar is 81 * complex. 82 */ 83 typedef std::complex<RealScalar> ComplexScalar; 84 85 /** \brief Type for vector of real scalar values eigenvalues as returned by betas(). 86 * 87 * This is a column vector with entries of type #Scalar. 88 * The length of the vector is the size of #MatrixType. 89 */ 90 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; 91 92 /** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). 93 * 94 * This is a column vector with entries of type #ComplexScalar. 95 * The length of the vector is the size of #MatrixType. 96 */ 97 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; 98 99 /** \brief Expression type for the eigenvalues as returned by eigenvalues(). 100 */ 101 typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; 102 103 /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 104 * 105 * This is a square matrix with entries of type #ComplexScalar. 106 * The size is the same as the size of #MatrixType. 107 */ 108 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 109 110 /** \brief Default constructor. 111 * 112 * The default constructor is useful in cases in which the user intends to 113 * perform decompositions via EigenSolver::compute(const MatrixType&, bool). 114 * 115 * \sa compute() for an example. 116 */ 117 GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} 118 119 /** \brief Default constructor with memory preallocation 120 * 121 * Like the default constructor but with preallocation of the internal data 122 * according to the specified problem \a size. 123 * \sa GeneralizedEigenSolver() 124 */ 125 GeneralizedEigenSolver(Index size) 126 : m_eivec(size, size), 127 m_alphas(size), 128 m_betas(size), 129 m_isInitialized(false), 130 m_eigenvectorsOk(false), 131 m_realQZ(size), 132 m_matS(size, size), 133 m_tmp(size) 134 {} 135 136 /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. 137 * 138 * \param[in] A Square matrix whose eigendecomposition is to be computed. 139 * \param[in] B Square matrix whose eigendecomposition is to be computed. 140 * \param[in] computeEigenvectors If true, both the eigenvectors and the 141 * eigenvalues are computed; if false, only the eigenvalues are computed. 142 * 143 * This constructor calls compute() to compute the generalized eigenvalues 144 * and eigenvectors. 145 * 146 * \sa compute() 147 */ 148 GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) 149 : m_eivec(A.rows(), A.cols()), 150 m_alphas(A.cols()), 151 m_betas(A.cols()), 152 m_isInitialized(false), 153 m_eigenvectorsOk(false), 154 m_realQZ(A.cols()), 155 m_matS(A.rows(), A.cols()), 156 m_tmp(A.cols()) 157 { 158 compute(A, B, computeEigenvectors); 159 } 160 161 /* \brief Returns the computed generalized eigenvectors. 162 * 163 * \returns %Matrix whose columns are the (possibly complex) eigenvectors. 164 * 165 * \pre Either the constructor 166 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function 167 * compute(const MatrixType&, const MatrixType& bool) has been called before, and 168 * \p computeEigenvectors was set to true (the default). 169 * 170 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 171 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 172 * eigenvectors are normalized to have (Euclidean) norm equal to one. The 173 * matrix returned by this function is the matrix \f$ V \f$ in the 174 * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. 175 * 176 * \sa eigenvalues() 177 */ 178 // EigenvectorsType eigenvectors() const; 179 180 /** \brief Returns an expression of the computed generalized eigenvalues. 181 * 182 * \returns An expression of the column vector containing the eigenvalues. 183 * 184 * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode 185 * Not that betas might contain zeros. It is therefore not recommended to use this function, 186 * but rather directly deal with the alphas and betas vectors. 187 * 188 * \pre Either the constructor 189 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function 190 * compute(const MatrixType&,const MatrixType&,bool) has been called before. 191 * 192 * The eigenvalues are repeated according to their algebraic multiplicity, 193 * so there are as many eigenvalues as rows in the matrix. The eigenvalues 194 * are not sorted in any particular order. 195 * 196 * \sa alphas(), betas(), eigenvectors() 197 */ 198 EigenvalueType eigenvalues() const 199 { 200 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 201 return EigenvalueType(m_alphas,m_betas); 202 } 203 204 /** \returns A const reference to the vectors containing the alpha values 205 * 206 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 207 * 208 * \sa betas(), eigenvalues() */ 209 ComplexVectorType alphas() const 210 { 211 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 212 return m_alphas; 213 } 214 215 /** \returns A const reference to the vectors containing the beta values 216 * 217 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 218 * 219 * \sa alphas(), eigenvalues() */ 220 VectorType betas() const 221 { 222 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 223 return m_betas; 224 } 225 226 /** \brief Computes generalized eigendecomposition of given matrix. 227 * 228 * \param[in] A Square matrix whose eigendecomposition is to be computed. 229 * \param[in] B Square matrix whose eigendecomposition is to be computed. 230 * \param[in] computeEigenvectors If true, both the eigenvectors and the 231 * eigenvalues are computed; if false, only the eigenvalues are 232 * computed. 233 * \returns Reference to \c *this 234 * 235 * This function computes the eigenvalues of the real matrix \p matrix. 236 * The eigenvalues() function can be used to retrieve them. If 237 * \p computeEigenvectors is true, then the eigenvectors are also computed 238 * and can be retrieved by calling eigenvectors(). 239 * 240 * The matrix is first reduced to real generalized Schur form using the RealQZ 241 * class. The generalized Schur decomposition is then used to compute the eigenvalues 242 * and eigenvectors. 243 * 244 * The cost of the computation is dominated by the cost of the 245 * generalized Schur decomposition. 246 * 247 * This method reuses of the allocated data in the GeneralizedEigenSolver object. 248 */ 249 GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); 250 251 ComputationInfo info() const 252 { 253 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 254 return m_realQZ.info(); 255 } 256 257 /** Sets the maximal number of iterations allowed. 258 */ 259 GeneralizedEigenSolver& setMaxIterations(Index maxIters) 260 { 261 m_realQZ.setMaxIterations(maxIters); 262 return *this; 263 } 264 265 protected: 266 MatrixType m_eivec; 267 ComplexVectorType m_alphas; 268 VectorType m_betas; 269 bool m_isInitialized; 270 bool m_eigenvectorsOk; 271 RealQZ<MatrixType> m_realQZ; 272 MatrixType m_matS; 273 274 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 275 ColumnVectorType m_tmp; 276 }; 277 278 //template<typename MatrixType> 279 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const 280 //{ 281 // eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 282 // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 283 // Index n = m_eivec.cols(); 284 // EigenvectorsType matV(n,n); 285 // // TODO 286 // return matV; 287 //} 288 289 template<typename MatrixType> 290 GeneralizedEigenSolver<MatrixType>& 291 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) 292 { 293 using std::sqrt; 294 using std::abs; 295 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); 296 297 // Reduce to generalized real Schur form: 298 // A = Q S Z and B = Q T Z 299 m_realQZ.compute(A, B, computeEigenvectors); 300 301 if (m_realQZ.info() == Success) 302 { 303 m_matS = m_realQZ.matrixS(); 304 if (computeEigenvectors) 305 m_eivec = m_realQZ.matrixZ().transpose(); 306 307 // Compute eigenvalues from matS 308 m_alphas.resize(A.cols()); 309 m_betas.resize(A.cols()); 310 Index i = 0; 311 while (i < A.cols()) 312 { 313 if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) 314 { 315 m_alphas.coeffRef(i) = m_matS.coeff(i, i); 316 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); 317 ++i; 318 } 319 else 320 { 321 Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); 322 Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); 323 m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); 324 m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); 325 326 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); 327 m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); 328 i += 2; 329 } 330 } 331 } 332 333 m_isInitialized = true; 334 m_eigenvectorsOk = false;//computeEigenvectors; 335 336 return *this; 337 } 338 339 } // end namespace Eigen 340 341 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H 342
