MatrixExponential.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
00005 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
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_MATRIX_EXPONENTIAL
00012 #define EIGEN_MATRIX_EXPONENTIAL
00013 
00014 #include "StemFunction.h"
00015 
00016 namespace Eigen {
00017 namespace internal {
00018 
00023 template <typename RealScalar>
00024 struct MatrixExponentialScalingOp
00025 {
00030   MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
00031 
00032 
00037   inline const RealScalar operator() (const RealScalar& x) const
00038   {
00039     using std::ldexp;
00040     return ldexp(x, -m_squarings);
00041   }
00042 
00043   typedef std::complex<RealScalar> ComplexScalar;
00044 
00049   inline const ComplexScalar operator() (const ComplexScalar& x) const
00050   {
00051     using std::ldexp;
00052     return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
00053   }
00054 
00055   private:
00056     int m_squarings;
00057 };
00058 
00064 template <typename MatA, typename MatU, typename MatV>
00065 void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V)
00066 {
00067   typedef typename MatA::PlainObject MatrixType;
00068   typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
00069   const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
00070   const MatrixType A2 = A * A;
00071   const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
00072   U.noalias() = A * tmp;
00073   V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
00074 }
00075 
00081 template <typename MatA, typename MatU, typename MatV>
00082 void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V)
00083 {
00084   typedef typename MatA::PlainObject MatrixType;
00085   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
00086   const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
00087   const MatrixType A2 = A * A;
00088   const MatrixType A4 = A2 * A2;
00089   const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
00090   U.noalias() = A * tmp;
00091   V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
00092 }
00093 
00099 template <typename MatA, typename MatU, typename MatV>
00100 void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V)
00101 {
00102   typedef typename MatA::PlainObject MatrixType;
00103   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
00104   const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
00105   const MatrixType A2 = A * A;
00106   const MatrixType A4 = A2 * A2;
00107   const MatrixType A6 = A4 * A2;
00108   const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 
00109     + b[1] * MatrixType::Identity(A.rows(), A.cols());
00110   U.noalias() = A * tmp;
00111   V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
00112 
00113 }
00114 
00120 template <typename MatA, typename MatU, typename MatV>
00121 void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V)
00122 {
00123   typedef typename MatA::PlainObject MatrixType;
00124   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
00125   const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
00126                           2162160.L, 110880.L, 3960.L, 90.L, 1.L};
00127   const MatrixType A2 = A * A;
00128   const MatrixType A4 = A2 * A2;
00129   const MatrixType A6 = A4 * A2;
00130   const MatrixType A8 = A6 * A2;
00131   const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 
00132     + b[1] * MatrixType::Identity(A.rows(), A.cols());
00133   U.noalias() = A * tmp;
00134   V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
00135 }
00136 
00142 template <typename MatA, typename MatU, typename MatV>
00143 void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V)
00144 {
00145   typedef typename MatA::PlainObject MatrixType;
00146   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
00147   const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
00148                           1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
00149                           33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
00150   const MatrixType A2 = A * A;
00151   const MatrixType A4 = A2 * A2;
00152   const MatrixType A6 = A4 * A2;
00153   V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
00154   MatrixType tmp = A6 * V;
00155   tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
00156   U.noalias() = A * tmp;
00157   tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
00158   V.noalias() = A6 * tmp;
00159   V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
00160 }
00161 
00169 #if LDBL_MANT_DIG > 64
00170 template <typename MatA, typename MatU, typename MatV>
00171 void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V)
00172 {
00173   typedef typename MatA::PlainObject MatrixType;
00174   typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
00175   const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
00176                           100610229646136770560000.L, 15720348382208870400000.L,
00177                           1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
00178                           595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
00179                           33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
00180                           46512.L, 306.L, 1.L};
00181   const MatrixType A2 = A * A;
00182   const MatrixType A4 = A2 * A2;
00183   const MatrixType A6 = A4 * A2;
00184   const MatrixType A8 = A4 * A4;
00185   V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
00186   MatrixType tmp = A8 * V;
00187   tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 
00188     + b[1] * MatrixType::Identity(A.rows(), A.cols());
00189   U.noalias() = A * tmp;
00190   tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
00191   V.noalias() = tmp * A8;
00192   V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 
00193     + b[0] * MatrixType::Identity(A.rows(), A.cols());
00194 }
00195 #endif
00196 
00197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
00198 struct matrix_exp_computeUV
00199 {
00207   static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
00208 };
00209 
00210 template <typename MatrixType>
00211 struct matrix_exp_computeUV<MatrixType, float>
00212 {
00213   template <typename ArgType>
00214   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
00215   {
00216     using std::frexp;
00217     using std::pow;
00218     const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
00219     squarings = 0;
00220     if (l1norm < 4.258730016922831e-001f) {
00221       matrix_exp_pade3(arg, U, V);
00222     } else if (l1norm < 1.880152677804762e+000f) {
00223       matrix_exp_pade5(arg, U, V);
00224     } else {
00225       const float maxnorm = 3.925724783138660f;
00226       frexp(l1norm / maxnorm, &squarings);
00227       if (squarings < 0) squarings = 0;
00228       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
00229       matrix_exp_pade7(A, U, V);
00230     }
00231   }
00232 };
00233 
00234 template <typename MatrixType>
00235 struct matrix_exp_computeUV<MatrixType, double>
00236 {
00237   template <typename ArgType>
00238   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
00239   {
00240     using std::frexp;
00241     using std::pow;
00242     const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
00243     squarings = 0;
00244     if (l1norm < 1.495585217958292e-002) {
00245       matrix_exp_pade3(arg, U, V);
00246     } else if (l1norm < 2.539398330063230e-001) {
00247       matrix_exp_pade5(arg, U, V);
00248     } else if (l1norm < 9.504178996162932e-001) {
00249       matrix_exp_pade7(arg, U, V);
00250     } else if (l1norm < 2.097847961257068e+000) {
00251       matrix_exp_pade9(arg, U, V);
00252     } else {
00253       const double maxnorm = 5.371920351148152;
00254       frexp(l1norm / maxnorm, &squarings);
00255       if (squarings < 0) squarings = 0;
00256       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<double>(squarings));
00257       matrix_exp_pade13(A, U, V);
00258     }
00259   }
00260 };
00261   
00262 template <typename MatrixType>
00263 struct matrix_exp_computeUV<MatrixType, long double>
00264 {
00265   template <typename ArgType>
00266   static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
00267   {
00268 #if   LDBL_MANT_DIG == 53   // double precision
00269     matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
00270   
00271 #else
00272   
00273     using std::frexp;
00274     using std::pow;
00275     const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
00276     squarings = 0;
00277   
00278 #if LDBL_MANT_DIG <= 64   // extended precision
00279   
00280     if (l1norm < 4.1968497232266989671e-003L) {
00281       matrix_exp_pade3(arg, U, V);
00282     } else if (l1norm < 1.1848116734693823091e-001L) {
00283       matrix_exp_pade5(arg, U, V);
00284     } else if (l1norm < 5.5170388480686700274e-001L) {
00285       matrix_exp_pade7(arg, U, V);
00286     } else if (l1norm < 1.3759868875587845383e+000L) {
00287       matrix_exp_pade9(arg, U, V);
00288     } else {
00289       const long double maxnorm = 4.0246098906697353063L;
00290       frexp(l1norm / maxnorm, &squarings);
00291       if (squarings < 0) squarings = 0;
00292       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
00293       matrix_exp_pade13(A, U, V);
00294     }
00295   
00296 #elif LDBL_MANT_DIG <= 106  // double-double
00297   
00298     if (l1norm < 3.2787892205607026992947488108213e-005L) {
00299       matrix_exp_pade3(arg, U, V);
00300     } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
00301       matrix_exp_pade5(arg, U, V);
00302     } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
00303       matrix_exp_pade7(arg, U, V);
00304     } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
00305       matrix_exp_pade9(arg, U, V);
00306     } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
00307       matrix_exp_pade13(arg, U, V);
00308     } else {
00309       const long double maxnorm = 3.2579440895405400856599663723517L;
00310       frexp(l1norm / maxnorm, &squarings);
00311       if (squarings < 0) squarings = 0;
00312       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
00313       matrix_exp_pade17(A, U, V);
00314     }
00315   
00316 #elif LDBL_MANT_DIG <= 112  // quadruple precison
00317   
00318     if (l1norm < 1.639394610288918690547467954466970e-005L) {
00319       matrix_exp_pade3(arg, U, V);
00320     } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
00321       matrix_exp_pade5(arg, U, V);
00322     } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
00323       matrix_exp_pade7(arg, U, V);
00324     } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
00325       matrix_exp_pade9(arg, U, V);
00326     } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
00327       matrix_exp_pade13(arg, U, V);
00328     } else {
00329       frexp(l1norm / maxnorm, &squarings);
00330       if (squarings < 0) squarings = 0;
00331       MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
00332       matrix_exp_pade17(A, U, V);
00333     }
00334   
00335 #else
00336   
00337     // this case should be handled in compute()
00338     eigen_assert(false && "Bug in MatrixExponential"); 
00339   
00340 #endif
00341 #endif  // LDBL_MANT_DIG
00342   }
00343 };
00344 
00345 
00346 /* Computes the matrix exponential
00347  *
00348  * \param arg    argument of matrix exponential (should be plain object)
00349  * \param result variable in which result will be stored
00350  */
00351 template <typename ArgType, typename ResultType>
00352 void matrix_exp_compute(const ArgType& arg, ResultType &result)
00353 {
00354   typedef typename ArgType::PlainObject MatrixType;
00355 #if LDBL_MANT_DIG > 112 // rarely happens
00356   typedef typename traits<MatrixType>::Scalar Scalar;
00357   typedef typename NumTraits<Scalar>::Real RealScalar;
00358   typedef typename std::complex<RealScalar> ComplexScalar;
00359   if (sizeof(RealScalar) > 14) {
00360     result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
00361     return;
00362   }
00363 #endif
00364   MatrixType U, V;
00365   int squarings; 
00366   matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
00367   MatrixType numer = U + V;
00368   MatrixType denom = -U + V;
00369   result = denom.partialPivLu().solve(numer);
00370   for (int i=0; i<squarings; i++)
00371     result *= result;   // undo scaling by repeated squaring
00372 }
00373 
00374 } // end namespace Eigen::internal
00375 
00386 template<typename Derived> struct MatrixExponentialReturnValue
00387 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
00388 {
00389     typedef typename Derived::Index Index;
00390   public:
00395     MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
00396 
00401     template <typename ResultType>
00402     inline void evalTo(ResultType& result) const
00403     {
00404       const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
00405       internal::matrix_exp_compute(tmp, result);
00406     }
00407 
00408     Index rows() const { return m_src.rows(); }
00409     Index cols() const { return m_src.cols(); }
00410 
00411   protected:
00412     const typename internal::ref_selector<Derived>::type m_src;
00413 };
00414 
00415 namespace internal {
00416 template<typename Derived>
00417 struct traits<MatrixExponentialReturnValue<Derived> >
00418 {
00419   typedef typename Derived::PlainObject ReturnType;
00420 };
00421 }
00422 
00423 template <typename Derived>
00424 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
00425 {
00426   eigen_assert(rows() == cols());
00427   return MatrixExponentialReturnValue<Derived>(derived());
00428 }
00429 
00430 } // end namespace Eigen
00431 
00432 #endif // EIGEN_MATRIX_EXPONENTIAL
 All Classes Functions Variables Typedefs Enumerator