Question about parallelization in BayesianTools

[florianhartig]

Question from a user

I tried the example with VSEM model from the BayesianTools manual with some minor edit changes (the code is attached). However, it seems to me that although the correct number of R process started (in Task manager), the other cores seem to be not busy at all (~zero % activity). I am trying that on the server with 60 cores. I tried with VSEM model as well as with more complex and time consuming model, however in all cases the cores seems to be not active.

I am using Windows 10 with R 3.5.0., however, I tried it with same result on Centos 7.

# Bayesian calibration

# Example from https://cran.r-project.org/web/packages/BayesianTools/BayesianTools.pdf

rm(list=ls())

library(BayesianTools)

## Example of how to use parallelization using the VSEM model

# Note that the parallelization produces overhead and is not always
# speeding things up. In this example, due to the small
# computational cost of the VSEM the parallelization is
# most likely to reduce the speed of the sampler.

# Creating reference data
PAR <- VSEMcreatePAR(1:1000)
refPars <- VSEMgetDefaults()
refPars[12,] <- c(0.2, 0.001, 1)
rownames(refPars)[12] <- "error-sd"
referenceData <- VSEM(refPars$best[1:11], PAR)
obs = apply(referenceData, 2, function(x) x + rnorm(length(x),
sd = abs(x) * refPars$best[12]))

# Selecting parameters
parSel = c(1:6, 12)

## Builidng the likelihood function
likelihood <- function(par, sum = TRUE){
x = refPars$best
x[parSel] = par
predicted <- VSEM(x[1:11], PAR)
diff = c(predicted[,1:3] - obs[,1:3])
llValues = dnorm(diff, sd = max(abs(c(predicted[,1:3])),0.0001) * x[12], log = TRUE)
if (sum == FALSE) return(llValues)
else return(sum(llValues))
}
# Prior
prior <- createUniformPrior(lower = refPars$lower[parSel], upper = refPars$upper[parSel])

####
## Definition of the packages and objects that are exported to the cluster.
# These are the objects that are used in the likelihood function.
opts <- list(packages = list("BayesianTools"), variables = list("refPars", "obs", "PAR", "parSel" ),
dlls = NULL)

# Create Bayesian Setup
BSVSEM <- createBayesianSetup(likelihood, prior, best = refPars$best[parSel],
names = rownames(refPars)[parSel], parallel = 59, parallelOptions = opts)

## The bayesianSetup can now be used in the runMCMC function.

# settings for the sampler, iterations should be increased for real applicatoin
settings <- list(iterations = 100000, nrChains = 2)

# Run Markov chain Monte Carlo with Differential-Evolution zs sampler
out <- runMCMC(bayesianSetup = BSVSEM, sampler = "DEzs", settings = settings)

# Remove the Bayesian Setup and close the cluster
stopParallel(BSVSEM)
rm(BSVSEM)

[florianhartig]

The general issue here is that MCMCs are not arbitrarily parallelizable. What you can do is to parallelize the independent MCMC chains, and you can usually partly parallelize samplers, but it's usually not efficient to use an arbitrary large number of cores.

For the case of BT, we currently only parallelize single MCMC samplers, but not the chains. How much speedup you can get per sampler depends on the sampler settings. Metropolis samplers can only use 1 core. The DE sampler family can use as many cores as you have internal chains. The default for the DEzs is 3 internal chains, thus you request 59 cores, but the sampler will only use 3 and the rest are idle.

We are working on parallelizing the chains as well, but with the current settings what you should do to get maximum parallelization is the following:

• start the chains as independent processes, see https://github.com/florianhartig/BayesianTools/issues/9

• for each process, request only as many cores as your MCMC can use. It might speed up the DEzs when you use more internal chains (see help of the sampler), but that might not always be the case (communication overhead, waiting for processes to finish). You have to try.