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.