Hi,

I have a spatiotemporal hierarchical model with a latent processes on which I place a Gaussian process prior. The model has been tested and validated and fits my data well, but sampling is rather slow. To make the model more computational tractable and to be able to apply it to larger data sets, I am implementing a nearest neighbor approximation for the Gaussian process. However, in doing this, I run into a puzzling issue. I will not show the entire model as it is rather involved (the model has been published here: https://www.pnas.org/content/117/4/1877) and will try to explain the issue in simple terms and with just a few lines of code.

The relevant code reads (note the distinction between the original model and the nearest neighbor model):

```
parameters {
vector[n] a;
...
}
transformed parameters {
vector[n] w; // this is the latent process
// in the original model I have:
L = cholesky_decompose(Q); // Q is the full covariance matrix
// in the nearest neighbor model I have:
L = nngp_chol(...); // nngp_chol is a user-defined function that computes the nearest neighbor Cholesky factor of the covariance matrix Q
w = L * a;
...
}
model{
a ~ normal(0,1);
...
}
```

By setting the number of neighbors equal to ‘n-1’, the nearest neighbor approximation should produce the exact full matrix ‘Q’. Indeed, printing variable ‘w’ in this case shows exactly the same output in both implementations of the model. That is,

L = cholesky_decompose(Q) and L = nngp_chol(…) produce exactly the same Cholesky factor. However, the model only runs when using the first. When I use L = nngp_chol(…) Stan throws the following error:

Log probability evaluates to log(0), i.e., negative infinity. ]

Stan can’t start sampling from this initial value.

This is puzzling to me. Apart from the line L = cholesky_decompose(Q) (or L = nngp_chol()) the two models are exactly the same, and both models produce the same ‘L’, but only one works. Any idea why this might be happening?

Thanks,

Louis