florianhartig / BayesianTools

General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics
https://cran.r-project.org/web/packages/BayesianTools/index.html
117 stars 29 forks source link

Possible to thin a BayesianOutput object prior to restart / -> memory problems in z-matrix #179

Open woodwards opened 5 years ago

woodwards commented 5 years ago

I wrote a loop around runMCMC to check convergence, but I am having trouble with memory errors when the BayesianOutput object grows too big after multiple calls. Is it possible to thin a BayesianOutput object prior to restart? Or some other way to avoid memory overflow? Thanks.

# run BT until stopping conditions met (these can be changed in the file BC_BASGRA_BT_stop.csv)
bt_chains <- nInternal * nChains 
bt_pars <- length(bt_names)
bt_conv <- NA
repeat{

  # run Bayesian Tools
  if (is.na(bt_conv)){ 
    # first run
    bt_out <- runMCMC(bayesianSetup = bt_setup, 
                      sampler = "DREAMzs", 
                      settings = bt_settings)
    cat(file=stderr(), " ", "\n")
  } else {
    # continuation
    if (nBurnin==0){
      cat(file=stderr(), paste("Greater than", conv_target, "so continuing..."), "\n")
    } else {
      cat(file=stderr(), "Restart doesn't work with nBurnin>0 so stopping...\n")
      # stop()
    }
    bt_out <- runMCMC(bayesianSetup = bt_out)
    cat(file=stderr(), " ", "\n")
  }

  # assess convergence  
  bt_length <- dim(bt_out[[1]]$chain[[1]])[[1]]
  cat(file=stderr(), paste("Total chains =", bt_chains), "\n")
  cat(file=stderr(), paste("Total samples per chain =", bt_length), "\n")
  bt_time <- sum(sapply(bt_out, function(x) x$settings$runtime[3]/60, simplify=TRUE))
  cat(file=stderr(), paste("Total time =", round(bt_time,2), "minutes"), "\n")
  bt_conv <- max(gelmanDiagnostics(bt_out)$psrf[,1])
  cat(file=stderr(), paste("Convergence max(psf) =", round(bt_conv,3)), "\n")

  # read stopping criterion from file (allows us to change it on the fly)
  conv_criteria <- read.csv("scripts/BC_BASGRA_BT_stop.csv", header=FALSE) 
  conv_target <- conv_criteria[1,1]
  conv_minutes <- conv_criteria[2,1]
  if ((bt_conv <= conv_target) || (bt_time >= conv_minutes)){
    break
  }
}
florianhartig commented 5 years ago

Is it necessary to thin after the sampler has run? Why don't you thin the sampler directly (option thin?)

woodwards commented 5 years ago

Every time I call runMCMC it adds another 10000 samples to my chain. The chain grows and eventually crashes R due to a memory problem.

I thought of solving the problem like this: when I get to 20000 and not yet converged, I throw away half of the old chain (either by truncation or thinning or some other way) before I call runMCMC again. So the chain will never grow bigger than 20000.

Thinning in runMCMC will thin the current samples., not the old ones, I think?

florianhartig commented 5 years ago

Yes, but if you thin them from the start, you shouldn't have the memory problems you describe.

I understand that an adaptive thinning might be a bit more convenient, but it will require a lot of changes on our side. Let's say you just set thin = 20 from the start. Wouldn't that already solve the problem?

marjoleinbruijning commented 3 years ago

Hi Florian,

Thanks for the great package!

I have a similar memory problem when fitting a large model with about 100 parameters, where I use a loop to continue adding samples until convergence is reached: eventually R crashes when the object gets too large.

I tried to increase thinning to solve the problem, but the issue remains, even when I set thin=100. I looked at the output object size as a function of thinning, and the effect of thinning on the object size seems to diminish at higher thinning values. For example, when thin=10, object size decreases by 56% compared to no thinning, while when thin=100, object size decreases by 62% compared to no thinning.

Is there any workaround here?

Many thanks,

Marjolein

florianhartig commented 3 years ago

Hi Marjolein,

100 parameters is usually not large enough to cause memory problems. How many iterations are you doing here, and what's the size of your MCMC sampler object? Can you confirm that thinning works as expected, i.e. that you get a thinned chain with only every 100th iteration recorded?

F

marjoleinbruijning commented 3 years ago

Thanks for the quick response.

I've been trying varying number of iterations, but around 1E5 iterations per run (thinning varying between 1-100).

As a quick example with only 1E4 iterations, this is the code I used (thinning works as expected):

settings <- list(iterations = 1E4, nrChains = 1, message=TRUE, thin=10)
out10 <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs",
                 settings = settings)
format(object.size(out10),units='Mb')
[1] "13.4 Mb"

nrow(getSample(out10))
[1] 1002
marjoleinbruijning commented 3 years ago

I have been looking into this a bit more, and maybe there is an issue with the thinning.

The number of rows in out$Z appears to equal the total number of iterations (without thinning), so this matrix can become very large. In one of my saved out objects, before it crashed, this matrix has 3000,000 rows, and makes up 90% of the total object size. (this while the output from getSample(out) has only 160,000 rows) So I guess this explains the large object size and the memory problems.

Is there something I can do to resolve this?

Thanks!

florianhartig commented 3 years ago

Hi @marjoleinbruijning,

OK, I see, this is the culprit. I have been meaning to look at the z matrix for a while, because I think in general it would be both better to thin it and to implement a bit of memory loss. I think this was https://github.com/florianhartig/BayesianTools/issues/8

I will try to implement this. For the moment, what I would recommend the following strategy between your updates:

## Generate a test likelihood function. 
ll <- generateTestDensityMultiNormal(sigma = "no correlation")

## Create a BayesianSetup object from the likelihood 
## is the recommended way of using the runMCMC() function.
bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

## Finally we can run the sampler and have a look
settings = list(iterations = 100000, thin = 10)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)

lenZ = nrow(out$Z)
sel = sample.int(lenZ, min(5000, lenZ),replace = T, prob = exp(-(lenZ:1)/lenZ))
out$Z = out$Z[sel,]

out2 = runMCMC(out)

If you do very frequent sampler updates, I might think about removing the prob option (which is the memory loss), because even so, you will get memory loss as the earlier updates will be sampled more often.

marjoleinbruijning commented 3 years ago

Thanks a lot, I'll try that!

Op vr 18 jun. 2021 om 10:00 schreef Florian Hartig @.***

:

Hi @marjoleinbruijning https://github.com/marjoleinbruijning,

OK, I see, this is the culprit. I have been meaning to look at the z matrix for a while, because I think in general it would be both better to thin it and to implement a bit of memory loss. I think this was #8 https://github.com/florianhartig/BayesianTools/issues/8

I will try to implement this. For the moment, what I would recommend the following strategy between your updates:

Generate a test likelihood function.

ll <- generateTestDensityMultiNormal(sigma = "no correlation")

Create a BayesianSetup object from the likelihood

is the recommended way of using the runMCMC() function.

bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

Finally we can run the sampler and have a look

settings = list(iterations = 100000, thin = 10) out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)

lenZ = nrow(out$Z) sel = sample.int(lenZ, min(5000, lenZ),replace = T, prob = exp(-(lenZ:1)/lenZ)) out$Z = out$Z[sel,]

out2 = runMCMC(out)

If you do very frequent sampler updates, I might think about removing the prob option (which is the memory loss), because even so, you will get memory loss as the earlier updates will be sampled more often.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/florianhartig/BayesianTools/issues/179#issuecomment-863841562, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABSPCRD3N7EY3IUEXPSQEFLTTL4JBANCNFSM4GX7HRMA .