Documentation for the action of the matrix exponential

I went through the documentation for solving system of linear ODEs, and saw that the prescribed method is to compute the full matrix exponential and then multiply by the initial condition vector y_0. A few years ago, @yizhang did implement a method which directly computes the action of the matrix exponential, see this PR.

Do we need to update the documentation or did the matrix exp action not make it in the final release? I haven’t dealt with linear ODEs in a while, so I want to check with people who might have been using this feature.

@yizhang @wds15


There was a bug in my implementation so it was turned off. I’m in the middle of fixing it, see


Excellent! Let me know if you need someone to review the upcoming PR.

1 Like

Is this faster? I’d assume so, but is it?

Depends on the nature & size of the matrix. In the PR there’s a unit test based on discretizing the Poisson equation, which gives a tridiagonal matrix with some off diagonal entries. Below is the benchmark with discretizing grid size = 5, 10, 20. (wall time in nanosec)

For matrices of dim < 10 I’d imagine there’d be not much difference.

1 Like

So what I’m seeing is what, a 4-16 times speedup for a single matrix exponential action with/out AD?

Again, the speedup(or the lack of it) depends on the matrix. You can find more benchmark results from the original paper:

It’s w/o AD. AD w.r.t matrix exponential is more involved and I need more time. Fixing this old bug is where I’d like to begin.

1 Like

The above PR fix matrix_exp_multiply accuracy bug by yizhang-yiz · Pull Request #2619 · stan-dev/math · GitHub adds new unit tests comparing the different matrix exponentiation methods on problematic matrices harvested from reported stan issues. looks great!

If time permits, it might also be worth stacking on additional unit tests comparing both methods on some other matrices that are known to cause problems to some methods when exponentiated.

The following 3x3 matrix is communicated by Moler in A Balancing Act for the Matrix Exponential » Cleve’s Corner: Cleve Moler on Mathematics and Computing - MATLAB & Simulink :

a = 2e10;
b = 4e8/6;
c = 200/3;
d = 3;
e = 1e-8;
A = [0 e 0; -(a+b) -d a; c 0 -c]

That matrix defeats one of the three matlab demo matrix exponentiation routines.

Another small example is given on page 8 of Moler, Van Loan - Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later :

[-49 24 ; -64 31]

That 2x2 matrix is apparently enough to defeat a Taylor series approximation to the matrix exponential – albeit one computed on an IBM 370 with only 7-exponent bits in its hexadecimal floating point format.

To turn these into unit tests for the action of the matrix exponential, we could sweep over a basis of vectors b and compute exp(A) b via each method.


Thanks @rfc . I just added the two matrices you brought up to unit tests. See math/matrix_exp_multiply_test.cpp at fac3d0c1b9f5b36998f76f6a788ef7bf25dcb814 · stan-dev/math · GitHub


@yizhang thank you! I see you also defined the unit tests cases to compare against known values of expm constant matrices, rather than whatever the full Stan matrix exponential computes, which is a better test design than what I was suggesting – more specific and also more robust against future changes to the code that might break or bias both Stan matrix exponential implementations in the same way.