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:exported at Tue Jan 25 21:56:31 UTC 2011 |