Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
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
Generated on Tue Jul 12 2022 17:46:54 by 1.7.2