Hello,

I am currently sampling a model with lots of multivariate normal data points, and I am interested in using a multivariate normal density based on sufficient statistics to speed up gradient computation. I am aware of several github issues discussing this problem, but since these features are still under development I decided to try and implement the multivariate normal density directly in the model, using a log-probability update of the form:

```
target += n*multi_normal_lpdf(mu_hat | mu, sigma) - 0.5*(n-1)*trace(mdivide_left_spd(sigma, sigma_hat));
```

In the univariate normal case this worked beautifully, and the model sampled quickly and accurately. However, in the multivariate case the model sampled slowly and I had some issues with divergent transitions and exceeding maximum tree-depth. I thought this might have been because `mdivide_left_spd`

is difficult to differentiate in when `sigma`

is an actual matrix.

I noticed that in `stan-dev/math`

many matrix-inverse computations involving covariance matrices use the Cholesky decomposition and perform computations on the LDLT factor. For instance, https://github.com/stan-dev/math/blob/ab9b5b6ceecdf4a8382ceeccfbde97aff9279b9b/stan/math/prim/mat/prob/multi_normal_sufficient.hpp uses `mdivide_left_ldlt(ldlt_Sigma, sampleSigma)`

to calculate the same density.

To test if this would help, I tried to expose `mdivide_left_ldlt`

to the Stan language and use it in my model. To accomplish this i did the following:

- Fork
`stan-dev/stan`

and add`mdivide_left_ldlt`

to`function_signatures.hpp`

- Fork
`stan-dev/rstan`

and add`mdivide_left_ldlt`

to`rstan/src/stan/lang/function_signatures.hpp`

- Change rstanâ€™s â€śupstreamâ€ť gitmodule to point to my own version of Stan
- Change rstanâ€™s StanHeaders/install-github.R script to point to my own versions of rstan and stan
- run StanHeaders/install-github.R
- build and install my modified version of rstan (
`make build`

,`make install`

)

The repos where I made these alterations are https://github.com/jproney/rstan and https://github.com/jproney/stan both on the `develop`

branch.

When I try to compile my model (attached below) in R, the parser recognizes the `mdivide_left_ldlt`

function signature, but I get the following error:

```
Error in compileCode(f, code, language = language, verbose = verbose) :
Compilation ERROR, function(s)/method(s) not created! In file included from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/RcppEigen/include/Eigen/Core:392,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/RcppEigen/include/Eigen/Dense:1,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/mat/meta/append_return_type.hpp:4,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/meta.hpp:9,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/memory/stack_alloc.hpp:7,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/rev/core/autodiffstackstorage.hpp:4,
from /home/petey/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math
```

This error only occurs when I try to compile a model using `mdivide_ldlt`

; all other models are fine.

My two general questions are as follows: 1.) Is exposing `mdivide_ldlt`

even a promising approach here? Is there another way to achieve the result I want or explain the issues Iâ€™m having? 2.) Is there a quick way to fix this compiler error?

My stan model (using `mdivide_left_ldlt`

) is attached. To alter the model so that it compiles correctly, use `mdivide_left_spd`

instead and omit the Cholesky decomposition. I have also attached a data file and set of initial conditions which I passed to the model. I attempted to sample the model with the statement:

```
stan_mod <- stan_model(file = "clusters_2type.stan", verbose = TRUE)
fit_data <- sampling(stan_mod, data = dat, control = list(adapt_delta = 0.8), chains = 4, refresh = 1, init =init)
```

(`dat`

and `init`

are lists which can be found in `standata.rda`

at https://filebin.net/n1hydf9dqrl1qqhl)

My R session info is as follows:

```
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.10
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] bindrcpp_0.2.2 bpinference_0.0.0.9000 testthat_2.0.1 rstan_2.19.1 ggplot2_3.1.0 StanHeaders_2.18.1-10
loaded via a namespace (and not attached):
[1] deSolve_1.21 tidyselect_0.2.5 remotes_2.0.2 purrr_0.2.5 lattice_0.20-38 colorspace_1.4-0 expm_0.999-4 usethis_1.5.0 stats4_3.5.2 loo_2.0.0
[11] yaml_2.2.0 rlang_0.3.1 pkgbuild_1.0.2 pillar_1.3.1 glue_1.3.0 withr_2.1.2 sessioninfo_1.1.1 matrixStats_0.54.0 bindr_0.1.1 plyr_1.8.4
[21] stringr_1.4.0 munsell_0.5.0 gtable_0.2.0 devtools_2.0.2 memoise_1.1.0 inline_0.3.15 callr_3.1.1 ps_1.3.0 parallel_3.5.2 Rcpp_1.0.0
[31] scales_1.0.0 backports_1.1.3 desc_1.2.0 pkgload_1.0.2 fs_1.2.6 gridExtra_2.3 digest_0.6.18 stringi_1.2.4 processx_3.2.1 dplyr_0.7.8
[41] grid_3.5.2 rprojroot_1.3-2 cli_1.0.1 tools_3.5.2 magrittr_1.5 lazyeval_0.2.1 tibble_2.0.1 crayon_1.3.4 pkgconfig_2.0.2 Matrix_1.2-15
[51] prettyunits_1.0.2 assertthat_0.2.0 rstudioapi_0.9.0 R6_2.3.0 compiler_3.5.2
```

Thank you for taking the time to read this long post. I would greatly appreciate any suggestions!

clusters_2type.stan (6.9 KB)

James