Automatic Differentiation
 
Loading...
Searching...
No Matches
matrix_exp_action_handler.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
2#define STAN_MATH_PRIM_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
3
7#include <unsupported/Eigen/MatrixFunctions>
8#include <vector>
9#include <cmath>
10
11namespace stan {
12namespace math {
13
26 const int _p_max = 8;
27 const int _m_max = 55;
28 const double _tol = 1.1e-16; // from the paper, double precision: 2^-53
29 const std::vector<double> _theta_m{
30 2.22044605e-16, 2.58095680e-08, 1.38634787e-05, 3.39716884e-04,
31 2.40087636e-03, 9.06565641e-03, 2.38445553e-02, 4.99122887e-02,
32 8.95776020e-02, 1.44182976e-01, 2.14235807e-01, 2.99615891e-01,
33 3.99777534e-01, 5.13914694e-01, 6.41083523e-01, 7.80287426e-01,
34 9.30532846e-01, 1.09086372e+00, 1.26038106e+00, 1.43825260e+00,
35 1.62371595e+00, 1.81607782e+00, 2.01471078e+00, 2.21904887e+00,
36 2.42858252e+00, 2.64285346e+00, 2.86144963e+00, 3.08400054e+00,
37 3.31017284e+00, 3.53966635e+00, 3.77221050e+00, 4.00756109e+00,
38 4.24549744e+00, 4.48581986e+00, 4.72834735e+00, 4.97291563e+00,
39 5.21937537e+00, 5.46759063e+00, 5.71743745e+00, 5.96880263e+00,
40 6.22158266e+00, 6.47568274e+00, 6.73101590e+00, 6.98750228e+00,
41 7.24506843e+00, 7.50364669e+00, 7.76317466e+00, 8.02359473e+00,
42 8.28485363e+00, 8.54690205e+00, 8.80969427e+00, 9.07318789e+00,
43 9.33734351e+00, 9.60212447e+00, 9.86749668e+00, 1.01334283e+01,
44 1.03998897e+01, 1.06668532e+01, 1.09342929e+01, 1.12021845e+01,
45 1.14705053e+01, 1.17392341e+01, 1.20083509e+01, 1.22778370e+01,
46 1.25476748e+01, 1.28178476e+01, 1.30883399e+01, 1.33591369e+01,
47 1.36302250e+01, 1.39015909e+01, 1.41732223e+01, 1.44451076e+01,
48 1.47172357e+01, 1.49895963e+01, 1.52621795e+01, 1.55349758e+01,
49 1.58079765e+01, 1.60811732e+01, 1.63545578e+01, 1.66281227e+01,
50 1.69018609e+01, 1.71757655e+01, 1.74498298e+01, 1.77240478e+01,
51 1.79984136e+01, 1.82729215e+01, 1.85475662e+01, 1.88223426e+01,
52 1.90972458e+01, 1.93722711e+01, 1.96474142e+01, 1.99226707e+01,
53 2.01980367e+01, 2.04735082e+01, 2.07490816e+01, 2.10247533e+01,
54 2.13005199e+01, 2.15763782e+01, 2.18523250e+01, 2.21283574e+01};
55
56 public:
64 template <typename EigMat1, typename EigMat2,
67 inline Eigen::MatrixXd action(const EigMat1& mat, const EigMat2& b,
68 const double& t = 1.0) {
69 Eigen::MatrixXd A = mat;
70 double mu = A.trace() / A.rows();
71 A.diagonal().array() -= mu;
72
73 int m, s;
74 set_approx_order(A, b, t, m, s);
75
76 double eta = exp(t * mu / s);
77
78 const auto& b_eval = b.eval();
79 Eigen::MatrixXd f = b_eval;
80 Eigen::MatrixXd bi = b_eval;
81
82 for (int i = 0; i < s; ++i) {
83 // maximum absolute row sum, aka inf-norm of matrix operator
84 double c1 = matrix_operator_inf_norm(bi);
85 for (int j = 0; j < m; ++j) {
86 bi = (t / (s * (j + 1))) * (A * bi);
87 f += bi;
88 double c2 = matrix_operator_inf_norm(bi);
89 if (c1 + c2 < _tol * matrix_operator_inf_norm(f))
90 break;
91 c1 = c2;
92 }
93 f *= eta;
94 bi = f;
95 }
96
97 return f;
98 }
99
105 double matrix_operator_inf_norm(Eigen::MatrixXd const& x) {
106 return x.cwiseAbs().rowwise().sum().maxCoeff();
107 }
108
126 template <typename EigMat1, require_all_eigen_t<EigMat1>* = nullptr,
127 require_all_st_same<double, EigMat1>* = nullptr>
128 double mat_power_1_norm(const EigMat1& mat, int m) {
129 if ((mat.array() > 0.0).all()) {
130 Eigen::VectorXd e = Eigen::VectorXd::Constant(mat.rows(), 1.0);
131 for (int j = 0; j < m; ++j) {
132 e = mat.transpose() * e;
133 }
134 return e.lpNorm<Eigen::Infinity>();
135 } else {
136 return mat.pow(m).cwiseAbs().colwise().sum().maxCoeff();
137 }
138 }
139
156 template <typename EigMat1, typename EigMat2,
159 inline void set_approx_order(const EigMat1& mat, const EigMat2& b,
160 const double& t, int& m, int& s) {
161 if (t < _tol) {
162 m = 0;
163 s = 1;
164 return;
165 }
166
167 // L1 norm
168 double normA = mat.colwise().template lpNorm<1>().maxCoeff();
169 if (normA < _tol) {
170 m = 0;
171 s = 1;
172 return;
173 }
174
175 Eigen::VectorXd alpha(_p_max - 1);
176 if (normA * (_m_max * b.cols())
177 < 4.0 * _theta_m[_m_max] * _p_max * (_p_max + 3)) {
178 alpha = Eigen::VectorXd::Constant(_p_max - 1, normA);
179 } else {
180 Eigen::VectorXd eta(_p_max);
181 for (int p = 0; p < _p_max; ++p) {
182 eta[p] = std::pow(mat_power_1_norm(mat, p + 2), 1.0 / (p + 2));
183 }
184 for (int p = 0; p < _p_max - 1; ++p) {
185 alpha[p] = std::max(eta[p], eta[p + 1]);
186 }
187 }
188
189 Eigen::MatrixXd mt = Eigen::MatrixXd::Zero(_p_max - 1, _m_max);
190 for (int p = 1; p < _p_max; ++p) {
191 for (int i = p * (p + 1) - 2; i < _m_max; ++i) {
192 mt(p - 1, i) = alpha[p - 1] / _theta_m[i];
193 }
194 }
195
196 Eigen::VectorXd u = Eigen::VectorXd::LinSpaced(_m_max, 1, _m_max);
197 Eigen::MatrixXd c = stan::math::ceil(mt) * u.asDiagonal();
198 for (Eigen::MatrixXd::Index i = 0; i < c.size(); ++i) {
199 if (c(i) <= 1.e-16) {
200 c(i) = std::numeric_limits<double>::infinity();
201 }
202 }
203 int cost = c.colwise().minCoeff().minCoeff(&m);
204 if (std::isinf(cost)) {
205 s = 1;
206 return;
207 }
208 s = std::max(cost / m, 1);
209 }
210};
211
212} // namespace math
213} // namespace stan
214
215#endif
void set_approx_order(const EigMat1 &mat, const EigMat2 &b, const double &t, int &m, int &s)
Approximation is based on parameter "m" and "s", proposed in CODE FRAGMENT 3.1 of the reference.
double matrix_operator_inf_norm(Eigen::MatrixXd const &x)
Eigen expression for matrix operator infinity norm.
double mat_power_1_norm(const EigMat1 &mat, int m)
Estimate the 1-norm of mat^m.
Eigen::MatrixXd action(const EigMat1 &mat, const EigMat2 &b, const double &t=1.0)
Perform the matrix exponential action exp(A*t)*B.
The implementation of the work by Awad H.
require_all_t< is_eigen< std::decay_t< Types > >... > require_all_eigen_t
Require all of the types satisfy is_eigen.
Definition is_eigen.hpp:65
require_all_t< std::is_same< scalar_type_t< std::decay_t< T > >, scalar_type_t< std::decay_t< Types > > >... > require_all_st_same
All scalar types of T and all of the Types satisfy std::is_same.
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
fvar< T > ceil(const fvar< T > &x)
Definition ceil.hpp:12
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
bool isinf(const stan::math::var &a)
Return 1 if the specified argument is positive infinity or negative infinity and 0 otherwise.
Definition std_isinf.hpp:16