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.
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)
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.
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.
@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.