MatrixExponential.h
Go to the documentation of this file.
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 // This file was edited for to the Stan math library to create
12 // the matrix exponential function (matrix_exp), 2016.
13 
14 #ifndef STAN_MATH_PRIM_MAT_FUN_MATRIXEXPONENTIAL_H
15 #define STAN_MATH_PRIM_MAT_FUN_MATRIXEXPONENTIAL_H
16 
18 
19 namespace Eigen {
20 
21  template <typename RealScalar>
23  {
24  /** \brief Constructor.
25  *
26  * \param[in] squarings The integer \f$ s \f$ in this document.
27  */
28  MatrixExponentialScalingOp(int squarings) : m_squarings(squarings) { }
29 
30 
31  /** \brief Scale a matrix coefficient.
32  *
33  * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
34  */
35  inline const RealScalar operator() (const RealScalar& x) const
36  {
37  using std::ldexp;
38  return ldexp(x, -m_squarings);
39  }
40 
41  typedef std::complex<RealScalar> ComplexScalar;
42 
43  /** \brief Scale a matrix coefficient.
44  *
45  * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$.
46  */
47  inline const ComplexScalar operator() (const ComplexScalar& x) const
48  {
49  using std::ldexp;
50  return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings));
51  }
52 
53  private:
55  };
56 
57 
58  /** \brief Compute the (5,5)-Pad&eacute; approximant to the exponential.
59  *
60  * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
61  * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
62  */
63  template <typename MatrixType>
64  void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V)
65  {
66  using Eigen::internal::traits;
67  typedef typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
68  const RealScalar b[] = {120.L, 60.L, 12.L, 1.L};
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 
75  /** \brief Compute the (5,5)-Pad&eacute; approximant to the exponential.
76  *
77  * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
78  * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
79  */
80  template <typename MatrixType>
81  void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V)
82  {
83  typedef typename Eigen::NumTraits<typename Eigen::internal::traits<MatrixType>::Scalar>::Real RealScalar;
84  const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L};
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 
92  /** \brief Compute the (7,7)-Pad&eacute; approximant to the exponential.
93  *
94  * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
95  * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
96  */
97  template <typename MatrixType>
98  void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V)
99  {
100  using Eigen::internal::traits;
101  typedef typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
102  const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L};
103  const MatrixType A2 = A * A;
104  const MatrixType A4 = A2 * A2;
105  const MatrixType A6 = A4 * A2;
106  const MatrixType tmp = b[7] * A6 + b[5] * A4 + b[3] * A2
107  + b[1] * MatrixType::Identity(A.rows(), A.cols());
108  U.noalias() = A * tmp;
109  V = b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
110  }
111 
112  /** \brief Compute the (9,9)-Pad&eacute; approximant to the exponential.
113  *
114  * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
115  * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
116  */
117  template <typename MatrixType>
118  void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V)
119  {
120  using Eigen::internal::traits;
121  typedef typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
122  const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L,
123  2162160.L, 110880.L, 3960.L, 90.L, 1.L};
124  const MatrixType A2 = A * A;
125  const MatrixType A4 = A2 * A2;
126  const MatrixType A6 = A4 * A2;
127  const MatrixType A8 = A6 * A2;
128  const MatrixType tmp = b[9] * A8 + b[7] * A6 + b[5] * A4 + b[3] * A2
129  + b[1] * MatrixType::Identity(A.rows(), A.cols());
130  U.noalias() = A * tmp;
131  V = b[8] * A8 + b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
132  }
133 
134  /** \brief Compute the (13,13)-Pad&eacute; approximant to the exponential.
135  *
136  * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
137  * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
138  */
139  template <typename MatrixType>
140  void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V)
141  {
142  using Eigen::internal::traits;
143  typedef typename Eigen::NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar;
144  const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L,
145  1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L,
146  33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L};
147  const MatrixType A2 = A * A;
148  const MatrixType A4 = A2 * A2;
149  const MatrixType A6 = A4 * A2;
150  V = b[13] * A6 + b[11] * A4 + b[9] * A2; // used for temporary storage
151  MatrixType tmp = A6 * V;
152  tmp += b[7] * A6 + b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols());
153  U.noalias() = A * tmp;
154  tmp = b[12] * A6 + b[10] * A4 + b[8] * A2;
155  V.noalias() = A6 * tmp;
156  V += b[6] * A6 + b[4] * A4 + b[2] * A2 + b[0] * MatrixType::Identity(A.rows(), A.cols());
157  }
158 
159  template <typename MatrixType, typename RealScalar = typename Eigen:: NumTraits<typename Eigen::internal::traits<MatrixType>::Scalar>::Real>
161  {
162  /** \brief Compute Pad&eacute; approximant to the exponential.
163  *
164  * Computes \c U, \c V and \c squarings such that \f$ (V+U)(V-U)^{-1} \f$ is a Pad&eacute;
165  * approximant of \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$, where \f$ M \f$
166  * denotes the matrix \c arg. The degree of the Pad&eacute; approximant and the value of squarings
167  * are chosen such that the approximation error is no more than the round-off error.
168  *
169  * <p> Edit for Stan: template ComputeUV::run so that it may used on
170  * autodiff variables (var and fvar). This required adding the scalar_type
171  * argument, which tells the function the type of elements used in the
172  * matrix.
173  */
174  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings);
175  };
176 
177  template <typename MatrixType>
178  struct matrix_exp_computeUV<MatrixType>
179  {
180  template <typename T>
181  static void run(const MatrixType& arg, MatrixType& U, MatrixType& V, int& squarings,
182  T scalar_type)
183  {
184  using std::frexp;
185  using std::pow;
187  const T l1norm = arg.cwiseAbs().colwise().sum().maxCoeff();
188  squarings = 0;
189  if (l1norm < 1.495585217958292e-002) {
190  matrix_exp_pade3(arg, U, V);
191  } else if (l1norm < 2.539398330063230e-001) {
192  matrix_exp_pade5(arg, U, V);
193  } else if (l1norm < 9.504178996162932e-001) {
194  matrix_exp_pade7(arg, U, V);
195  } else if (l1norm < 2.097847961257068e+000) {
196  matrix_exp_pade9(arg, U, V);
197  } else {
198  const double maxnorm = 5.371920351148152;
199  frexp(value_of_rec(l1norm) / value_of_rec(maxnorm), &squarings);
200  if (squarings < 0) squarings = 0;
201  MatrixType A = arg.unaryExpr(MatrixExponentialScalingOp<T>(squarings));
202  matrix_exp_pade13(A, U, V);
203  }
204  }
205  };
206 
207 }
208 
209 #endif
void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V)
Compute the (9,9)-Padé approximant to the exponential.
std::complex< RealScalar > ComplexScalar
double value_of_rec(const fvar< T > &v)
T ldexp(const T &a, int b)
Definition: ldexp.hpp:19
static void run(const MatrixType &arg, MatrixType &U, MatrixType &V, int &squarings, T scalar_type)
constexpr T pow(T x)
Definition: pow.h:75
Definition: StanUtils.h:6
Float_t tmp
Definition: plot.C:36
void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V)
Compute the (5,5)-Padé approximant to the exponential.
MatrixExponentialScalingOp(int squarings)
Constructor.
const RealScalar operator()(const RealScalar &x) const
Scale a matrix coefficient.
void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V)
Compute the (13,13)-Padé approximant to the exponential.
void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V)
Compute the (7,7)-Padé approximant to the exponential.
static const double A
Definition: Units.h:82
const hit & b
Definition: hits.cxx:21
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V)
Compute the (5,5)-Padé approximant to the exponential.