Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
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
16namespace Eigen {
17namespace internal {
18
23template <typename RealScalar>
24struct 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
64template <typename MatA, typename MatU, typename MatV>
65void 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
81template <typename MatA, typename MatU, typename MatV>
82void 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
99template <typename MatA, typename MatU, typename MatV>
100void 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
120template <typename MatA, typename MatU, typename MatV>
121void 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
142template <typename MatA, typename MatU, typename MatV>
143void 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
170template <typename MatA, typename MatU, typename MatV>
171void 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
197template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
198struct matrix_exp_computeUV
199{
207 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
208};
209
210template <typename MatrixType>
211struct 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
234template <typename MatrixType>
235struct 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
263template <typename MatrixType>
264struct 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 <= 113 // quadruple precision
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
347template<typename T> struct is_exp_known_type : false_type {};
348template<> struct is_exp_known_type<float> : true_type {};
349template<> struct is_exp_known_type<double> : true_type {};
350#if LDBL_MANT_DIG <= 113
351template<> struct is_exp_known_type<long double> : true_type {};
352#endif
353
354template <typename ArgType, typename ResultType>
355void 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 */
374template <typename ArgType, typename ResultType>
375void 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
396template<typename Derived> struct MatrixExponentialReturnValue
397: public ReturnByValue<MatrixExponentialReturnValue<Derived> >
398{
399 public:
404 MatrixExponentialReturnValue(const Derived& src) : m_src(src) { }
405
410 template <typename ResultType>
411 inline void evalTo(ResultType& result) const
412 {
413 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
414 internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
415 }
416
417 Index rows() const { return m_src.rows(); }
418 Index cols() const { return m_src.cols(); }
419
420 protected:
421 const typename internal::ref_selector<Derived>::type m_src;
422};
423
424namespace internal {
425template<typename Derived>
426struct traits<MatrixExponentialReturnValue<Derived> >
427{
428 typedef typename Derived::PlainObject ReturnType;
429};
430}
431
432template <typename Derived>
433const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
434{
435 eigen_assert(rows() == cols());
436 return MatrixExponentialReturnValue<Derived>(derived());
437}
438
439} // end namespace Eigen
440
441#endif // EIGEN_MATRIX_EXPONENTIAL
const MatrixExponentialReturnValue< Derived > exp() const
Definition: MatrixExponential.h:433
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:398
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition: MatrixExponential.h:411
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition: MatrixExponential.h:404