Unusually slow GPU (OpenCL) vs CPU sampling times in blavaan/cmdstanr on MareNostrum 5

Interfaces CmdStan
, ,
1

I am encountering a severe performance degradation when using the GPU backend compared to the CPU in Bayesian SEM simulations via blavaan and cmdstanr. The GPU OpenCL sampling times are significantly slower than the CPU across almost all conditions. For small models (N=200), the GPU is up to 10x slower. The performance gap narrows as complexity increases, but even at N=5000 with 20 indicators, the GPU merely breaks even with the CPU rather than outperforming it. This contradicts my prior experience on other HPC clusters.

Code to simulate data and run the model

R

library(blavaan)
library(cmdstanr)

# Simulation parameters
nobs <- 200
nind <- 3
ncat <- 2

# Data generation function
makeord <- function(Data, vars = NULL, ncat = 2) {
  if (length(vars) == 0) vars <- 1:NCOL(Data)
  Data <- rbind(Data, NA)
  Data[nrow(Data), vars] <- ncat
  Data[, vars] <- apply(Data[, vars, drop = FALSE], 2, function(x) {
    nc <- tail(x, 1)
    tmpp <- (1 / nc) + runif(1, -.1 / nc, .1 / nc)
    brks <- c(
      min(x, na.rm = TRUE) - .1,
      seq(
        quantile(x, tmpp, na.rm = TRUE),
        quantile(x, 1 - tmpp, na.rm = TRUE),
        length.out = (nc - 1)
      ),
      max(x, na.rm = TRUE) + .1
    )
    cut(x, breaks = brks, labels = FALSE)
  })
  Data[-nrow(Data), ]
}

# Generate data
HS.model <- paste0(
  ' visual  =~ ', paste0('x', 1:nind, collapse = ' + '), "\n",
  ' textual =~ ', paste0('x', (nind + 1):(2 * nind), collapse = ' + '), "\n",
  ' speed   =~ ', paste0('x', (2 * nind + 1):(3 * nind), collapse = ' + ')
)

set.seed(123)
dcont <- lavaan::simulateData(HS.model, sample.nobs = nobs)
Data <- makeord(dcont, ncat = ncat)

# 1. Generate syntax via blavaan
export_dir <- "lavExport"
bfit2 <- bcfa(
  model = HS.model, data = Data, ordered = TRUE, do.fit = FALSE, 
  mcmcfile = export_dir, test = 'none', bcontrol = list(backend = "cmdstanr")
)

load(file.path(export_dir, 'semstan.rda'))

# 2. CPU Execution
mod_cpu <- cmdstan_model(
  stan_file = file.path(export_dir, 'sem.stan'),
  cpp_options = list(stan_threads = TRUE), stanc_options = list("O1")
)

rjarg_cpu <- list(
  data = stantrans$data, init = stantrans$inits, iter_warmup = 1000, 
  iter_sampling = 2000, chains = 3, parallel_chains = 3, threads_per_chain = 1
)
res_cpu <- do.call(mod_cpu$sample, rjarg_cpu)

# 3. GPU Execution (OpenCL)
mod_gpu <- cmdstan_model(
  stan_file = file.path(export_dir, 'sem.stan'),
  cpp_options = list(stan_opencl = TRUE), stanc_options = list("O1")
)

rjarg_gpu <- list(
  data = stantrans$data, init = stantrans$inits, iter_warmup = 1000, 
  iter_sampling = 2000, chains = 3, parallel_chains = 3, opencl_ids = c(0, 0)
)
res_gpu <- do.call(mod_gpu$sample, rjarg_gpu)

Benchmark Results Summary:

(Note: Values represent the mean pure sampling time in seconds across simulation replicates. Some CPU conditions at N=5000 are marked N/A as those jobs had not yet completed in this batch).

N Categories Indicators Mean GPU OpenCL Time (s) Mean CPU Time (s)
200 2 3 2264.50 173.42
200 2 10 1424.88 504.31
200 2 20 3977.25 3596.15
200 5 3 987.98 89.79
200 5 10 1133.79 503.65
200 5 20 3354.69 2474.70
200 7 3 772.99 81.00
200 7 10 1294.91 331.78
200 7 20 3489.41 2147.39
1000 2 3 1453.45 423.86
1000 2 10 2752.43 1878.29
1000 2 20 7764.20 7078.44
1000 5 3 1367.43 430.94
1000 5 10 2830.29 1942.51
1000 5 20 7834.80 6446.41
1000 7 3 1231.51 402.08
1000 7 10 2889.64 1950.95
1000 7 20 7714.20 6526.29
5000 2 3 3585.96 2669.57
5000 2 10 18808.70 17080.56
5000 2 20 45143.42 45608.28
5000 5 3 3337.66 N/A
5000 5 10 16497.83 N/A
5000 5 20 43381.96 N/A
5000 7 3 3436.32 N/A
5000 7 10 14119.93 N/A
5000 7 20 43013.03 N/A

Operating System: Linux (MareNostrum 5 HPC Cluster)

CmdStan Version: 2.35.0

Compiler/Toolkit: GCC (R/4.5.0-gcc, openblas/0.3.27-gcc), CUDA

Hardware Specifications:

  • CPU: Intel Sapphire Rapids (Platinum 8480+)

  • GPU: NVIDIA Hopper H100

Could this bottleneck on H100 GPUs be caused by data transfer overheads for these specific dimension sizes, the way blavaan generates the .stan code, or a known OpenCL driver configuration issue with this specific hardware architecture?

Thanks,

Sinval

2

I have limited experience trying to run cmdstan via OpenCL but there are a few key details that limit which models benefit from the GPU.

Importantly, cmdstan can only offload certain parts of the execution to the GPU. These are mostly matrix operations. So if you have very large matrix multiplications or your model uses _glm functions, these are good candidates for acceleration. If your stan code lacks these or the matrix operations aren’t sufficiently large, then you won’t see any benefits. Now because only certain operations can execute on the GPU, data will have to move back and forth between the CPU and GPU, this introduces overhead that can result in pretty noticeable slowdowns when the computations being done on the GPU aren’t worth the overhead.

I suspect that’s what happening here. Looking at the Stan code that blavaan uses, it does not very look opencl friendly to me (although it’s somewhat complex so I could be missing something). At low N in your testing, you see major slowdowns, which is likely due to GPU operations giving no performance benefit at a low scale and the GPU-transfer overhead dominates leading to slowdowns. At high N, you get closer to parity, my guess is that the non-GPU parts of the process are actually dominating runtime here and the previously seen overhead (and any differences in perf on the GPU) get washed out by other more expensive parts of the execution.

It may be the case that you could speed this up a bit, but you would likely need to write some custom Stan code to make it as GPU friendly as possible. Not sure how it would pan out, but that may be an option.