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 template<typename _MatrixType>
00033 class GeneralizedSelfAdjointEigenSolver;
00034 
00078 template<typename _MatrixType> class SelfAdjointEigenSolver
00079 {
00080   public:
00081 
00082     typedef _MatrixType MatrixType;
00083     enum {
00084       Size = MatrixType::RowsAtCompileTime,
00085       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00086       Options = MatrixType::Options,
00087       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00088     };
00089 
00091     typedef typename MatrixType::Scalar Scalar;
00092     typedef typename MatrixType::Index Index;
00093 
00100     typedef typename NumTraits<Scalar>::Real RealScalar;
00101 
00107     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
00108     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
00109 
00120     SelfAdjointEigenSolver()
00121         : m_eivec(),
00122           m_eivalues(),
00123           m_subdiag(),
00124           m_isInitialized(false)
00125     { }
00126 
00139     SelfAdjointEigenSolver(Index size)
00140         : m_eivec(size, size),
00141           m_eivalues(size),
00142           m_subdiag(size > 1 ? size - 1 : 1),
00143           m_isInitialized(false)
00144     {}
00145 
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 
00200     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
00201 
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 
00242     const RealVectorType& eigenvalues() const
00243     {
00244       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00245       return m_eivalues;
00246     }
00247 
00266     MatrixType operatorSqrt() const
00267     {
00268       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00269       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00270       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00271     }
00272 
00291     MatrixType operatorInverseSqrt() const
00292     {
00293       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00294       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00295       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00296     }
00297 
00302     ComputationInfo info() const
00303     {
00304       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00305       return m_info;
00306     }
00307 
00312     static const int m_maxIterations = 30;
00313 
00314     #ifdef EIGEN2_SUPPORT
00315     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
00316       : m_eivec(matrix.rows(), matrix.cols()),
00317         m_eivalues(matrix.cols()),
00318         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00319         m_isInitialized(false)
00320     {
00321       compute(matrix, computeEigenvectors);
00322     }
00323     
00324     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00325         : m_eivec(matA.cols(), matA.cols()),
00326           m_eivalues(matA.cols()),
00327           m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
00328           m_isInitialized(false)
00329     {
00330       static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00331     }
00332     
00333     void compute(const MatrixType& matrix, bool computeEigenvectors)
00334     {
00335       compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00336     }
00337 
00338     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00339     {
00340       compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00341     }
00342     #endif // EIGEN2_SUPPORT
00343 
00344   protected:
00345     MatrixType m_eivec;
00346     RealVectorType m_eivalues;
00347     typename TridiagonalizationType::SubDiagonalType m_subdiag;
00348     ComputationInfo m_info;
00349     bool m_isInitialized;
00350     bool m_eigenvectorsOk;
00351 };
00352 
00369 namespace internal {
00370 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00371 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
00372 }
00373 
00374 template<typename MatrixType>
00375 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00376 ::compute(const MatrixType& matrix, int options)
00377 {
00378   eigen_assert(matrix.cols() == matrix.rows());
00379   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00380           && (options&EigVecMask)!=EigVecMask
00381           && "invalid option parameter");
00382   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00383   Index n = matrix.cols();
00384   m_eivalues.resize(n,1);
00385 
00386   if(n==1)
00387   {
00388     m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
00389     if(computeEigenvectors)
00390       m_eivec.setOnes(n,n);
00391     m_info = Success;
00392     m_isInitialized = true;
00393     m_eigenvectorsOk = computeEigenvectors;
00394     return *this;
00395   }
00396 
00397   // declare some aliases
00398   RealVectorType& diag = m_eivalues;
00399   MatrixType& mat = m_eivec;
00400 
00401   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00402   RealScalar scale = matrix.cwiseAbs().maxCoeff();
00403   if(scale==Scalar(0)) scale = 1;
00404   mat = matrix / scale;
00405   m_subdiag.resize(n-1);
00406   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
00407   
00408   Index end = n-1;
00409   Index start = 0;
00410   Index iter = 0; // number of iterations we are working on one element
00411 
00412   while (end>0)
00413   {
00414     for (Index i = start; i<end; ++i)
00415       if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
00416         m_subdiag[i] = 0;
00417 
00418     // find the largest unreduced block
00419     while (end>0 && m_subdiag[end-1]==0)
00420     {
00421       iter = 0;
00422       end--;
00423     }
00424     if (end<=0)
00425       break;
00426 
00427     // if we spent too many iterations on the current element, we give up
00428     iter++;
00429     if(iter > m_maxIterations) break;
00430 
00431     start = end - 1;
00432     while (start>0 && m_subdiag[start-1]!=0)
00433       start--;
00434 
00435     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
00436   }
00437 
00438   if (iter <= m_maxIterations)
00439     m_info = Success;
00440   else
00441     m_info = NoConvergence;
00442 
00443   // Sort eigenvalues and corresponding vectors.
00444   // TODO make the sort optional ?
00445   // TODO use a better sort algorithm !!
00446   if (m_info == Success)
00447   {
00448     for (Index i = 0; i < n-1; ++i)
00449     {
00450       Index k;
00451       m_eivalues.segment(i,n-i).minCoeff(&k);
00452       if (k > 0)
00453       {
00454         std::swap(m_eivalues[i], m_eivalues[k+i]);
00455         if(computeEigenvectors)
00456           m_eivec.col(i).swap(m_eivec.col(k+i));
00457       }
00458     }
00459   }
00460   
00461   // scale back the eigen values
00462   m_eivalues *= scale;
00463 
00464   m_isInitialized = true;
00465   m_eigenvectorsOk = computeEigenvectors;
00466   return *this;
00467 }
00468 
00469 namespace internal {
00470 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00471 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
00472 {
00473   // NOTE this version avoids over & underflow, however since the matrix is prescaled, overflow cannot occur,
00474   // and underflows should be meaningless anyway. So I don't any reason to enable this version, but I keep
00475   // it here for reference:
00476 //   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
00477 //   RealScalar e = subdiag[end-1];
00478 //   RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
00479   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
00480   RealScalar e2 = abs2(subdiag[end-1]);
00481   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
00482   RealScalar x = diag[start] - mu;
00483   RealScalar z = subdiag[start];
00484   for (Index k = start; k < end; ++k)
00485   {
00486     JacobiRotation<RealScalar> rot;
00487     rot.makeGivens(x, z);
00488 
00489     // do T = G' T G
00490     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
00491     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
00492 
00493     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
00494     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
00495     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
00496     
00497 
00498     if (k > start)
00499       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
00500 
00501     x = subdiag[k];
00502 
00503     if (k < end - 1)
00504     {
00505       z = -rot.s() * subdiag[k+1];
00506       subdiag[k + 1] = rot.c() * subdiag[k+1];
00507     }
00508     
00509     // apply the givens rotation to the unit matrix Q = Q * G
00510     if (matrixQ)
00511     {
00512       // FIXME if StorageOrder == RowMajor this operation is not very efficient
00513       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
00514       q.applyOnTheRight(k,k+1,rot);
00515     }
00516   }
00517 }
00518 } // end namespace internal
00519 
00520 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:20