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 MatrixType>
65 void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V)
66 {
67  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
68  const RealScalar b[] = {120., 60., 12., 1.};
69  const MatrixType A2 = A * A;
70  const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
71  U.noalias() = A * tmp;
72  V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
73 }
74 
80 template <typename MatrixType>
81 void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V)
82 {
83  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
84  const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.};
85  const MatrixType A2 = A * A;
86  const MatrixType A4 = A2 * A2;
87  const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
88  U.noalias() = A * tmp;
89  V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
90 }
91 
97 template <typename MatrixType>
98 void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V)
99 {
100  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
101  const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
102  const MatrixType A2 = A * A;
103  const MatrixType A4 = A2 * A2;
104  const MatrixType A6 = A4 * A2;
105  const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2
106  + b[1] * MatrixType::Identity(A.rows(), A.cols());
107  U.noalias() = A * tmp;
108  V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
109 
110 }
111 
117 template <typename MatrixType>
118 void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V)
119 {
120  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
121  const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
122  2162160., 110880., 3960., 90., 1.};
123  const MatrixType A2 = A * A;
124  const MatrixType A4 = A2 * A2;
125  const MatrixType A6 = A4 * A2;
126  const MatrixType A8 = A6 * A2;
127  const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
128  + b[1] * MatrixType::Identity(A.rows(), A.cols());
129  U.noalias() = A * tmp;
130  V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
131 }
132 
138 template <typename MatrixType>
139 void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V)
140 {
141  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
142  const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600.,
143  1187353796428800., 129060195264000., 10559470521600., 670442572800.,
144  33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
145  const MatrixType A2 = A * A;
146  const MatrixType A4 = A2 * A2;
147  const MatrixType A6 = A4 * A2;
148  V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
149  MatrixType tmp = A6 * V;
150  tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
151  U.noalias() = A * tmp;
152  tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
153  V.noalias() = A6 * tmp;
154  V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
155 }
156 
164 #if LDBL_MANT_DIG > 64
165 template <typename MatrixType>
166 void matrix_exp_pade17(const MatrixType &A, MatrixType &U, MatrixType &V)
167 {
168  typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
169  const RealScalar b[] = {830034394580628357120000.L, 415017197290314178560000.L,
170  100610229646136770560000.L, 15720348382208870400000.L,
171  1774878043152614400000.L, 153822763739893248000.L, 10608466464820224000.L,
172  595373117923584000.L, 27563570274240000.L, 1060137318240000.L,
173  33924394183680.L, 899510451840.L, 19554575040.L, 341863200.L, 4651200.L,
174  46512.L, 306.L, 1.L};
175  const MatrixType A2 = A * A;
176  const MatrixType A4 = A2 * A2;
177  const MatrixType A6 = A4 * A2;
178  const MatrixType A8 = A4 * A4;
179  V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
180  MatrixType tmp = A8 * V;
181  tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
182  + b[1] * MatrixType::Identity(A.rows(), A.cols());
183  U.noalias() = A * tmp;
184  tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
185  V.noalias() = tmp * A8;
186  V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2
187  + b[0] * MatrixType::Identity(A.rows(), A.cols());
188 }
189 #endif
190 
191 template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
192 struct matrix_exp_computeUV
193 {
201  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
202 };
203 
204 template <typename MatrixType>
205 struct matrix_exp_computeUV<MatrixType, float>
206 {
207  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings)
208  {
209  using std::frexp;
210  using std::pow;
211  const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
212  squarings = 0;
213  if (l1norm < 4.258730016922831e-001) {
214  matrix_exp_pade3(arg, U, V);
215  } else if (l1norm < 1.880152677804762e+000) {
216  matrix_exp_pade5(arg, U, V);
217  } else {
218  const float maxnorm = 3.925724783138660f;
219  frexp(l1norm / maxnorm, &squarings);
220  if (squarings < 0) squarings = 0;
221  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<float>(squarings));
222  matrix_exp_pade7(A, U, V);
223  }
224  }
225 };
226 
227 template <typename MatrixType>
228 struct matrix_exp_computeUV<MatrixType, double>
229 {
230  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings)
231  {
232  using std::frexp;
233  using std::pow;
234  const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
235  squarings = 0;
236  if (l1norm < 1.495585217958292e-002) {
237  matrix_exp_pade3(arg, U, V);
238  } else if (l1norm < 2.539398330063230e-001) {
239  matrix_exp_pade5(arg, U, V);
240  } else if (l1norm < 9.504178996162932e-001) {
241  matrix_exp_pade7(arg, U, V);
242  } else if (l1norm < 2.097847961257068e+000) {
243  matrix_exp_pade9(arg, U, V);
244  } else {
245  const double maxnorm = 5.371920351148152;
246  frexp(l1norm / maxnorm, &squarings);
247  if (squarings < 0) squarings = 0;
248  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<double>(squarings));
249  matrix_exp_pade13(A, U, V);
250  }
251  }
252 };
253 
254 template <typename MatrixType>
255 struct matrix_exp_computeUV<MatrixType, long double>
256 {
257  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings)
258  {
259 #if LDBL_MANT_DIG == 53 // double precision
260 
261  matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
262 
263 #else
264 
265  using std::frexp;
266  using std::pow;
267  const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
268  squarings = 0;
269 
270 #if LDBL_MANT_DIG <= 64 // extended precision
271 
272  if (l1norm < 4.1968497232266989671e-003L) {
273  matrix_exp_pade3(arg, U, V);
274  } else if (l1norm < 1.1848116734693823091e-001L) {
275  matrix_exp_pade5(arg, U, V);
276  } else if (l1norm < 5.5170388480686700274e-001L) {
277  matrix_exp_pade7(arg, U, V);
278  } else if (l1norm < 1.3759868875587845383e+000L) {
279  matrix_exp_pade9(arg, U, V);
280  } else {
281  const long double maxnorm = 4.0246098906697353063L;
282  frexp(l1norm / maxnorm, &squarings);
283  if (squarings < 0) squarings = 0;
284  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
285  matrix_exp_pade13(A, U, V);
286  }
287 
288 #elif LDBL_MANT_DIG <= 106 // double-double
289 
290  if (l1norm < 3.2787892205607026992947488108213e-005L) {
291  matrix_exp_pade3(arg, U, V);
292  } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
293  matrix_exp_pade5(arg, U, V);
294  } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
295  matrix_exp_pade7(arg, U, V);
296  } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
297  matrix_exp_pade9(arg, U, V);
298  } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
299  matrix_exp_pade13(arg, U, V);
300  } else {
301  const long double maxnorm = 3.2579440895405400856599663723517L;
302  frexp(l1norm / maxnorm, &squarings);
303  if (squarings < 0) squarings = 0;
304  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
305  matrix_exp_pade17(A, U, V);
306  }
307 
308 #elif LDBL_MANT_DIG <= 112 // quadruple precison
309 
310  if (l1norm < 1.639394610288918690547467954466970e-005L) {
311  matrix_exp_pade3(arg, U, V);
312  } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
313  matrix_exp_pade5(arg, U, V);
314  } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
315  matrix_exp_pade7(arg, U, V);
316  } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
317  matrix_exp_pade9(arg, U, V);
318  } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
319  matrix_exp_pade13(arg, U, V);
320  } else {
321  frexp(l1norm / maxnorm, &squarings);
322  if (squarings < 0) squarings = 0;
323  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<long double>(squarings));
324  matrix_exp_pade17(A, U, V);
325  }
326 
327 #else
328 
329  // this case should be handled in compute()
330  eigen_assert(false && "Bug in MatrixExponential");
331 
332 #endif
333 #endif // LDBL_MANT_DIG
334  }
335 };
336 
337 
338 /* Computes the matrix exponential
339  *
340  * \param arg argument of matrix exponential (should be plain object)
341  * \param result variable in which result will be stored
342  */
343 template <typename MatrixType, typename ResultType>
344 void matrix_exp_compute(const MatrixType& arg, ResultType &result)
345 {
346 #if LDBL_MANT_DIG > 112 // rarely happens
347  typedef typename traits<MatrixType>::Scalar Scalar;
348  typedef typename NumTraits<Scalar>::Real RealScalar;
349  typedef typename std::complex<RealScalar> ComplexScalar;
350  if (sizeof(RealScalar) > 14) {
351  result = arg.matrixFunction(StdStemFunctions<ComplexScalar>::exp);
352  return;
353  }
354 #endif
355  MatrixType U, V;
356  int squarings;
357  matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
358  MatrixType numer = U + V;
359  MatrixType denom = -U + V;
360  result = denom.partialPivLu().solve(numer);
361  for (int i=0; i<squarings; i++)
362  result *= result; // undo scaling by repeated squaring
363 }
364 
365 } // end namespace Eigen::internal
366 
377 template<typename Derived> struct MatrixExponentialReturnValue
378 : public ReturnByValue<MatrixExponentialReturnValue<Derived> >
379 {
380  typedef typename Derived::Index Index;
381  public:
386  MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
387 
392  template <typename ResultType>
393  inline void evalTo(ResultType& result) const
394  {
395  const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
396  internal::matrix_exp_compute(tmp, result);
397  }
398 
399  Index rows() const { return m_src.rows(); }
400  Index cols() const { return m_src.cols(); }
401 
402  protected:
403  const typename internal::ref_selector<Derived>::type m_src;
404 };
405 
406 namespace internal {
407 template<typename Derived>
408 struct traits<MatrixExponentialReturnValue<Derived> >
409 {
410  typedef typename Derived::PlainObject ReturnType;
411 };
412 }
413 
414 template <typename Derived>
415 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
416 {
417  eigen_assert(rows() == cols());
418  return MatrixExponentialReturnValue<Derived>(derived());
419 }
420 
421 } // end namespace Eigen
422 
423 #endif // EIGEN_MATRIX_EXPONENTIAL
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:393
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
Proxy for the matrix exponential of some matrix (expression).
Definition: MatrixExponential.h:377
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:386