![]() |
Eigen
3.3.3
|
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,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_EIGENSOLVER_H 00012 #define EIGEN_EIGENSOLVER_H 00013 00014 #include "./RealSchur.h" 00015 00016 namespace Eigen { 00017 00064 template<typename _MatrixType> class EigenSolver 00065 { 00066 public: 00067 00069 typedef _MatrixType MatrixType; 00070 00071 enum { 00072 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00073 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00074 Options = MatrixType::Options, 00075 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00076 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00077 }; 00078 00080 typedef typename MatrixType::Scalar Scalar; 00081 typedef typename NumTraits<Scalar>::Real RealScalar; 00082 typedef Eigen::Index Index; 00083 00090 typedef std::complex<RealScalar> ComplexScalar; 00091 00097 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 00098 00104 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 00105 00113 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {} 00114 00121 explicit EigenSolver(Index size) 00122 : m_eivec(size, size), 00123 m_eivalues(size), 00124 m_isInitialized(false), 00125 m_eigenvectorsOk(false), 00126 m_realSchur(size), 00127 m_matT(size, size), 00128 m_tmp(size) 00129 {} 00130 00146 template<typename InputType> 00147 explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) 00148 : m_eivec(matrix.rows(), matrix.cols()), 00149 m_eivalues(matrix.cols()), 00150 m_isInitialized(false), 00151 m_eigenvectorsOk(false), 00152 m_realSchur(matrix.cols()), 00153 m_matT(matrix.rows(), matrix.cols()), 00154 m_tmp(matrix.cols()) 00155 { 00156 compute(matrix.derived(), computeEigenvectors); 00157 } 00158 00179 EigenvectorsType eigenvectors() const; 00180 00199 const MatrixType& pseudoEigenvectors() const 00200 { 00201 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00202 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00203 return m_eivec; 00204 } 00205 00224 MatrixType pseudoEigenvalueMatrix() const; 00225 00244 const EigenvalueType& eigenvalues() const 00245 { 00246 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00247 return m_eivalues; 00248 } 00249 00277 template<typename InputType> 00278 EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); 00279 00281 ComputationInfo info() const 00282 { 00283 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00284 return m_info; 00285 } 00286 00288 EigenSolver& setMaxIterations(Index maxIters) 00289 { 00290 m_realSchur.setMaxIterations(maxIters); 00291 return *this; 00292 } 00293 00295 Index getMaxIterations() 00296 { 00297 return m_realSchur.getMaxIterations(); 00298 } 00299 00300 private: 00301 void doComputeEigenvectors(); 00302 00303 protected: 00304 00305 static void check_template_parameters() 00306 { 00307 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00308 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 00309 } 00310 00311 MatrixType m_eivec; 00312 EigenvalueType m_eivalues; 00313 bool m_isInitialized; 00314 bool m_eigenvectorsOk; 00315 ComputationInfo m_info; 00316 RealSchur<MatrixType> m_realSchur; 00317 MatrixType m_matT; 00318 00319 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 00320 ColumnVectorType m_tmp; 00321 }; 00322 00323 template<typename MatrixType> 00324 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const 00325 { 00326 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00327 const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); 00328 Index n = m_eivalues.rows(); 00329 MatrixType matD = MatrixType::Zero(n,n); 00330 for (Index i=0; i<n; ++i) 00331 { 00332 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision)) 00333 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i)); 00334 else 00335 { 00336 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)), 00337 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)); 00338 ++i; 00339 } 00340 } 00341 return matD; 00342 } 00343 00344 template<typename MatrixType> 00345 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const 00346 { 00347 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00348 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00349 const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); 00350 Index n = m_eivec.cols(); 00351 EigenvectorsType matV(n,n); 00352 for (Index j=0; j<n; ++j) 00353 { 00354 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n) 00355 { 00356 // we have a real eigen value 00357 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); 00358 matV.col(j).normalize(); 00359 } 00360 else 00361 { 00362 // we have a pair of complex eigen values 00363 for (Index i=0; i<n; ++i) 00364 { 00365 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); 00366 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); 00367 } 00368 matV.col(j).normalize(); 00369 matV.col(j+1).normalize(); 00370 ++j; 00371 } 00372 } 00373 return matV; 00374 } 00375 00376 template<typename MatrixType> 00377 template<typename InputType> 00378 EigenSolver<MatrixType>& 00379 EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) 00380 { 00381 check_template_parameters(); 00382 00383 using std::sqrt; 00384 using std::abs; 00385 using numext::isfinite; 00386 eigen_assert(matrix.cols() == matrix.rows()); 00387 00388 // Reduce to real Schur form. 00389 m_realSchur.compute(matrix.derived(), computeEigenvectors); 00390 00391 m_info = m_realSchur.info(); 00392 00393 if (m_info == Success) 00394 { 00395 m_matT = m_realSchur.matrixT(); 00396 if (computeEigenvectors) 00397 m_eivec = m_realSchur.matrixU(); 00398 00399 // Compute eigenvalues from matT 00400 m_eivalues.resize(matrix.cols()); 00401 Index i = 0; 00402 while (i < matrix.cols()) 00403 { 00404 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 00405 { 00406 m_eivalues.coeffRef(i) = m_matT.coeff(i, i); 00407 if(!(isfinite)(m_eivalues.coeffRef(i))) 00408 { 00409 m_isInitialized = true; 00410 m_eigenvectorsOk = false; 00411 m_info = NumericalIssue; 00412 return *this; 00413 } 00414 ++i; 00415 } 00416 else 00417 { 00418 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); 00419 Scalar z; 00420 // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); 00421 // without overflow 00422 { 00423 Scalar t0 = m_matT.coeff(i+1, i); 00424 Scalar t1 = m_matT.coeff(i, i+1); 00425 Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1))); 00426 t0 /= maxval; 00427 t1 /= maxval; 00428 Scalar p0 = p/maxval; 00429 z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); 00430 } 00431 00432 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); 00433 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); 00434 if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1)))) 00435 { 00436 m_isInitialized = true; 00437 m_eigenvectorsOk = false; 00438 m_info = NumericalIssue; 00439 return *this; 00440 } 00441 i += 2; 00442 } 00443 } 00444 00445 // Compute eigenvectors. 00446 if (computeEigenvectors) 00447 doComputeEigenvectors(); 00448 } 00449 00450 m_isInitialized = true; 00451 m_eigenvectorsOk = computeEigenvectors; 00452 00453 return *this; 00454 } 00455 00456 00457 template<typename MatrixType> 00458 void EigenSolver<MatrixType>::doComputeEigenvectors() 00459 { 00460 using std::abs; 00461 const Index size = m_eivec.cols(); 00462 const Scalar eps = NumTraits<Scalar>::epsilon(); 00463 00464 // inefficient! this is already computed in RealSchur 00465 Scalar norm(0); 00466 for (Index j = 0; j < size; ++j) 00467 { 00468 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 00469 } 00470 00471 // Backsubstitute to find vectors of upper triangular form 00472 if (norm == Scalar(0)) 00473 { 00474 return; 00475 } 00476 00477 for (Index n = size-1; n >= 0; n--) 00478 { 00479 Scalar p = m_eivalues.coeff(n).real(); 00480 Scalar q = m_eivalues.coeff(n).imag(); 00481 00482 // Scalar vector 00483 if (q == Scalar(0)) 00484 { 00485 Scalar lastr(0), lastw(0); 00486 Index l = n; 00487 00488 m_matT.coeffRef(n,n) = Scalar(1); 00489 for (Index i = n-1; i >= 0; i--) 00490 { 00491 Scalar w = m_matT.coeff(i,i) - p; 00492 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 00493 00494 if (m_eivalues.coeff(i).imag() < Scalar(0)) 00495 { 00496 lastw = w; 00497 lastr = r; 00498 } 00499 else 00500 { 00501 l = i; 00502 if (m_eivalues.coeff(i).imag() == Scalar(0)) 00503 { 00504 if (w != Scalar(0)) 00505 m_matT.coeffRef(i,n) = -r / w; 00506 else 00507 m_matT.coeffRef(i,n) = -r / (eps * norm); 00508 } 00509 else // Solve real equations 00510 { 00511 Scalar x = m_matT.coeff(i,i+1); 00512 Scalar y = m_matT.coeff(i+1,i); 00513 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); 00514 Scalar t = (x * lastr - lastw * r) / denom; 00515 m_matT.coeffRef(i,n) = t; 00516 if (abs(x) > abs(lastw)) 00517 m_matT.coeffRef(i+1,n) = (-r - w * t) / x; 00518 else 00519 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; 00520 } 00521 00522 // Overflow control 00523 Scalar t = abs(m_matT.coeff(i,n)); 00524 if ((eps * t) * t > Scalar(1)) 00525 m_matT.col(n).tail(size-i) /= t; 00526 } 00527 } 00528 } 00529 else if (q < Scalar(0) && n > 0) // Complex vector 00530 { 00531 Scalar lastra(0), lastsa(0), lastw(0); 00532 Index l = n-1; 00533 00534 // Last vector component imaginary so matrix is triangular 00535 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n))) 00536 { 00537 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); 00538 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); 00539 } 00540 else 00541 { 00542 ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q); 00543 m_matT.coeffRef(n-1,n-1) = numext::real(cc); 00544 m_matT.coeffRef(n-1,n) = numext::imag(cc); 00545 } 00546 m_matT.coeffRef(n,n-1) = Scalar(0); 00547 m_matT.coeffRef(n,n) = Scalar(1); 00548 for (Index i = n-2; i >= 0; i--) 00549 { 00550 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); 00551 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 00552 Scalar w = m_matT.coeff(i,i) - p; 00553 00554 if (m_eivalues.coeff(i).imag() < Scalar(0)) 00555 { 00556 lastw = w; 00557 lastra = ra; 00558 lastsa = sa; 00559 } 00560 else 00561 { 00562 l = i; 00563 if (m_eivalues.coeff(i).imag() == RealScalar(0)) 00564 { 00565 ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q); 00566 m_matT.coeffRef(i,n-1) = numext::real(cc); 00567 m_matT.coeffRef(i,n) = numext::imag(cc); 00568 } 00569 else 00570 { 00571 // Solve complex equations 00572 Scalar x = m_matT.coeff(i,i+1); 00573 Scalar y = m_matT.coeff(i+1,i); 00574 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; 00575 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 00576 if ((vr == Scalar(0)) && (vi == Scalar(0))) 00577 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); 00578 00579 ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi); 00580 m_matT.coeffRef(i,n-1) = numext::real(cc); 00581 m_matT.coeffRef(i,n) = numext::imag(cc); 00582 if (abs(x) > (abs(lastw) + abs(q))) 00583 { 00584 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; 00585 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; 00586 } 00587 else 00588 { 00589 cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q); 00590 m_matT.coeffRef(i+1,n-1) = numext::real(cc); 00591 m_matT.coeffRef(i+1,n) = numext::imag(cc); 00592 } 00593 } 00594 00595 // Overflow control 00596 Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); 00597 if ((eps * t) * t > Scalar(1)) 00598 m_matT.block(i, n-1, size-i, 2) /= t; 00599 00600 } 00601 } 00602 00603 // We handled a pair of complex conjugate eigenvalues, so need to skip them both 00604 n--; 00605 } 00606 else 00607 { 00608 eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen 00609 } 00610 } 00611 00612 // Back transformation to get eigenvectors of original matrix 00613 for (Index j = size-1; j >= 0; j--) 00614 { 00615 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1); 00616 m_eivec.col(j) = m_tmp; 00617 } 00618 } 00619 00620 } // end namespace Eigen 00621 00622 #endif // EIGEN_EIGENSOLVER_H