Multivariate normal sufficient statistics + exposing mdivide_ldlt

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

Hmm, I’m not sure what that error is.

The issue I see with exposing mdivide_left_ldlt is that the first argument (the LDLT_Factor) has no equivalent Stan type (this function: https://github.com/stan-dev/math/blob/develop/stan/math/rev/mat/fun/mdivide_left_ldlt.hpp#L296).

The easiest way to do something like this (and take advantage of whatever trace_inv_quad_form_ldlt has to offer) would probably be to use the Rstan external C++ stuff: https://cran.r-project.org/web/packages/rstan/vignettes/external.html . That ends up looking sloppier, but it might be easier to work with in terms of trying to get your code working.

Thanks for the reply. The problem actually had to do with the way I was attempting to compress the data into summary statistics, and not with the gradient computation. The summary statistics which I thought would be sufficient to replicate the posterior of the entire dataset actually did not contain enough information, so the posterior ended up being flat along certain dimensions.

1 Like