MatrixPower.h
Go to the documentation of this file.
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 
00017 template<typename MatrixType>
00018 class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval<MatrixType> >
00019 {
00020   public:
00021     typedef typename MatrixType::RealScalar RealScalar;
00022     typedef typename MatrixType::Index Index;
00023 
00024     MatrixPowerRetval(MatrixPower<MatrixType>& pow, RealScalar p) : m_pow(pow), m_p(p)
00025     { }
00026 
00027     template<typename ResultType>
00028     inline void evalTo(ResultType& res) const
00029     { m_pow.compute(res, m_p); }
00030 
00031     Index rows() const { return m_pow.rows(); }
00032     Index cols() const { return m_pow.cols(); }
00033 
00034   private:
00035     MatrixPower<MatrixType>& m_pow;
00036     const RealScalar m_p;
00037     MatrixPowerRetval& operator=(const MatrixPowerRetval&);
00038 };
00039 
00040 template<typename MatrixType>
00041 class MatrixPowerAtomic
00042 {
00043   private:
00044     enum {
00045       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00046       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
00047     };
00048     typedef typename MatrixType::Scalar Scalar;
00049     typedef typename MatrixType::RealScalar RealScalar;
00050     typedef std::complex<RealScalar> ComplexScalar;
00051     typedef typename MatrixType::Index Index;
00052     typedef Array<Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> ArrayType;
00053 
00054     const MatrixType& m_A;
00055     RealScalar m_p;
00056 
00057     void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const;
00058     void compute2x2(MatrixType& res, RealScalar p) const;
00059     void computeBig(MatrixType& res) const;
00060     static int getPadeDegree(float normIminusT);
00061     static int getPadeDegree(double normIminusT);
00062     static int getPadeDegree(long double normIminusT);
00063     static ComplexScalar computeSuperDiag(const ComplexScalar&, const ComplexScalar&, RealScalar p);
00064     static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p);
00065 
00066   public:
00067     MatrixPowerAtomic(const MatrixType& T, RealScalar p);
00068     void compute(MatrixType& res) const;
00069 };
00070 
00071 template<typename MatrixType>
00072 MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) :
00073   m_A(T), m_p(p)
00074 { eigen_assert(T.rows() == T.cols()); }
00075 
00076 template<typename MatrixType>
00077 void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res) const
00078 {
00079   res.resizeLike(m_A);
00080   switch (m_A.rows()) {
00081     case 0:
00082       break;
00083     case 1:
00084       res(0,0) = std::pow(m_A(0,0), m_p);
00085       break;
00086     case 2:
00087       compute2x2(res, m_p);
00088       break;
00089     default:
00090       computeBig(res);
00091   }
00092 }
00093 
00094 template<typename MatrixType>
00095 void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const
00096 {
00097   int i = degree<<1;
00098   res = (m_p-degree) / ((i-1)<<1) * IminusT;
00099   for (--i; i; --i) {
00100     res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
00101         .solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval();
00102   }
00103   res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
00104 }
00105 
00106 // This function assumes that res has the correct size (see bug 614)
00107 template<typename MatrixType>
00108 void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
00109 {
00110   using std::abs;
00111   using std::pow;
00112   
00113   ArrayType logTdiag = m_A.diagonal().array().log();
00114   res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
00115 
00116   for (Index i=1; i < m_A.cols(); ++i) {
00117     res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
00118     if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
00119       res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
00120     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)))
00121       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));
00122     else
00123       res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
00124     res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
00125   }
00126 }
00127 
00128 template<typename MatrixType>
00129 void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) const
00130 {
00131   const int digits = std::numeric_limits<RealScalar>::digits;
00132   const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
00133                                     digits <=  53? 2.789358995219730e-1:                    // double precision
00134                                     digits <=  64? 2.4471944416607995472e-1L:               // extended precision
00135                                     digits <= 106? 1.1016843812851143391275867258512e-1L:   // double-double
00136                                                    9.134603732914548552537150753385375e-2L; // quadruple precision
00137   MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
00138   RealScalar normIminusT;
00139   int degree, degree2, numberOfSquareRoots = 0;
00140   bool hasExtraSquareRoot = false;
00141 
00142   /* FIXME
00143    * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
00144    * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
00145    * worry about tiny values (e.g. 1e-300) because they will reach 1 if
00146    * repetitively sqrt'ed.
00147    *
00148    * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the
00149    * bottom right corner.
00150    *
00151    * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
00152    * [      ]   = [                   ]
00153    * [ 0  0 ]     [  0         0      ]
00154    */
00155   for (Index i=0; i < m_A.cols(); ++i)
00156     eigen_assert(m_A(i,i) != RealScalar(0));
00157 
00158   while (true) {
00159     IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
00160     normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
00161     if (normIminusT < maxNormForPade) {
00162       degree = getPadeDegree(normIminusT);
00163       degree2 = getPadeDegree(normIminusT/2);
00164       if (degree - degree2 <= 1 || hasExtraSquareRoot)
00165         break;
00166       hasExtraSquareRoot = true;
00167     }
00168     MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
00169     T = sqrtT.template triangularView<Upper>();
00170     ++numberOfSquareRoots;
00171   }
00172   computePade(degree, IminusT, res);
00173 
00174   for (; numberOfSquareRoots; --numberOfSquareRoots) {
00175     compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
00176     res = res.template triangularView<Upper>() * res;
00177   }
00178   compute2x2(res, m_p);
00179 }
00180   
00181 template<typename MatrixType>
00182 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT)
00183 {
00184   const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
00185   int degree = 3;
00186   for (; degree <= 4; ++degree)
00187     if (normIminusT <= maxNormForPade[degree - 3])
00188       break;
00189   return degree;
00190 }
00191 
00192 template<typename MatrixType>
00193 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT)
00194 {
00195   const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
00196       1.999045567181744e-1, 2.789358995219730e-1 };
00197   int degree = 3;
00198   for (; degree <= 7; ++degree)
00199     if (normIminusT <= maxNormForPade[degree - 3])
00200       break;
00201   return degree;
00202 }
00203 
00204 template<typename MatrixType>
00205 inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT)
00206 {
00207 #if   LDBL_MANT_DIG == 53
00208   const int maxPadeDegree = 7;
00209   const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
00210       1.999045567181744e-1L, 2.789358995219730e-1L };
00211 #elif LDBL_MANT_DIG <= 64
00212   const int maxPadeDegree = 8;
00213   const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
00214       6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
00215 #elif LDBL_MANT_DIG <= 106
00216   const int maxPadeDegree = 10;
00217   const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
00218       1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
00219       2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
00220       1.1016843812851143391275867258512e-1L };
00221 #else
00222   const int maxPadeDegree = 10;
00223   const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
00224       6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
00225       9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
00226       3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
00227       9.134603732914548552537150753385375e-2L };
00228 #endif
00229   int degree = 3;
00230   for (; degree <= maxPadeDegree; ++degree)
00231     if (normIminusT <= maxNormForPade[degree - 3])
00232       break;
00233   return degree;
00234 }
00235 
00236 template<typename MatrixType>
00237 inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
00238 MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p)
00239 {
00240   ComplexScalar logCurr = std::log(curr);
00241   ComplexScalar logPrev = std::log(prev);
00242   int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI));
00243   ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber);
00244   return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev);
00245 }
00246 
00247 template<typename MatrixType>
00248 inline typename MatrixPowerAtomic<MatrixType>::RealScalar
00249 MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
00250 {
00251   RealScalar w = numext::atanh2(curr - prev, curr + prev);
00252   return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev);
00253 }
00254 
00274 template<typename MatrixType>
00275 class MatrixPower
00276 {
00277   private:
00278     enum {
00279       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00280       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00281       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00282       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00283     };
00284     typedef typename MatrixType::Scalar Scalar;
00285     typedef typename MatrixType::RealScalar RealScalar;
00286     typedef typename MatrixType::Index Index;
00287 
00288   public:
00297     explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0)
00298     { eigen_assert(A.rows() == A.cols()); }
00299 
00307     const MatrixPowerRetval<MatrixType> operator()(RealScalar p)
00308     { return MatrixPowerRetval<MatrixType>(*this, p); }
00309 
00317     template<typename ResultType>
00318     void compute(ResultType& res, RealScalar p);
00319     
00320     Index rows() const { return m_A.rows(); }
00321     Index cols() const { return m_A.cols(); }
00322 
00323   private:
00324     typedef std::complex<RealScalar> ComplexScalar;
00325     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, MatrixType::Options,
00326               MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrix;
00327 
00328     typename MatrixType::Nested m_A;
00329     MatrixType m_tmp;
00330     ComplexMatrix m_T, m_U, m_fT;
00331     RealScalar m_conditionNumber;
00332 
00333     RealScalar modfAndInit(RealScalar, RealScalar*);
00334 
00335     template<typename ResultType>
00336     void computeIntPower(ResultType&, RealScalar);
00337 
00338     template<typename ResultType>
00339     void computeFracPower(ResultType&, RealScalar);
00340 
00341     template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
00342     static void revertSchur(
00343         Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
00344         const ComplexMatrix& T,
00345         const ComplexMatrix& U);
00346 
00347     template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
00348     static void revertSchur(
00349         Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
00350         const ComplexMatrix& T,
00351         const ComplexMatrix& U);
00352 };
00353 
00354 template<typename MatrixType>
00355 template<typename ResultType>
00356 void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p)
00357 {
00358   switch (cols()) {
00359     case 0:
00360       break;
00361     case 1:
00362       res(0,0) = std::pow(m_A.coeff(0,0), p);
00363       break;
00364     default:
00365       RealScalar intpart, x = modfAndInit(p, &intpart);
00366       computeIntPower(res, intpart);
00367       computeFracPower(res, x);
00368   }
00369 }
00370 
00371 template<typename MatrixType>
00372 typename MatrixPower<MatrixType>::RealScalar
00373 MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
00374 {
00375   typedef Array<RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime> RealArray;
00376 
00377   *intpart = std::floor(x);
00378   RealScalar res = x - *intpart;
00379 
00380   if (!m_conditionNumber && res) {
00381     const ComplexSchur<MatrixType> schurOfA(m_A);
00382     m_T = schurOfA.matrixT();
00383     m_U = schurOfA.matrixU();
00384     
00385     const RealArray absTdiag = m_T.diagonal().array().abs();
00386     m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
00387   }
00388 
00389   if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
00390     --res;
00391     ++*intpart;
00392   }
00393   return res;
00394 }
00395 
00396 template<typename MatrixType>
00397 template<typename ResultType>
00398 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
00399 {
00400   RealScalar pp = std::abs(p);
00401 
00402   if (p<0)  m_tmp = m_A.inverse();
00403   else      m_tmp = m_A;
00404 
00405   res = MatrixType::Identity(rows(), cols());
00406   while (pp >= 1) {
00407     if (std::fmod(pp, 2) >= 1)
00408       res = m_tmp * res;
00409     m_tmp *= m_tmp;
00410     pp /= 2;
00411   }
00412 }
00413 
00414 template<typename MatrixType>
00415 template<typename ResultType>
00416 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
00417 {
00418   if (p) {
00419     eigen_assert(m_conditionNumber);
00420     MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
00421     revertSchur(m_tmp, m_fT, m_U);
00422     res = m_tmp * res;
00423   }
00424 }
00425 
00426 template<typename MatrixType>
00427 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
00428 inline void MatrixPower<MatrixType>::revertSchur(
00429     Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
00430     const ComplexMatrix& T,
00431     const ComplexMatrix& U)
00432 { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
00433 
00434 template<typename MatrixType>
00435 template<int Rows, int Cols, int Options, int MaxRows, int MaxCols>
00436 inline void MatrixPower<MatrixType>::revertSchur(
00437     Matrix<RealScalar, Rows, Cols, Options, MaxRows, MaxCols>& res,
00438     const ComplexMatrix& T,
00439     const ComplexMatrix& U)
00440 { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
00441 
00455 template<typename Derived>
00456 class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> >
00457 {
00458   public:
00459     typedef typename Derived::PlainObject PlainObject;
00460     typedef typename Derived::RealScalar RealScalar;
00461     typedef typename Derived::Index Index;
00462 
00469     MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p)
00470     { }
00471 
00478     template<typename ResultType>
00479     inline void evalTo(ResultType& res) const
00480     { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); }
00481 
00482     Index rows() const { return m_A.rows(); }
00483     Index cols() const { return m_A.cols(); }
00484 
00485   private:
00486     const Derived& m_A;
00487     const RealScalar m_p;
00488     MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
00489 };
00490 
00491 namespace internal {
00492 
00493 template<typename MatrixPowerType>
00494 struct traits< MatrixPowerRetval<MatrixPowerType> >
00495 { typedef typename MatrixPowerType::PlainObject ReturnType; };
00496 
00497 template<typename Derived>
00498 struct traits< MatrixPowerReturnValue<Derived> >
00499 { typedef typename Derived::PlainObject ReturnType; };
00500 
00501 }
00502 
00503 template<typename Derived>
00504 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const
00505 { return MatrixPowerReturnValue<Derived>(derived(), p); }
00506 
00507 } // namespace Eigen
00508 
00509 #endif // EIGEN_MATRIX_POWER


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:11