Eigen-Contrib  5.0.1
 
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// SPDX-License-Identifier: MPL-2.0
11
12#ifndef EIGEN_MATRIX_EXPONENTIAL
13#define EIGEN_MATRIX_EXPONENTIAL
14
15#include "StemFunction.h"
16
17// IWYU pragma: private
18#include "./InternalHeaderCheck.h"
19
20namespace Eigen {
21namespace internal {
22
24template <typename MatrixType,
25 bool HasWritableRealView = !NumTraits<typename traits<MatrixType>::Scalar>::IsComplex ||
26 complex_array_access<typename traits<MatrixType>::Scalar>::value>
27struct matrix_exp_scale_impl {
28 template <typename ArgType>
29 static MatrixType run(const ArgType& arg, int squarings) {
30 MatrixType result(arg.rows(), arg.cols());
31 result.realView().array() = arg.realView().array().ldexp(-squarings);
32 return result;
33 }
34};
35
36template <typename MatrixType>
37struct matrix_exp_scale_impl<MatrixType, false> {
38 template <typename ArgType>
39 static MatrixType run(const ArgType& arg, int squarings) {
40 using Scalar = typename traits<MatrixType>::Scalar;
41 return arg.unaryExpr([squarings](const Scalar& x) {
42 using std::ldexp;
43 return Scalar(ldexp(numext::real(x), -squarings), ldexp(numext::imag(x), -squarings));
44 });
45 }
46};
47
48template <typename MatrixType, typename ArgType>
49MatrixType matrix_exp_scale(const ArgType& arg, int squarings) {
50 return matrix_exp_scale_impl<MatrixType>::run(arg, squarings);
51}
52
58template <typename MatA, typename MatU, typename MatV>
59void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V) {
60 typedef typename MatA::PlainObject MatrixType;
61 typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
62 const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
63 const MatrixType A2 = A * A;
64 const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
65 U.noalias() = A * tmp;
66 V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
67}
68
74template <typename MatA, typename MatU, typename MatV>
75void matrix_exp_pade5(const MatA& A, MatU& U, MatV& V) {
76 typedef typename MatA::PlainObject MatrixType;
77 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
78 const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
79 const MatrixType A2 = A * A;
80 const MatrixType A4 = A2 * A2;
81 const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
82 U.noalias() = A * tmp;
83 V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
84}
85
91template <typename MatA, typename MatU, typename MatV>
92void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V) {
93 typedef typename MatA::PlainObject MatrixType;
94 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
95 const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
96 const MatrixType A2 = A * A;
97 const MatrixType A4 = A2 * A2;
98 const MatrixType A6 = A4 * A2;
99 const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
100 U.noalias() = A * tmp;
101 V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
102}
103
109template <typename MatA, typename MatU, typename MatV>
110void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V) {
111 typedef typename MatA::PlainObject MatrixType;
112 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
113 const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
114 2162160.L, 110880.L, 3960.L, 90.L, 1.L};
115 const MatrixType A2 = A * A;
116 const MatrixType A4 = A2 * A2;
117 const MatrixType A6 = A4 * A2;
118 const MatrixType A8 = A6 * A2;
119 const MatrixType tmp =
120 b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
121 U.noalias() = A * tmp;
122 V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
123}
124
130template <typename MatA, typename MatU, typename MatV>
131void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V) {
132 typedef typename MatA::PlainObject MatrixType;
133 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
134 const RealScalar b[] = {64764752532480000.L,
135 32382376266240000.L,
136 7771770303897600.L,
137 1187353796428800.L,
138 129060195264000.L,
139 10559470521600.L,
140 670442572800.L,
141 33522128640.L,
142 1323241920.L,
143 40840800.L,
144 960960.L,
145 16380.L,
146 182.L,
147 1.L};
148 const MatrixType A2 = A * A;
149 const MatrixType A4 = A2 * A2;
150 const MatrixType A6 = A4 * A2;
151 V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
152 MatrixType tmp = A6 * V;
153 tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
154 U.noalias() = A * tmp;
155 tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
156 V.noalias() = A6 * tmp;
157 V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
158}
159
167#if LDBL_MANT_DIG > 64
168template <typename MatA, typename MatU, typename MatV>
169void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V) {
170 typedef typename MatA::PlainObject MatrixType;
171 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
172 const RealScalar b[] = {830034394580628357120000.L,
173 415017197290314178560000.L,
174 100610229646136770560000.L,
175 15720348382208870400000.L,
176 1774878043152614400000.L,
177 153822763739893248000.L,
178 10608466464820224000.L,
179 595373117923584000.L,
180 27563570274240000.L,
181 1060137318240000.L,
182 33924394183680.L,
183 899510451840.L,
184 19554575040.L,
185 341863200.L,
186 4651200.L,
187 46512.L,
188 306.L,
189 1.L};
190 const MatrixType A2 = A * A;
191 const MatrixType A4 = A2 * A2;
192 const MatrixType A6 = A4 * A2;
193 const MatrixType A8 = A4 * A4;
194 V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
195 MatrixType tmp = A8 * V;
196 tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
197 U.noalias() = A * tmp;
198 tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
199 V.noalias() = tmp * A8;
200 V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
201}
202#endif
203
204template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
205struct matrix_exp_computeUV {
213 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
214};
215
216template <typename MatrixType>
217struct matrix_exp_computeUV<MatrixType, float> {
218 template <typename ArgType>
219 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
220 using std::frexp;
221 using std::pow;
222 const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
223 squarings = 0;
224 if (l1norm < 4.258730016922831e-001f) {
225 matrix_exp_pade3(arg, U, V);
226 } else if (l1norm < 1.880152677804762e+000f) {
227 matrix_exp_pade5(arg, U, V);
228 } else {
229 const float maxnorm = 3.925724783138660f;
230 frexp(l1norm / maxnorm, &squarings);
231 if (squarings < 0) squarings = 0;
232 MatrixType A = matrix_exp_scale<MatrixType>(arg, squarings);
233 matrix_exp_pade7(A, U, V);
234 }
235 }
236};
237
238template <typename MatrixType>
239struct matrix_exp_computeUV<MatrixType, double> {
240 template <typename ArgType>
241 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
242 using std::frexp;
243 using std::pow;
244 const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
245 squarings = 0;
246 if (l1norm < 1.495585217958292e-002) {
247 matrix_exp_pade3(arg, U, V);
248 } else if (l1norm < 2.539398330063230e-001) {
249 matrix_exp_pade5(arg, U, V);
250 } else if (l1norm < 9.504178996162932e-001) {
251 matrix_exp_pade7(arg, U, V);
252 } else if (l1norm < 2.097847961257068e+000) {
253 matrix_exp_pade9(arg, U, V);
254 } else {
255 const double maxnorm = 5.371920351148152;
256 frexp(l1norm / maxnorm, &squarings);
257 if (squarings < 0) squarings = 0;
258 MatrixType A = matrix_exp_scale<MatrixType>(arg, squarings);
259 matrix_exp_pade13(A, U, V);
260 }
261 }
262};
263
264template <typename MatrixType>
265struct matrix_exp_computeUV<MatrixType, long double> {
266 template <typename ArgType>
267 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
268#if LDBL_MANT_DIG == 53 // double precision
269 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
270
271#else
272
273 using std::frexp;
274 using std::pow;
275 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
276 squarings = 0;
277
278#if LDBL_MANT_DIG <= 64 // extended precision
279
280 if (l1norm < 4.1968497232266989671e-003L) {
281 matrix_exp_pade3(arg, U, V);
282 } else if (l1norm < 1.1848116734693823091e-001L) {
283 matrix_exp_pade5(arg, U, V);
284 } else if (l1norm < 5.5170388480686700274e-001L) {
285 matrix_exp_pade7(arg, U, V);
286 } else if (l1norm < 1.3759868875587845383e+000L) {
287 matrix_exp_pade9(arg, U, V);
288 } else {
289 const long double maxnorm = 4.0246098906697353063L;
290 frexp(l1norm / maxnorm, &squarings);
291 if (squarings < 0) squarings = 0;
292 MatrixType A = matrix_exp_scale<MatrixType>(arg, squarings);
293 matrix_exp_pade13(A, U, V);
294 }
295
296#elif LDBL_MANT_DIG <= 106 // double-double
297
298 if (l1norm < 3.2787892205607026992947488108213e-005L) {
299 matrix_exp_pade3(arg, U, V);
300 } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
301 matrix_exp_pade5(arg, U, V);
302 } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
303 matrix_exp_pade7(arg, U, V);
304 } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
305 matrix_exp_pade9(arg, U, V);
306 } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
307 matrix_exp_pade13(arg, U, V);
308 } else {
309 const long double maxnorm = 3.2579440895405400856599663723517L;
310 frexp(l1norm / maxnorm, &squarings);
311 if (squarings < 0) squarings = 0;
312 MatrixType A = matrix_exp_scale<MatrixType>(arg, squarings);
313 matrix_exp_pade17(A, U, V);
314 }
315
316#elif LDBL_MANT_DIG <= 113 // quadruple precision
317
318 if (l1norm < 1.639394610288918690547467954466970e-005L) {
319 matrix_exp_pade3(arg, U, V);
320 } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
321 matrix_exp_pade5(arg, U, V);
322 } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
323 matrix_exp_pade7(arg, U, V);
324 } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
325 matrix_exp_pade9(arg, U, V);
326 } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
327 matrix_exp_pade13(arg, U, V);
328 } else {
329 const long double maxnorm = 2.884233277829519311757165057717815L;
330 frexp(l1norm / maxnorm, &squarings);
331 if (squarings < 0) squarings = 0;
332 MatrixType A = matrix_exp_scale<MatrixType>(arg, squarings);
333 matrix_exp_pade17(A, U, V);
334 }
335
336#else
337
338 // this case should be handled in compute()
339 eigen_assert(false && "Bug in MatrixExponential");
340
341#endif
342#endif // LDBL_MANT_DIG
343 }
344};
345
346template <typename T>
347using is_exp_known_type = bool_constant<std::is_same<T, float>::value || std::is_same<T, double>::value
348#if LDBL_MANT_DIG <= 113
349 || std::is_same<T, long double>::value
350#endif
351 >;
352
353template <typename ArgType, typename ResultType>
354void matrix_exp_compute(const ArgType& arg, ResultType& result, std::true_type) // natively supported scalar type
355{
356 typedef typename ArgType::PlainObject MatrixType;
357 MatrixType U, V;
358 int squarings;
359 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
360 MatrixType numer = U + V;
361 MatrixType denom = -U + V;
362 result = denom.partialPivLu().solve(numer);
363 for (int i = 0; i < squarings; i++) result *= result; // undo scaling by repeated squaring
364}
365
366/* Computes the matrix exponential
367 *
368 * \param arg argument of matrix exponential (should be plain object)
369 * \param result variable in which result will be stored
370 */
371template <typename ArgType, typename ResultType>
372void matrix_exp_compute(const ArgType& arg, ResultType& result, std::false_type) // default
373{
374 typedef typename ArgType::PlainObject MatrixType;
375 typedef make_complex_t<typename traits<MatrixType>::Scalar> ComplexScalar;
376 result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
377}
378
379} // namespace internal
380
391template <typename Derived>
392struct MatrixExponentialReturnValue : public ReturnByValue<MatrixExponentialReturnValue<Derived> > {
393 public:
398 MatrixExponentialReturnValue(const Derived& src) : m_src(src) {}
399
404 template <typename ResultType>
405 inline void evalTo(ResultType& result) const {
406 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
407 internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
408 }
409
410 Index rows() const { return m_src.rows(); }
411 Index cols() const { return m_src.cols(); }
412
413 protected:
414 const typename internal::ref_selector<Derived>::type m_src;
415};
416
417namespace internal {
418template <typename Derived>
419struct traits<MatrixExponentialReturnValue<Derived> > {
420 typedef typename Derived::PlainObject ReturnType;
421};
422} // namespace internal
423
424template <typename Derived>
425const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const {
426 eigen_assert(rows() == cols());
427 return MatrixExponentialReturnValue<Derived>(derived());
428}
429
430} // end namespace Eigen
431
432#endif // EIGEN_MATRIX_EXPONENTIAL
Eigen::MatrixBase::exp
const MatrixExponentialReturnValue< Derived > exp() const
Definition MatrixExponential.h:425
Eigen
Namespace containing all symbols from the Eigen library.
Eigen::MatrixExponentialReturnValue
Proxy for the matrix exponential of some matrix (expression).
Definition MatrixExponential.h:392
Eigen::MatrixExponentialReturnValue::evalTo
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition MatrixExponential.h:405
Eigen::MatrixExponentialReturnValue::MatrixExponentialReturnValue
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition MatrixExponential.h:398
Eigen::internal::matrix_exp_computeUV
Compute the (17,17)-Padé approximant to the exponential.
Definition MatrixExponential.h:205
Eigen::internal::matrix_exp_computeUV::run
static void run(const MatrixType &arg, MatrixType &U, MatrixType &V, int &squarings)
Compute Padé approximant to the exponential.
Eigen::internal::matrix_exp_scale_impl
Scale a real or complex matrix by .
Definition MatrixExponential.h:27