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