Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
MatrixExponential.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009, 2010, 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
5 // Copyright (C) 2011, 2013 Chen-Pang He <jdh8@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 namespace internal {
18 
23 template <typename RealScalar>
24 struct MatrixExponentialScalingOp
25 {
30  MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
31 
32 
37  inline const RealScalar operator() (const RealScalar& x) const
38  {
39  using std::ldexp;
40  return ldexp(x, -m_squarings);
41  }
42 
43  typedef std::complex<RealScalar> ComplexScalar;
44 
49  inline const ComplexScalar operator() (const ComplexScalar& x) const
50  {
51  using std::ldexp;
52  return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
53  }
54 
55  private:
56  int m_squarings;
57 };
58 
64 template <typename MatA, typename MatU, typename MatV>
65 void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V)
66 {
67  typedef typename MatA::PlainObject MatrixType;
68  typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
69  const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
70  const MatrixType A2 = A * A;
71  const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
72  U.noalias() = A * tmp;
73  V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
74 }
75 
81 template <typename MatA, typename MatU, typename MatV>
82 void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V)
83 {
84  typedef typename MatA::PlainObject MatrixType;
85  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
86  const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
87  const MatrixType A2 = A * A;
88  const MatrixType A4 = A2 * A2;
89  const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
90  U.noalias() = A * tmp;
91  V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
92 }
93 
99 template <typename MatA, typename MatU, typename MatV>
100 void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V)
101 {
102  typedef typename MatA::PlainObject MatrixType;
103  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
104  const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
105  const MatrixType A2 = A * A;
106  const MatrixType A4 = A2 * A2;
107  const MatrixType A6 = A4 * A2;
108  const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2
109  + b[1] * MatrixType::Identity(A.rows(), A.cols());
110  U.noalias() = A * tmp;
111  V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
112 
113 }
114 
120 template <typename MatA, typename MatU, typename MatV>
121 void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V)
122 {
123  typedef typename MatA::PlainObject MatrixType;
124  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
125  const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
126  2162160.L, 110880.L, 3960.L, 90.L, 1.L};
127  const MatrixType A2 = A * A;
128  const MatrixType A4 = A2 * A2;
129  const MatrixType A6 = A4 * A2;
130  const MatrixType A8 = A6 * A2;
131  const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
132  + b[1] * MatrixType::Identity(A.rows(), A.cols());
133  U.noalias() = A * tmp;
134  V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
135 }
136 
142 template <typename MatA, typename MatU, typename MatV>
143 void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V)
144 {
145  typedef typename MatA::PlainObject MatrixType;
146  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
147  const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
148  1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
149  33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
150  const MatrixType A2 = A * A;
151  const MatrixType A4 = A2 * A2;
152  const MatrixType A6 = A4 * A2;
153  V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
154  MatrixType tmp = A6 * V;
155  tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
156  U.noalias() = A * tmp;
157  tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
158  V.noalias() = A6 * tmp;
159  V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
160 }
161 
169 #if LDBL_MANT_DIG > 64
170 template <typename MatA, typename MatU, typename MatV>
171 void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V)
172 {
173  typedef typename MatA::PlainObject MatrixType;
174  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
175  const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
176  100610229646136770560000.L, 15720348382208870400000.L,
177  1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
178  595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
179  33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
180  46512.L, 306.L, 1.L};
181  const MatrixType A2 = A * A;
182  const MatrixType A4 = A2 * A2;
183  const MatrixType A6 = A4 * A2;
184  const MatrixType A8 = A4 * A4;
185  V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
186  MatrixType tmp = A8 * V;
187  tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
188  + b[1] * MatrixType::Identity(A.rows(), A.cols());
189  U.noalias() = A * tmp;
190  tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
191  V.noalias() = tmp * A8;
192  V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2
193  + b[0] * MatrixType::Identity(A.rows(), A.cols());
194 }
195 #endif
196 
197 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
198 struct matrix_exp_computeUV
199 {
207  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
208 };
209 
210 template <typename MatrixType>
211 struct matrix_exp_computeUV<MatrixType, float>
212 {
213  template <typename ArgType>
214  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
215  {
216  using std::frexp;
217  using std::pow;
218  const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
219  squarings = 0;
220  if (l1norm < 4.258730016922831e-001f) {
221  matrix_exp_pade3(arg, U, V);
222  } else if (l1norm < 1.880152677804762e+000f) {
223  matrix_exp_pade5(arg, U, V);
224  } else {
225  const float maxnorm = 3.925724783138660f;
226  frexp(l1norm / maxnorm, &squarings);
227  if (squarings < 0) squarings = 0;
228  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
229  matrix_exp_pade7(A, U, V);
230  }
231  }
232 };
233 
234 template <typename MatrixType>
235 struct matrix_exp_computeUV<MatrixType, double>
236 {
237  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
238  template <typename ArgType>
239  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
240  {
241  using std::frexp;
242  using std::pow;
243  const RealScalar l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
244  squarings = 0;
245  if (l1norm < 1.495585217958292e-002) {
246  matrix_exp_pade3(arg, U, V);
247  } else if (l1norm < 2.539398330063230e-001) {
248  matrix_exp_pade5(arg, U, V);
249  } else if (l1norm < 9.504178996162932e-001) {
250  matrix_exp_pade7(arg, U, V);
251  } else if (l1norm < 2.097847961257068e+000) {
252  matrix_exp_pade9(arg, U, V);
253  } else {
254  const RealScalar maxnorm = 5.371920351148152;
255  frexp(l1norm / maxnorm, &squarings);
256  if (squarings < 0) squarings = 0;
257  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<RealScalar>(squarings));
258  matrix_exp_pade13(A, U, V);
259  }
260  }
261 };
262 
263 template <typename MatrixType>
264 struct matrix_exp_computeUV<MatrixType, long double>
265 {
266  template <typename ArgType>
267  static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings)
268  {
269 #if LDBL_MANT_DIG == 53 // double precision
270  matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
271 
272 #else
273 
274  using std::frexp;
275  using std::pow;
276  const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
277  squarings = 0;
278 
279 #if LDBL_MANT_DIG <= 64 // extended precision
280 
281  if (l1norm < 4.1968497232266989671e-003L) {
282  matrix_exp_pade3(arg, U, V);
283  } else if (l1norm < 1.1848116734693823091e-001L) {
284  matrix_exp_pade5(arg, U, V);
285  } else if (l1norm < 5.5170388480686700274e-001L) {
286  matrix_exp_pade7(arg, U, V);
287  } else if (l1norm < 1.3759868875587845383e+000L) {
288  matrix_exp_pade9(arg, U, V);
289  } else {
290  const long double maxnorm = 4.0246098906697353063L;
291  frexp(l1norm / maxnorm, &squarings);
292  if (squarings < 0) squarings = 0;
293  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
294  matrix_exp_pade13(A, U, V);
295  }
296 
297 #elif LDBL_MANT_DIG <= 106 // double-double
298 
299  if (l1norm < 3.2787892205607026992947488108213e-005L) {
300  matrix_exp_pade3(arg, U, V);
301  } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
302  matrix_exp_pade5(arg, U, V);
303  } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
304  matrix_exp_pade7(arg, U, V);
305  } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
306  matrix_exp_pade9(arg, U, V);
307  } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
308  matrix_exp_pade13(arg, U, V);
309  } else {
310  const long double maxnorm = 3.2579440895405400856599663723517L;
311  frexp(l1norm / maxnorm, &squarings);
312  if (squarings < 0) squarings = 0;
313  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
314  matrix_exp_pade17(A, U, V);
315  }
316 
317 #elif LDBL_MANT_DIG <= 112 // quadruple precison
318 
319  if (l1norm < 1.639394610288918690547467954466970e-005L) {
320  matrix_exp_pade3(arg, U, V);
321  } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
322  matrix_exp_pade5(arg, U, V);
323  } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
324  matrix_exp_pade7(arg, U, V);
325  } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
326  matrix_exp_pade9(arg, U, V);
327  } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
328  matrix_exp_pade13(arg, U, V);
329  } else {
330  const long double maxnorm = 2.884233277829519311757165057717815L;
331  frexp(l1norm / maxnorm, &squarings);
332  if (squarings < 0) squarings = 0;
333  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
334  matrix_exp_pade17(A, U, V);
335  }
336 
337 #else
338 
339  // this case should be handled in compute()
340  eigen_assert(false && "Bug in MatrixExponential");
341 
342 #endif
343 #endif // LDBL_MANT_DIG
344  }
345 };
346 
347 template<typename T> struct is_exp_known_type : false_type {};
348 template<> struct is_exp_known_type<float> : true_type {};
349 template<> struct is_exp_known_type<double> : true_type {};
350 #if LDBL_MANT_DIG <= 112
351 template<> struct is_exp_known_type<long double> : true_type {};
352 #endif
353 
354 template <typename ArgType, typename ResultType>
355 void matrix_exp_compute(const ArgType& arg, ResultType &result, true_type) // natively supported scalar type
356 {
357  typedef typename ArgType::PlainObject MatrixType;
358  MatrixType U, V;
359  int squarings;
360  matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
361  MatrixType numer = U + V;
362  MatrixType denom = -U + V;
363  result = denom.partialPivLu().solve(numer);
364  for (int i=0; i<squarings; i++)
365  result *= result; // undo scaling by repeated squaring
366 }
367 
368 
369 /* Computes the matrix exponential
370  *
371  * \param arg argument of matrix exponential (should be plain object)
372  * \param result variable in which result will be stored
373  */
374 template <typename ArgType, typename ResultType>
375 void matrix_exp_compute(const ArgType& arg, ResultType &result, false_type) // default
376 {
377  typedef typename ArgType::PlainObject MatrixType;
378  typedef typename traits<MatrixType>::Scalar Scalar;
379  typedef typename NumTraits<Scalar>::Real RealScalar;
380  typedef typename std::complex<RealScalar> ComplexScalar;
381  result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
382 }
383 
384 } // end namespace Eigen::internal
385 
396 template<typename Derived> struct MatrixExponentialReturnValue
397 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
398 {
399  typedef typename Derived::Index Index;
400  public:
405  MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
406 
411  template <typename ResultType>
412  inline void evalTo(ResultType& result) const
413  {
414  const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
415  internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
416  }
417 
418  Index rows() const { return m_src.rows(); }
419  Index cols() const { return m_src.cols(); }
420 
421  protected:
422  const typename internal::ref_selector<Derived>::type m_src;
423 };
424 
425 namespace internal {
426 template<typename Derived>
427 struct traits<MatrixExponentialReturnValue<Derived> >
428 {
429  typedef typename Derived::PlainObject ReturnType;
430 };
431 }
432 
433 template <typename Derived>
434 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
435 {
436  eigen_assert(rows() == cols());
437  return MatrixExponentialReturnValue<Derived>(derived());
438 }
439 
440 } // end namespace Eigen
441 
442 #endif // EIGEN_MATRIX_EXPONENTIAL
Eigen::MatrixExponentialReturnValue::evalTo
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:412
Eigen
Namespace containing all symbols from the Eigen library.
Eigen::MatrixBase::exp
const MatrixExponentialReturnValue< Derived > exp() const
Definition: MatrixExponential.h:434
Eigen::MatrixExponentialReturnValue
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:396
Eigen::MatrixExponentialReturnValue::MatrixExponentialReturnValue
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:405
Eigen::arg
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_arg_op< typename Derived::Scalar >, const Derived > arg(const Eigen::ArrayBase< Derived > &x)
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index