EigenSolver.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_EIGENSOLVER_H
00027 #define EIGEN_EIGENSOLVER_H
00028 
00029 #include "./RealSchur.h"
00030 
00031 namespace Eigen { 
00032 
00079 template<typename _MatrixType> class EigenSolver
00080 {
00081   public:
00082 
00084     typedef _MatrixType MatrixType;
00085 
00086     enum {
00087       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00088       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00089       Options = MatrixType::Options,
00090       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00091       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00092     };
00093 
00095     typedef typename MatrixType::Scalar Scalar;
00096     typedef typename NumTraits<Scalar>::Real RealScalar;
00097     typedef typename MatrixType::Index Index;
00098 
00105     typedef std::complex<RealScalar> ComplexScalar;
00106 
00112     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
00113 
00119     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
00120 
00128  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
00129 
00136     EigenSolver(Index size)
00137       : m_eivec(size, size),
00138         m_eivalues(size),
00139         m_isInitialized(false),
00140         m_eigenvectorsOk(false),
00141         m_realSchur(size),
00142         m_matT(size, size), 
00143         m_tmp(size)
00144     {}
00145 
00161     EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00162       : m_eivec(matrix.rows(), matrix.cols()),
00163         m_eivalues(matrix.cols()),
00164         m_isInitialized(false),
00165         m_eigenvectorsOk(false),
00166         m_realSchur(matrix.cols()),
00167         m_matT(matrix.rows(), matrix.cols()), 
00168         m_tmp(matrix.cols())
00169     {
00170       compute(matrix, computeEigenvectors);
00171     }
00172 
00193     EigenvectorsType eigenvectors() const;
00194 
00213     const MatrixType& pseudoEigenvectors() const
00214     {
00215       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00216       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00217       return m_eivec;
00218     }
00219 
00238     MatrixType pseudoEigenvalueMatrix() const;
00239 
00258     const EigenvalueType& eigenvalues() const
00259     {
00260       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00261       return m_eivalues;
00262     }
00263 
00291     EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00292 
00293     ComputationInfo info() const
00294     {
00295       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00296       return m_realSchur.info();
00297     }
00298 
00299   private:
00300     void doComputeEigenvectors();
00301 
00302   protected:
00303     MatrixType m_eivec;
00304     EigenvalueType m_eivalues;
00305     bool m_isInitialized;
00306     bool m_eigenvectorsOk;
00307     RealSchur<MatrixType> m_realSchur;
00308     MatrixType m_matT;
00309 
00310     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00311     ColumnVectorType m_tmp;
00312 };
00313 
00314 template<typename MatrixType>
00315 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
00316 {
00317   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00318   Index n = m_eivalues.rows();
00319   MatrixType matD = MatrixType::Zero(n,n);
00320   for (Index i=0; i<n; ++i)
00321   {
00322     if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
00323       matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
00324     else
00325     {
00326       matD.template block<2,2>(i,i) <<  internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
00327                                        -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
00328       ++i;
00329     }
00330   }
00331   return matD;
00332 }
00333 
00334 template<typename MatrixType>
00335 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
00336 {
00337   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00338   eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00339   Index n = m_eivec.cols();
00340   EigenvectorsType matV(n,n);
00341   for (Index j=0; j<n; ++j)
00342   {
00343     if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
00344     {
00345       // we have a real eigen value
00346       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
00347       matV.col(j).normalize();
00348     }
00349     else
00350     {
00351       // we have a pair of complex eigen values
00352       for (Index i=0; i<n; ++i)
00353       {
00354         matV.coeffRef(i,j)   = ComplexScalar(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
00355         matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
00356       }
00357       matV.col(j).normalize();
00358       matV.col(j+1).normalize();
00359       ++j;
00360     }
00361   }
00362   return matV;
00363 }
00364 
00365 template<typename MatrixType>
00366 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00367 {
00368   assert(matrix.cols() == matrix.rows());
00369 
00370   // Reduce to real Schur form.
00371   m_realSchur.compute(matrix, computeEigenvectors);
00372   if (m_realSchur.info() == Success)
00373   {
00374     m_matT = m_realSchur.matrixT();
00375     if (computeEigenvectors)
00376       m_eivec = m_realSchur.matrixU();
00377   
00378     // Compute eigenvalues from matT
00379     m_eivalues.resize(matrix.cols());
00380     Index i = 0;
00381     while (i < matrix.cols()) 
00382     {
00383       if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
00384       {
00385         m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
00386         ++i;
00387       }
00388       else
00389       {
00390         Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
00391         Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
00392         m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
00393         m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
00394         i += 2;
00395       }
00396     }
00397     
00398     // Compute eigenvectors.
00399     if (computeEigenvectors)
00400       doComputeEigenvectors();
00401   }
00402 
00403   m_isInitialized = true;
00404   m_eigenvectorsOk = computeEigenvectors;
00405 
00406   return *this;
00407 }
00408 
00409 // Complex scalar division.
00410 template<typename Scalar>
00411 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
00412 {
00413   Scalar r,d;
00414   if (internal::abs(yr) > internal::abs(yi))
00415   {
00416       r = yi/yr;
00417       d = yr + r*yi;
00418       return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
00419   }
00420   else
00421   {
00422       r = yr/yi;
00423       d = yi + r*yr;
00424       return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
00425   }
00426 }
00427 
00428 
00429 template<typename MatrixType>
00430 void EigenSolver<MatrixType>::doComputeEigenvectors()
00431 {
00432   const Index size = m_eivec.cols();
00433   const Scalar eps = NumTraits<Scalar>::epsilon();
00434 
00435   // inefficient! this is already computed in RealSchur
00436   Scalar norm(0);
00437   for (Index j = 0; j < size; ++j)
00438   {
00439     norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
00440   }
00441   
00442   // Backsubstitute to find vectors of upper triangular form
00443   if (norm == 0.0)
00444   {
00445     return;
00446   }
00447 
00448   for (Index n = size-1; n >= 0; n--)
00449   {
00450     Scalar p = m_eivalues.coeff(n).real();
00451     Scalar q = m_eivalues.coeff(n).imag();
00452 
00453     // Scalar vector
00454     if (q == Scalar(0))
00455     {
00456       Scalar lastr(0), lastw(0);
00457       Index l = n;
00458 
00459       m_matT.coeffRef(n,n) = 1.0;
00460       for (Index i = n-1; i >= 0; i--)
00461       {
00462         Scalar w = m_matT.coeff(i,i) - p;
00463         Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00464 
00465         if (m_eivalues.coeff(i).imag() < 0.0)
00466         {
00467           lastw = w;
00468           lastr = r;
00469         }
00470         else
00471         {
00472           l = i;
00473           if (m_eivalues.coeff(i).imag() == 0.0)
00474           {
00475             if (w != 0.0)
00476               m_matT.coeffRef(i,n) = -r / w;
00477             else
00478               m_matT.coeffRef(i,n) = -r / (eps * norm);
00479           }
00480           else // Solve real equations
00481           {
00482             Scalar x = m_matT.coeff(i,i+1);
00483             Scalar y = m_matT.coeff(i+1,i);
00484             Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
00485             Scalar t = (x * lastr - lastw * r) / denom;
00486             m_matT.coeffRef(i,n) = t;
00487             if (internal::abs(x) > internal::abs(lastw))
00488               m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
00489             else
00490               m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
00491           }
00492 
00493           // Overflow control
00494           Scalar t = internal::abs(m_matT.coeff(i,n));
00495           if ((eps * t) * t > Scalar(1))
00496             m_matT.col(n).tail(size-i) /= t;
00497         }
00498       }
00499     }
00500     else if (q < Scalar(0) && n > 0) // Complex vector
00501     {
00502       Scalar lastra(0), lastsa(0), lastw(0);
00503       Index l = n-1;
00504 
00505       // Last vector component imaginary so matrix is triangular
00506       if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
00507       {
00508         m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
00509         m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
00510       }
00511       else
00512       {
00513         std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
00514         m_matT.coeffRef(n-1,n-1) = internal::real(cc);
00515         m_matT.coeffRef(n-1,n) = internal::imag(cc);
00516       }
00517       m_matT.coeffRef(n,n-1) = 0.0;
00518       m_matT.coeffRef(n,n) = 1.0;
00519       for (Index i = n-2; i >= 0; i--)
00520       {
00521         Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
00522         Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00523         Scalar w = m_matT.coeff(i,i) - p;
00524 
00525         if (m_eivalues.coeff(i).imag() < 0.0)
00526         {
00527           lastw = w;
00528           lastra = ra;
00529           lastsa = sa;
00530         }
00531         else
00532         {
00533           l = i;
00534           if (m_eivalues.coeff(i).imag() == RealScalar(0))
00535           {
00536             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
00537             m_matT.coeffRef(i,n-1) = internal::real(cc);
00538             m_matT.coeffRef(i,n) = internal::imag(cc);
00539           }
00540           else
00541           {
00542             // Solve complex equations
00543             Scalar x = m_matT.coeff(i,i+1);
00544             Scalar y = m_matT.coeff(i+1,i);
00545             Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
00546             Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
00547             if ((vr == 0.0) && (vi == 0.0))
00548               vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
00549 
00550             std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
00551             m_matT.coeffRef(i,n-1) = internal::real(cc);
00552             m_matT.coeffRef(i,n) = internal::imag(cc);
00553             if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
00554             {
00555               m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
00556               m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
00557             }
00558             else
00559             {
00560               cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
00561               m_matT.coeffRef(i+1,n-1) = internal::real(cc);
00562               m_matT.coeffRef(i+1,n) = internal::imag(cc);
00563             }
00564           }
00565 
00566           // Overflow control
00567           using std::max;
00568           Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
00569           if ((eps * t) * t > Scalar(1))
00570             m_matT.block(i, n-1, size-i, 2) /= t;
00571 
00572         }
00573       }
00574       
00575       // We handled a pair of complex conjugate eigenvalues, so need to skip them both
00576       n--;
00577     }
00578     else
00579     {
00580       eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
00581     }
00582   }
00583 
00584   // Back transformation to get eigenvectors of original matrix
00585   for (Index j = size-1; j >= 0; j--)
00586   {
00587     m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
00588     m_eivec.col(j) = m_tmp;
00589   }
00590 }
00591 
00592 } // end namespace Eigen
00593 
00594 #endif // EIGEN_EIGENSOLVER_H