I have read, here and there across the months, about several models in Stan that involved a system of linear ODEs, and where the use of a numerical integrator was faster than using a matrix exponential.
Is there a good theory on when this might happen and perhaps examples? My thinking is cases where we require the solutions at many closeby (relatively speaking) time points.
No example sorry, but at some point I was working on a model with a sparse transition matrix (mostly diagonal) and found good speedups by writing a wrapper around the matrix_exp function to check for zeroes in the relevant off diagonals and compute accordingly. I just ditched it from my own code as part of ‘chuck everything overboard it won’t compile on 32 bit and it’s cran submission time’, but here it is ;)
int[] checkoffdiagzero(matrix M){
int z[rows(M)];
for(i in 1:rows(M)){
z[i] = 0;
for(j in 1:cols(M)){
if(i!=j){
if(M[i,j] != 0){
z[i] = 1;
break;
}
}
}
if(z[i]==0){ //check rows
for(j in 1:rows(M)){
if(i!=j){
if(M[j,i] != 0){
z[i] = 1;
break;
}
}
}
}
}
return z;
}
matrix matrix_exp(matrix M){
matrix[rows(M),rows(M)] out;
int z[rows(out)] = checkoffdiagzero(M);
int z1[sum(z)];
int z0[rows(M)-sum(z)];
int cz1 = 1;
int cz0 = 1;
for(i in 1:rows(M)){
if(z[i] == 1){
z1[cz1] = i;
cz1 += 1;
}
if(z[i] == 0){
z0[cz0] = i;
cz0 += 1;
}
}
if(size(z1) > 0) out[z1,z1] = matrix_exp(M[z1,z1]);
out[z0,] = rep_matrix(0,size(z0),rows(M));
out[,z0] = rep_matrix(0,rows(M),size(z0));
for(i in 1:size(z0)) out[z0[i],z0[i]] = exp(M[z0[i],z0[i]]);
return out;
}Depends what the matrix exponential function is using under the hood. A lot of the times they’re just using a Runge-Kutta method or a Pade approximation. So depending on the properties of the ODE, e.g. stiffness, something like BDF could be more appropriate than whatever the matrix exponential is using. As another example, if the problem is Hamiltonian, then I bet a symplectic integrator would do better than a matrix exponential implementation.
The matrix exp we have uses a Pade approximation, coupled with scaling and squaring. The matrix exp action function uses the method proposed by Al-Mohy & Higam (2011).
Do you think it would have a hard time dealing with stiff systems?
p.s.: love the paper you referred to.
A-priori, I’m not sure if it’d do worse or better on stiff systems than say, BDF, but it wouldn’t be too hard to try. Just experiment with a 2D system with one eigenvalue that’s like -1e6 and another one that’s -1e0. You can just use the matrix exponential and DeSolve in R.
I’d also try the Hamiltonian system
and see how the Pade approximation compares to leapfrog. Then also try it with a stiff Hamiltonian system, by replacing the -1 with say -1e6. If you tried forward Euler on this system you’d notice that over time you’re gaining energy so that if you try to integrate for a long time the solutions will be pretty bad compared to a symplectic integrator like leapfrog. If you tried backward Euler you’d also get bad solutions but you’d be losing energy over time. I’m not sure what the Pade would do, but I’d be curious.
These are all good suggestions, I’ll start playing with it.
It might be that uneven scales mess up the pade approximation because of arithmetic precision error (the scaling and squaring may not work as well). It indeed shouldn’t be too hard to test. I’ll use RStan’s function to extract Stan functions and run some rest in R – or maybe I’ll directly write some tests in C++.
Yeah that’s a good point. It could be that an eigenvalue of 1e-6 isn’t bad all on its own, but the uneven scales are what give the Pade trouble. I don’t know enough about the matrix exp algorithms enough to say, but uneven scales are always trouble numerically.
This seems like an excellent situation for speculative evaluation with multithreading.