# Matrix exponential function

Dear all,

I am attempting to solve a system of ordinary differental equations, using matrix exponentials, at a number of different timepoints. I have used the expose_stan_functions code to access my stan function (using RStan) and compared the Stan matrix exponential, at time 10, to the matrix exponential produced by the R ‘expm’ function, with an example matrix of the differential equation parameter values.

When using the Stan function, with these parameter values at this time (and a few other times but not all), the second to last row of the matrix exponential produces negative values. Whilst the ‘expm’ function in R does not.

Other than these differences the only difference between the matrix exponentials are small floating point errors.

Any help would be much appreciated.

Many thanks,
Isaac

“expm_example.stan”

``````functions{
matrix expm_stan(matrix A, real time_point){
matrix[5,5] A_test;
A_test = matrix_exp(A * time_point);
return A_test;
}
}
``````
``````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

expose_stan_functions(stanc("expm_example.stan"))

results_R <- expm(A * 10)
results_stan <- expm_stan(A, 10)
``````

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;

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
}

'

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``````
1 Like

It’d be good to figure out which one of `matrix_exp` or `expm` is returning an incorrect result. @IJS91 could you solve the problem with the numerical integrator (maybe bdf, to be conservative), and report the result?

I can also do it, though I’ll probably only get to it later this week.

expm() has a lot of different methods that all generate positive values, just fyi…

@Charles_Driver Thanks a lot for doing this, knowing that has saved me a lot of time. @charlesm93 Yes - happy to do that. @Charles_Driver yes I was thinking they should be positive.