slow buildMCMC()

[KK Sasa]

Hi Nimble team,

The complied R code amazingly runs fast :), but the building and compiling time is  daunting :(.

An example is below.

1. Any tips to shorten the time?
2.  I did control of monitor, but "gamma" is still monitored. Is it a bug?

library(nimble)
code <- quote(
{
    for(i in 1:I) {
        d[i] ~ dnorm(0, sd = 100)
    }
    for(i in 1:block){
        for(c in 2:pairedItem){
            weight[i,c-1] ~ dunif(0,1)
            tau[i,c-1] ~ dnorm(0, sd = 100)           
        }       
    }   
    s ~ dunif(0,100)
    for(n in 1:N) {
        for(i in 1:I) {           
            nu[n,i,1] <- 1
            for(c in 2:numP){
                nu[n,i,c] <- exp(theta[n] - d[i])
            }           
            de[n,i] <- sum(nu[n,i,1:numP])
            for(c in 1:numP){
                p[n,i,c] <- nu[n,i,c]/de[n,i]
            }           
            pcm[n,i] ~ dcat(p[n,i,1:numP])
        }
        theta[n] ~ dnorm(0,sd = s)
        for(i in 1:block){           
            nu2[n,i,1] <- 1
            for(c in 2:pairedItem){
                nu2[n,i,c] <- exp(weight[i,c-1]*theta[n]+(1-weight[i,c-1])*gamma[n] + tau[i,c-1])
            }
            de2[n,i] <- sum(nu2[n,i,1:pairedItem])
            for(c in 1:pairedItem){
                pr[n,i,c] <- nu2[n,i,c]/de2[n,i]
            }                   
            nrm[n,i] ~ dcat(pr[n,i,1:pairedItem])
        }       
        gamma[n] ~ dnorm(0,sd = 1)       
    }
}   
)   

load("data.R")
constants <- list(N = 1000, I = 10,numP=2,block=5,pairedItem=2)
data <- list(pcm = data$pcm,nrm=data$nrm)
inits <- list(d=rep(0,10))

ptm <- proc.time()
Rmodel <- nimbleModel(code=code, name='Rmodel',constants=constants, inits=inits, check = FALSE)
Rmodel$resetData() # remove data
Rmodel$setData(data) # put data in
RmodelSpec <- configureMCMC(Rmodel, print=TRUE)
RmodelSpec$addMonitors(c('d','s','tau','weight'))
Rmcmc <- buildMCMC(RmodelSpec) # uncomplied R code
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
Cmcmc$run(2000)
samples <- as.matrix(Cmcmc$mvSamples)
require("coda")
burnin = 1000
samples = samples[-(1:burnin),]
coda_samples <- mcmc(samples)
summary(coda_samples)
print(proc.time() - ptm)

[Perry de Valpine]

Dear KK,

Thanks.  I’m glad to see this use case.

Unfortunately for large models we do see slow building and compiling times.  It’s on our agenda to improve, but here are some suggestions for what you can do.

Make sure you have the latest version of NIMBLE (building and compiling is still slow, but faster than it once was).

You can replace your two lines

  RmodelSpec <- configureMCMC(Rmodel, print=TRUE)

  RmodelSpec$addMonitors(c('d','s','tau','weight'))

with

  RmodelSpec <- configureMCMC(Rmodel, print=TRUE, useConjugacy = FALSE, monitors = c('d','s','tau','weight'))

The conjugacy checking is a slow part of MCMC configuration, and your model has no conjugacy anyway, so skipping it (useConjugacy = FALSE) helps speed, although that’s only a small bit of the total time.  You can set the monitors as you like in the call to configureMCMC.  It also works to directly modify RmodelSpec$monitors, which is simply a character vector of variable names to monitor.

The most important factor in building and compiling time is the number of nodes in the model.  One way to reduce this is to use nimbleFunctions in your BUGS model to skip over intermediate variables that don’t really need to be in the model.  For example, it appears to me that “nu" and “de” really only exist as intermediate steps to calculate elements of p.  In that case, you could replace  

for(i in 1:I) {            
            nu[n,i,1] <- 1
            for(c in 2:numP){
                nu[n,i,c] <- exp(theta[n] - d[i])
            }            
            de[n,i] <- sum(nu[n,i,1:numP])
            for(c in 1:numP){
                p[n,i,c] <- nu[n,i,c]/de[n,i]
            }            
            pcm[n,i] ~ dcat(p[n,i,1:numP])

        }

by something like this 

for(i in 1:I) {

            p[n, i, 1:numP] <- calculateP(theta = theta[n], d = d[i])           

            pcm[n,i] ~ dcat(p[n,i,1:numP])

        }

where you provide the nimbleFunction

calculateP <- nimbleFunction(

  run = function(theta = double(0), d = double(0)) {

    returnType(double(1))

    expthetaminusd <- exp(theta-d)

    denom <- 1 + expthetaminusd ## did I follow that correctly?

    declare(ans, double(1))

    setSize(ans, 2)

    ans[1] <- 1/denom

    ans[2] <- expthetaminusd/denom

    return(ans)

 }

)

(Or you could use this kind of idea for each scalar element of p, rather than for the vector p[n, i, 1:numP])

By cutting nu and de out of the model, it looks like you would remove 30000 nodes.  By calculating p with the last index vectorized, you would be cutting p from having 20000 scalar nodes to having 10000 vector nodes, each of length 2.  A similar approach could be used for pr, I think.

I haven’t tested this with your code, and you’ll have to validate if I converted the calculations correctly in the calculateP example above, but I’m just trying to illustrate the idea.

Please let us know if this helps.

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/9edfece0-0b6c-4dff-b852-f168c5c429fe%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<data.R>