Running MCMC chains on multiple cores.

[Sarah Power]

I was wondering if it is possible to run MCMC chains on more than one core per chain. I tried:

options(mc.cores = parallel::detectCores())

 I am updated to R version 4.3.2. Sean Anderson has given me these reproducible examples:

library(RTMB)
set.seed(123)
data <- list(Y = rnorm(10) + 1:10, x=1:10)
parameters <- list(a=0, b=0, logSigma=0)

f <- function(parms) {
Y <- data$Y
x <- data$x
a <- parms$a
b <- parms$b
logSigma <- parms$logSigma
ADREPORT(exp(2*logSigma));
nll = -sum(dnorm(Y, a+b*x, exp(logSigma), TRUE))
nll
}
obj <- MakeADFun(f, parameters)

library(tmbstan)
s <- tmbstan(obj, iter = 100, chains = 2, cores = 1) # works
s <- tmbstan(obj, iter = 100, chains = 2, cores = 2) # doesn’t work

  s <- tmbstan(obj, iter = 100, chains = 2, cores = 3) # doesn’t work

Also this ran for Sean, but not for me:
TMB::runExample("simple")
fit <- tmbstan(obj, chains=2, iter = 50, cores = 2) # works

[jcavieresg]

Hi!

Well, it is weird because I have the same problem running your code. My output from that is:

[[1]] Stan model 'tmb_generic' does not contain samples. [[2]]

Stan model 'tmb_generic' does not contain samples.

However, using TMB (not RTMB) and then running tmbstan I don't have any problem to use the cores of the computer. Maybe (I'm thinking aloud), it could be related to the priors used as default on RTMB (it does?). Since you are not considering priors for your parameters in your function, I assume that some default prior could be causing the problem in the sampling.  

[Kasper Kristensen]

Sarah, I can reproduce your RTMB example on a Windows machine, however the 'simple' TMB example seems to work fine (?). Note, that parallelization works fundamentally different on Windows than other OSes, so the following diagnostics are Windows specific.

1) First problem is that library(RTMB) is not automatically loaded on the nodes. We can fix this by adding 'require(RTMB)' to the body of the function.
2) Second problem is that 'data' is not exported to the nodes. We can fix this using closures, i.e. by explicitly assigning data to evironment(f).

Ideally, these problems should be fixed in the tmbstan package, but for now one must do something like this:

library(RTMB)
set.seed(123)
data <- list(Y = rnorm(10) + 1:10, x=1:10)
parameters <- list(a=0, b=0, logSigma=0)
f <- function(parms) {

  require(RTMB)

  Y <- data$Y
  x <- data$x
  a <- parms$a
  b <- parms$b
  logSigma <- parms$logSigma
  ADREPORT(exp(2*logSigma));
  nll = -sum(dnorm(Y, a+b*x, exp(logSigma), TRUE))
  nll
}

environment(f) <- new.env()
environment(f)$data <- data

obj <- MakeADFun(f, parameters)
library(tmbstan)

s <- tmbstan(obj, iter = 100, chains = 2, cores = 2) ## Works !

[Sarah Power]

Can we have more than one chain per core? Or multithreads per chain? Sean passed this along to me: Reduce Sum: A Minimal Example (mc-stan.org)

[Kasper Kristensen]

Yes, multiple threads per chain has worked in TMB for a long time via the 'parallel_accumulator'. In RTMB the parallel split can be fully automated. However, currently the parallel version of RTMB is shipped under a different name 'RTMBp'.

You'll need a fairly large dataset to get a significant speedup, something like:

library(RTMBp)
n <- 1e6
set.seed(123)
x <- rnorm(n)
data <- list(Y = rnorm(n) + x, x=x)

parameters <- list(a=0, b=0, logSigma=0)
f <- function(parms) {
  Y <- data$Y
  x <- data$x
  a <- parms$a
  b <- parms$b
  logSigma <- parms$logSigma
  ADREPORT(exp(2*logSigma));
  nll = -sum(dnorm(Y, a+b*x, exp(logSigma), TRUE))
  nll
}

obj <- MakeADFun(f, parameters) ## NOTE 'Autopar work split' info on how much work is assigned to each thread
library(tmbstan)
system.time(s <- tmbstan(obj, iter = 1e2, chains = 1,seed=1))