Main MRPT website > C++ reference
MRPT logo

ComplexEigenSolver.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) 2009 Claire Maurice
00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00007 //
00008 // Eigen is free software; you can redistribute it and/or
00009 // modify it under the terms of the GNU Lesser General Public
00010 // License as published by the Free Software Foundation; either
00011 // version 3 of the License, or (at your option) any later version.
00012 //
00013 // Alternatively, you can redistribute it and/or
00014 // modify it under the terms of the GNU General Public License as
00015 // published by the Free Software Foundation; either version 2 of
00016 // the License, or (at your option) any later version.
00017 //
00018 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00019 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00020 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00021 // GNU General Public License for more details.
00022 //
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License and a copy of the GNU General Public License along with
00025 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
00028 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
00029 
00030 #include "./EigenvaluesCommon.h"
00031 #include "./ComplexSchur.h"
00032 
00033 /** \eigenvalues_module \ingroup Eigenvalues_Module
00034   *
00035   *
00036   * \class ComplexEigenSolver
00037   *
00038   * \brief Computes eigenvalues and eigenvectors of general complex matrices
00039   *
00040   * \tparam _MatrixType the type of the matrix of which we are
00041   * computing the eigendecomposition; this is expected to be an
00042   * instantiation of the Matrix class template.
00043   *
00044   * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
00045   * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
00046   * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
00047   * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
00048   * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
00049   * almost always invertible, in which case we have \f$ A = V D V^{-1}
00050   * \f$. This is called the eigendecomposition.
00051   *
00052   * The main function in this class is compute(), which computes the
00053   * eigenvalues and eigenvectors of a given function. The
00054   * documentation for that function contains an example showing the
00055   * main features of the class.
00056   *
00057   * \sa class EigenSolver, class SelfAdjointEigenSolver
00058   */
00059 template<typename _MatrixType> class ComplexEigenSolver
00060 {
00061   public:
00062 
00063     /** \brief Synonym for the template parameter \p _MatrixType. */
00064     typedef _MatrixType MatrixType;
00065 
00066     enum {
00067       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00068       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00069       Options = MatrixType::Options,
00070       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00071       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00072     };
00073 
00074     /** \brief Scalar type for matrices of type #MatrixType. */
00075     typedef typename MatrixType::Scalar Scalar;
00076     typedef typename NumTraits<Scalar>::Real RealScalar;
00077     typedef typename MatrixType::Index Index;
00078 
00079     /** \brief Complex scalar type for #MatrixType.
00080       *
00081       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
00082       * \c float or \c double) and just \c Scalar if #Scalar is
00083       * complex.
00084       */
00085     typedef std::complex<RealScalar> ComplexScalar;
00086 
00087     /** \brief Type for vector of eigenvalues as returned by eigenvalues().
00088       *
00089       * This is a column vector with entries of type #ComplexScalar.
00090       * The length of the vector is the size of #MatrixType.
00091       */
00092     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
00093 
00094     /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
00095       *
00096       * This is a square matrix with entries of type #ComplexScalar.
00097       * The size is the same as the size of #MatrixType.
00098       */
00099     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, ColsAtCompileTime> EigenvectorType;
00100 
00101     /** \brief Default constructor.
00102       *
00103       * The default constructor is useful in cases in which the user intends to
00104       * perform decompositions via compute().
00105       */
00106     ComplexEigenSolver()
00107             : m_eivec(),
00108               m_eivalues(),
00109               m_schur(),
00110               m_isInitialized(false),
00111               m_eigenvectorsOk(false),
00112               m_matX()
00113     {}
00114 
00115     /** \brief Default Constructor with memory preallocation
00116       *
00117       * Like the default constructor but with preallocation of the internal data
00118       * according to the specified problem \a size.
00119       * \sa ComplexEigenSolver()
00120       */
00121     ComplexEigenSolver(Index size)
00122             : m_eivec(size, size),
00123               m_eivalues(size),
00124               m_schur(size),
00125               m_isInitialized(false),
00126               m_eigenvectorsOk(false),
00127               m_matX(size, size)
00128     {}
00129 
00130     /** \brief Constructor; computes eigendecomposition of given matrix.
00131       *
00132       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00133       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00134       *    eigenvalues are computed; if false, only the eigenvalues are
00135       *    computed.
00136       *
00137       * This constructor calls compute() to compute the eigendecomposition.
00138       */
00139       ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00140             : m_eivec(matrix.rows(),matrix.cols()),
00141               m_eivalues(matrix.cols()),
00142               m_schur(matrix.rows()),
00143               m_isInitialized(false),
00144               m_eigenvectorsOk(false),
00145               m_matX(matrix.rows(),matrix.cols())
00146     {
00147       compute(matrix, computeEigenvectors);
00148     }
00149 
00150     /** \brief Returns the eigenvectors of given matrix.
00151       *
00152       * \returns  A const reference to the matrix whose columns are the eigenvectors.
00153       *
00154       * \pre Either the constructor
00155       * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
00156       * function compute(const MatrixType& matrix, bool) has been called before
00157       * to compute the eigendecomposition of a matrix, and
00158       * \p computeEigenvectors was set to true (the default).
00159       *
00160       * This function returns a matrix whose columns are the eigenvectors. Column
00161       * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
00162       * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
00163       * have (Euclidean) norm equal to one. The matrix returned by this
00164       * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
00165       * V^{-1} \f$, if it exists.
00166       *
00167       * Example: \include ComplexEigenSolver_eigenvectors.cpp
00168       * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
00169       */
00170     const EigenvectorType& eigenvectors() const
00171     {
00172       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00173       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00174       return m_eivec;
00175     }
00176 
00177     /** \brief Returns the eigenvalues of given matrix.
00178       *
00179       * \returns A const reference to the column vector containing the eigenvalues.
00180       *
00181       * \pre Either the constructor
00182       * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
00183       * function compute(const MatrixType& matrix, bool) has been called before
00184       * to compute the eigendecomposition of a matrix.
00185       *
00186       * This function returns a column vector containing the
00187       * eigenvalues. Eigenvalues are repeated according to their
00188       * algebraic multiplicity, so there are as many eigenvalues as
00189       * rows in the matrix.
00190       *
00191       * Example: \include ComplexEigenSolver_eigenvalues.cpp
00192       * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
00193       */
00194     const EigenvalueType& eigenvalues() const
00195     {
00196       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00197       return m_eivalues;
00198     }
00199 
00200     /** \brief Computes eigendecomposition of given matrix.
00201       *
00202       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00203       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00204       *    eigenvalues are computed; if false, only the eigenvalues are
00205       *    computed.
00206       * \returns    Reference to \c *this
00207       *
00208       * This function computes the eigenvalues of the complex matrix \p matrix.
00209       * The eigenvalues() function can be used to retrieve them.  If
00210       * \p computeEigenvectors is true, then the eigenvectors are also computed
00211       * and can be retrieved by calling eigenvectors().
00212       *
00213       * The matrix is first reduced to Schur form using the
00214       * ComplexSchur class. The Schur decomposition is then used to
00215       * compute the eigenvalues and eigenvectors.
00216       *
00217       * The cost of the computation is dominated by the cost of the
00218       * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
00219       * is the size of the matrix.
00220       *
00221       * Example: \include ComplexEigenSolver_compute.cpp
00222       * Output: \verbinclude ComplexEigenSolver_compute.out
00223       */
00224     ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00225 
00226     /** \brief Reports whether previous computation was successful.
00227       *
00228       * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
00229       */
00230     ComputationInfo info() const
00231     {
00232       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00233       return m_schur.info();
00234     }
00235 
00236   protected:
00237     EigenvectorType m_eivec;
00238     EigenvalueType m_eivalues;
00239     ComplexSchur<MatrixType> m_schur;
00240     bool m_isInitialized;
00241     bool m_eigenvectorsOk;
00242     EigenvectorType m_matX;
00243 
00244   private:
00245     void doComputeEigenvectors(RealScalar matrixnorm);
00246     void sortEigenvalues(bool computeEigenvectors);
00247 };
00248 
00249 
00250 template<typename MatrixType>
00251 ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00252 {
00253   // this code is inspired from Jampack
00254   assert(matrix.cols() == matrix.rows());
00255 
00256   // Do a complex Schur decomposition, A = U T U^*
00257   // The eigenvalues are on the diagonal of T.
00258   m_schur.compute(matrix, computeEigenvectors);
00259 
00260   if(m_schur.info() == Success)
00261   {
00262     m_eivalues = m_schur.matrixT().diagonal();
00263     if(computeEigenvectors)
00264       doComputeEigenvectors(matrix.norm());
00265     sortEigenvalues(computeEigenvectors);
00266   }
00267 
00268   m_isInitialized = true;
00269   m_eigenvectorsOk = computeEigenvectors;
00270   return *this;
00271 }
00272 
00273 
00274 template<typename MatrixType>
00275 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00276 {
00277   const Index n = m_eivalues.size();
00278 
00279   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00280   // The matrix X is unit triangular.
00281   m_matX = EigenvectorType::Zero(n, n);
00282   for(Index k=n-1 ; k>=0 ; k--)
00283   {
00284     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00285     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00286     for(Index i=k-1 ; i>=0 ; i--)
00287     {
00288       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00289       if(k-i-1>0)
00290         m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
00291       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00292       if(z==ComplexScalar(0))
00293       {
00294         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00295         // Use a small value instead, to prevent division by zero.
00296         internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00297       }
00298       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00299     }
00300   }
00301 
00302   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00303   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00304   // .. and normalize the eigenvectors
00305   for(Index k=0 ; k<n ; k++)
00306   {
00307     m_eivec.col(k).normalize();
00308   }
00309 }
00310 
00311 
00312 template<typename MatrixType>
00313 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00314 {
00315   const Index n =  m_eivalues.size();
00316   for (Index i=0; i<n; i++)
00317   {
00318     Index k;
00319     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00320     if (k != 0)
00321     {
00322       k += i;
00323       std::swap(m_eivalues[k],m_eivalues[i]);
00324       if(computeEigenvectors)
00325         m_eivec.col(i).swap(m_eivec.col(k));
00326     }
00327   }
00328 }
00329 
00330 
00331 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011