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/stanand addmdivide_left_ldlttofunction_signatures.hpp - Fork
stan-dev/rstanand addmdivide_left_ldlttorstan/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