The ctsem software (hierarchical state space modelling in continuous / discrete time, model coded in stan) requires the evaluation of two matrix exponentials (as well as a bunch of other, less costly, matrix operations) for each timepoint of data, one for the state matrix (governing how the system behaves with respect to itself) and one for the jacobian of the state matrix. This allows an extended Kalman filtering approach to integrate out the latent states.

Last week I was asked about a large problem where the system was running out of memory, and ended up addressing it by modifying ctsem to pre-compute independent subsets of the state and jacobian matrices, then loop over these subsets and compose the full matrices. This gave about a 300% improvement in performance and memory usage, which is wonderful. I guess the same approach could be applied to other matrix operations so thought mentioning it might be helpful, though donâ€™t know if it would have the same impact. Since there was such an improvement tucked away for relatively little difficulty implementing, Iâ€™m wondering if anyone has other ideas I might be missing out on :)

subset matrix exponential function for the interestedâ€¦ where the subsets integer array is of size n by nrows, n is the number of subsets, nrows is the number of rows of the square matrix m. the first elements of each subsets[i,] array denote the indices of the subset, the rest should be zero.

```
matrix expmSubsets(matrix m, int[,] subsets){
int nr = rows(m);
matrix[nr,nr] e = rep_matrix(0,nr,nr);
for(si in 1:size(subsets)){
int n=0;
for(j in 1:nr){
if(subsets[si,j]!=0) n+=1;
}
if(n > 1) e[subsets[si,1:n],subsets[si,1:n] ] = matrix_exp(m[subsets[si,1:n],subsets[si,1:n]]);
if(n == 1) e[subsets[si,1],subsets[si,1] ] = exp(m[subsets[si,1],subsets[si,1]]);
}
return e;
}
```