Matrix exponential function

Modeling
1

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)
1 Like
2

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

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.

4

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

5

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