1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009, 2010 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 5 // Copyright (C) 2011 Chen-Pang He <jdh8 (at) ms63.hinet.net> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_MATRIX_EXPONENTIAL 12 #define EIGEN_MATRIX_EXPONENTIAL 13 14 #include "StemFunction.h" 15 16 namespace Eigen { 17 18 #if defined(_MSC_VER) || defined(__FreeBSD__) 19 template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); } 20 #endif 21 22 23 /** \ingroup MatrixFunctions_Module 24 * \brief Class for computing the matrix exponential. 25 * \tparam MatrixType type of the argument of the exponential, 26 * expected to be an instantiation of the Matrix class template. 27 */ 28 template <typename MatrixType> 29 class MatrixExponential { 30 31 public: 32 33 /** \brief Constructor. 34 * 35 * The class stores a reference to \p M, so it should not be 36 * changed (or destroyed) before compute() is called. 37 * 38 * \param[in] M matrix whose exponential is to be computed. 39 */ 40 MatrixExponential(const MatrixType &M); 41 42 /** \brief Computes the matrix exponential. 43 * 44 * \param[out] result the matrix exponential of \p M in the constructor. 45 */ 46 template <typename ResultType> 47 void compute(ResultType &result); 48 49 private: 50 51 // Prevent copying 52 MatrixExponential(const MatrixExponential&); 53 MatrixExponential& operator=(const MatrixExponential&); 54 55 /** \brief Compute the (3,3)-Padé approximant to the exponential. 56 * 57 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 58 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 59 * 60 * \param[in] A Argument of matrix exponential 61 */ 62 void pade3(const MatrixType &A); 63 64 /** \brief Compute the (5,5)-Padé approximant to the exponential. 65 * 66 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 67 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 68 * 69 * \param[in] A Argument of matrix exponential 70 */ 71 void pade5(const MatrixType &A); 72 73 /** \brief Compute the (7,7)-Padé approximant to the exponential. 74 * 75 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 76 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 77 * 78 * \param[in] A Argument of matrix exponential 79 */ 80 void pade7(const MatrixType &A); 81 82 /** \brief Compute the (9,9)-Padé approximant to the exponential. 83 * 84 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 85 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 86 * 87 * \param[in] A Argument of matrix exponential 88 */ 89 void pade9(const MatrixType &A); 90 91 /** \brief Compute the (13,13)-Padé approximant to the exponential. 92 * 93 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 94 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 95 * 96 * \param[in] A Argument of matrix exponential 97 */ 98 void pade13(const MatrixType &A); 99 100 /** \brief Compute the (17,17)-Padé approximant to the exponential. 101 * 102 * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé 103 * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. 104 * 105 * This function activates only if your long double is double-double or quadruple. 106 * 107 * \param[in] A Argument of matrix exponential 108 */ 109 void pade17(const MatrixType &A); 110 111 /** \brief Compute Padé approximant to the exponential. 112 * 113 * Computes \c m_U, \c m_V and \c m_squarings such that 114 * \f$ (V+U)(V-U)^{-1} \f$ is a Padé of 115 * \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$. The 116 * degree of the Padé approximant and the value of 117 * squarings are chosen such that the approximation error is no 118 * more than the round-off error. 119 * 120 * The argument of this function should correspond with the (real 121 * part of) the entries of \c m_M. It is used to select the 122 * correct implementation using overloading. 123 */ 124 void computeUV(double); 125 126 /** \brief Compute Padé approximant to the exponential. 127 * 128 * \sa computeUV(double); 129 */ 130 void computeUV(float); 131 132 /** \brief Compute Padé approximant to the exponential. 133 * 134 * \sa computeUV(double); 135 */ 136 void computeUV(long double); 137 138 typedef typename internal::traits<MatrixType>::Scalar Scalar; 139 typedef typename NumTraits<Scalar>::Real RealScalar; 140 typedef typename std::complex<RealScalar> ComplexScalar; 141 142 /** \brief Reference to matrix whose exponential is to be computed. */ 143 typename internal::nested<MatrixType>::type m_M; 144 145 /** \brief Odd-degree terms in numerator of Padé approximant. */ 146 MatrixType m_U; 147 148 /** \brief Even-degree terms in numerator of Padé approximant. */ 149 MatrixType m_V; 150 151 /** \brief Used for temporary storage. */ 152 MatrixType m_tmp1; 153 154 /** \brief Used for temporary storage. */ 155 MatrixType m_tmp2; 156 157 /** \brief Identity matrix of the same size as \c m_M. */ 158 MatrixType m_Id; 159 160 /** \brief Number of squarings required in the last step. */ 161 int m_squarings; 162 163 /** \brief L1 norm of m_M. */ 164 RealScalar m_l1norm; 165 }; 166 167 template <typename MatrixType> 168 MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M) : 169 m_M(M), 170 m_U(M.rows(),M.cols()), 171 m_V(M.rows(),M.cols()), 172 m_tmp1(M.rows(),M.cols()), 173 m_tmp2(M.rows(),M.cols()), 174 m_Id(MatrixType::Identity(M.rows(), M.cols())), 175 m_squarings(0), 176 m_l1norm(M.cwiseAbs().colwise().sum().maxCoeff()) 177 { 178 /* empty body */ 179 } 180 181 template <typename MatrixType> 182 template <typename ResultType> 183 void MatrixExponential<MatrixType>::compute(ResultType &result) 184 { 185 #if LDBL_MANT_DIG > 112 // rarely happens 186 if(sizeof(RealScalar) > 14) { 187 result = m_M.matrixFunction(StdStemFunctions<ComplexScalar>::exp); 188 return; 189 } 190 #endif 191 computeUV(RealScalar()); 192 m_tmp1 = m_U + m_V; // numerator of Pade approximant 193 m_tmp2 = -m_U + m_V; // denominator of Pade approximant 194 result = m_tmp2.partialPivLu().solve(m_tmp1); 195 for (int i=0; i<m_squarings; i++) 196 result *= result; // undo scaling by repeated squaring 197 } 198 199 template <typename MatrixType> 200 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A) 201 { 202 const RealScalar b[] = {120., 60., 12., 1.}; 203 m_tmp1.noalias() = A * A; 204 m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id; 205 m_U.noalias() = A * m_tmp2; 206 m_V = b[2]*m_tmp1 + b[0]*m_Id; 207 } 208 209 template <typename MatrixType> 210 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A) 211 { 212 const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.}; 213 MatrixType A2 = A * A; 214 m_tmp1.noalias() = A2 * A2; 215 m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id; 216 m_U.noalias() = A * m_tmp2; 217 m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id; 218 } 219 220 template <typename MatrixType> 221 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A) 222 { 223 const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; 224 MatrixType A2 = A * A; 225 MatrixType A4 = A2 * A2; 226 m_tmp1.noalias() = A4 * A2; 227 m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 228 m_U.noalias() = A * m_tmp2; 229 m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 230 } 231 232 template <typename MatrixType> 233 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A) 234 { 235 const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., 236 2162160., 110880., 3960., 90., 1.}; 237 MatrixType A2 = A * A; 238 MatrixType A4 = A2 * A2; 239 MatrixType A6 = A4 * A2; 240 m_tmp1.noalias() = A6 * A2; 241 m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 242 m_U.noalias() = A * m_tmp2; 243 m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 244 } 245 246 template <typename MatrixType> 247 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A) 248 { 249 const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., 250 1187353796428800., 129060195264000., 10559470521600., 670442572800., 251 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; 252 MatrixType A2 = A * A; 253 MatrixType A4 = A2 * A2; 254 m_tmp1.noalias() = A4 * A2; 255 m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage 256 m_tmp2.noalias() = m_tmp1 * m_V; 257 m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 258 m_U.noalias() = A * m_tmp2; 259 m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2; 260 m_V.noalias() = m_tmp1 * m_tmp2; 261 m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 262 } 263 264 #if LDBL_MANT_DIG > 64 265 template <typename MatrixType> 266 EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade17(const MatrixType &A) 267 { 268 const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L, 269 100610229646136770560000.L, 15720348382208870400000.L, 270 1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L, 271 595373117923584000.L, 27563570274240000.L, 1060137318240000.L, 272 33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L, 273 46512.L, 306.L, 1.L}; 274 MatrixType A2 = A * A; 275 MatrixType A4 = A2 * A2; 276 MatrixType A6 = A4 * A2; 277 m_tmp1.noalias() = A4 * A4; 278 m_V = b[17]*m_tmp1 + b[15]*A6 + b[13]*A4 + b[11]*A2; // used for temporary storage 279 m_tmp2.noalias() = m_tmp1 * m_V; 280 m_tmp2 += b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; 281 m_U.noalias() = A * m_tmp2; 282 m_tmp2 = b[16]*m_tmp1 + b[14]*A6 + b[12]*A4 + b[10]*A2; 283 m_V.noalias() = m_tmp1 * m_tmp2; 284 m_V += b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; 285 } 286 #endif 287 288 template <typename MatrixType> 289 void MatrixExponential<MatrixType>::computeUV(float) 290 { 291 using std::max; 292 using std::pow; 293 using std::ceil; 294 if (m_l1norm < 4.258730016922831e-001) { 295 pade3(m_M); 296 } else if (m_l1norm < 1.880152677804762e+000) { 297 pade5(m_M); 298 } else { 299 const float maxnorm = 3.925724783138660f; 300 m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); 301 MatrixType A = m_M / pow(Scalar(2), m_squarings); 302 pade7(A); 303 } 304 } 305 306 template <typename MatrixType> 307 void MatrixExponential<MatrixType>::computeUV(double) 308 { 309 using std::max; 310 using std::pow; 311 using std::ceil; 312 if (m_l1norm < 1.495585217958292e-002) { 313 pade3(m_M); 314 } else if (m_l1norm < 2.539398330063230e-001) { 315 pade5(m_M); 316 } else if (m_l1norm < 9.504178996162932e-001) { 317 pade7(m_M); 318 } else if (m_l1norm < 2.097847961257068e+000) { 319 pade9(m_M); 320 } else { 321 const double maxnorm = 5.371920351148152; 322 m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); 323 MatrixType A = m_M / pow(Scalar(2), m_squarings); 324 pade13(A); 325 } 326 } 327 328 template <typename MatrixType> 329 void MatrixExponential<MatrixType>::computeUV(long double) 330 { 331 using std::max; 332 using std::pow; 333 using std::ceil; 334 #if LDBL_MANT_DIG == 53 // double precision 335 computeUV(double()); 336 #elif LDBL_MANT_DIG <= 64 // extended precision 337 if (m_l1norm < 4.1968497232266989671e-003L) { 338 pade3(m_M); 339 } else if (m_l1norm < 1.1848116734693823091e-001L) { 340 pade5(m_M); 341 } else if (m_l1norm < 5.5170388480686700274e-001L) { 342 pade7(m_M); 343 } else if (m_l1norm < 1.3759868875587845383e+000L) { 344 pade9(m_M); 345 } else { 346 const long double maxnorm = 4.0246098906697353063L; 347 m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); 348 MatrixType A = m_M / pow(Scalar(2), m_squarings); 349 pade13(A); 350 } 351 #elif LDBL_MANT_DIG <= 106 // double-double 352 if (m_l1norm < 3.2787892205607026992947488108213e-005L) { 353 pade3(m_M); 354 } else if (m_l1norm < 6.4467025060072760084130906076332e-003L) { 355 pade5(m_M); 356 } else if (m_l1norm < 6.8988028496595374751374122881143e-002L) { 357 pade7(m_M); 358 } else if (m_l1norm < 2.7339737518502231741495857201670e-001L) { 359 pade9(m_M); 360 } else if (m_l1norm < 1.3203382096514474905666448850278e+000L) { 361 pade13(m_M); 362 } else { 363 const long double maxnorm = 3.2579440895405400856599663723517L; 364 m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); 365 MatrixType A = m_M / pow(Scalar(2), m_squarings); 366 pade17(A); 367 } 368 #elif LDBL_MANT_DIG <= 112 // quadruple precison 369 if (m_l1norm < 1.639394610288918690547467954466970e-005L) { 370 pade3(m_M); 371 } else if (m_l1norm < 4.253237712165275566025884344433009e-003L) { 372 pade5(m_M); 373 } else if (m_l1norm < 5.125804063165764409885122032933142e-002L) { 374 pade7(m_M); 375 } else if (m_l1norm < 2.170000765161155195453205651889853e-001L) { 376 pade9(m_M); 377 } else if (m_l1norm < 1.125358383453143065081397882891878e+000L) { 378 pade13(m_M); 379 } else { 380 const long double maxnorm = 2.884233277829519311757165057717815L; 381 m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); 382 MatrixType A = m_M / pow(Scalar(2), m_squarings); 383 pade17(A); 384 } 385 #else 386 // this case should be handled in compute() 387 eigen_assert(false && "Bug in MatrixExponential"); 388 #endif // LDBL_MANT_DIG 389 } 390 391 /** \ingroup MatrixFunctions_Module 392 * 393 * \brief Proxy for the matrix exponential of some matrix (expression). 394 * 395 * \tparam Derived Type of the argument to the matrix exponential. 396 * 397 * This class holds the argument to the matrix exponential until it 398 * is assigned or evaluated for some other reason (so the argument 399 * should not be changed in the meantime). It is the return type of 400 * MatrixBase::exp() and most of the time this is the only way it is 401 * used. 402 */ 403 template<typename Derived> struct MatrixExponentialReturnValue 404 : public ReturnByValue<MatrixExponentialReturnValue<Derived> > 405 { 406 typedef typename Derived::Index Index; 407 public: 408 /** \brief Constructor. 409 * 410 * \param[in] src %Matrix (expression) forming the argument of the 411 * matrix exponential. 412 */ 413 MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } 414 415 /** \brief Compute the matrix exponential. 416 * 417 * \param[out] result the matrix exponential of \p src in the 418 * constructor. 419 */ 420 template <typename ResultType> 421 inline void evalTo(ResultType& result) const 422 { 423 const typename Derived::PlainObject srcEvaluated = m_src.eval(); 424 MatrixExponential<typename Derived::PlainObject> me(srcEvaluated); 425 me.compute(result); 426 } 427 428 Index rows() const { return m_src.rows(); } 429 Index cols() const { return m_src.cols(); } 430 431 protected: 432 const Derived& m_src; 433 private: 434 MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&); 435 }; 436 437 namespace internal { 438 template<typename Derived> 439 struct traits<MatrixExponentialReturnValue<Derived> > 440 { 441 typedef typename Derived::PlainObject ReturnType; 442 }; 443 } 444 445 template <typename Derived> 446 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const 447 { 448 eigen_assert(rows() == cols()); 449 return MatrixExponentialReturnValue<Derived>(derived()); 450 } 451 452 } // end namespace Eigen 453 454 #endif // EIGEN_MATRIX_EXPONENTIAL 455
