Since I also rely on the matrix exponential, I was curious about this. The following code matches the above but includes my own re-write of the Higham08 implementation into a stan function. The same difference occurs, suggesting either that a) the c++ stan function of CharlesM has some minor failing in exactly the same way as my own code, or b) lower level stan functions / numeric issues are generating the difference and it’s nothing to do with the matrix exponential coding. I’m assuming the latter!

```
st='functions{
matrix expm_stan(matrix A, real time_point){
matrix[5,5] A_test;
A_test = matrix_exp(A * time_point);
return A_test;
}
}
'
st2='functions{
matrix expmp(matrix A, matrix padeC, vector padeCbig){
int n;
real nA;
real colsum;
int l;
matrix[4,10] C;
vector[4] t;
matrix[rows(A),rows(A)] I;
matrix[rows(A),rows(A)] P;
matrix[rows(A),rows(A)] U;
matrix[rows(A),rows(A)] V;
matrix[rows(A),rows(A)] X;
vector[14] Cbig;
real s;
real si;
matrix[rows(A),rows(A)] B;
matrix[rows(A),rows(A)] B2;
matrix[rows(A),rows(A)] B4;
matrix[rows(A),rows(A)] B6;
matrix[rows(A),rows(A)] A2;
si = 0;
C = padeC;
Cbig = padeCbig;
n =rows(A);
if(n != cols(A)) print("expmp: Matrix not square!")
// if (n <= 1) X = exp(A);
if(sum(A)==sum(diagonal(A))) X = diag_matrix(exp(diagonal(A)));
else{
// nA = Matrix::norm(A, "1")
nA = 0;
for(coli in 1:n){
colsum=0;
for(rowi in 1:n){
colsum=colsum+fabs(A[rowi,coli]);
}
if(colsum > nA) nA = colsum;
}
I = diag_matrix(rep_vector(1,n));
if (nA <= 2.1) {
t[1] = 0.015; t[2]= 0.25; t[3]= 0.95; t[4]= 2.1;
//l = which.max(nA <= t)
for(ti in 1:4){
if(l==0){
if(nA <= t[ti]) l = ti;
}
}
A2 = A * A;
P = I;
U = C[l, 2] * I;
V = C[l, 1] * I;
for (k in 1:l) {
P = P * A2;
U = U + C[l, (2 * k) + 2] * P;
V = V + C[l, (2 * k) + 1] * P;
}
U = A * U;
X = inverse(V - U) * (V + U);
}
else {
s = log2(nA/5.4);
B = A;
if (s > 0) {
s = ceil(s);
B = B/(2^s);
}
B2 = B * B;
B4 = B2 * B2;
B6 = B2*B4;
U = B*(B6*(Cbig[14] * B6 + Cbig[12] * B4 + Cbig[10] * B2) + Cbig[8] * B6 + Cbig[6] * B4 + Cbig[4] * B2 + Cbig[2] * I);
V = B6*(Cbig[13] * B6 + Cbig[11] * B4 + Cbig[9] * B2) + Cbig[7] * B6 + Cbig[5] * B4 + Cbig[3] * B2 + Cbig[1] * I;
X = inverse(V - U) * (V + U);
if (s > 0) {
while (si < s){
si = si + 1;
X = X * X;
}
}
}
}
return X;
}
}
data {
matrix[4,10] padeC; // for matrix exponential
vector[14] padeCbig;
}
'
library(rstan)
library(expm)
rate_one <- 10
rate_two <- 5
scale_one <- 1000
A <- matrix(0, nrow = 5, ncol = 5)
A[1,1] <- - rate_one
A[2,1] <- rate_one
A[2,2] <- -rate_one
A[3,2] <- rate_one
A[3,3] <- -rate_two
A[4,3] <- rate_two
A[4,4] <- -rate_two
A[5,4] <- rate_two * scale_one
padeC=rbind(c(120, 60, 12, 1, 0, 0, 0, 0, 0, 0), c(30240,
15120, 3360, 420, 30, 1, 0, 0, 0, 0), c(17297280,
8648640, 1995840, 277200, 25200, 1512, 56, 1, 0,
0), c(17643225600, 8821612800, 2075673600, 302702400,
30270240, 2162160, 110880, 3960, 90, 1))
padeCbig= matrix(c(64764752532480000, 32382376266240000, 7771770303897600,
1187353796428800, 129060195264000, 10559470521600,
670442572800, 33522128640, 1323241920, 40840800,
960960, 16380, 182, 1))
expose_stan_functions(stanc(model_code = st2))
results_R <- expm(A * 10,method=c('Higham08'))
results_stan <- expm_stan(A, 10) #stan approach
results_stan2 <- expmp(A * 10,padeC,padeCbig) #higham08 approach written in stan
```