Stan Math Library  2.20.0
reverse mode automatic differentiation
matrix_exp_action_handler.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
3 
5 #include <vector>
6 
7 namespace stan {
8 namespace math {
9 
10 /*
11  * The implemention of the work by Awad H. Al-Mohy and Nicholas J. Higham
12  * "Computing the Action of the Matrix Exponential,
13  * with an Application to Exponential Integrators"
14  * Read More: https://epubs.siam.org/doi/abs/10.1137/100788860
15  *
16  * Calculates exp(mat*t)*b, where mat & b are matrices,
17  * and t is double.
18  */
20  static constexpr int p_max = 8;
21  static constexpr int m_max = 55;
22  static constexpr double tol = 1.1e-16;
23 
24  // table 3.1 in the reference
25  const std::vector<double> theta_m_single_precision{
26  1.3e-1, 1.0e0, 2.2e0, 3.6e0, 4.9e0, 6.3e0,
27  7.7e0, 9.1e0, 1.1e1, 1.2e1, 1.3e1};
28  const std::vector<double> theta_m_double_precision{
29  2.4e-3, 1.4e-1, 6.4e-1, 1.4e0, 2.4e0, 3.5e0,
30  4.7e0, 6.0e0, 7.2e0, 8.5e0, 9.9e0};
31 
32  double l1norm(const Eigen::MatrixXd& m) {
33  return m.colwise().lpNorm<1>().maxCoeff();
34  }
35 
36  public:
37  /* Constructor
38  */
40 
41  /* Perform the matrix exponential action exp(A*t)*B
42  * @param [in] mat matrix A
43  * @param [in] b matrix B
44  * @param [in] t double t, e.g. time.
45  * @return matrix exp(A*t)*B
46  */
47  inline Eigen::MatrixXd action(const Eigen::MatrixXd& mat,
48  const Eigen::MatrixXd& b,
49  const double& t = 1.0) {
50  Eigen::MatrixXd A = mat;
51  double mu = A.trace() / A.rows();
52  for (int i = 0; i < A.rows(); ++i) {
53  A(i, i) -= mu;
54  }
55 
56  int m{0}, s{0};
57  set_approximation_parameter(A, t, m, s);
58 
59  Eigen::MatrixXd res(A.rows(), b.cols());
60 
61  for (int col = 0; col < b.cols(); ++col) {
62  bool conv = false;
63  Eigen::VectorXd B = b.col(col);
64  Eigen::VectorXd F = B;
65  const auto eta = std::exp(t * mu / s);
66  for (int i = 1; i < s + 1; ++i) {
67  auto c1 = B.template lpNorm<Eigen::Infinity>();
68  if (m > 0) {
69  for (int j = 1; j < m + 1; ++j) {
70  B = t * A * B / (s * j);
71  auto c2 = B.template lpNorm<Eigen::Infinity>();
72  F += B;
73  if (c1 + c2 < tol * F.template lpNorm<Eigen::Infinity>()) {
74  conv = true;
75  break;
76  }
77  c1 = c2;
78  }
79  }
80  F *= eta;
81  B = F;
82  if (conv)
83  break;
84  }
85  res.col(col) = F;
86  } // loop b columns
87  return res;
88  }
89 
90  /* Approximation is based on parameter "m" and "s",
91  * proposed in CODE FRAGMENT 3.1 of the reference. The
92  * parameters depend on matrix A, as well as limits
93  * "m_max" & "p_max", defined in table 3.1 and eq. 3.11,
94  * respectively. Ideally one can choose "m_max" and
95  * "p_max" to suite the specific computing precision needs,
96  * here we just use the one suggested in the paper
97  * (paragraph between eq. 3.11 & eq. 3.12).
98  *
99  * @param [in] mat matrix A
100  * @param [in] t double t in exp(A*t)*B
101  * @param [out] m int parameter m
102  * @param [out] s int parameter s
103  */
104  inline void set_approximation_parameter(const Eigen::MatrixXd& mat,
105  const double& t, int& m, int& s) {
106  if (l1norm(mat) < tol || t < tol) {
107  m = 0;
108  s = 1;
109  } else {
110  m = m_max;
111 
112  Eigen::MatrixXd a = mat * t;
113  const double theta_m = theta_m_double_precision.back();
114  for (auto i = 0; i < std::ceil(std::log2(p_max)); ++i) {
115  a *= a;
116  }
117  double ap = std::pow(l1norm(a), 1.0 / p_max);
118  int c = std::ceil(ap / theta_m);
119  s = (c < 1 ? 1 : c);
120  }
121  }
122 };
123 
124 } // namespace math
125 } // namespace stan
126 
127 #endif
Eigen::MatrixXd action(const Eigen::MatrixXd &mat, const Eigen::MatrixXd &b, const double &t=1.0)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
Eigen::Matrix< T, Eigen::Dynamic, 1 > col(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, size_t j)
Return the specified column of the specified matrix using start-at-1 indexing.
Definition: col.hpp:23
void set_approximation_parameter(const Eigen::MatrixXd &mat, const double &t, int &m, int &s)
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
fvar< T > log2(const fvar< T > &x)
Return the base two logarithm of the specified argument.
Definition: log2.hpp:20
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:12

     [ Stan Home Page ] © 2011–2018, Stan Development Team.