I was working on my sde / state space code today, and at some point became frustrated with Stan occasionally complaining about non spd matrices during Cholesky decomposition, even though the printed matrices could be decomposed in R. So I wrote a user defined Cholesky function, and was surprised at the fairly substantial speed up this gave to my models – given that the Cholesky decomp is just a modest piece of the puzzle. This is with smallish (2-6 dimensional) matrices, I didn’t test any further. I also don’t know why this is occurring – whether comp time or precision / gradients, or if something peculiar to the cases I tried – but thought I should mention it. This R script compares output of my function to the stan function.

```
sm <- '
functions{
matrix chol(matrix a){
matrix[rows(a),rows(a)] l;
for(j in 1:cols(a)){
for(i in 1:rows(a)){
if(j==i) {
if(j == 1){
if(a[j,j] <=0) {
l[j,j] = 1e-8;
print("Negative variance ", a[j,j], " set to 1e-8 for Cholesky decomp");
}
if(a[j,j] > 0) l[j,j] = sqrt(a[j,j]);
}
if(j > 1) l[j,j] = sqrt(a[j,j] - dot_self(l[i, 1 : j-1]));
}
if(i > j) l[i,j] = ( a[i,j] - dot_product( l[ i, 1:(j-1) ], l[j, 1:(j-1)]) ) / l[j,j];
if(j > i) l[i,j] = 0;
}
}
return l;
}
}
data{
int dim;
}
parameters{
cov_matrix[dim] cmat;
}
model{
print("cov = ", cmat);
print("Chol diff = ",chol(cmat)-cholesky_decompose(cmat));
}'
stan(model_code = sm,data = list(dim=4),iter=1)
```