Main MRPT website > C++ reference
MRPT logo

SelfAdjointEigenSolver.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-2010 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_SELFADJOINTEIGENSOLVER_H
00027 #define EIGEN_SELFADJOINTEIGENSOLVER_H
00028 
00029 #include "./EigenvaluesCommon.h"
00030 #include "./Tridiagonalization.h"
00031 
00032 /** \eigenvalues_module \ingroup Eigenvalues_Module
00033   *
00034   *
00035   * \class SelfAdjointEigenSolver
00036   *
00037   * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices
00038   *
00039   * \tparam _MatrixType the type of the matrix of which we are computing the
00040   * eigendecomposition; this is expected to be an instantiation of the Matrix
00041   * class template.
00042   *
00043   * A matrix \f$ A \f$ is selfadjoint if it equals its adjoint. For real
00044   * matrices, this means that the matrix is symmetric: it equals its
00045   * transpose. This class computes the eigenvalues and eigenvectors of a
00046   * selfadjoint matrix. These are the scalars \f$ \lambda \f$ and vectors
00047   * \f$ v \f$ such that \f$ Av = \lambda v \f$.  The eigenvalues of a
00048   * selfadjoint matrix are always real. If \f$ D \f$ is a diagonal matrix with
00049   * the eigenvalues on the diagonal, and \f$ V \f$ is a matrix with the
00050   * eigenvectors as its columns, then \f$ A = V D V^{-1} \f$ (for selfadjoint
00051   * matrices, the matrix \f$ V \f$ is always invertible). This is called the
00052   * eigendecomposition.
00053   *
00054   * The algorithm exploits the fact that the matrix is selfadjoint, making it
00055   * faster and more accurate than the general purpose eigenvalue algorithms
00056   * implemented in EigenSolver and ComplexEigenSolver.
00057   *
00058   * Only the \b lower \b triangular \b part of the input matrix is referenced.
00059   *
00060   * Call the function compute() to compute the eigenvalues and eigenvectors of
00061   * a given matrix. Alternatively, you can use the
00062   * SelfAdjointEigenSolver(const MatrixType&, bool) constructor which computes
00063   * the eigenvalues and eigenvectors at construction time. Once the eigenvalue
00064   * and eigenvectors are computed, they can be retrieved with the eigenvalues()
00065   * and eigenvectors() functions.
00066   *
00067   * The documentation for SelfAdjointEigenSolver(const MatrixType&, bool)
00068   * contains an example of the typical use of this class.
00069   *
00070   * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and
00071   * the likes, see the class GeneralizedSelfAdjointEigenSolver.
00072   *
00073   * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver
00074   */
00075 template<typename _MatrixType> class SelfAdjointEigenSolver
00076 {
00077   public:
00078 
00079     typedef _MatrixType MatrixType;
00080     enum {
00081       Size = MatrixType::RowsAtCompileTime,
00082       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00083       Options = MatrixType::Options,
00084       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00085     };
00086 
00087     /** \brief Scalar type for matrices of type \p _MatrixType. */
00088     typedef typename MatrixType::Scalar Scalar;
00089     typedef typename MatrixType::Index Index;
00090 
00091     /** \brief Real scalar type for \p _MatrixType.
00092       *
00093       * This is just \c Scalar if #Scalar is real (e.g., \c float or
00094       * \c double), and the type of the real part of \c Scalar if #Scalar is
00095       * complex.
00096       */
00097     typedef typename NumTraits<Scalar>::Real RealScalar;
00098 
00099     /** \brief Type for vector of eigenvalues as returned by eigenvalues().
00100       *
00101       * This is a column vector with entries of type #RealScalar.
00102       * The length of the vector is the size of \p _MatrixType.
00103       */
00104     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
00105     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
00106 
00107     /** \brief Default constructor for fixed-size matrices.
00108       *
00109       * The default constructor is useful in cases in which the user intends to
00110       * perform decompositions via compute(const MatrixType&, bool) or
00111       * compute(const MatrixType&, const MatrixType&, bool). This constructor
00112       * can only be used if \p _MatrixType is a fixed-size matrix; use
00113       * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
00114       *
00115       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
00116       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
00117       */
00118     SelfAdjointEigenSolver()
00119         : m_eivec(),
00120           m_eivalues(),
00121           m_subdiag(),
00122           m_isInitialized(false)
00123     { }
00124 
00125     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
00126       *
00127       * \param [in]  size  Positive integer, size of the matrix whose
00128       * eigenvalues and eigenvectors will be computed.
00129       *
00130       * This constructor is useful for dynamic-size matrices, when the user
00131       * intends to perform decompositions via compute(const MatrixType&, bool)
00132       * or compute(const MatrixType&, const MatrixType&, bool). The \p size
00133       * parameter is only used as a hint. It is not an error to give a wrong
00134       * \p size, but it may impair performance.
00135       *
00136       * \sa compute(const MatrixType&, bool) for an example
00137       */
00138     SelfAdjointEigenSolver(Index size)
00139         : m_eivec(size, size),
00140           m_eivalues(size),
00141           m_subdiag(size > 1 ? size - 1 : 1),
00142           m_isInitialized(false)
00143     {}
00144 
00145     /** \brief Constructor; computes eigendecomposition of given matrix.
00146       *
00147       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
00148       *    be computed. Only the lower triangular part of the matrix is referenced.
00149       * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
00150       *
00151       * This constructor calls compute(const MatrixType&, bool) to compute the
00152       * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
00153       * \p options equals ComputeEigenvectors.
00154       *
00155       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
00156       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out
00157       *
00158       * \sa compute(const MatrixType&, bool),
00159       *     SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, bool)
00160       */
00161     SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
00162       : m_eivec(matrix.rows(), matrix.cols()),
00163         m_eivalues(matrix.cols()),
00164         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00165         m_isInitialized(false)
00166     {
00167       compute(matrix, options);
00168     }
00169 
00170     /** \brief Computes eigendecomposition of given matrix.
00171       *
00172       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
00173       *    be computed. Only the lower triangular part of the matrix is referenced.
00174       * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
00175       * \returns    Reference to \c *this
00176       *
00177       * This function computes the eigenvalues of \p matrix.  The eigenvalues()
00178       * function can be used to retrieve them.  If \p options equals ComputeEigenvectors,
00179       * then the eigenvectors are also computed and can be retrieved by
00180       * calling eigenvectors().
00181       *
00182       * This implementation uses a symmetric QR algorithm. The matrix is first
00183       * reduced to tridiagonal form using the Tridiagonalization class. The
00184       * tridiagonal matrix is then brought to diagonal form with implicit
00185       * symmetric QR steps with Wilkinson shift. Details can be found in
00186       * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>.
00187       *
00188       * The cost of the computation is about \f$ 9n^3 \f$ if the eigenvectors
00189       * are required and \f$ 4n^3/3 \f$ if they are not required.
00190       *
00191       * This method reuses the memory in the SelfAdjointEigenSolver object that
00192       * was allocated when the object was constructed, if the size of the
00193       * matrix does not change.
00194       *
00195       * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp
00196       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out
00197       *
00198       * \sa SelfAdjointEigenSolver(const MatrixType&, bool)
00199       */
00200     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
00201 
00202     /** \brief Returns the eigenvectors of given matrix (pencil).
00203       *
00204       * \returns  A const reference to the matrix whose columns are the eigenvectors.
00205       *
00206       * \pre The eigenvectors have been computed before.
00207       *
00208       * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
00209       * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
00210       * eigenvectors are normalized to have (Euclidean) norm equal to one. If
00211       * this object was used to solve the eigenproblem for the selfadjoint
00212       * matrix \f$ A \f$, then the matrix returned by this function is the
00213       * matrix \f$ V \f$ in the eigendecomposition \f$ A = V D V^{-1} \f$.
00214       *
00215       * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp
00216       * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out
00217       *
00218       * \sa eigenvalues()
00219       */
00220     const MatrixType& eigenvectors() const
00221     {
00222       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00223       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00224       return m_eivec;
00225     }
00226 
00227     /** \brief Returns the eigenvalues of given matrix (pencil).
00228       *
00229       * \returns A const reference to the column vector containing the eigenvalues.
00230       *
00231       * \pre The eigenvalues have been computed before.
00232       *
00233       * The eigenvalues are repeated according to their algebraic multiplicity,
00234       * so there are as many eigenvalues as rows in the matrix.
00235       *
00236       * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
00237       * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
00238       *
00239       * \sa eigenvectors(), MatrixBase::eigenvalues()
00240       */
00241     const RealVectorType& eigenvalues() const
00242     {
00243       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00244       return m_eivalues;
00245     }
00246 
00247     /** \brief Computes the positive-definite square root of the matrix.
00248       *
00249       * \returns the positive-definite square root of the matrix
00250       *
00251       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
00252       * have been computed before.
00253       *
00254       * The square root of a positive-definite matrix \f$ A \f$ is the
00255       * positive-definite matrix whose square equals \f$ A \f$. This function
00256       * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
00257       * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
00258       *
00259       * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp
00260       * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out
00261       *
00262       * \sa operatorInverseSqrt(),
00263       *     \ref MatrixFunctions_Module "MatrixFunctions Module"
00264       */
00265     MatrixType operatorSqrt() const
00266     {
00267       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00268       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00269       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00270     }
00271 
00272     /** \brief Computes the inverse square root of the matrix.
00273       *
00274       * \returns the inverse positive-definite square root of the matrix
00275       *
00276       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
00277       * have been computed before.
00278       *
00279       * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
00280       * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
00281       * cheaper than first computing the square root with operatorSqrt() and
00282       * then its inverse with MatrixBase::inverse().
00283       *
00284       * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp
00285       * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out
00286       *
00287       * \sa operatorSqrt(), MatrixBase::inverse(),
00288       *     \ref MatrixFunctions_Module "MatrixFunctions Module"
00289       */
00290     MatrixType operatorInverseSqrt() const
00291     {
00292       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00293       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00294       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00295     }
00296 
00297     /** \brief Reports whether previous computation was successful.
00298       *
00299       * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
00300       */
00301     ComputationInfo info() const
00302     {
00303       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00304       return m_info;
00305     }
00306 
00307     /** \brief Maximum number of iterations.
00308       *
00309       * Maximum number of iterations allowed for an eigenvalue to converge.
00310       */
00311     static const int m_maxIterations = 30;
00312 
00313   protected:
00314     MatrixType m_eivec;
00315     RealVectorType m_eivalues;
00316     typename TridiagonalizationType::SubDiagonalType m_subdiag;
00317     ComputationInfo m_info;
00318     bool m_isInitialized;
00319     bool m_eigenvectorsOk;
00320 };
00321 
00322 /** \internal
00323   *
00324   * \eigenvalues_module \ingroup Eigenvalues_Module
00325   *
00326   * Performs a QR step on a tridiagonal symmetric matrix represented as a
00327   * pair of two vectors \a diag and \a subdiag.
00328   *
00329   * \param matA the input selfadjoint matrix
00330   * \param hCoeffs returned Householder coefficients
00331   *
00332   * For compilation efficiency reasons, this procedure does not use eigen expression
00333   * for its arguments.
00334   *
00335   * Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
00336   * "implicit symmetric QR step with Wilkinson shift"
00337   */
00338 namespace internal {
00339 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00340 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
00341 }
00342 
00343 template<typename MatrixType>
00344 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00345 ::compute(const MatrixType& matrix, int options)
00346 {
00347   eigen_assert(matrix.cols() == matrix.rows());
00348   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00349           && (options&EigVecMask)!=EigVecMask
00350           && "invalid option parameter");
00351   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00352   Index n = matrix.cols();
00353   m_eivalues.resize(n,1);
00354 
00355   if(n==1)
00356   {
00357     m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
00358     if(computeEigenvectors)
00359       m_eivec.setOnes();
00360     m_info = Success;
00361     m_isInitialized = true;
00362     m_eigenvectorsOk = computeEigenvectors;
00363     return *this;
00364   }
00365 
00366   // declare some aliases
00367   RealVectorType& diag = m_eivalues;
00368   MatrixType& mat = m_eivec;
00369 
00370   mat = matrix;
00371   m_subdiag.resize(n-1);
00372   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
00373 
00374   Index end = n-1;
00375   Index start = 0;
00376   Index iter = 0; // number of iterations we are working on one element
00377 
00378   while (end>0)
00379   {
00380     for (Index i = start; i<end; ++i)
00381       if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
00382         m_subdiag[i] = 0;
00383 
00384     // find the largest unreduced block
00385     while (end>0 && m_subdiag[end-1]==0)
00386     {
00387       iter = 0;
00388       end--;
00389     }
00390     if (end<=0)
00391       break;
00392 
00393     // if we spent too many iterations on the current element, we give up
00394     iter++;
00395     if(iter > m_maxIterations) break;
00396 
00397     start = end - 1;
00398     while (start>0 && m_subdiag[start-1]!=0)
00399       start--;
00400 
00401     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
00402   }
00403 
00404   if (iter <= m_maxIterations)
00405     m_info = Success;
00406   else
00407     m_info = NoConvergence;
00408 
00409   // Sort eigenvalues and corresponding vectors.
00410   // TODO make the sort optional ?
00411   // TODO use a better sort algorithm !!
00412   if (m_info == Success)
00413   {
00414     for (Index i = 0; i < n-1; ++i)
00415     {
00416       Index k;
00417       m_eivalues.segment(i,n-i).minCoeff(&k);
00418       if (k > 0)
00419       {
00420         std::swap(m_eivalues[i], m_eivalues[k+i]);
00421         if(computeEigenvectors)
00422           m_eivec.col(i).swap(m_eivec.col(k+i));
00423       }
00424     }
00425   }
00426 
00427   m_isInitialized = true;
00428   m_eigenvectorsOk = computeEigenvectors;
00429   return *this;
00430 }
00431 
00432 namespace internal {
00433 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00434 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
00435 {
00436   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
00437   RealScalar e2 = abs2(subdiag[end-1]);
00438   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
00439   RealScalar x = diag[start] - mu;
00440   RealScalar z = subdiag[start];
00441 
00442   for (Index k = start; k < end; ++k)
00443   {
00444     JacobiRotation<RealScalar> rot;
00445     rot.makeGivens(x, z);
00446 
00447     // do T = G' T G
00448     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
00449     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
00450 
00451     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
00452     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
00453     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
00454 
00455     if (k > start)
00456       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
00457 
00458     x = subdiag[k];
00459 
00460     if (k < end - 1)
00461     {
00462       z = -rot.s() * subdiag[k+1];
00463       subdiag[k + 1] = rot.c() * subdiag[k+1];
00464     }
00465 
00466     // apply the givens rotation to the unit matrix Q = Q * G
00467     if (matrixQ)
00468     {
00469       // FIXME if StorageOrder == RowMajor this operation is not very efficient
00470       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
00471       q.applyOnTheRight(k,k+1,rot);
00472     }
00473   }
00474 }
00475 } // end namespace internal
00476 
00477 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011