1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2010,2012 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 6 // Copyright (C) 2016 Tobias Wood <tobias (at) spinicist.org.uk> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H 13 #define EIGEN_GENERALIZEDEIGENSOLVER_H 14 15 #include "./RealQZ.h" 16 17 namespace Eigen { 18 19 /** \eigenvalues_module \ingroup Eigenvalues_Module 20 * 21 * 22 * \class GeneralizedEigenSolver 23 * 24 * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices 25 * 26 * \tparam _MatrixType the type of the matrices of which we are computing the 27 * eigen-decomposition; this is expected to be an instantiation of the Matrix 28 * class template. Currently, only real matrices are supported. 29 * 30 * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars 31 * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If 32 * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and 33 * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = 34 * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we 35 * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. 36 * 37 * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the 38 * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is 39 * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ 40 * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, 41 * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: 42 * \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 43 * called the left eigenvector. 44 * 45 * Call the function compute() to compute the generalized eigenvalues and eigenvectors of 46 * a given matrix pair. Alternatively, you can use the 47 * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the 48 * eigenvalues and eigenvectors at construction time. Once the eigenvalue and 49 * eigenvectors are computed, they can be retrieved with the eigenvalues() and 50 * eigenvectors() functions. 51 * 52 * Here is an usage example of this class: 53 * Example: \include GeneralizedEigenSolver.cpp 54 * Output: \verbinclude GeneralizedEigenSolver.out 55 * 56 * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver 57 */ 58 template<typename _MatrixType> class GeneralizedEigenSolver 59 { 60 public: 61 62 /** \brief Synonym for the template parameter \p _MatrixType. */ 63 typedef _MatrixType MatrixType; 64 65 enum { 66 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 67 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 68 Options = MatrixType::Options, 69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 71 }; 72 73 /** \brief Scalar type for matrices of type #MatrixType. */ 74 typedef typename MatrixType::Scalar Scalar; 75 typedef typename NumTraits<Scalar>::Real RealScalar; 76 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 77 78 /** \brief Complex scalar type for #MatrixType. 79 * 80 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 81 * \c float or \c double) and just \c Scalar if #Scalar is 82 * complex. 83 */ 84 typedef std::complex<RealScalar> ComplexScalar; 85 86 /** \brief Type for vector of real scalar values eigenvalues as returned by betas(). 87 * 88 * This is a column vector with entries of type #Scalar. 89 * The length of the vector is the size of #MatrixType. 90 */ 91 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; 92 93 /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). 94 * 95 * This is a column vector with entries of type #ComplexScalar. 96 * The length of the vector is the size of #MatrixType. 97 */ 98 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; 99 100 /** \brief Expression type for the eigenvalues as returned by eigenvalues(). 101 */ 102 typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; 103 104 /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 105 * 106 * This is a square matrix with entries of type #ComplexScalar. 107 * The size is the same as the size of #MatrixType. 108 */ 109 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 110 111 /** \brief Default constructor. 112 * 113 * The default constructor is useful in cases in which the user intends to 114 * perform decompositions via EigenSolver::compute(const MatrixType&, bool). 115 * 116 * \sa compute() for an example. 117 */ 118 GeneralizedEigenSolver() 119 : m_eivec(), 120 m_alphas(), 121 m_betas(), 122 m_valuesOkay(false), 123 m_vectorsOkay(false), 124 m_realQZ() 125 {} 126 127 /** \brief Default constructor with memory preallocation 128 * 129 * Like the default constructor but with preallocation of the internal data 130 * according to the specified problem \a size. 131 * \sa GeneralizedEigenSolver() 132 */ 133 explicit GeneralizedEigenSolver(Index size) 134 : m_eivec(size, size), 135 m_alphas(size), 136 m_betas(size), 137 m_valuesOkay(false), 138 m_vectorsOkay(false), 139 m_realQZ(size), 140 m_tmp(size) 141 {} 142 143 /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. 144 * 145 * \param[in] A Square matrix whose eigendecomposition is to be computed. 146 * \param[in] B Square matrix whose eigendecomposition is to be computed. 147 * \param[in] computeEigenvectors If true, both the eigenvectors and the 148 * eigenvalues are computed; if false, only the eigenvalues are computed. 149 * 150 * This constructor calls compute() to compute the generalized eigenvalues 151 * and eigenvectors. 152 * 153 * \sa compute() 154 */ 155 GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) 156 : m_eivec(A.rows(), A.cols()), 157 m_alphas(A.cols()), 158 m_betas(A.cols()), 159 m_valuesOkay(false), 160 m_vectorsOkay(false), 161 m_realQZ(A.cols()), 162 m_tmp(A.cols()) 163 { 164 compute(A, B, computeEigenvectors); 165 } 166 167 /* \brief Returns the computed generalized eigenvectors. 168 * 169 * \returns %Matrix whose columns are the (possibly complex) right eigenvectors. 170 * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. 171 * 172 * \pre Either the constructor 173 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function 174 * compute(const MatrixType&, const MatrixType& bool) has been called before, and 175 * \p computeEigenvectors was set to true (the default). 176 * 177 * \sa eigenvalues() 178 */ 179 EigenvectorsType eigenvectors() const { 180 eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated."); 181 return m_eivec; 182 } 183 184 /** \brief Returns an expression of the computed generalized eigenvalues. 185 * 186 * \returns An expression of the column vector containing the eigenvalues. 187 * 188 * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode 189 * Not that betas might contain zeros. It is therefore not recommended to use this function, 190 * but rather directly deal with the alphas and betas vectors. 191 * 192 * \pre Either the constructor 193 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function 194 * compute(const MatrixType&,const MatrixType&,bool) has been called before. 195 * 196 * The eigenvalues are repeated according to their algebraic multiplicity, 197 * so there are as many eigenvalues as rows in the matrix. The eigenvalues 198 * are not sorted in any particular order. 199 * 200 * \sa alphas(), betas(), eigenvectors() 201 */ 202 EigenvalueType eigenvalues() const 203 { 204 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 205 return EigenvalueType(m_alphas,m_betas); 206 } 207 208 /** \returns A const reference to the vectors containing the alpha values 209 * 210 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 211 * 212 * \sa betas(), eigenvalues() */ 213 ComplexVectorType alphas() const 214 { 215 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 216 return m_alphas; 217 } 218 219 /** \returns A const reference to the vectors containing the beta values 220 * 221 * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 222 * 223 * \sa alphas(), eigenvalues() */ 224 VectorType betas() const 225 { 226 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 227 return m_betas; 228 } 229 230 /** \brief Computes generalized eigendecomposition of given matrix. 231 * 232 * \param[in] A Square matrix whose eigendecomposition is to be computed. 233 * \param[in] B Square matrix whose eigendecomposition is to be computed. 234 * \param[in] computeEigenvectors If true, both the eigenvectors and the 235 * eigenvalues are computed; if false, only the eigenvalues are 236 * computed. 237 * \returns Reference to \c *this 238 * 239 * This function computes the eigenvalues of the real matrix \p matrix. 240 * The eigenvalues() function can be used to retrieve them. If 241 * \p computeEigenvectors is true, then the eigenvectors are also computed 242 * and can be retrieved by calling eigenvectors(). 243 * 244 * The matrix is first reduced to real generalized Schur form using the RealQZ 245 * class. The generalized Schur decomposition is then used to compute the eigenvalues 246 * and eigenvectors. 247 * 248 * The cost of the computation is dominated by the cost of the 249 * generalized Schur decomposition. 250 * 251 * This method reuses of the allocated data in the GeneralizedEigenSolver object. 252 */ 253 GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); 254 255 ComputationInfo info() const 256 { 257 eigen_assert(m_valuesOkay && "EigenSolver is not initialized."); 258 return m_realQZ.info(); 259 } 260 261 /** Sets the maximal number of iterations allowed. 262 */ 263 GeneralizedEigenSolver& setMaxIterations(Index maxIters) 264 { 265 m_realQZ.setMaxIterations(maxIters); 266 return *this; 267 } 268 269 protected: 270 271 static void check_template_parameters() 272 { 273 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 274 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 275 } 276 277 EigenvectorsType m_eivec; 278 ComplexVectorType m_alphas; 279 VectorType m_betas; 280 bool m_valuesOkay, m_vectorsOkay; 281 RealQZ<MatrixType> m_realQZ; 282 ComplexVectorType m_tmp; 283 }; 284 285 template<typename MatrixType> 286 GeneralizedEigenSolver<MatrixType>& 287 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) 288 { 289 check_template_parameters(); 290 291 using std::sqrt; 292 using std::abs; 293 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); 294 Index size = A.cols(); 295 m_valuesOkay = false; 296 m_vectorsOkay = false; 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 if (m_realQZ.info() == Success) 301 { 302 // Resize storage 303 m_alphas.resize(size); 304 m_betas.resize(size); 305 if (computeEigenvectors) 306 { 307 m_eivec.resize(size,size); 308 m_tmp.resize(size); 309 } 310 311 // Aliases: 312 Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size); 313 ComplexVectorType &cv = m_tmp; 314 const MatrixType &mZ = m_realQZ.matrixZ(); 315 const MatrixType &mS = m_realQZ.matrixS(); 316 const MatrixType &mT = m_realQZ.matrixT(); 317 318 Index i = 0; 319 while (i < size) 320 { 321 if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0)) 322 { 323 // Real eigenvalue 324 m_alphas.coeffRef(i) = mS.diagonal().coeff(i); 325 m_betas.coeffRef(i) = mT.diagonal().coeff(i); 326 if (computeEigenvectors) 327 { 328 v.setConstant(Scalar(0.0)); 329 v.coeffRef(i) = Scalar(1.0); 330 // For singular eigenvalues do nothing more 331 if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)()) 332 { 333 // Non-singular eigenvalue 334 const Scalar alpha = real(m_alphas.coeffRef(i)); 335 const Scalar beta = m_betas.coeffRef(i); 336 for (Index j = i-1; j >= 0; j--) 337 { 338 const Index st = j+1; 339 const Index sz = i-j; 340 if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) 341 { 342 // 2x2 block 343 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) ); 344 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); 345 v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); 346 j--; 347 } 348 else 349 { 350 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)); 351 } 352 } 353 } 354 m_eivec.col(i).real().noalias() = mZ.transpose() * v; 355 m_eivec.col(i).real().normalize(); 356 m_eivec.col(i).imag().setConstant(0); 357 } 358 ++i; 359 } 360 else 361 { 362 // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T 363 // 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): 364 365 // T = [a 0] 366 // [0 b] 367 RealScalar a = mT.diagonal().coeff(i), 368 b = mT.diagonal().coeff(i+1); 369 const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b; 370 371 // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. 372 Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal(); 373 374 Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); 375 Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); 376 const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z); 377 m_alphas.coeffRef(i) = conj(alpha); 378 m_alphas.coeffRef(i+1) = alpha; 379 380 if (computeEigenvectors) { 381 // Compute eigenvector in position (i+1) and then position (i) is just the conjugate 382 cv.setZero(); 383 cv.coeffRef(i+1) = Scalar(1.0); 384 // here, the "static_cast" workaound expression template issues. 385 cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1)) 386 / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i)); 387 for (Index j = i-1; j >= 0; j--) 388 { 389 const Index st = j+1; 390 const Index sz = i+1-j; 391 if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) 392 { 393 // 2x2 block 394 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) ); 395 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); 396 cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); 397 j--; 398 } else { 399 cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() 400 / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j))); 401 } 402 } 403 m_eivec.col(i+1).noalias() = (mZ.transpose() * cv); 404 m_eivec.col(i+1).normalize(); 405 m_eivec.col(i) = m_eivec.col(i+1).conjugate(); 406 } 407 i += 2; 408 } 409 } 410 411 m_valuesOkay = true; 412 m_vectorsOkay = computeEigenvectors; 413 } 414 return *this; 415 } 416 417 } // end namespace Eigen 418 419 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H 420