Reducing loo memory utilisation

[mattmoo]

I'm running out of memory to perform LOO-CV using loo.

From the documentation, I gather that defining loo.function might lower memory utilisation.

Some of it seems straightforward, but there is maybe some difference in terminology where I can't quite figure out what 'draws' is.

Would there be a general function for mcp.fit objects compatible with loo.function? (i.e., such that loo.mat and loo.fn are essentially identical in the example below.

set.seed(024)

library(mcp)
options(mc.cores = 3)  # Speed up sampling
set.seed(42)  # Make the script deterministic

model = list(
  recalled | trials(items) ~ 1,
  ~ 0 + items
)
prior = list(
  int_1 = "dnorm(2, 1)",  # high recall for easy trials
  cp_1 = "dnorm(4, 1)",  # performance dicontinuity around 4
  items_2 = "dnorm(-0.4, 1) T( , -0.2)"  # performance deteriorates considerably
)


# Get just the model - no fit
empty = mcp(model, family = binomial(), sample = FALSE)

# Simulate data
df = data.frame(items = rep(1:9, each = 40))
df$recalled = empty$simulate(
  df$items,
  cp_1 = 4,
  int_1 = 2.5,
  items_2 = -0.4,
  sigma = 0.5)

fit_info = mcp(model, data = df, prior = prior, family = binomial(), sample = "both", iter = 10000)

# Do LOO using conventional matrix.
loo.mat = loo::loo(fit_info)


#  A function that takes arguments i, data, and draws and returns a vector
#  containing the log-likelihood for the ith observation evaluated at each
#  posterior draw.
llfun <- function(i, data, draws) {
  return(NULL)
}

# The args argument must also be specified and should be a named list with the
# following components.
args = nlist(
  # An object containing data (e.g. observed outcome and predictors) needed to
  # compute the pointwise log-likelihood. data should be in an appropriate form
  # so that (i=i, data=data[i,,drop=FALSE], draws=draws) returns the S-vector of
  # log-likelihoods for the ith observation.
  data = NULL,
  # An object containing the posterior draws for any parameters needed to
  # compute the pointwise log-likelihood.
  draws = NULL,
  N = dim(fit_info$mcmc_loglik[[1]])[2],
  S = dim(fit_info$mcmc_loglik[[1]])[1])

# Do LOO using loo.function.
loo.fn = loo::loo(llfun,
                  args = args,
                  cores = 1)

print(loo.mat)
print(loo.fn)

[lindeloev]

This is timely; the next version of mcp did already switch to post-sampling computation of log-likelihoods. So it should be possible to do it quite efficiently. I'll look into it.

But maybe the problem is your memory limit. For this fit, my laptop peaked at 2 GB of RAM usage over and above base R/Rstudio. Running with 3xdata, it maxed at 6GB, including everything. If you're on Windows, have you checked memory.limit()? You can increase it using, e.g., memory.limit(2^13) for 8GB. There are other ways on Mac and Linux, I think.

[lindeloev]

OK, I'll try and get this into the next release! If all goes will, I hope to submit to CRAN in a few weaks.

[mattmoo]

Thanks!

[mattmoo]

Somewhat related to this, is there any thinning of chains in mcp?

[lindeloev]

mcp doesn't have thinning because it just throws away information. But I can see a use case for memory-limited operations, so will consider adding it. stan has a similar approach. Perhaps I'll just add an nsamples argument to loo.mcpfit.

[mattmoo]

That could be an option, in the interim, would it be as simple as removing however many rows from the mcmc tables in the mcpfit object? Or maybe there are other fields that need adjusting?

[lindeloev]

Yes, try removing rows from fit$mcmc_loglik. I don't think loo.mcpfit takes any other input in mcp v0.3. (will from v0.4 onwards)

…

On Mon, Mar 15, 2021, 22:43 mattmoo ***@***.***> wrote: That could be an option, in the interim, would it be as simple as removing however many rows from the mcmc tables in the mcpfit object? Or maybe there are other fields that need adjusting? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#109 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAK7NOJAD6AA4K4CZHW5U7LTDZ5OVANCNFSM4ZBHRQ7Q> .

[mattmoo]

I have used these functions, which are butchered versions of yours (pp_eval, loo.mcpfit), without much generalisability.

llfun = function(data_i, draws = NA, fit)  {
  # Not used
  draws = NA
  
  loglik = mcp:::tidy_samples(
    fit,
    population = TRUE,
    varying = FALSE,
    prior = FALSE,
    nsamples = NULL
  ) %>%
    tidyr::expand_grid(mcp:::add_rhs_predictors(newdata = data_i %>%
                                                  dplyr::select(-!!fit$pars$weights),
                                                fit)) %>%
    dplyr::mutate(
      loglik := rlang::exec(
        mcp:::simulate_vectorized,
        fit,!!!.,
        .type = 'loglik',
        .rate = TRUE,
        .which_y = 'mu',
        .arma = TRUE,
        .scale = 'response'
      )
    ) %>%
    dplyr::select(loglik) %>%
    as.matrix()
  
  return(loglik)
}



calculate.loo = function(fit, n.cores = 15) {
  
  chain_id = rep(seq_along(fit$mcmc_post), 
                 each = nrow(fit$mcmc_post[[1]]))
  r_eff = loo::relative_eff(x = llfun,
                            chain_id = chain_id,
                            cores = n.cores,
                            data = fit$data,
                            draws = NA,
                            fit = fit)
  
  loo.result = loo::loo(x = llfun, 
                        data = fit$data,
                        fit = fit,
                        draws = NA,
                        r_eff = r_eff,
                        cores = n.cores)
  
  return(loo.result)
}

[lindeloev]

I just committed support for pointwise loo to the v0.4 branch. @mattmoo, could you try it out and see if it works as expected? loo(fit, pointwise = TRUE, nsamples = 1000). For AR(N) models there will be a very slight deviation with pointwise because the log-likelihood is evaluated differently for the first N data points.

In addition to being slow in it's own right, mcp adds an overhead of ~0.2 seconds for each row of data. This can be optimized down to basically zero upon popular demand. But for now, I prioritize code simplicity.

[mattmoo]

Just using pointwise = TRUE, it still seems to fail, giving the error:

Error in vec_slice(x, idx) : 
  long vectors not supported yet: ../../src/include/Rinlinedfuns.h:535
Error in system(paste(which, shQuote(names[i])), intern = TRUE, ignore.stderr = TRUE) : 
  cannot popen '/usr/bin/which 'pdflatex' 2>/dev/null', probable reason 'Cannot allocate memory'
Error in system(paste(which, shQuote(names[i])), intern = TRUE, ignore.stderr = TRUE) : 
  cannot popen '/usr/bin/which 'pdflatex' 2>/dev/null', probable reason 'Cannot allocate memory'
Error in system(paste(which, shQuote(names[i])), intern = TRUE, ignore.stderr = TRUE) : 
  cannot popen '/usr/bin/which 'pdflatex' 2>/dev/null', probable reason 'Cannot allocate memory'
Error in system(paste(which, shQuote(names[i])), intern = TRUE, ignore.stderr = TRUE) : 
  cannot popen '/usr/bin/which 'pdflatex' 2>/dev/null', probable reason 'Cannot allocate memory'

Also, just for a sanity check, might be a good idea to increment the package version in the v0.4 branch. I assume that if I was actually using 0.3.0.9000, it would hit me with a warning about the pointwise argument.

[lindeloev]

Did your "butchered" version did run OK? The only real difference I can see is that you used 15 cores.

[mattmoo]

Yes, my version runs fine. I would have thought that more cores would tend to exacerbate any memory shortages, but I must say that I don't have a deep understanding of the parallel libraries!

The current v0.4 seems to behave pretty similarly to what I recall from v0.3.0.9000. I do see pointwise in ?mcp::loo, so I assume that I'm using the right version (?). Below is system utilisation over about 10 minutes before the failure on vec_slice().

Utilisation for my version looks is stable at a much lower level (again 10 minutes, it does take about 3 hours total)

[lindeloev]

Cool, thanks! Yes, more core should exacerbate the problem. I'll check my implementation later; it's likely a simple error somewhere. It's great that you got this working, so I know that there is light at the end of the tunnel :-)

How many rows of data do you have?

[mattmoo]

This is for 57,077 rows of Bernoulli data, 250,000 iterations in 10 chains.