Eigen-unsupported  5.0.1-dev+284dcc12
 
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
16// IWYU pragma: private
17#include "./InternalHeaderCheck.h"
18
19namespace Eigen {
20namespace internal {
21
26template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
27struct MatrixExponentialScalingOp {
28 using RealScalar = typename NumTraits<Scalar>::Real;
29
34 MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) {}
35
40 inline const Scalar operator()(const Scalar& x) const {
41 using std::ldexp;
42 return Scalar(ldexp(Eigen::numext::real(x), -m_squarings), ldexp(Eigen::numext::imag(x), -m_squarings));
43 }
44
45 private:
46 int m_squarings;
47};
48
49template <typename Scalar>
50struct MatrixExponentialScalingOp<Scalar, /*IsComplex=*/false> {
55 MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) {}
56
61 inline const Scalar operator()(const Scalar& x) const {
62 using std::ldexp;
63 return ldexp(x, -m_squarings);
64 }
65
66 private:
67 int m_squarings;
68};
69
75template <typename MatA, typename MatU, typename MatV>
76void matrix_exp_pade3(const MatA& A, MatU& U, MatV& V) {
77 typedef typename MatA::PlainObject MatrixType;
78 typedef typename NumTraits<typename traits<MatA>::Scalar>::Real RealScalar;
79 const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
80 const MatrixType A2 = A * A;
81 const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
82 U.noalias() = A * tmp;
83 V = b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
84}
85
91template <typename MatA, typename MatU, typename MatV>
92void matrix_exp_pade5(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[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
96 const MatrixType A2 = A * A;
97 const MatrixType A4 = A2 * A2;
98 const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
99 U.noalias() = A * tmp;
100 V = b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
101}
102
108template <typename MatA, typename MatU, typename MatV>
109void matrix_exp_pade7(const MatA& A, MatU& U, MatV& V) {
110 typedef typename MatA::PlainObject MatrixType;
111 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
112 const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
113 const MatrixType A2 = A * A;
114 const MatrixType A4 = A2 * A2;
115 const MatrixType A6 = A4 * A2;
116 const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
117 U.noalias() = A * tmp;
118 V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
119}
120
126template <typename MatA, typename MatU, typename MatV>
127void matrix_exp_pade9(const MatA& A, MatU& U, MatV& V) {
128 typedef typename MatA::PlainObject MatrixType;
129 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
130 const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
131 2162160.L, 110880.L, 3960.L, 90.L, 1.L};
132 const MatrixType A2 = A * A;
133 const MatrixType A4 = A2 * A2;
134 const MatrixType A6 = A4 * A2;
135 const MatrixType A8 = A6 * A2;
136 const MatrixType tmp =
137 b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
138 U.noalias() = A * tmp;
139 V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
140}
141
147template <typename MatA, typename MatU, typename MatV>
148void matrix_exp_pade13(const MatA& A, MatU& U, MatV& V) {
149 typedef typename MatA::PlainObject MatrixType;
150 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
151 const RealScalar b[] = {64764752532480000.L,
152 32382376266240000.L,
153 7771770303897600.L,
154 1187353796428800.L,
155 129060195264000.L,
156 10559470521600.L,
157 670442572800.L,
158 33522128640.L,
159 1323241920.L,
160 40840800.L,
161 960960.L,
162 16380.L,
163 182.L,
164 1.L};
165 const MatrixType A2 = A * A;
166 const MatrixType A4 = A2 * A2;
167 const MatrixType A6 = A4 * A2;
168 V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
169 MatrixType tmp = A6 * V;
170 tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
171 U.noalias() = A * tmp;
172 tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
173 V.noalias() = A6 * tmp;
174 V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
175}
176
184#if LDBL_MANT_DIG > 64
185template <typename MatA, typename MatU, typename MatV>
186void matrix_exp_pade17(const MatA& A, MatU& U, MatV& V) {
187 typedef typename MatA::PlainObject MatrixType;
188 typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
189 const RealScalar b[] = {830034394580628357120000.L,
190 415017197290314178560000.L,
191 100610229646136770560000.L,
192 15720348382208870400000.L,
193 1774878043152614400000.L,
194 153822763739893248000.L,
195 10608466464820224000.L,
196 595373117923584000.L,
197 27563570274240000.L,
198 1060137318240000.L,
199 33924394183680.L,
200 899510451840.L,
201 19554575040.L,
202 341863200.L,
203 4651200.L,
204 46512.L,
205 306.L,
206 1.L};
207 const MatrixType A2 = A * A;
208 const MatrixType A4 = A2 * A2;
209 const MatrixType A6 = A4 * A2;
210 const MatrixType A8 = A4 * A4;
211 V = b[17] * A8 + b[15] * A6 + b[13] * A4 + b[11] * A2; // used for temporary storage
212 MatrixType tmp = A8 * V;
213 tmp += b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
214 U.noalias() = A * tmp;
215 tmp = b[16] * A8 + b[14] * A6 + b[12] * A4 + b[10] * A2;
216 V.noalias() = tmp * A8;
217 V += b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
218}
219#endif
220
221template <typename MatrixType, typename RealScalar = typename NumTraits<typename traits<MatrixType>::Scalar>::Real>
222struct matrix_exp_computeUV {
230 static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
231};
232
233template <typename MatrixType>
234struct matrix_exp_computeUV<MatrixType, float> {
235 using Scalar = typename traits<MatrixType>::Scalar;
236 template <typename ArgType>
237 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
238 using std::frexp;
239 using std::pow;
240 const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
241 squarings = 0;
242 if (l1norm < 4.258730016922831e-001f) {
243 matrix_exp_pade3(arg, U, V);
244 } else if (l1norm < 1.880152677804762e+000f) {
245 matrix_exp_pade5(arg, U, V);
246 } else {
247 const float maxnorm = 3.925724783138660f;
248 frexp(l1norm / maxnorm, &squarings);
249 if (squarings < 0) squarings = 0;
250 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
251 matrix_exp_pade7(A, U, V);
252 }
253 }
254};
255
256template <typename MatrixType>
257struct matrix_exp_computeUV<MatrixType, double> {
258 using Scalar = typename traits<MatrixType>::Scalar;
259 template <typename ArgType>
260 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
261 using std::frexp;
262 using std::pow;
263 const double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
264 squarings = 0;
265 if (l1norm < 1.495585217958292e-002) {
266 matrix_exp_pade3(arg, U, V);
267 } else if (l1norm < 2.539398330063230e-001) {
268 matrix_exp_pade5(arg, U, V);
269 } else if (l1norm < 9.504178996162932e-001) {
270 matrix_exp_pade7(arg, U, V);
271 } else if (l1norm < 2.097847961257068e+000) {
272 matrix_exp_pade9(arg, U, V);
273 } else {
274 const double maxnorm = 5.371920351148152;
275 frexp(l1norm / maxnorm, &squarings);
276 if (squarings < 0) squarings = 0;
277 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
278 matrix_exp_pade13(A, U, V);
279 }
280 }
281};
282
283template <typename MatrixType>
284struct matrix_exp_computeUV<MatrixType, long double> {
285 template <typename ArgType>
286 static void run(const ArgType& arg, MatrixType& U, MatrixType& V, int& squarings) {
287#if LDBL_MANT_DIG == 53 // double precision
288 matrix_exp_computeUV<MatrixType, double>::run(arg, U, V, squarings);
289
290#else
291
292 using Scalar = typename traits<MatrixType>::Scalar;
293
294 using std::frexp;
295 using std::pow;
296 const long double l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
297 squarings = 0;
298
299#if LDBL_MANT_DIG <= 64 // extended precision
300
301 if (l1norm < 4.1968497232266989671e-003L) {
302 matrix_exp_pade3(arg, U, V);
303 } else if (l1norm < 1.1848116734693823091e-001L) {
304 matrix_exp_pade5(arg, U, V);
305 } else if (l1norm < 5.5170388480686700274e-001L) {
306 matrix_exp_pade7(arg, U, V);
307 } else if (l1norm < 1.3759868875587845383e+000L) {
308 matrix_exp_pade9(arg, U, V);
309 } else {
310 const long double maxnorm = 4.0246098906697353063L;
311 frexp(l1norm / maxnorm, &squarings);
312 if (squarings < 0) squarings = 0;
313 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
314 matrix_exp_pade13(A, U, V);
315 }
316
317#elif LDBL_MANT_DIG <= 106 // double-double
318
319 if (l1norm < 3.2787892205607026992947488108213e-005L) {
320 matrix_exp_pade3(arg, U, V);
321 } else if (l1norm < 6.4467025060072760084130906076332e-003L) {
322 matrix_exp_pade5(arg, U, V);
323 } else if (l1norm < 6.8988028496595374751374122881143e-002L) {
324 matrix_exp_pade7(arg, U, V);
325 } else if (l1norm < 2.7339737518502231741495857201670e-001L) {
326 matrix_exp_pade9(arg, U, V);
327 } else if (l1norm < 1.3203382096514474905666448850278e+000L) {
328 matrix_exp_pade13(arg, U, V);
329 } else {
330 const long double maxnorm = 3.2579440895405400856599663723517L;
331 frexp(l1norm / maxnorm, &squarings);
332 if (squarings < 0) squarings = 0;
333 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
334 matrix_exp_pade17(A, U, V);
335 }
336
337#elif LDBL_MANT_DIG <= 113 // quadruple precision
338
339 if (l1norm < 1.639394610288918690547467954466970e-005L) {
340 matrix_exp_pade3(arg, U, V);
341 } else if (l1norm < 4.253237712165275566025884344433009e-003L) {
342 matrix_exp_pade5(arg, U, V);
343 } else if (l1norm < 5.125804063165764409885122032933142e-002L) {
344 matrix_exp_pade7(arg, U, V);
345 } else if (l1norm < 2.170000765161155195453205651889853e-001L) {
346 matrix_exp_pade9(arg, U, V);
347 } else if (l1norm < 1.125358383453143065081397882891878e+000L) {
348 matrix_exp_pade13(arg, U, V);
349 } else {
350 const long double maxnorm = 2.884233277829519311757165057717815L;
351 frexp(l1norm / maxnorm, &squarings);
352 if (squarings < 0) squarings = 0;
353 MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<Scalar>(squarings));
354 matrix_exp_pade17(A, U, V);
355 }
356
357#else
358
359 // this case should be handled in compute()
360 eigen_assert(false && "Bug in MatrixExponential");
361
362#endif
363#endif // LDBL_MANT_DIG
364 }
365};
366
367template <typename T>
368struct is_exp_known_type : false_type {};
369template <>
370struct is_exp_known_type<float> : true_type {};
371template <>
372struct is_exp_known_type<double> : true_type {};
373#if LDBL_MANT_DIG <= 113
374template <>
375struct is_exp_known_type<long double> : true_type {};
376#endif
377
378template <typename ArgType, typename ResultType>
379void matrix_exp_compute(const ArgType& arg, ResultType& result, true_type) // natively supported scalar type
380{
381 typedef typename ArgType::PlainObject MatrixType;
382 MatrixType U, V;
383 int squarings;
384 matrix_exp_computeUV<MatrixType>::run(arg, U, V, squarings); // Pade approximant is (U+V) / (-U+V)
385 MatrixType numer = U + V;
386 MatrixType denom = -U + V;
387 result = denom.partialPivLu().solve(numer);
388 for (int i = 0; i < squarings; i++) result *= result; // undo scaling by repeated squaring
389}
390
391/* Computes the matrix exponential
392 *
393 * \param arg argument of matrix exponential (should be plain object)
394 * \param result variable in which result will be stored
395 */
396template <typename ArgType, typename ResultType>
397void matrix_exp_compute(const ArgType& arg, ResultType& result, false_type) // default
398{
399 typedef typename ArgType::PlainObject MatrixType;
400 typedef make_complex_t<typename traits<MatrixType>::Scalar> ComplexScalar;
401 result = arg.matrixFunction(internal::stem_function_exp<ComplexScalar>);
402}
403
404} // namespace internal
405
416template <typename Derived>
417struct MatrixExponentialReturnValue : public ReturnByValue<MatrixExponentialReturnValue<Derived> > {
418 public:
423 MatrixExponentialReturnValue(const Derived& src) : m_src(src) {}
424
429 template <typename ResultType>
430 inline void evalTo(ResultType& result) const {
431 const typename internal::nested_eval<Derived, 10>::type tmp(m_src);
432 internal::matrix_exp_compute(tmp, result, internal::is_exp_known_type<typename Derived::RealScalar>());
433 }
434
435 Index rows() const { return m_src.rows(); }
436 Index cols() const { return m_src.cols(); }
437
438 protected:
439 const typename internal::ref_selector<Derived>::type m_src;
440};
441
442namespace internal {
443template <typename Derived>
444struct traits<MatrixExponentialReturnValue<Derived> > {
445 typedef typename Derived::PlainObject ReturnType;
446};
447} // namespace internal
448
449template <typename Derived>
450const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const {
451 eigen_assert(rows() == cols());
452 return MatrixExponentialReturnValue<Derived>(derived());
453}
454
455} // end namespace Eigen
456
457#endif // EIGEN_MATRIX_EXPONENTIAL
Eigen::MatrixBase::exp
const MatrixExponentialReturnValue< Derived > exp() const
Definition MatrixExponential.h:450
Eigen
Namespace containing all symbols from the Eigen library.
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
Eigen::arg
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_arg_op< typename Derived::Scalar >, const Derived > arg(const Eigen::ArrayBase< Derived > &x)
Eigen::MatrixExponentialReturnValue
Proxy for the matrix exponential of some matrix (expression).
Definition MatrixExponential.h:417
Eigen::MatrixExponentialReturnValue::evalTo
void evalTo(ResultType &result) const
Compute the matrix exponential.
Definition MatrixExponential.h:430
Eigen::MatrixExponentialReturnValue::MatrixExponentialReturnValue
MatrixExponentialReturnValue(const Derived &src)
Constructor.
Definition MatrixExponential.h:423
Eigen::internal::MatrixExponentialScalingOp::operator()
const Scalar operator()(const Scalar &x) const
Scale a matrix coefficient.
Definition MatrixExponential.h:40
Eigen::internal::MatrixExponentialScalingOp::MatrixExponentialScalingOp
MatrixExponentialScalingOp(int squarings)
Constructor.
Definition MatrixExponential.h:34
Eigen::internal::matrix_exp_computeUV
Compute the (17,17)-Padé approximant to the exponential.
Definition MatrixExponential.h:222
Eigen::internal::matrix_exp_computeUV::run
static void run(const MatrixType &arg, MatrixType &U, MatrixType &V, int &squarings)
Compute Padé approximant to the exponential.