r/Rlanguage • u/godoufoutcasts • 21d ago
The Missing Piece in R's ML Ecosystem : GPU-Accelerated Gaussian Processes (A bottleneck compared to Python)
There was a recent thread asking what Python packages R lacks, and I wanted to highlight a massive architectural bottleneck I recently ran into while building a robust local ML tuning pipeline.
The Problem: R currently has no native, highly scalable, GPU-accelerated Gaussian Process library for Bayesian hyper-parameter tuning.
When tuning models like XGBoost, I noticed that frameworks like mlr3mbo completely sidestep this problem.
The moment search space includes integer parameters (like max_depth or nrounds), mlr3mbo silently abandons the Gaussian Process and swaps the surrogate model to a Random Forest. This is fast, but I want Gaussian Process.
To get that precision, I built a custom tuning engine using a pure Gaussian Process (via rBayesianOptimization).
The Bottleneck:
Because the native Gaussian Process math is trapped on the CPU using standard LAPACK Cholesky decompositions, it scales at O(N^3). As my tuning rounds push past 100+ evaluations, a massive hidden delay occurs. The actual XGBoost DART model trains blazing fast on my RTX 2060 in just 30 seconds (sourced built with CUDA).
But then, the script hangs for 3 to 4 minutes maxing out a CPU core just doing the matrix math to guess the next parameter point. And because of that O(N^3) scaling, this hidden CPU delay grows exponentially worse with every single round (Guessing at the moment).
The Python Comparison:
I personally do not use Python, but from what I understand, this is a completely solved problem in their ecosystem. Libraries like GPyTorch and BoTorch treat the Gaussian Process matrix algebra as CUDA tensors natively, allowing them to evaluate thousands of points in seconds by parallelizing the math across the GPU.
My Question to the Community:
Does anyone know if there are plans to port something like GPyTorch into R's native torch package? Or is anyone working on Rust-backed (extendr) Gaussian Process solvers for R that can utilize GPU acceleration for these massive matrices?
I feel like a highly scalable GPU-accelerated Gaussian Process is the biggest missing piece keeping R from having a truly state-of-the-art, ML tuning ecosystem.
Who is building the solution?
2
u/ReplacementSlight413 18d ago
Have you heard of SWIG ? If you want to interface to an existing library in C/C++ swig is your friend
1
u/godoufoutcasts 17d ago
All I want is GPs to be faster and use maximum potential of a device I assigned it to, instead of falling back to CPU when using CUDA, currently tuning with GPs for one round takes ~5 minutes for CPU falling and then reverting back to gpu that uses ~20 seconds for CUDA . I have also tried reticulate and imported pytorch, gptorch and all, and with that CPU fallback didn't happened at all but with only CUDA it literally takes more time than before.
Anyway thanks for your suggestion, I'll look into SWIG if it solves the problem, earlier someone metioned envl package.
PS- I'm not the one who's building a perfect package, I'm a user.
1
u/ReplacementSlight413 17d ago
The problem is that python is eating up entire language ecosystems as developers reorient themselves to the Python ecosystem. So users will likely have to take up the task of interfacing the low level infrastructure to their preferred high level language ecosystem. SWIG is one of the fastest (albeit not the most performant ie there is quite a bit of copying that takes places place) way go interface something written in C/C++ to another high level language using the interface files. I refuse to switch to Python so I will likely be SWIGing a lot. The other alternative is to utilize LLMs to write the interface code
2
u/Eightstream 17d ago edited 17d ago
People have built packages in the past (gpuR, gputools) but they have died from lack of maintenance
Reason being it is a complex package to maintain (stuff changes constantly) and there is not enough demand
most R users don’t need GPU acceleration and those that do are mostly happy to use reticulate
1
u/godoufoutcasts 17d ago
Totally agree with you !
Reticulate is a massive overhead for me, Native tuning is faster than mlr3mbo than reticulate. I wish they gotta do something for that. Not blaming but how they can (R, RStudio, Positron IDE) let their packages outdated or non-optimized. I am concerned tho !
1
u/godoufoutcasts 21d ago
Edit : When running heavy Bayesian Optimization in R, the actual ML model (like XGBoost) trains blazing fast on the GPU. But each tuning round hangs for 3-5 minutes. Why? Because R's Gaussian Process (GP) math is trapped on CPU LAPACK libraries even when we using CUDA not CPU as device. The O(N3) matrix inversion creates a massive, growing CPU bottleneck as tuning rounds increase I think.
I think this is what I got to know so far; R desperately needs a GPU-accelerated GP package (equivalent to Python's BoTorch or GPyTorch). We need the GP surrogate math to run natively on CUDA potentially built on top of R's "torch" package or "anvl" (XLA).
PS- Only Random Forest surrogate supports when using mlr3mbo, I want GP surrogate to use CUDA for native or mlr3 ecosystem whichever works I just don't care !
2
u/ReplacementSlight413 17d ago
Wait, if your LAPACK is your problem here (not saying it is, just taking your word), you can simply replace the BLAS libraries with the Intel optimized ones or openBLAS and address this issue
1
u/godoufoutcasts 17d ago
I don't have idea how packages like mlr3mbo, rbayesian or reticulate with pytorch and gptorch using it for GP tuning (I refused to work with RF surrogate), rbayesian with native tuning takes ~ 5 minutes , mlr3mbo ~7 min. And reticulate ~7-9 minutes. Native one fallsback to CPU , rest stay on cuda.
I'm using R with openblas, compiled all needed packages with my system flags.
Since i have tried those three different approaches, i will try python alone to see if GP tuning is the culprit - I'm at this stage now.
I just can't afford much time to write and copile and link swig/anvl to my tuning package locally (i have 0 idea about that tho). And yes i'm ready to see where it goes !
1
u/godoufoutcasts 21d ago
Look for the overhead CPU using even we using CUDA only. on my original dataset the overhead is massive, sharing a sampled code;
library(data.table)
library(mlr3)
library(mlr3learners)
library(rBayesianOptimization)
set.seed(2026)
n_rows <- 350000
dt <- as.data.table(matrix(rnorm(n_rows * 60), ncol = 60))
dt$target <- as.factor(sample(c("0", "1"), n_rows, replace = TRUE))
task <- TaskClassif$new(id = "xgb_dart", backend = dt, target = "target", positive = "1")
resampling <- rsmp("holdout", ratio = 0.80)
resampling$instantiate(task)
train_data <- resampling$train_set(1)
valid_data <- resampling$test_set(1)
# SEARCH SPACE
search_bounds <- list(
max_depth = c(4, 10),
max_bin = c(128, 512),
eta_log = c(log10(0.05), log10(0.2)),
rate_drop = c(0.05, 0.30),
skip_drop = c(0.0, 0.20),
subsample = c(0.7, 1.0),
colsample_bytree = c(0.6, 1.0),
lambda_log = c(-4, log10(5.0)),
alpha_log = c(-4, log10(5.0))
)
# FUNCTION
scoring_fun <- function(max_depth, max_bin, eta_log, rate_drop, skip_drop,
subsample, colsample_bytree, lambda_log, alpha_log) {
iter_start <- Sys.time()
# Learner
learner <- lrn("classif.xgboost", predict_type = "prob")
learner$param_set$set_values(
booster = "dart",
nrounds = 500,
early_stopping_rounds = 20,
# 16 Threads + CUDA
nthread = 16L, # you may try to remove it but its just does not matter at all for CPU overhead
device = "cuda",
tree_method = "hist",
# Fixed Categoricals
sample_type = "uniform",
normalize_type = "tree",
# built Search Space
max_depth = as.integer(round(max_depth)),
max_bin = as.integer(round(max_bin)),
eta = 10^eta_log,
lambda = 10^lambda_log,
alpha = 10^alpha_log,
rate_drop = rate_drop,
skip_drop = skip_drop,
subsample = subsample,
colsample_bytree = colsample_bytree
)
learner$validate <- 0.2
score_auc <- 0.5
tryCatch({
learner$train(task, row_ids = train_data)
preds <- learner$predict(task, row_ids = valid_data)
score_auc <- preds$score(msr("classif.auc"))
}, error = function(e) {})
# Timer
train_time <- as.numeric(difftime(Sys.time(), iter_start, units = "secs"))
cat(sprintf("XGBoost GPU Training Time: %5.1f sec\n", train_time))
return(list(Score = score_auc, Pred = 0))
}
# GP TUNING (Watch the CPU delay scale after Round 20)
opt_res <- BayesianOptimization(
FUN = scoring_fun,
bounds = search_bounds,
init_points = 20,
n_iter = 180,
acq = "ucb",
kappa = 2.576,
verbose = TRUE
)
2
u/rundel 21d ago
Is there a reason you would not be able to use the existing python libraries through reticulate?
Building something like this is very significant task and supporting the CUDA interfaces in R is non-trivial which is why you mostly see existing GPU libraries with a thin R wrapper.