1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012, 2013 Chen-Pang He <jdh8 (at) ms63.hinet.net> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_MATRIX_POWER 11 #define EIGEN_MATRIX_POWER 12 13 namespace Eigen { 14 15 template<typename MatrixType> class MatrixPower; 16 17 template<typename MatrixType> 18 class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval<MatrixType> > 19 { 20 public: 21 typedef typename MatrixType::RealScalar RealScalar; 22 typedef typename MatrixType::Index Index; 23 24 MatrixPowerRetval(MatrixPower<MatrixType>& pow, RealScalar p) : m_pow(pow), m_p(p) 25 { } 26 27 template<typename ResultType> 28 inline void evalTo(ResultType& res) const 29 { m_pow.compute(res, m_p); } 30 31 Index rows() const { return m_pow.rows(); } 32 Index cols() const { return m_pow.cols(); } 33 34 private: 35 MatrixPower<MatrixType>& m_pow; 36 const RealScalar m_p; 37 MatrixPowerRetval& operator=(const MatrixPowerRetval&); 38 }; 39 40 template<typename MatrixType> 41 class MatrixPowerAtomic 42 { 43 private: 44 enum { 45 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 46 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 47 }; 48 typedef typename MatrixType::Scalar Scalar; 49 typedef typename MatrixType::RealScalar RealScalar; 50 typedef std::complex<RealScalar> ComplexScalar; 51 typedef typename MatrixType::Index Index; 52 typedef Array<Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> ArrayType; 53 54 const MatrixType& m_A; 55 RealScalar m_p; 56 57 void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const; 58 void compute2x2(MatrixType& res, RealScalar p) const; 59 void computeBig(MatrixType& res) const; 60 static int getPadeDegree(float normIminusT); 61 static int getPadeDegree(double normIminusT); 62 static int getPadeDegree(long double normIminusT); 63 static ComplexScalar computeSuperDiag(const ComplexScalar&, const ComplexScalar&, RealScalar p); 64 static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p); 65 66 public: 67 MatrixPowerAtomic(const MatrixType& T, RealScalar p); 68 void compute(MatrixType& res) const; 69 }; 70 71 template<typename MatrixType> 72 MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) : 73 m_A(T), m_p(p) 74 { eigen_assert(T.rows() == T.cols()); } 75 76 template<typename MatrixType> 77 void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res) const 78 { 79 res.resizeLike(m_A); 80 switch (m_A.rows()) { 81 case 0: 82 break; 83 case 1: 84 res(0,0) = std::pow(m_A(0,0), m_p); 85 break; 86 case 2: 87 compute2x2(res, m_p); 88 break; 89 default: 90 computeBig(res); 91 } 92 } 93 94 template<typename MatrixType> 95 void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const 96 { 97 int i = degree<<1; 98 res = (m_p-degree) / ((i-1)<<1) * IminusT; 99 for (--i; i; --i) { 100 res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>() 101 .solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval(); 102 } 103 res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); 104 } 105 106 // This function assumes that res has the correct size (see bug 614) 107 template<typename MatrixType> 108 void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const 109 { 110 using std::abs; 111 using std::pow; 112 113 ArrayType logTdiag = m_A.diagonal().array().log(); 114 res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); 115 116 for (Index i=1; i < m_A.cols(); ++i) { 117 res.coeffRef(i,i) = pow(m_A.coeff(i,i), p); 118 if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) 119 res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1); 120 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))) 121 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)); 122 else 123 res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p); 124 res.coeffRef(i-1,i) *= m_A.coeff(i-1,i); 125 } 126 } 127 128 template<typename MatrixType> 129 void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) const 130 { 131 const int digits = std::numeric_limits<RealScalar>::digits; 132 const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision 133 digits <= 53? 2.789358995219730e-1: // double precision 134 digits <= 64? 2.4471944416607995472e-1L: // extended precision 135 digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double 136 9.134603732914548552537150753385375e-2L; // quadruple precision 137 MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>(); 138 RealScalar normIminusT; 139 int degree, degree2, numberOfSquareRoots = 0; 140 bool hasExtraSquareRoot = false; 141 142 /* FIXME 143 * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite 144 * loop. We should move 0 eigenvalues to bottom right corner. We need not 145 * worry about tiny values (e.g. 1e-300) because they will reach 1 if 146 * repetitively sqrt'ed. 147 * 148 * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the 149 * bottom right corner. 150 * 151 * [ T A ]^p [ T^p (T^-1 T^p A) ] 152 * [ ] = [ ] 153 * [ 0 0 ] [ 0 0 ] 154 */ 155 for (Index i=0; i < m_A.cols(); ++i) 156 eigen_assert(m_A(i,i) != RealScalar(0)); 157 158 while (true) { 159 IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T; 160 normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); 161 if (normIminusT < maxNormForPade) { 162 degree = getPadeDegree(normIminusT); 163 degree2 = getPadeDegree(normIminusT/2); 164 if (degree - degree2 <= 1 || hasExtraSquareRoot) 165 break; 166 hasExtraSquareRoot = true; 167 } 168 MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); 169 T = sqrtT.template triangularView<Upper>(); 170 ++numberOfSquareRoots; 171 } 172 computePade(degree, IminusT, res); 173 174 for (; numberOfSquareRoots; --numberOfSquareRoots) { 175 compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots)); 176 res = res.template triangularView<Upper>() * res; 177 } 178 compute2x2(res, m_p); 179 } 180 181 template<typename MatrixType> 182 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT) 183 { 184 const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f }; 185 int degree = 3; 186 for (; degree <= 4; ++degree) 187 if (normIminusT <= maxNormForPade[degree - 3]) 188 break; 189 return degree; 190 } 191 192 template<typename MatrixType> 193 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT) 194 { 195 const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1, 196 1.999045567181744e-1, 2.789358995219730e-1 }; 197 int degree = 3; 198 for (; degree <= 7; ++degree) 199 if (normIminusT <= maxNormForPade[degree - 3]) 200 break; 201 return degree; 202 } 203 204 template<typename MatrixType> 205 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT) 206 { 207 #if LDBL_MANT_DIG == 53 208 const int maxPadeDegree = 7; 209 const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L, 210 1.999045567181744e-1L, 2.789358995219730e-1L }; 211 #elif LDBL_MANT_DIG <= 64 212 const int maxPadeDegree = 8; 213 const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L, 214 6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L }; 215 #elif LDBL_MANT_DIG <= 106 216 const int maxPadeDegree = 10; 217 const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ , 218 1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L, 219 2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L, 220 1.1016843812851143391275867258512e-1L }; 221 #else 222 const int maxPadeDegree = 10; 223 const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ , 224 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L, 225 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L, 226 3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L, 227 9.134603732914548552537150753385375e-2L }; 228 #endif 229 int degree = 3; 230 for (; degree <= maxPadeDegree; ++degree) 231 if (normIminusT <= maxNormForPade[degree - 3]) 232 break; 233 return degree; 234 } 235 236 template<typename MatrixType> 237 inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar 238 MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p) 239 { 240 ComplexScalar logCurr = std::log(curr); 241 ComplexScalar logPrev = std::log(prev); 242 int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI)); 243 ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber); 244 return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev); 245 } 246 247 template<typename MatrixType> 248 inline typename MatrixPowerAtomic<MatrixType>::RealScalar 249 MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p) 250 { 251 RealScalar w = numext::atanh2(curr - prev, curr + prev); 252 return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev); 253 } 254 255 /** 256 * \ingroup MatrixFunctions_Module 257 * 258 * \brief Class for computing matrix powers. 259 * 260 * \tparam MatrixType type of the base, expected to be an instantiation 261 * of the Matrix class template. 262 * 263 * This class is capable of computing real/complex matrices raised to 264 * an arbitrary real power. Meanwhile, it saves the result of Schur 265 * decomposition if an non-integral power has even been calculated. 266 * Therefore, if you want to compute multiple (>= 2) matrix powers 267 * for the same matrix, using the class directly is more efficient than 268 * calling MatrixBase::pow(). 269 * 270 * Example: 271 * \include MatrixPower_optimal.cpp 272 * Output: \verbinclude MatrixPower_optimal.out 273 */ 274 template<typename MatrixType> 275 class MatrixPower 276 { 277 private: 278 enum { 279 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 280 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 281 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 282 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 283 }; 284 typedef typename MatrixType::Scalar Scalar; 285 typedef typename MatrixType::RealScalar RealScalar; 286 typedef typename MatrixType::Index Index; 287 288 public: 289 /** 290 * \brief Constructor. 291 * 292 * \param[in] A the base of the matrix power. 293 * 294 * The class stores a reference to A, so it should not be changed 295 * (or destroyed) before evaluation. 296 */ 297 explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0) 298 { eigen_assert(A.rows() == A.cols()); } 299 300 /** 301 * \brief Returns the matrix power. 302 * 303 * \param[in] p exponent, a real scalar. 304 * \return The expression \f$ A^p \f$, where A is specified in the 305 * constructor. 306 */ 307 const MatrixPowerRetval<MatrixType> operator()(RealScalar p) 308 { return MatrixPowerRetval<MatrixType>(*this, p); } 309 310 /** 311 * \brief Compute the matrix power. 312 * 313 * \param[in] p exponent, a real scalar. 314 * \param[out] res \f$ A^p \f$ where A is specified in the 315 * constructor. 316 */ 317 template<typename ResultType> 318 void compute(ResultType& res, RealScalar p); 319 320 Index rows() const { return m_A.rows(); } 321 Index cols() const { return m_A.cols(); } 322 323 private: 324 typedef std::complex<RealScalar> ComplexScalar; 325 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, MatrixType::Options, 326 MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrix; 327 328 typename MatrixType::Nested m_A; 329 MatrixType m_tmp; 330 ComplexMatrix m_T, m_U, m_fT; 331 RealScalar m_conditionNumber; 332 333 RealScalar modfAndInit(RealScalar, RealScalar*); 334 335 template<typename ResultType> 336 void computeIntPower(ResultType&, RealScalar); 337 338 template<typename ResultType> 339 void computeFracPower(ResultType&, RealScalar); 340 341 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 342 static void revertSchur( 343 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 344 const ComplexMatrix& T, 345 const ComplexMatrix& U); 346 347 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 348 static void revertSchur( 349 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 350 const ComplexMatrix& T, 351 const ComplexMatrix& U); 352 }; 353 354 template<typename MatrixType> 355 template<typename ResultType> 356 void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p) 357 { 358 switch (cols()) { 359 case 0: 360 break; 361 case 1: 362 res(0,0) = std::pow(m_A.coeff(0,0), p); 363 break; 364 default: 365 RealScalar intpart, x = modfAndInit(p, &intpart); 366 computeIntPower(res, intpart); 367 computeFracPower(res, x); 368 } 369 } 370 371 template<typename MatrixType> 372 typename MatrixPower<MatrixType>::RealScalar 373 MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart) 374 { 375 typedef Array<RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> RealArray; 376 377 *intpart = std::floor(x); 378 RealScalar res = x - *intpart; 379 380 if (!m_conditionNumber && res) { 381 const ComplexSchur<MatrixType> schurOfA(m_A); 382 m_T = schurOfA.matrixT(); 383 m_U = schurOfA.matrixU(); 384 385 const RealArray absTdiag = m_T.diagonal().array().abs(); 386 m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); 387 } 388 389 if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) { 390 --res; 391 ++*intpart; 392 } 393 return res; 394 } 395 396 template<typename MatrixType> 397 template<typename ResultType> 398 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p) 399 { 400 RealScalar pp = std::abs(p); 401 402 if (p<0) m_tmp = m_A.inverse(); 403 else m_tmp = m_A; 404 405 res = MatrixType::Identity(rows(), cols()); 406 while (pp >= 1) { 407 if (std::fmod(pp, 2) >= 1) 408 res = m_tmp * res; 409 m_tmp *= m_tmp; 410 pp /= 2; 411 } 412 } 413 414 template<typename MatrixType> 415 template<typename ResultType> 416 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p) 417 { 418 if (p) { 419 eigen_assert(m_conditionNumber); 420 MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT); 421 revertSchur(m_tmp, m_fT, m_U); 422 res = m_tmp * res; 423 } 424 } 425 426 template<typename MatrixType> 427 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 428 inline void MatrixPower<MatrixType>::revertSchur( 429 Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 430 const ComplexMatrix& T, 431 const ComplexMatrix& U) 432 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); } 433 434 template<typename MatrixType> 435 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols> 436 inline void MatrixPower<MatrixType>::revertSchur( 437 Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res, 438 const ComplexMatrix& T, 439 const ComplexMatrix& U) 440 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); } 441 442 /** 443 * \ingroup MatrixFunctions_Module 444 * 445 * \brief Proxy for the matrix power of some matrix (expression). 446 * 447 * \tparam Derived type of the base, a matrix (expression). 448 * 449 * This class holds the arguments to the matrix power until it is 450 * assigned or evaluated for some other reason (so the argument 451 * should not be changed in the meantime). It is the return type of 452 * MatrixBase::pow() and related functions and most of the 453 * time this is the only way it is used. 454 */ 455 template<typename Derived> 456 class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> > 457 { 458 public: 459 typedef typename Derived::PlainObject PlainObject; 460 typedef typename Derived::RealScalar RealScalar; 461 typedef typename Derived::Index Index; 462 463 /** 464 * \brief Constructor. 465 * 466 * \param[in] A %Matrix (expression), the base of the matrix power. 467 * \param[in] p scalar, the exponent of the matrix power. 468 */ 469 MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p) 470 { } 471 472 /** 473 * \brief Compute the matrix power. 474 * 475 * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the 476 * constructor. 477 */ 478 template<typename ResultType> 479 inline void evalTo(ResultType& res) const 480 { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); } 481 482 Index rows() const { return m_A.rows(); } 483 Index cols() const { return m_A.cols(); } 484 485 private: 486 const Derived& m_A; 487 const RealScalar m_p; 488 MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); 489 }; 490 491 namespace internal { 492 493 template<typename MatrixPowerType> 494 struct traits< MatrixPowerRetval<MatrixPowerType> > 495 { typedef typename MatrixPowerType::PlainObject ReturnType; }; 496 497 template<typename Derived> 498 struct traits< MatrixPowerReturnValue<Derived> > 499 { typedef typename Derived::PlainObject ReturnType; }; 500 501 } 502 503 template<typename Derived> 504 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const 505 { return MatrixPowerReturnValue<Derived>(derived(), p); } 506 507 } // namespace Eigen 508 509 #endif // EIGEN_MATRIX_POWER 510
