MatrixPower.h
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
 All Classes Functions Variables Typedefs Enumerator