From 8e5d2b6fc49f58a62d8c90f61157296de2861a04 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Sun, 21 Mar 2010 21:57:34 +0000 Subject: [PATCH] Rename Complex in ComplexSchur and ComplexEigenSolver to ComplexScalar for consistency with the RealScalar type; correct ComplexEigenSolver docs to take non-diagonalizable matrices into account; refactor ComplexEigenSolver::compute(). --- Eigen/src/Eigenvalues/ComplexEigenSolver.h | 110 +++++++++++---------- Eigen/src/Eigenvalues/ComplexSchur.h | 20 ++-- 2 files changed, 68 insertions(+), 62 deletions(-) diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 7ca87b8cc..c1e65cbfd 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -38,11 +38,12 @@ * instantiation of the Matrix class template. * * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars - * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$. - * The eigendecomposition of a matrix is \f$ A = V D V^{-1} \f$, - * where \f$ D \f$ is a diagonal matrix. The entries on the diagonal - * of \f$ D \f$ are the eigenvalues and the columns of \f$ V \f$ are - * the eigenvectors. + * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v + * \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 = V D \f$. The matrix \f$ V \f$ is + * almost always invertible, in which case we have \f$ A = V D V^{-1} + * \f$. This is called the eigendecomposition. * * The main function in this class is compute(), which computes the * eigenvalues and eigenvectors of a given function. The @@ -73,21 +74,21 @@ template class ComplexEigenSolver * \c float or \c double) and just \c Scalar if #Scalar is * complex. */ - typedef std::complex Complex; + typedef std::complex ComplexScalar; /** \brief Type for vector of eigenvalues as returned by eigenvalues(). * - * This is a column vector with entries of type #Complex. + * This is a column vector with entries of type #ComplexScalar. * The length of the vector is the size of \p _MatrixType. */ - typedef Matrix EigenvalueType; + typedef Matrix EigenvalueType; /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). * - * This is a square matrix with entries of type #Complex. + * This is a square matrix with entries of type #ComplexScalar. * The size is the same as the size of \p _MatrixType. */ - typedef Matrix EigenvectorType; + typedef Matrix EigenvectorType; /** \brief Default constructor. * @@ -99,7 +100,7 @@ template class ComplexEigenSolver /** \brief Constructor; computes eigendecomposition of given matrix. * - * \param[in] matrix %Matrix whose eigendecomposition is to be computed. + * \param[in] matrix Sqarae matrix whose eigendecomposition is to be computed. * * This constructor calls compute() to compute the eigendecomposition. */ @@ -117,11 +118,13 @@ template class ComplexEigenSolver * ComplexEigenSolver(const MatrixType& matrix) or the member * function compute(const MatrixType& matrix) has been called * before to compute the eigendecomposition of a matrix. This - * function returns the matrix \f$ V \f$ in the - * eigendecomposition \f$ A = V D V^{-1} \f$. The columns of \f$ - * V \f$ are the eigenvectors. The eigenvectors are normalized to - * have (Euclidean) norm equal to one, and are in the same order - * as the eigenvalues as returned by eigenvalues(). + * function returns a matrix whose columns are the + * eigenvectors. Column \f$ k \f$ 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 eigendecomposition \f$ + * A = V D V^{-1} \f$, if it exists. * * Example: \include ComplexEigenSolver_eigenvectors.cpp * Output: \verbinclude ComplexEigenSolver_eigenvectors.out @@ -138,7 +141,10 @@ template class ComplexEigenSolver * ComplexEigenSolver(const MatrixType& matrix) or the member * function compute(const MatrixType& matrix) has been called * before to compute the eigendecomposition of a matrix. This - * function returns a column vector containing the eigenvalues. + * function returns a column vector containing the + * eigenvalues. Eigenvalues are repeated according to their + * algebraic multiplicity, so there are as many eigenvalues as + * rows in the matrix. * * Example: \include ComplexEigenSolver_eigenvalues.cpp * Output: \verbinclude ComplexEigenSolver_eigenvalues.out @@ -151,7 +157,7 @@ template class ComplexEigenSolver /** \brief Computes eigendecomposition of given matrix. * - * \param[in] matrix %Matrix whose eigendecomposition is to be computed. + * \param[in] matrix Square matrix whose eigendecomposition is to be computed. * * This function computes the eigenvalues and eigenvectors of \p * matrix. The eigenvalues() and eigenvectors() functions can be @@ -182,56 +188,56 @@ void ComplexEigenSolver::compute(const MatrixType& matrix) { // this code is inspired from Jampack assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - m_eivec.resize(n,n); + const int n = matrix.cols(); + const RealScalar matrixnorm = matrix.norm(); - RealScalar eps = NumTraits::epsilon(); - - // Reduce to complex Schur form + // Step 1: Do a complex Schur decomposition, A = U T U^* + // The eigenvalues are on the diagonal of T. ComplexSchur schur(matrix); - m_eivalues = schur.matrixT().diagonal(); - m_eivec.setZero(); - - Complex d2, z; - RealScalar norm = matrix.norm(); - - // compute the (normalized) eigenvectors + // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T. + // The matrix X is unit triangular. + EigenvectorType X = EigenvectorType::Zero(n, n); for(int k=n-1 ; k>=0 ; k--) { - d2 = schur.matrixT().coeff(k,k); - m_eivec.coeffRef(k,k) = Complex(1.0,0.0); + X.coeffRef(k,k) = ComplexScalar(1.0,0.0); + // Compute X(i,k) using the (i,k) entry of the equation X T = D X for(int i=k-1 ; i>=0 ; i--) { - m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); + X.coeffRef(i,k) = -schur.matrixT().coeff(i,k); if(k-i-1>0) - m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); - z = schur.matrixT().coeff(i,i) - d2; - if(z==Complex(0)) - ei_real_ref(z) = eps * norm; - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; - + X.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value(); + ComplexScalar z = schur.matrixT().coeff(i,i) - schur.matrixT().coeff(k,k); + if(z==ComplexScalar(0)) + { + // If the i-th and k-th eigenvalue are equal, then z equals 0. + // Use a small value instead, to prevent division by zero. + ei_real_ref(z) = NumTraits::epsilon() * matrixnorm; + } + X.coeffRef(i,k) = X.coeff(i,k) / z; } - m_eivec.col(k).normalize(); } - m_eivec = schur.matrixU() * m_eivec; + // Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) + m_eivec = schur.matrixU() * X; + // .. and normalize the eigenvectors + for(int k=0 ; k class ComplexSchur }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix ComplexMatrixType; + typedef std::complex ComplexScalar; + typedef Matrix ComplexMatrixType; enum { Size = MatrixType::RowsAtCompileTime }; @@ -158,8 +158,8 @@ void ComplexSchur::compute(const MatrixType& matrix, bool skipU) // TODO skip Q if skipU = true HessenbergDecomposition hess(matrix); - m_matT = hess.matrixH().template cast(); - if(!skipU) m_matU = hess.matrixQ().template cast(); + m_matT = hess.matrixH().template cast(); + if(!skipU) m_matU = hess.matrixQ().template cast(); // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. @@ -172,7 +172,7 @@ void ComplexSchur::compute(const MatrixType& matrix, bool skipU) int iu = m_matT.cols() - 1; int il; RealScalar d,sd,sf; - Complex c,b,disc,r1,r2,kappa; + ComplexScalar c,b,disc,r1,r2,kappa; RealScalar eps = NumTraits::epsilon(); @@ -188,7 +188,7 @@ void ComplexSchur::compute(const MatrixType& matrix, bool skipU) if(!ei_isMuchSmallerThan(sd,d,eps)) break; - m_matT.coeffRef(iu,iu-1) = Complex(0); + m_matT.coeffRef(iu,iu-1) = ComplexScalar(0); iter = 0; --iu; } @@ -216,12 +216,12 @@ void ComplexSchur::compute(const MatrixType& matrix, bool skipU) --il; } - if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); + if( il != 0 ) m_matT.coeffRef(il,il-1) = ComplexScalar(0); // compute the shift kappa as one of the eigenvalues of the 2x2 // diagonal block on the bottom of the active submatrix - Matrix t = m_matT.template block<2,2>(iu-1,iu-1); + Matrix t = m_matT.template block<2,2>(iu-1,iu-1); sf = t.cwiseAbs().sum(); t /= sf; // the normalization by sf is to avoid under/overflow @@ -251,7 +251,7 @@ void ComplexSchur::compute(const MatrixType& matrix, bool skipU) } // perform the QR step using Givens rotations - PlanarRotation rot; + PlanarRotation rot; rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); for(int i=il ; i::compute(const MatrixType& matrix, bool skipU) int i2 = i+2; rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i)); - m_matT.coeffRef(i2,i) = Complex(0); + m_matT.coeffRef(i2,i) = ComplexScalar(0); } } }