Hi everyone,

I’ve got a big model with about half million observations. I’m modeling environmental variables associated with species occurrence records with multivariate normals, and I’m estimating separate MVNs for each species and decade hierarchically. Because `multi_normal_cholesky`

requires a matrix for the Cholesky factor of the covariance matrix, I believe I have to call the likelihood inside for-loops. Even though all of the records for each species and decade can be vectorised within the for-loops, I’m still stuck with `n_species * n_decade`

likelihood calls, leading to the following model template:

```
data {
// ..
int<lower=1> n_obs, n_species, n_decade, n_dim;
array[n_obs] vector[n_dim] y;
// ..
}
parameters {
// ..
array[n_species, n_decade] vector[n_dim] mu;
array[n_species, n_decade] vector<lower=0>[n_dim] sigma;
array[n_species, n_decade] cholesky_factor_corr[n_dim] Omega_L;
// ..
}
model {
// ..
for (i in 1:n_species) {
for (j in 1:n_decade) {
y[i, j] ~ multi_normal_cholesky(mu[i, j], diag_pre_multiply(sigma[i, j], Omega_L[i, j])
}
}
// ..
}
```

At the moment, I’ve written a partial_sum function to parallelise the likelihood, but this happens for each species/decade combination separately. The model has become a bit more complex with `first`

and `n_obs`

giving the index of the first observation and the number of records for each species/decade, respectively, and there are also observation-level random effects leading to the new array of mean vectors `mu_f`

.

```
functions {
real partial_sum_lpdf(array[] vector y_slice, int start, int end, array[] vector mu, matrix Sigma_L) {
return multi_normal_cholesky_lupdf(y_slice | mu[start:end], Sigma_L);
}
}
model {
// ..
for (i in 1:n_species) {
for (j in 1:n_decade) {
target += reduce_sum(partial_sum_lupdf,
segment(y, first[i, j], n_obs[i, j]),
grainsize[i, j],
mu_f[i, j, 1:n_obs[i, j]],
diag_pre_multiply(sigma[i, j], Omega_L[i, j]));
}
}
// ..
}
```

My question is, is there some obvious way to parallelise the model further or more efficiently? TL;DR, my main concern is that the likelihood requires a matrix for the covariance for the Cholesky factor, and I can’t think of an alternative approach to what I’m currently using, i.e. an array of matrices.

Thanks!

Matt