![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012, 2013 Chen-Pang He <jdh8@ms63.hinet.net> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_MATRIX_POWER 00011 #define EIGEN_MATRIX_POWER 00012 00013 namespace Eigen { 00014 00015 template<typename MatrixType> class MatrixPower; 00016 00030 /* TODO This class is only used by MatrixPower, so it should be nested 00031 * into MatrixPower, like MatrixPower::ReturnValue. However, my 00032 * compiler complained about unused template parameter in the 00033 * following declaration in namespace internal. 00034 * 00035 * template<typename MatrixType> 00036 * struct traits<MatrixPower<MatrixType>::ReturnValue>; 00037 */ 00038 template<typename MatrixType> 00039 class MatrixPowerParenthesesReturnValue : public ReturnByValue< MatrixPowerParenthesesReturnValue<MatrixType> > 00040 { 00041 public: 00042 typedef typename MatrixType::RealScalar RealScalar; 00043 typedef typename MatrixType::Index Index; 00044 00051 MatrixPowerParenthesesReturnValue(MatrixPower<MatrixType>& pow, RealScalar p) : m_pow(pow), m_p(p) 00052 { } 00053 00059 template<typename ResultType> 00060 inline void evalTo(ResultType& res) const 00061 { m_pow.compute(res, m_p); } 00062 00063 Index rows() const { return m_pow.rows(); } 00064 Index cols() const { return m_pow.cols(); } 00065 00066 private: 00067 MatrixPower<MatrixType>& m_pow; 00068 const RealScalar m_p; 00069 }; 00070 00086 template<typename MatrixType> 00087 class MatrixPowerAtomic : internal::noncopyable 00088 { 00089 private: 00090 enum { 00091 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00092 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 00093 }; 00094 typedef typename MatrixType::Scalar Scalar; 00095 typedef typename MatrixType::RealScalar RealScalar; 00096 typedef std::complex<RealScalar> ComplexScalar; 00097 typedef typename MatrixType::Index Index; 00098 typedef Block<MatrixType,Dynamic,Dynamic> ResultType; 00099 00100 const MatrixType& m_A; 00101 RealScalar m_p; 00102 00103 void computePade(int degree, const MatrixType& IminusT, ResultType& res) const; 00104 void compute2x2(ResultType& res, RealScalar p) const; 00105 void computeBig(ResultType& res) const; 00106 static int getPadeDegree(float normIminusT); 00107 static int getPadeDegree(double normIminusT); 00108 static int getPadeDegree(long double normIminusT); 00109 static ComplexScalar computeSuperDiag(const ComplexScalar&, const ComplexScalar&, RealScalar p); 00110 static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p); 00111 00112 public: 00124 MatrixPowerAtomic(const MatrixType& T, RealScalar p); 00125 00132 void compute(ResultType& res) const; 00133 }; 00134 00135 template<typename MatrixType> 00136 MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) : 00137 m_A(T), m_p(p) 00138 { 00139 eigen_assert(T.rows() == T.cols()); 00140 eigen_assert(p > -1 && p < 1); 00141 } 00142 00143 template<typename MatrixType> 00144 void MatrixPowerAtomic<MatrixType>::compute(ResultType& res) const 00145 { 00146 using std::pow; 00147 switch (m_A.rows()) { 00148 case 0: 00149 break; 00150 case 1: 00151 res(0,0) = pow(m_A(0,0), m_p); 00152 break; 00153 case 2: 00154 compute2x2(res, m_p); 00155 break; 00156 default: 00157 computeBig(res); 00158 } 00159 } 00160 00161 template<typename MatrixType> 00162 void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, ResultType& res) const 00163 { 00164 int i = 2*degree; 00165 res = (m_p-degree) / (2*i-2) * IminusT; 00166 00167 for (--i; i; --i) { 00168 res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>() 00169 .solve((i==1 ? -m_p : i&1 ? (-m_p-i/2)/(2*i) : (m_p-i/2)/(2*i-2)) * IminusT).eval(); 00170 } 00171 res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); 00172 } 00173 00174 // This function assumes that res has the correct size (see bug 614) 00175 template<typename MatrixType> 00176 void MatrixPowerAtomic<MatrixType>::compute2x2(ResultType& res, RealScalar p) const 00177 { 00178 using std::abs; 00179 using std::pow; 00180 res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); 00181 00182 for (Index i=1; i < m_A.cols(); ++i) { 00183 res.coeffRef(i,i) = pow(m_A.coeff(i,i), p); 00184 if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) 00185 res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1); 00186 else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) 00187 res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1)); 00188 else 00189 res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p); 00190 res.coeffRef(i-1,i) *= m_A.coeff(i-1,i); 00191 } 00192 } 00193 00194 template<typename MatrixType> 00195 void MatrixPowerAtomic<MatrixType>::computeBig(ResultType& res) const 00196 { 00197 using std::ldexp; 00198 const int digits = std::numeric_limits<RealScalar>::digits; 00199 const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1L // single precision 00200 : digits <= 53? 2.789358995219730e-1L // double precision 00201 : digits <= 64? 2.4471944416607995472e-1L // extended precision 00202 : digits <= 106? 1.1016843812851143391275867258512e-1L // double-double 00203 : 9.134603732914548552537150753385375e-2L; // quadruple precision 00204 MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>(); 00205 RealScalar normIminusT; 00206 int degree, degree2, numberOfSquareRoots = 0; 00207 bool hasExtraSquareRoot = false; 00208 00209 for (Index i=0; i < m_A.cols(); ++i) 00210 eigen_assert(m_A(i,i) != RealScalar(0)); 00211 00212 while (true) { 00213 IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T; 00214 normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); 00215 if (normIminusT < maxNormForPade) { 00216 degree = getPadeDegree(normIminusT); 00217 degree2 = getPadeDegree(normIminusT/2); 00218 if (degree - degree2 <= 1 || hasExtraSquareRoot) 00219 break; 00220 hasExtraSquareRoot = true; 00221 } 00222 matrix_sqrt_triangular(T, sqrtT); 00223 T = sqrtT.template triangularView<Upper>(); 00224 ++numberOfSquareRoots; 00225 } 00226 computePade(degree, IminusT, res); 00227 00228 for (; numberOfSquareRoots; --numberOfSquareRoots) { 00229 compute2x2(res, ldexp(m_p, -numberOfSquareRoots)); 00230 res = res.template triangularView<Upper>() * res; 00231 } 00232 compute2x2(res, m_p); 00233 } 00234 00235 template<typename MatrixType> 00236 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT) 00237 { 00238 const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f }; 00239 int degree = 3; 00240 for (; degree <= 4; ++degree) 00241 if (normIminusT <= maxNormForPade[degree - 3]) 00242 break; 00243 return degree; 00244 } 00245 00246 template<typename MatrixType> 00247 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT) 00248 { 00249 const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1, 00250 1.999045567181744e-1, 2.789358995219730e-1 }; 00251 int degree = 3; 00252 for (; degree <= 7; ++degree) 00253 if (normIminusT <= maxNormForPade[degree - 3]) 00254 break; 00255 return degree; 00256 } 00257 00258 template<typename MatrixType> 00259 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT) 00260 { 00261 #if LDBL_MANT_DIG == 53 00262 const int maxPadeDegree = 7; 00263 const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L, 00264 1.999045567181744e-1L, 2.789358995219730e-1L }; 00265 #elif LDBL_MANT_DIG <= 64 00266 const int maxPadeDegree = 8; 00267 const long double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L, 00268 6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L }; 00269 #elif LDBL_MANT_DIG <= 106 00270 const int maxPadeDegree = 10; 00271 const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ , 00272 1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L, 00273 2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L, 00274 1.1016843812851143391275867258512e-1L }; 00275 #else 00276 const int maxPadeDegree = 10; 00277 const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ , 00278 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L, 00279 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L, 00280 3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L, 00281 9.134603732914548552537150753385375e-2L }; 00282 #endif 00283 int degree = 3; 00284 for (; degree <= maxPadeDegree; ++degree) 00285 if (normIminusT <= maxNormForPade[degree - 3]) 00286 break; 00287 return degree; 00288 } 00289 00290 template<typename MatrixType> 00291 inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar 00292 MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p) 00293 { 00294 using std::ceil; 00295 using std::exp; 00296 using std::log; 00297 using std::sinh; 00298 00299 ComplexScalar logCurr = log(curr); 00300 ComplexScalar logPrev = log(prev); 00301 int unwindingNumber = ceil((numext::imag(logCurr - logPrev) - RealScalar(EIGEN_PI)) / RealScalar(2*EIGEN_PI)); 00302 ComplexScalar w = numext::log1p((curr-prev)/prev)/RealScalar(2) + ComplexScalar(0, EIGEN_PI*unwindingNumber); 00303 return RealScalar(2) * exp(RealScalar(0.5) * p * (logCurr + logPrev)) * sinh(p * w) / (curr - prev); 00304 } 00305 00306 template<typename MatrixType> 00307 inline typename MatrixPowerAtomic<MatrixType>::RealScalar 00308 MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p) 00309 { 00310 using std::exp; 00311 using std::log; 00312 using std::sinh; 00313 00314 RealScalar w = numext::log1p((curr-prev)/prev)/RealScalar(2); 00315 return 2 * exp(p * (log(curr) + log(prev)) / 2) * sinh(p * w) / (curr - prev); 00316 } 00317 00337 template<typename MatrixType> 00338 class MatrixPower : internal::noncopyable 00339 { 00340 private: 00341 typedef typename MatrixType::Scalar Scalar; 00342 typedef typename MatrixType::RealScalar RealScalar; 00343 typedef typename MatrixType::Index Index; 00344 00345 public: 00354 explicit MatrixPower(const MatrixType& A) : 00355 m_A(A), 00356 m_conditionNumber(0), 00357 m_rank(A.cols()), 00358 m_nulls(0) 00359 { eigen_assert(A.rows() == A.cols()); } 00360 00368 const MatrixPowerParenthesesReturnValue<MatrixType> operator()(RealScalar p) 00369 { return MatrixPowerParenthesesReturnValue<MatrixType>(*this, p); } 00370 00378 template<typename ResultType> 00379 void compute(ResultType& res, RealScalar p); 00380 00381 Index rows() const { return m_A.rows(); } 00382 Index cols() const { return m_A.cols(); } 00383 00384 private: 00385 typedef std::complex<RealScalar> ComplexScalar; 00386 typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, 00387 MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> ComplexMatrix; 00388 00390 typename MatrixType::Nested m_A; 00391 00393 MatrixType m_tmp; 00394 00396 ComplexMatrix m_T, m_U; 00397 00399 ComplexMatrix m_fT; 00400 00407 RealScalar m_conditionNumber; 00408 00410 Index m_rank; 00411 00413 Index m_nulls; 00414 00424 void split(RealScalar& p, RealScalar& intpart); 00425 00427 void initialize(); 00428 00429 template<typename ResultType> 00430 void computeIntPower(ResultType& res, RealScalar p); 00431 00432 template<typename ResultType> 00433 void computeFracPower(ResultType& res, RealScalar p); 00434 00435 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 00436 static void revertSchur( 00437 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 00438 const ComplexMatrix& T, 00439 const ComplexMatrix& U); 00440 00441 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 00442 static void revertSchur( 00443 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 00444 const ComplexMatrix& T, 00445 const ComplexMatrix& U); 00446 }; 00447 00448 template<typename MatrixType> 00449 template<typename ResultType> 00450 void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p) 00451 { 00452 using std::pow; 00453 switch (cols()) { 00454 case 0: 00455 break; 00456 case 1: 00457 res(0,0) = pow(m_A.coeff(0,0), p); 00458 break; 00459 default: 00460 RealScalar intpart; 00461 split(p, intpart); 00462 00463 res = MatrixType::Identity(rows(), cols()); 00464 computeIntPower(res, intpart); 00465 if (p) computeFracPower(res, p); 00466 } 00467 } 00468 00469 template<typename MatrixType> 00470 void MatrixPower<MatrixType>::split(RealScalar& p, RealScalar& intpart) 00471 { 00472 using std::floor; 00473 using std::pow; 00474 00475 intpart = floor(p); 00476 p -= intpart; 00477 00478 // Perform Schur decomposition if it is not yet performed and the power is 00479 // not an integer. 00480 if (!m_conditionNumber && p) 00481 initialize(); 00482 00483 // Choose the more stable of intpart = floor(p) and intpart = ceil(p). 00484 if (p > RealScalar(0.5) && p > (1-p) * pow(m_conditionNumber, p)) { 00485 --p; 00486 ++intpart; 00487 } 00488 } 00489 00490 template<typename MatrixType> 00491 void MatrixPower<MatrixType>::initialize() 00492 { 00493 const ComplexSchur<MatrixType> schurOfA(m_A); 00494 JacobiRotation<ComplexScalar> rot; 00495 ComplexScalar eigenvalue; 00496 00497 m_fT.resizeLike(m_A); 00498 m_T = schurOfA.matrixT(); 00499 m_U = schurOfA.matrixU(); 00500 m_conditionNumber = m_T.diagonal().array().abs().maxCoeff() / m_T.diagonal().array().abs().minCoeff(); 00501 00502 // Move zero eigenvalues to the bottom right corner. 00503 for (Index i = cols()-1; i>=0; --i) { 00504 if (m_rank <= 2) 00505 return; 00506 if (m_T.coeff(i,i) == RealScalar(0)) { 00507 for (Index j=i+1; j < m_rank; ++j) { 00508 eigenvalue = m_T.coeff(j,j); 00509 rot.makeGivens(m_T.coeff(j-1,j), eigenvalue); 00510 m_T.applyOnTheRight(j-1, j, rot); 00511 m_T.applyOnTheLeft(j-1, j, rot.adjoint()); 00512 m_T.coeffRef(j-1,j-1) = eigenvalue; 00513 m_T.coeffRef(j,j) = RealScalar(0); 00514 m_U.applyOnTheRight(j-1, j, rot); 00515 } 00516 --m_rank; 00517 } 00518 } 00519 00520 m_nulls = rows() - m_rank; 00521 if (m_nulls) { 00522 eigen_assert(m_T.bottomRightCorner(m_nulls, m_nulls).isZero() 00523 && "Base of matrix power should be invertible or with a semisimple zero eigenvalue."); 00524 m_fT.bottomRows(m_nulls).fill(RealScalar(0)); 00525 } 00526 } 00527 00528 template<typename MatrixType> 00529 template<typename ResultType> 00530 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p) 00531 { 00532 using std::abs; 00533 using std::fmod; 00534 RealScalar pp = abs(p); 00535 00536 if (p<0) 00537 m_tmp = m_A.inverse(); 00538 else 00539 m_tmp = m_A; 00540 00541 while (true) { 00542 if (fmod(pp, 2) >= 1) 00543 res = m_tmp * res; 00544 pp /= 2; 00545 if (pp < 1) 00546 break; 00547 m_tmp *= m_tmp; 00548 } 00549 } 00550 00551 template<typename MatrixType> 00552 template<typename ResultType> 00553 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p) 00554 { 00555 Block<ComplexMatrix,Dynamic,Dynamic> blockTp(m_fT, 0, 0, m_rank, m_rank); 00556 eigen_assert(m_conditionNumber); 00557 eigen_assert(m_rank + m_nulls == rows()); 00558 00559 MatrixPowerAtomic<ComplexMatrix>(m_T.topLeftCorner(m_rank, m_rank), p).compute(blockTp); 00560 if (m_nulls) { 00561 m_fT.topRightCorner(m_rank, m_nulls) = m_T.topLeftCorner(m_rank, m_rank).template triangularView<Upper>() 00562 .solve(blockTp * m_T.topRightCorner(m_rank, m_nulls)); 00563 } 00564 revertSchur(m_tmp, m_fT, m_U); 00565 res = m_tmp * res; 00566 } 00567 00568 template<typename MatrixType> 00569 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 00570 inline void MatrixPower<MatrixType>::revertSchur( 00571 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 00572 const ComplexMatrix& T, 00573 const ComplexMatrix& U) 00574 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); } 00575 00576 template<typename MatrixType> 00577 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 00578 inline void MatrixPower<MatrixType>::revertSchur( 00579 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 00580 const ComplexMatrix& T, 00581 const ComplexMatrix& U) 00582 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); } 00583 00597 template<typename Derived> 00598 class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> > 00599 { 00600 public: 00601 typedef typename Derived::PlainObject PlainObject; 00602 typedef typename Derived::RealScalar RealScalar; 00603 typedef typename Derived::Index Index; 00604 00611 MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p) 00612 { } 00613 00620 template<typename ResultType> 00621 inline void evalTo(ResultType& res) const 00622 { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); } 00623 00624 Index rows() const { return m_A.rows(); } 00625 Index cols() const { return m_A.cols(); } 00626 00627 private: 00628 const Derived& m_A; 00629 const RealScalar m_p; 00630 }; 00631 00645 template<typename Derived> 00646 class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerReturnValue<Derived> > 00647 { 00648 public: 00649 typedef typename Derived::PlainObject PlainObject; 00650 typedef typename std::complex<typename Derived::RealScalar> ComplexScalar; 00651 typedef typename Derived::Index Index; 00652 00659 MatrixComplexPowerReturnValue(const Derived& A, const ComplexScalar& p) : m_A(A), m_p(p) 00660 { } 00661 00671 template<typename ResultType> 00672 inline void evalTo(ResultType& res) const 00673 { res = (m_p * m_A.log()).exp(); } 00674 00675 Index rows() const { return m_A.rows(); } 00676 Index cols() const { return m_A.cols(); } 00677 00678 private: 00679 const Derived& m_A; 00680 const ComplexScalar m_p; 00681 }; 00682 00683 namespace internal { 00684 00685 template<typename MatrixPowerType> 00686 struct traits< MatrixPowerParenthesesReturnValue<MatrixPowerType> > 00687 { typedef typename MatrixPowerType::PlainObject ReturnType; }; 00688 00689 template<typename Derived> 00690 struct traits< MatrixPowerReturnValue<Derived> > 00691 { typedef typename Derived::PlainObject ReturnType; }; 00692 00693 template<typename Derived> 00694 struct traits< MatrixComplexPowerReturnValue<Derived> > 00695 { typedef typename Derived::PlainObject ReturnType; }; 00696 00697 } 00698 00699 template<typename Derived> 00700 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const 00701 { return MatrixPowerReturnValue<Derived>(derived(), p); } 00702 00703 template<typename Derived> 00704 const MatrixComplexPowerReturnValue<Derived> MatrixBase<Derived>::pow(const std::complex<RealScalar>& p) const 00705 { return MatrixComplexPowerReturnValue<Derived>(derived(), p); } 00706 00707 } // namespace Eigen 00708 00709 #endif // EIGEN_MATRIX_POWER