A native R package for building, training, and deploying neural networks. Backed by the ggml C library, designed primarily for Vulkan GPU acceleration with full CPU fallback — no Python, no TensorFlow, everything runs inside your R session.
GPU-first design: when a Vulkan-capable GPU is available (NVIDIA, AMD, Intel, ARM Mali, Qualcomm Adreno), all operations run on GPU automatically. On machines without a GPU the package falls back to CPU transparently — no code changes needed.
Two complementary APIs:
| API | Style | When to use |
|---|---|---|
| Sequential / Functional | Keras-like, static graph | Production models, CRAN-standard workflow |
Dynamic autograd (ag_*) |
PyTorch-like, eager | Research, custom architectures, Transformers |
Also serves as the backend engine for llamaR (LLM inference) and sd2R (Stable Diffusion).
install.packages("ggmlR")GPU (Vulkan) support is auto-detected at build time.
Ubuntu / Debian — to enable GPU:
sudo apt install libvulkan-dev glslcWindows — install Rtools and optionally the Vulkan SDK for GPU support.
Force-enable or disable Vulkan GPU backend:
install.packages("ggmlR", configure.args = "--with-vulkan")
install.packages("ggmlR", configure.args = "--without-vulkan")Enable CPU SIMD acceleration (AVX2, SSE4, etc.) for faster inference on your machine:
install.packages("ggmlR", configure.args = "--with-simd")Enable the hard-exit path used by multi-GPU standalone scripts
(ggml_vulkan_shutdown(hard = TRUE); see Clean shutdown below and
vignette("multi-gpu")). Off by default — the released package must not call
_exit(), as CRAN policy forbids a package terminating the R session:
install.packages("ggmlR", configure.args = "--enable-hard-exit")Options can be combined:
install.packages("ggmlR", configure.args = "--with-vulkan --with-simd --enable-hard-exit")Full step-by-step setup on Ubuntu/Debian, from a clean system to a working GPU build.
1. R and the Vulkan loader/tools:
sudo apt install -y r-base
sudo apt install vulkan-tools libvulkan-dev2. The glslc shader compiler (needed to build ggmlR's Vulkan backend):
# Ubuntu 24.04 (Noble)
sudo add-apt-repository universe
sudo apt update
sudo apt install glslc
# Ubuntu 22.04 (Jammy) — install the LunarG Vulkan SDK instead
wget -qO- https://packages.lunarg.com/lunarg-signing-key-pub.asc | \
sudo tee /etc/apt/trusted.gpg.d/lunarg.asc
sudo wget -qO /etc/apt/sources.list.d/lunarg-vulkan-jammy.list \
https://packages.lunarg.com/vulkan/lunarg-vulkan-jammy.list
sudo apt update
sudo apt install -y vulkan-sdk3. Verify the GPU is visible to Vulkan:
vulkaninfo --summary4. Install ggmlR with CPU SIMD acceleration:
sudo Rscript -e 'install.packages("ggmlR", configure.args = "--with-simd")'5. Confirm GPU support from R:
Rscript -e 'library(ggmlR)
ggml_vulkan_status()'Important: on Windows, R ignores
configure.args/--configure-args(they are only honoured by the Unix./configurepath). Use environment variables instead, set in the same R session before installing:
Sys.setenv(GGML_USE_SIMD = "1") # enable CPU SIMD (AVX2/SSE4/FMA/F16C)
Sys.setenv(GGML_USE_VULKAN = "1") # force-enable Vulkan (or "0" to disable)
Sys.setenv(GGML_VK_HARD_EXIT = "1") # enable ggml_vulkan_shutdown(hard = TRUE)The flags are read by configure.win, which only runs when the package is
built from source. The CRAN binary for Windows is pre-built without them,
so pass type = "source" to force a source build:
# From CRAN (force a source build so the flags take effect):
install.packages("ggmlR", type = "source")
# From GitHub (always builds from source):
remotes::install_github("Zabis13/ggmlR", force = TRUE)Vulkan is still auto-detected when the Vulkan SDK is present, so
GGML_USE_VULKAN is only needed to force it on/off. Accepted values:
1 / yes / true / on (and 0 / no / false / off for Vulkan).
OpenMP (multi-threaded CPU executor) is detected automatically on Rtools.
GGMLR_LOG_DEVICE — when set to a non-empty value other than 0, enables
Vulkan telemetry logging. Two messages are emitted:
-
The selected device and its capabilities, once a backend context is created:
ggml_vulkan: device #0 'Vulkan0' selected | fp16=1 bf16=1 coopmat=1 coopmat2=0 bda=1 max_buffer=4095 MB suballoc_block=1024 MB sysmem_fallback=0 -
A per-graph summary of any ops the Vulkan backend could not run and fell back to CPU (e.g.
GET_ROWS_BACK, the backward of an embedding table):ggml_vulkan: 12 op(s) not supported on GPU during graph compute, ran on CPU (per-type: GET_ROWS_BACK=12)This is worth turning on before trusting a training benchmark: a fallback changes only the speed, never the result, so it is otherwise invisible.
inst/examples/backward_gpu_demo.Rchecks the backward ops one at a time. -
Where device buffers are placed, when Resizable BAR is off:
ggml_vulkan: no Resizable BAR (host-visible VRAM heap 256 MiB of 16128 MiB) -- device buffers use device-local memory only; GGML_VK_FORCE_HOST_VISIBLE_VIDMEM=1 overridesThe choice is made from the device's memory heap sizes, not from the driver name or version. Without Resizable BAR the host-visible VRAM window is small, and placing device buffers there sent part of them over PCIe (measured: SAGEConv 127 ms instead of 10). Two variables override the automatic choice:
GGML_VK_FORCE_HOST_VISIBLE_VIDMEM=1keeps the window (the previous behaviour),GGML_VK_DISABLE_HOST_VISIBLE_VIDMEM=1never uses it.
All three are off by default so they do not clutter output (notably during tests, where many contexts and training graphs are created). Enable them to diagnose which GPU was picked, which acceleration paths (fp16 / bf16 / coopmat) are active, and which ops are silently running on CPU:
Sys.setenv(GGMLR_LOG_DEVICE = "1")ggmlR runs the heavy linear algebra on the GPU, but some steps stay on the CPU
through R's BLAS/LAPACK — most notably the eigendecomposition in PCA
(RunGGML(op = "embed")) and any eigen() / prcomp() fallback. R ships with
the reference BLAS/LAPACK, which is single-threaded: on a multi-core machine
those CPU phases use just one core. Switching to OpenBLAS parallelises them
(e.g. PCA on a 3000-gene object dropped from ~24 s to ~1.6 s here).
This is a system-level change (outside the package, needs sudo). On
Debian/Ubuntu:
# 1. Install OpenBLAS. Prefer the OpenMP variant over the default pthreads one
# (the R-admin manual recommends the openmp build): its threads are governed
# by OMP_NUM_THREADS and it is better behaved under tools like valgrind.
sudo apt update
sudo apt install libopenblas-openmp-dev # or libopenblas-dev for the pthreads build
# 2. Point BLAS at OpenBLAS (pick the openblas entry from the list):
sudo update-alternatives --config libblas.so.3-x86_64-linux-gnu
# 3. Point LAPACK at OpenBLAS too — eigen() goes through LAPACK, not BLAS:
sudo update-alternatives --config liblapack.so.3-x86_64-linux-gnu# 4. Confirm R picked it up — both lines should name openblas, not
# /usr/lib/.../blas/libblas.so or /usr/lib/.../lapack/liblapack.so:
si <- sessionInfo(); cat("BLAS: ", si$BLAS, "\nLAPACK:", si$LAPACK, "\n")OpenBLAS grabs all cores by default. If that oversubscribes alongside parallel
resampling (tidymodels / mlr3) or trips CRAN's 2-thread limit in checks, cap it:
Sys.setenv(OPENBLAS_NUM_THREADS = 2).
Note (valgrind). A multithreaded OpenBLAS spawns worker threads, and valgrind reports each worker's thread-local storage as a one-off "possibly lost" block (
definitely lost: 0). This is a known false positive, not a leak. Run checks withOPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1for a clean valgrind report; CRAN's own machines use single-threaded reference BLAS, so it does not appear there.
Accelerate an ordinary R matrix multiply on the GPU without rewriting your code — plain matrices in, plain matrices out, with a transparent CPU fallback.
library(ggmlR)
A <- matrix(rnorm(2000 * 1500), 2000, 1500)
B <- matrix(rnorm(1500 * 1000), 1500, 1000)
C <- ggml_matmul(A, B) # A %*% B on the GPU, returns a plain matrix
G <- ggml_crossprod(A) # t(A) %*% A
H <- ggml_tcrossprod(A) # A %*% t(A)Prefer the operators? Wrap one operand and %*% / crossprod() / tcrossprod()
dispatch to the GPU — nothing else in your code changes:
Ag <- as_gpu_matrix(A)
C <- Ag %*% B # GPU
G <- crossprod(Ag) # GPU- Dispatch —
device = "auto"(default) uses the GPU when one is present and the multiply is large enough to amortise the host↔VRAM transfer; small multiplies and machines without a GPU stay on the CPU. Force it withdevice = "gpu"/"cpu". - Precision — R multiplies in double precision (f64); a GPU offers f32 at
best, so the GPU path is a fast approximate multiply, not a bit-for-bit
replacement.
prec = "f32"(default) requests f32 accumulation; how close the result lands depends on the Vulkan driver (some, e.g. RADV/Mesa, accumulate in f16 regardless, giving ~1e-3relative error either way).prec = "f16"only lowers precision further, for speed. - Full double precision — need bit-accurate results?
ggml_matmul_f64(A, B)runs the multiply in fp64 on the GPU, matching R's%*%to ~1e-15. It only pays off on data-centre cards with fast fp64 (Tesla P100/V100, Instinct); consumer GPUs cripple fp64, so it falls back to (and is usually slower than) the CPU there.
The fastest way to get a model running — stack layers with the pipe, compile, train.
library(ggmlR)
model <- ggml_model_sequential() |>
ggml_layer_dense(128L, activation = "relu", input_shape = 784L) |>
ggml_layer_dropout(rate = 0.3) |>
ggml_layer_dense(10L, activation = "softmax")
model <- ggml_compile(model,
optimizer = "adam",
loss = "categorical_crossentropy",
metrics = "accuracy")
model <- ggml_fit(model, x_train, y_train,
epochs = 10L,
batch_size = 32L,
validation_split = 0.1,
verbose = 1L)
# Important: always assign the return value back to model —
# ggml_fit() returns the model with updated weights.
plot(model$history)
ggml_evaluate(model, x_test, y_test)
preds <- ggml_predict(model, x_new)| Layer | Function |
|---|---|
| Dense | ggml_layer_dense(units, activation) |
| Conv1D | ggml_layer_conv_1d(filters, kernel_size) |
| Conv2D | ggml_layer_conv_2d(filters, kernel_size, padding) |
| MaxPooling2D | ggml_layer_max_pooling_2d(pool_size) |
| GlobalAvgPool2D | ggml_layer_global_average_pooling_2d() |
| BatchNorm | ggml_layer_batch_norm() (RMS-normalizes, then scales/shifts) |
| RMSNorm | ggml_layer_rms_norm() |
| LayerNorm | ggml_layer_layer_norm() (centres as well as scales) |
| Flatten | ggml_layer_flatten() |
| Permute | ggml_layer_permute(dims) — reorder a sample's axes, aperm() order |
| Reshape | ggml_layer_reshape(shape) — one axis may be -1 |
| SequencePooling | ggml_layer_sequence_pooling(mode) — c(seq, d_model) to d_model |
| Dropout | ggml_layer_dropout(rate) |
| Embedding | ggml_layer_embedding(vocab_size, dim) |
| PositionalEmbedding | ggml_layer_positional_embedding() |
| LSTM | ggml_layer_lstm(units, return_sequences) |
| GRU | ggml_layer_gru(units, return_sequences) |
| Attention | ggml_layer_attention(d_model, n_heads, causal, mask, rope, dropout) (functional API) |
| TransformerBlock | ggml_layer_transformer_block(d_model, n_heads, ...) (functional API) |
loss = |
Notes |
|---|---|
"categorical_crossentropy" |
Multi-class; expects a softmax output |
"mse" |
Regression, squared error |
"mae" |
Regression, absolute error — robust to outliers |
"huber" |
Quadratic near zero, linear beyond delta = 1 |
"binary_crossentropy" |
Independent per-output labels; expects a sigmoid output |
"mean_squared_error", "mean_absolute_error", "huber_loss" and
"binary_cross_entropy" are accepted as aliases. In a multi-output model each
head takes its own loss — see Multi-output model.
ggml_fit(sample_weight = ) takes a vector of one weight per sample, or — for
"mse" — an nrow(x) by ncol(y) matrix acting as a per-output mask. A mask
trains only the outputs whose weight is non-zero, which is what a Q-model
needs: 1 on the action that has a target, 0 elsewhere.
model <- ggml_model_sequential() |>
ggml_layer_conv_2d(32L, kernel_size = c(3L, 3L), activation = "relu",
input_shape = c(28L, 28L, 1L)) |>
ggml_layer_max_pooling_2d(pool_size = c(2L, 2L)) |>
ggml_layer_conv_2d(64L, kernel_size = c(3L, 3L), activation = "relu") |>
ggml_layer_global_average_pooling_2d() |>
ggml_layer_dense(10L, activation = "softmax")x and y have to fit in memory twice — once as an R matrix, once as the ggml
copy. Pass a generator instead and only the current batch is ever resident.
It is a function of no arguments returning list(x, y), or NULL when the data
runs out:
# Read batches off disk instead of loading the whole set
gen <- function() {
f <- next_shard() # your own bookkeeping
if (is.null(f)) {
rewind_shards() # reset, so the next epoch starts over
return(NULL) # NULL ends the epoch
}
d <- readRDS(f)
list(d$x, d$y) # exactly batch_size rows
}
model <- ggml_fit(model, generator = gen, epochs = 10L, batch_size = 32L)Note the reset: training runs several epochs over the same source, and a bare
closure gives the loop no handle to rewind it, so a finite generator has to
rewind itself before returning NULL. Forget it and the second epoch gets
nothing — which stops training with a warning rather than looping over an empty
source.
A generator that never returns NULL is fine too — augmentation, simulators, RL
rollouts — but then steps_per_epoch says where an epoch ends:
augmenting <- function() {
i <- sample(nrow(x_train), 32L)
list(random_flip(x_train[i, , drop = FALSE]), y_train[i, , drop = FALSE])
}
model <- ggml_fit(model, generator = augmenting,
steps_per_epoch = 500L, epochs = 10L, batch_size = 32L,
validation_generator = val_gen, validation_steps = 20L)Batches are in ordinary R layout, same as x/y, and must contain exactly
batch_size samples — a short batch is an error rather than being padded.
initial_epoch shifts the epoch numbering seen by callbacks and LR schedules
(as in Keras); it does not restore optimizer state, so Adam starts with zeroed
moments. The arguments that need a whole dataset to be meaningful —
validation_split, validation_data, class_weight, sample_weight — are
rejected with a generator rather than silently ignored.
For a training loop you drive yourself (RL, or anything where "epoch" is not a
pass over a dataset), ggml_fit_opt_gen() is the same loop one level down.
When progress is not measured in passes over data — reinforcement learning,
curricula, anything where the next batch depends on what the model just did —
ggml_trainer() hands back a context that lives until you free it, and you step
it yourself:
tr <- ggml_trainer(model, batch_size = 32L)
on.exit(tr$free())
target <- tr$model() # target network for bootstrapping
for (i in seq_len(10000L)) {
rollout(env, target) # your agent fills a replay buffer
b <- sample_batch(replay, 32L)
loss <- tr$step(b$x, b$y) # one forward + backward, returns the loss
if (i %% 500L == 0L) target <- tr$model() # refresh every 500 steps
}
model <- tr$model() # weights back into a model object
tr$free()Weights, Adam moments and the compute graph stay alive between steps, so a step
costs one forward and backward pass and nothing else. $eval(x, y) measures
without training, $set_lr() changes the rate mid-run, and gradient
accumulation follows nbatch_logical — the optimizer updates only on the period
boundary, and the trainer tracks that internally. The context holds GPU buffers,
so free it when done.
Wire layers into arbitrary graphs — residual connections, multi-input/output, shared weights.
inp <- ggml_input(shape = 64L)
x <- inp |> ggml_layer_dense(64L, activation = "relu")
res <- ggml_layer_add(list(inp, x)) # element-wise add
out <- res |> ggml_layer_dense(10L, activation = "softmax")
m <- ggml_model(inputs = inp, outputs = out)
m <- ggml_compile(m, optimizer = "adam", loss = "categorical_crossentropy")
m <- ggml_fit(m, x_train, y_train, epochs = 5L, batch_size = 32L)inp <- ggml_input(shape = 30L, dtype = "int32", name = "tokens")
emb <- inp |> ggml_layer_embedding(vocab_size = 500L, dim = 32L)
# Branch A: GRU path
proj_a <- emb |>
ggml_layer_gru(32L, return_sequences = FALSE) |>
ggml_layer_dense(32L)
# Branch B: flatten + projection
proj_b <- emb |>
ggml_layer_flatten() |>
ggml_layer_dense(32L, activation = "relu") |>
ggml_layer_dense(32L)
# Residual merge
out <- ggml_layer_add(list(proj_a, proj_b)) |>
ggml_layer_dropout(rate = 0.3) |>
ggml_layer_dense(2L, activation = "softmax")
m <- ggml_model(inputs = inp, outputs = out)Token values must be 0-based integers in
[0, vocab_size - 1].
inp1 <- ggml_input(shape = 20L, name = "timeseries")
inp2 <- ggml_input(shape = 3L, name = "metadata")
br1 <- inp1 |> ggml_layer_dense(16L, activation = "relu")
br2 <- inp2 |> ggml_layer_dense(8L, activation = "relu")
out <- ggml_layer_concatenate(list(br1, br2), axis = 0L) |>
ggml_layer_dense(2L, activation = "softmax")
m <- ggml_model(inputs = list(inp1, inp2), outputs = out)
m <- ggml_compile(m, optimizer = "adam", loss = "categorical_crossentropy")
# Pass x as a list — one matrix per input
m <- ggml_fit(m, x = list(x_ts, x_meta), y = y,
epochs = 10L, batch_size = 32L)
preds <- ggml_predict(m, list(x_ts, x_meta))inp <- ggml_input(shape = 64L)
hidden <- inp |> ggml_layer_dense(64L, activation = "relu")
out <- hidden |> ggml_layer_dense(10L, activation = "softmax")
m <- ggml_model(inputs = inp, outputs = list(hidden, out))
preds <- ggml_predict(m, x)
# preds[[1]] — hidden activations [n × 64]
# preds[[2]] — class probabilities [n × 10]Each output gets its own loss and weight; the optimized total is
sum(loss_weights[i] * loss_i). Losses may be mixed — a cross-entropy head and
a regression head train in the same model.
inp <- ggml_input(shape = 64L)
trunk <- inp |> ggml_layer_dense(64L, activation = "relu")
policy <- trunk |> ggml_layer_dense(10L, activation = "softmax", name = "policy")
value <- trunk |> ggml_layer_dense(1L, name = "value")
m <- ggml_model(inputs = inp, outputs = list(policy, value))
m <- ggml_compile(m, optimizer = "adam",
loss = list(policy = "categorical_crossentropy",
value = "mse"),
loss_weights = c(policy = 1.0, value = 0.5))
# y is a list — one matrix per output, matched by name
m <- ggml_fit(m, x, list(policy = y_policy, value = y_value),
epochs = 10L, batch_size = 32L)
m$history$policy_loss # per-head losses, so a stalled head is visible
m$history$value_lossloss also accepts a single name (applied to every output) or an unnamed list
matched by position. ggml_evaluate() likewise reports <output>_loss per head
alongside the weighted total.
Attention takes and returns a sequence node c(seq_len, d_model), so blocks
stack and a residual merge needs no reshaping. time_distributed = TRUE applies
one dense kernel per position — the position-wise feed-forward sublayer.
inp <- ggml_input(shape = c(64L, 128L)) # c(seq_len, d_model)
h <- inp |> ggml_layer_positional_embedding()
for (i in 1:4) h <- h |> ggml_layer_transformer_block(128L, n_heads = 8L)
out <- h |> ggml_layer_sequence_pooling() |>
ggml_layer_dense(2L, activation = "softmax")
m <- ggml_model(inputs = inp, outputs = out)ggml_layer_sequence_pooling() collapses the sequence into one vector, so the
head's width is d_model however long the sequence is — unlike
ggml_layer_flatten(), whose output grows with it. mode = "first" takes
position 1 instead of the mean, the CLS-token convention.
ggml_layer_transformer_block() is the pre-LN block — normalize, attend,
residual, normalize, feed-forward, residual — assembled from the layers below,
so a block that needs to differ can still be wired by hand. norm = "layer"
swaps RMSNorm for ggml_layer_layer_norm(), which centres as well as scales:
attn <- inp |> ggml_layer_rms_norm() |>
ggml_layer_attention(d_model = 128L, n_heads = 8L)
h <- ggml_layer_add(list(inp, attn)) # residual 1
ff <- h |> ggml_layer_rms_norm() |>
ggml_layer_dense(512L, activation = "silu", time_distributed = TRUE)
ff <- ff |> ggml_layer_dense(128L, time_distributed = TRUE)
h2 <- ggml_layer_add(list(h, ff)) # residual 2Note that activation = "gelu" cannot be trained: ggml has no backward rule
for it, so the graph builds and ggml_fit() then aborts. Use "silu".
causal = TRUE masks keys after the query (GPT-style decoder); mask takes an
explicit additive mask of shape c(seq_q, seq_kv) — 0 where a key may be
attended to, a large negative where it may not — which is how padding is
excluded from a batch of unequal lengths. rope = TRUE applies rotary position
embedding to queries and keys, encoding relative position without weights.
dropout drops query-key links while training and is the identity at inference.
Cross-attention takes queries from one node and keys/values from another, which
may have a different length:
dec <- inp |> ggml_layer_attention(128L, n_heads = 8L, causal = TRUE, rope = TRUE)
pad <- ggml_input(shape = c(64L, 64L))
enc <- inp |> ggml_layer_attention(128L, n_heads = 8L, mask = pad)
ctx <- ggml_input(shape = c(96L, 128L))
x <- dec |> ggml_layer_attention(128L, n_heads = 8L, context = ctx)ggml_attention() returns a reusable layer object — applying it to several
nodes shares one set of projections, which is how an encoder block is reused
across a stack. All heads are computed in one batched pass, so the graph does
not grow with n_heads.
residual_block <- function(x, filters, name) {
main <- x |> ggml_layer_conv_2d(filters, c(3L, 3L), padding = "same",
name = paste0(name, "_conv"))
shortcut <- x |> ggml_layer_conv_2d(filters, c(1L, 1L), padding = "same",
name = paste0(name, "_proj"))
ggml_layer_add(list(main, shortcut), name = paste0(name, "_add"))
}
inp <- ggml_input(shape = c(32L, 32L, 3L))
x <- inp |> ggml_layer_conv_2d(16L, c(3L, 3L), activation = "relu",
padding = "same")
x <- residual_block(x, 16L, "res1")
x <- residual_block(x, 32L, "res2")
out <- x |>
ggml_layer_global_average_pooling_2d() |>
ggml_layer_dropout(rate = 0.4) |>
ggml_layer_dense(3L, activation = "softmax")
m <- ggml_model(inputs = inp, outputs = out)enc <- ggml_dense(32L, activation = "relu", name = "encoder")
x1 <- ggml_input(shape = 16L, name = "left")
x2 <- ggml_input(shape = 16L, name = "right")
h1 <- ggml_apply(x1, enc) # identical weights
h2 <- ggml_apply(x2, enc)
out <- ggml_layer_add(list(h1, h2)) |>
ggml_layer_dense(2L, activation = "softmax")
m <- ggml_model(inputs = list(x1, x2), outputs = out)| Feature | Keras (Python) | ggmlR |
|---|---|---|
| Batch dimension | part of input_shape |
excluded from shape |
| Merge layers | add([a, b]) |
ggml_layer_add(list(a, b)) |
| Shared layers | reuse layer object | ggml_dense() + ggml_apply() |
| Multi-input data | list of arrays | list() of R matrices |
| Multi-output predict | list of numpy arrays | R list of matrices |
| Backend | TensorFlow / JAX / PyTorch | ggml (Vulkan GPU, CPU fallback) |
compile() and fit() are re-exported from the generics package, so they
are shared with keras3 and other packages that use the same generics.
evaluate() is not: ggmlR declares its own.
generics::evaluate names its first argument x, so S3 dispatch happens on
whatever is bound to x. With the keras3 argument order that works positionally
— evaluate(model, x, y) — but naming the data explicitly,
evaluate(model, x = x, y = y), binds the data to the dispatch argument and
fails with "no applicable method". Declaring the generic here as
evaluate(object, ...) makes both forms work.
The trade-off: attaching both ggmlR and keras3 makes the two evaluate generics
mask one another, and whichever package was attached later wins. If you use both,
qualify the call:
ggmlR::evaluate(model, x_test, y_test)Build and train arbitrary architectures with eager execution and automatic differentiation.
library(ggmlR)
# Define parameters
W <- ag_param(matrix(rnorm(4 * 8) * 0.1, 8, 4))
b <- ag_param(matrix(0, 8, 1))
# Forward + backward
with_grad_tape({
h <- ag_add(ag_matmul(W, x_batch), b)
h <- ag_relu(h)
loss <- ag_mse_loss(h, y_batch)
})
grads <- backward(loss)
opt <- optimizer_adam(list(W = W, b = b), lr = 1e-3)
opt$step(grads)
opt$zero_grad()model <- ag_sequential(
ag_linear(64L, 128L, activation = "relu"),
ag_layer_norm(128L),
ag_dropout(0.1),
ag_linear(128L, 10L)
)
params <- model$parameters()
opt <- optimizer_adam(params, lr = 1e-3)
sch <- lr_scheduler_cosine(opt, T_max = 50L, lr_min = 1e-5)
dl <- ag_dataloader(x_train, y_train, batch_size = 32L, shuffle = TRUE)
for (epoch in 1:50) {
for (batch in dl$epoch()) {
with_grad_tape({
out <- model$forward(batch$x)
loss <- ag_softmax_cross_entropy_loss(out, batch$y)
})
grads <- backward(loss)
clip_grad_norm(params, grads, max_norm = 1.0)
opt$step(grads)
opt$zero_grad()
}
sch$step()
}dp_train() splits data across N replicas, averages gradients, and keeps weights in sync automatically.
make_model <- function() {
W <- ag_param(matrix(rnorm(4 * 2) * 0.1, 2, 4))
b <- ag_param(matrix(0, 2, 1))
list(
forward = function(x) ag_add(ag_matmul(W, x), b),
parameters = function() list(W = W, b = b)
)
}
result <- dp_train(
make_model = make_model,
data = my_dataset, # list of samples
loss_fn = function(out, tgt) ag_mse_loss(out, tgt),
forward_fn = function(model, s) model$forward(s$x),
target_fn = function(s) s$y,
n_gpu = 2L, # number of replicas
n_iter = 100L,
lr = 1e-3,
max_norm = 5.0
)
result$loss_history # numeric vector, one value per iteration
result$model # trained replica 0For machines with several Vulkan GPUs, ggmlR ships a native tensor-parallel and pipeline-parallel path (not in upstream ggml, which does tensor-split only for CUDA/SYCL). The cross-device transport defaults to portable host staging (device→host→device), correct on every driver.
Tensor parallelism — split a weight matrix's rows across GPUs, compute each slice on its own device, gather the result:
W <- matrix(rnorm(4096 * 256), nrow = 4096) # 4096 outputs x 256 inputs
X <- matrix(rnorm(8 * 256), nrow = 8) # batch of 8
Y <- ggml_vulkan_split_mul_mat(W, X, n_devices = 2) # == X %*% t(W)TP × DP hybrid — data-parallel over the batch across replicas of tensor-parallel device groups (e.g. 2 replicas × TP=2 on a 4-GPU box):
# replica A = GPUs {0,1}, replica B = GPUs {2,3}; batch split across replicas,
# weights tensor-split within each pair. No cross-replica traffic at inference.
Y <- ggml_tp_dp_forward(W, X, replicas = list(c(0, 1), c(2, 3)))Pipeline parallelism — split a model by layers across GPUs; the activation tensor is handed between stages just once per pass (a single cross-device copy). Suits models too large for one card:
# stage 1 (layers 1..k) on GPU 0 -> stage 2 (layers k+1..n) on GPU 1
mk <- function(dev, Wt) list(
device = dev, in_shape = c(K, M),
build = function(ctx, input) {
w <- ggml_new_tensor_2d(ctx, GGML_TYPE_F32, K, K)
list(output = ggml_relu(ctx, ggml_mul_mat(ctx, w, input)),
set_weights = function() ggml_backend_tensor_set_data(w, as.numeric(Wt)))
})
y <- ggml_pp_forward(list(mk(0L, W1), mk(1L, W2)), x = as.numeric(X), out_shape = c(K, M))See inst/examples/tp_dp_hybrid.R and inst/examples/pp_pipeline.R for complete runnable demos.
Measured with llamaR (which links libggml.a statically) driving Qwen2.5-1.5B-Instruct Q4_K_M (Vulkan, host-staging cross-device transport) on two 4-GPU hosts, decode throughput, median of 3 runs of 128 tokens:
- P100 — 4× Tesla P100-SXM2-16GB
- V100 — 4× Tesla V100-32GB, 2× Xeon E5-2698 v4, 256 GB RAM
| Strategy | GPUs | Split | Decode t/s (P100) | Decode t/s (V100) | Notes |
|---|---|---|---|---|---|
| Baseline | 1 | none | 419.7 | 516.9 | model fits in one card — fastest |
| Pipeline (PP) | 2 | layer | 150.4 | 221.5 | layers spread across 2 GPUs |
| Tensor (TP) | 2 | row | 150.4 | 223.2 | rows split, all-reduce per layer |
| Pipeline (PP) | 4 | layer | 133.3 | 176.3 | more hops → slower |
| Tensor (TP) | 4 | row | 130.0 | 176.6 | more hops → slower |
| TP=2 × DP=2 | 4 | row + replicas | 306 | 446 | 2 replicas × TP=2, run concurrently |
| DP=4 | 4 | replicas | 975 | 1300 | 4 single-GPU replicas, run concurrently |
Same model on an 8× Tesla V100-32GB host (2× Xeon E5-2698 v4, 256 GB RAM), showing how the pattern holds as GPU count doubles:
| Strategy | GPUs | Split | Decode t/s | Notes |
|---|---|---|---|---|
| Baseline | 1 | none | 684.9 | model fits in one card — fastest |
| Pipeline (PP) | 2 | layer | 229.7 | layers spread across 2 GPUs |
| Tensor (TP) | 2 | row | 231.3 | rows split, all-reduce per layer |
| Pipeline (PP) | 4 | layer | 187.5 | more hops → slower |
| Tensor (TP) | 4 | row | 191.0 | more hops → slower |
| Pipeline (PP) | 8 | layer | 142.8 | 8-way split — slowest single-context |
| Tensor (TP) | 8 | row | 137.1 | per-layer all-reduce across 8 cards |
| TP=2 × DP=4 | 8 | row + replicas | 897 | 4 replicas × TP=2, run concurrently |
| DP=8 | 8 | replicas | 2290 | 8 single-GPU replicas, run concurrently |
Takeaway: when a model fits in one GPU, data parallelism (DP — independent replicas) wins by a wide margin: DP=4 delivers ~2.3× the single-card throughput and ~7.5× any split mode. Splitting such a model across cards (PP/TP) only adds cross-device overhead — the ~1 GB/s host-staging transport dominates. PP and TP earn their keep only when the model does not fit in one card (e.g. a 30B+ model on 16 GB cards): there a split is the only way to run it at all, and PP minimizes cross-device hops (one activation copy per pass) while TP maximizes per-token parallelism at the cost of a per-layer all-reduce. Reproduce with llamaR's inst/examples/bench_pp_tp_dp.sh.
Clean shutdown: when a standalone script uses several GPUs, make
ggml_vulkan_shutdown(hard = TRUE)its last statement. This tears down Vulkan and then calls_exit(0), skipping the exit-time loader-static-destruction phase that can otherwise flakily segfault after your results are printed (the results are already computed by then, so the crash is harmless-but-noisy). Use plainggml_vulkan_shutdown()(nohard) mid-session — it releases the devices and is safe to call repeatedly, but does not guarantee a clean process exit on its own.The
_exit(0)path is compiled out by default, because CRAN policy forbids a package terminating the R session. In a default buildhard = TRUEdoes the normal teardown and warns instead of exiting — never silently — so the flaky exit-time segfault can still fire. Compile it in with--configure-args="--enable-hard-exit"(Windows:Sys.setenv(GGML_VK_HARD_EXIT = "1")before installing), and check the current build withggml_vulkan_hard_exit_available().When to use it — and when not.
hard = TRUEonly pays off for multi-GPU scripts (tensor/pipeline parallelism,split_mul_mat/pp_forward), where late device finalizers race the loader teardown on exit. A single-GPU script never triggers that race, so it does not needhard = TRUE. Do not use it under a Jupyter / Kaggle notebook: there_exit(0)kills the R kernel, and the notebook runner (papermill/nbclient) reports it asDeadKernelError: Kernel diedeven though all results were already computed and written. In a notebook, either drop the call or use plainggml_vulkan_shutdown()(soft teardown, kernel stays alive). A portable one-liner that keeps the hard exit only outside notebooks:is_notebook <- nzchar(Sys.getenv("KAGGLE_KERNEL_RUN_TYPE")) || nzchar(Sys.getenv("JPY_PARENT_PID")) ggml_vulkan_shutdown(hard = !is_notebook)
| Category | Functions |
|---|---|
| Linear | ag_matmul, ag_add, ag_sub, ag_mul, ag_scale |
| Activations | ag_relu, ag_sigmoid, ag_tanh, ag_softmax |
| Reductions | ag_sum, ag_mean (with dim, keepdim) |
| Math | ag_log, ag_exp, ag_pow, ag_clamp |
| Shape | ag_reshape, ag_transpose, ag_slice_rows, ag_concat_rows |
| Device buffers | ag_buffer, ag_buffer_write, ag_buffer_read, ag_buffer_free, ag_get_rows (see Device buffers) |
| Attention | ag_multihead_attention, ag_flash_attention (all heads in one fused op) |
| Loss | ag_mse_loss, ag_cross_entropy_loss, ag_softmax_cross_entropy_loss |
| Layers | ag_linear, ag_batch_norm, ag_layer_norm, ag_dropout, ag_embedding |
| Containers | ag_sequential |
| Optimizers | optimizer_sgd, optimizer_adam |
| Schedulers | lr_scheduler_step, lr_scheduler_cosine (SGDR via T_mult), lr_scheduler_onecycle, lr_scheduler_cyclic, lr_scheduler_warmup_cosine |
| Memory | ag_checkpoint (recompute a segment instead of storing its activations) |
| Utilities | clip_grad_norm, clip_grad_value, check_grad_anomaly, ag_gradcheck, dp_train |
ag_multihead_attention() loops over heads. Because ag_* has no slice
operation, pulling one head out of a [d_model, seq] matrix is a matmul
against a selector matrix; each head then costs two more matmuls, two
transposes and a softmax. A 4-head block records 52 tape nodes for the
attention core alone.
ag_flash_attention() does all heads in a single ggml_flash_attn_ext() call,
and its gradient in a single ggml_flash_attn_back() — one tape node.
d_model <- 64L; seq_len <- 32L; n_heads <- 4L
q <- ag_param(matrix(runif(d_model * seq_len, -1, 1), d_model, seq_len))
k <- ag_param(matrix(runif(d_model * seq_len, -1, 1), d_model, seq_len))
v <- ag_param(matrix(runif(d_model * seq_len, -1, 1), d_model, seq_len))
with_grad_tape({
out <- ag_flash_attention(q, k, v, n_heads)
loss <- ag_mse_loss(out, target)
})
backward(loss)Rows are split into n_heads contiguous blocks — the same head layout
ag_multihead_attention() uses. Projections are not included: apply W_q,
W_k, W_v before the call and W_o after it. k and v may be a different
length from q (cross-attention) as long as they match each other.
Attention core only, forward and backward, CPU backend
(inst/scripts/measure_ag_flash_attn.R):
| tape nodes | head loop | fused | ||
|---|---|---|---|---|
| d32 h2 seq16 | 26 | 1.47 ms | 0.04 ms | 36.5× |
| d64 h4 seq32 | 52 | 3.52 ms | 0.15 ms | 23.2× |
| d128 h8 seq64 | 104 | 10.72 ms | 0.70 ms | 15.4× |
| d256 h8 seq128 | 104 | 41.69 ms | 2.16 ms | 19.3× |
Those figures compare the same work on both sides — the projections and the loss are excluded from each. They are an upper bound: the fused column pays no R-level tape bookkeeping at all.
Masking is supported in both directions — the same mask is applied to the forward pass and to the gradient:
# Decoder-style: query j attends to keys 1..j
out <- ag_flash_attention(q, k, v, n_heads, causal = TRUE)
# Or supply one. NOTE the orientation: [seq_kv, seq_q] -- keys index the ROWS,
# queries the columns, which is the transpose of the usual "row = query"
# convention. Logical (TRUE = attend) or numeric (0 / -Inf) both work.
allow <- matrix(TRUE, seq_kv, seq_q)
allow[padding_positions, ] <- FALSE
out <- ag_flash_attention(q, k, v, n_heads, mask = allow)A [seq_q, seq_kv] matrix is rejected rather than accepted quietly, because
for a square case it would mask the transposed entries and nothing would say
so.
Not supported: ALiBi (max_bias), logit softcap and attention sinks —
ggml_flash_attn_back() asserts on all three.
ag_multihead_attention() stays as the fallback for what the fused op does not
cover: attention dropout, and the packaged projections — it owns W_q,
W_k, W_v, W_o and b_o, where ag_flash_attention() takes
already-projected Q/K/V and leaves the output projection to you. For plain and
causal attention without dropout, prefer the fused op.
Every autograd op records what its backward rule will need — ag_matmul keeps
the activations flowing through it, and on a deep stack those snapshots are
most of what the tape holds. ag_checkpoint() runs a segment of the forward
pass without recording it, keeping only the segment's input, and re-runs
that segment during backward() to rebuild the values the rules need. Memory
drops by whatever the segment would have stored; the price is running its
forward pass twice.
block <- ag_linear(256L, 256L, activation = "relu")
with_grad_tape({
h <- ag_checkpoint(function(inp) block$forward(inp), x)
loss <- ag_mse_loss(h, y)
})
backward(loss)Typically you checkpoint every n-th block of a deep stack. Note the
local(): the segment is called again later, during backward(), so it must
not close over the loop variable — otherwise every segment ends up sharing one
binding, all of them see the last layer at replay time, and the earlier layers
silently receive no gradient while the loss still looks normal.
with_grad_tape({
h <- x
for (i in seq_along(layers)) {
seg <- local({
lyr <- layers[[i]] # captured per iteration
function(inp) lyr$forward(inp)
})
h <- if (i %% 2L == 0L) ag_checkpoint(seg, h) else seg(h)
}
loss <- ag_mse_loss(h, y)
})
backward(loss)On a 12-layer 256×256 stack at batch 64, measured on the tape itself:
| tape held | saved | |
|---|---|---|
| no checkpointing | 10.13 MB | — |
| every 3rd block | 6.88 MB | 32% |
| every 2nd block | 5.23 MB | 48% |
| every block | 0.25 MB | 98% |
Gradients are identical to a non-checkpointed run in each case.
fn must be reproducible: given the same inputs and RNG state it has to
compute the same thing. Randomness inside a segment is handled — the RNG state
is saved at the forward pass and restored before the replay, so an
ag_dropout() inside a checkpoint replays with the same mask rather than
drawing a fresh one (which would apply a mask the forward never used). Anything
non-deterministic beyond R's RNG is not.
The output re-enters the tape as one opaque node: detached from the segment's internals, so gradients cannot arrive twice, but still linked to the inputs you passed in, so earlier layers train as usual. A tape containing a checkpoint always uses the closure backward — the segment is R code, not something the graph path below can turn into nodes.
backward() builds the entire backward pass as one ggml graph and computes
it in a single call, so intermediate gradients never leave the device. The
alternative — one R closure per tape node, each computing its gradient with
%*% and t() on the host — is still there as a fallback.
The graph path is the default. Turn it off per session or per run:
ggmlR:::ag_backward_graph(FALSE) # or: GGMLR_AG_BACKWARD_GRAPH=0 Rscript ...
ggmlR:::ag_backward_path() # "graph", "closures", or "closures (<why>)"The path is GPU-only and all-or-nothing: it emits ag_matmul,
ag_add (both broadcasts), the three losses, relu/sigmoid/tanh,
ag_transpose, ag_softmax, ag_scale and ag_mul — enough for dense
stacks, classifiers, ag_multihead_attention (causal included) and
ag_dropout. A tape holding any other operation — ag_sub, ag_sum,
ag_layer_norm, ag_embedding, ag_batch_norm and the rest — runs the
closure path in full rather than splitting the pass between the two.
ag_backward_path() names the operation that caused the fallback. Gradients
are identical either way; the test suite checks the graph path against the
closures on every covered shape.
It became the faster path once weights stopped moving. Measured across five
shapes it runs 1.2–2.7× the closures, with only the smallest (64×64, depth 2)
at parity. It was a measured slowdown before resident weights: stage profiling
put the actual GPU compute at 12–25% of backward(), with the rest spread over
uploading forward snapshots, downloading gradients, building nodes in R and
writing $grad back. Two of those four are now gone — a snapshot that is a
weight is already on the device, and gradients stay there — which is what turned
the ratio around.
Measure before trusting either path on your own shapes:
Rscript inst/scripts/measure_ag_backward_real.R # ratios per model
GGMLR_AG_BENCH_PROF=1 Rscript inst/scripts/measure_ag_backward_real.R # + stage breakdownGGMLR_AG_BWD_PROF=1 turns the same per-stage timings on inside any script;
read them with ggmlR:::ag_backward_profile_report().
Short answer: for the ag_* ops on their own, on the hardware measured so far,
it is not — and the reason is worth knowing before reaching for ag_device("gpu").
ag_* operation uploaded its operands and downloaded
its result, which is no longer what happens on the ops listed in "Keeping
tensors on the device" — a step now costs 8 crossings rather than 10, and the
forward's share of them no longer grows with depth. Treat these numbers as the
shape of the old problem, and re-run
inst/scripts/measure_ag_gpu_threshold.R before deciding anything from them.
Every ag_* operation is executed on its own: .ag_run_op creates tensors for
the inputs, uploads their data, builds a one-node graph, computes it and reads
the result back. The arithmetic is a small part of that. A trivial ag_scale
on a 4 MB operand costs about 16 ms, and the cost tracks the bytes moved at
roughly 4 ms/MB — more than a 1024×1024 matmul takes on the CPU altogether.
The consequence is that the usual intuition, "it will pay off once the matrices
are big enough", does not hold here. Measured with
inst/scripts/measure_ag_gpu_threshold.R (radv, multi-threaded BLAS on the CPU
side):
| what grows | CPU vs GPU |
|---|---|
layer size d, batch 64 |
CPU ahead throughout, and the gap widens: 1.5× at d=1024, 12.7× at d=2048 |
| batch, layer fixed at 512 | gap narrows 14× → 1.35× and then flattens; it does not cross over |
| 4-layer MLP forward | steady 4.2–4.9× for the CPU at every size — per-op cost compounds with depth |
ag_flash_attention |
1.26× for the GPU at d_model=64, a tie at 256 and 512 |
Two things follow. Growing the layer makes it worse, not better: the upload
grows with d² while the CPU side has BLAS threading to fall back on. Growing
the batch is the direction that helps — the weight upload is paid once per call
however wide the batch is — but it converges to about 1.35× rather than crossing
zero.
The last row is the useful one. ag_flash_attention uploads Q, K and V once and
runs the whole attention inside a single call, which is the best case for this
execution model: no repeated uploads at all. It reaches parity with the CPU, not
a win. So an upload cache for unchanged operands — the obvious next optimisation,
and a real one — would narrow the gap rather than reverse it.
Where that leaves the GPU path:
- Use it when the alternative is not the CPU: models whose weights do not fit
in host memory comfortably, work already resident on the device, or the
nn_*sequential/functional path, which builds one graph for the whole model and does not pay per-op upload at all. - Do not reach for it to speed up an
ag_*chain on this class of hardware. Setag_device("cpu")— the default — and the same code runs faster. - Re-measure on different hardware before generalising. These numbers come from radv, which announces itself as non-conformant, against a multi-threaded CPU BLAS. A proprietary driver, or a card with much greater memory bandwidth, is a different experiment.
Rscript inst/scripts/measure_ag_gpu_threshold.R # the table above
Rscript inst/scripts/measure_ag_upload_cost.R # where the time inside one op goesThis is about ag_* specifically. The nn_* path compiles a model into a single
ggml graph, so it is not subject to the per-op upload described here, and the
Vulkan backend behaves quite differently there.
The per-op round trip above has an obvious remedy — leave values on the device — and it now applies to the whole training step. Four kinds of tensor stay in backend buffers:
- weights —
ag_param()uploads once into a residency pool that survives the tape resetwith_grad_tape()performs at every step; - optimizer moments — Adam's
mandvlive beside the weights, and the step updates all three on the device; - gradients — the graph backward puts a device handle in each leaf's
$grad, and the optimizer consumes it there; - forward activations — the ops take handles instead of materialising their operands.
Measured on one training step: 4 host/device crossings and 0.047 MB, against 10 and 0.188 before, split as forward 3 / backward 1 / step 0. What remains is the batch going up, the loss scalar coming back, and one backward snapshot — the optimizer step touches the host not at all. Forward traffic no longer scales with depth either: a four-layer network with four different activations pays the same three crossings as a single layer, because nothing between the input and the loss goes through R.
Both defaults can be turned off:
ggmlR:::ag_backward_graph(FALSE) # GGMLR_AG_BACKWARD_GRAPH=0
ggmlR:::ag_backward_resident(FALSE) # GGMLR_AG_RESIDENT_GRADS=0Numerically nothing changes: this is transport, not arithmetic, and the test suite checks gradients against the closure path and against finite differences.
One thing to know if you write code against $grad. With resident gradients
it holds a device handle, not a matrix. Read it through the supported accessors
and nothing changes; do arithmetic on it directly and you get a handle, which
has no arithmetic methods and so fails loudly rather than computing something
wrong. The same applies to holding a gradient across a tape boundary — its
buffer belongs to the step that produced it. dp_train() and ag_checkpoint()
both do this and both materialise first.
A note on how this was measured, because the subsystem figure misleads. The forward is more transfer-bound than the backward (46–85%, upload alone up to 72%), so making it resident speeds the forward up by 1.9–8.2x — but the forward is only 39–58% of a step, so the step gains 1.2–1.4x. Measure what a change does to the step, not to the part it touches: the two differ by a factor of six here, and the subsystem number is the noisier of the two.
Two reports come with it:
ggmlR:::ag_forward_profile(TRUE) # then run a forward pass
ggmlR:::ag_forward_profile_report() # ctx/create/flush/upload/graph/compute/download
ag_tape_memory() # what the tape is holding, right nowag_tape_memory() splits the tape into operands, which outlive it, and
activations, which clearing it would release — counting storage shared between
the two only once, and charging it to the operand. A weight's matrix and the
snapshot a backward closure captured of it are the same object, so reporting it
as an activation would claim the tape can free memory the parameter keeps alive
regardless.
Rscript inst/scripts/measure_ag_resident_gain.R # what residency delivered
Rscript inst/scripts/measure_ag_forward_profile.R # where the forward goes
Rscript inst/scripts/measure_ag_tape_memory.R # tape composition by shapeA large table that is written a few columns at a time and read in random
batches — a replay buffer in off-policy reinforcement learning is the typical
case — would cost a full upload per write as an ag_tensor, or a batch upload
per read if kept in R. ag_buffer() keeps it on the device instead: writes send
only the new columns, and ag_get_rows() sends only the indices.
ag_device("gpu")
buf <- ag_buffer(rows = 6, capacity = 100000) # one column per transition
ag_buffer_write(buf, new_columns, col_offset = pos) # in place, synchronous
# a batch, split into fields on the device
f <- ag_capture(function(idx) {
x <- ag_get_rows(buf, idx) # rows x B, no gradient
list(obs = ag_slice_rows(x, 0, 3), act = ag_slice_rows(x, 3, 2))
}, params = list(buf))
batch <- f(idx = matrix(sample.int(100000, 256) - 1, ncol = 1))
ag_buffer_read(buf, 0, 1000) # download a range (checkpoints)
ag_buffer_free(buf) # release now, not at gc()ag_get_rows()selects columns of the Rrows x capacityview; the name followsggml_get_rows, which selects along ggml's second dimension. Indices are 0-based. Host indices are checked; indices already on the device (inside a recording) are clamped to0..capacity-1.- Inside
ag_capture()pass the indices as an argument and list the buffer inparams: one recording serves every batch of the same size, and freeing or re-uploading the buffer records anew. ag_slice_rows()cuts a range of rows on the device (a view made contiguous); its gradient is zero outside the range, so it can feed a network.- Creation fails when the buffer needs more than 85% of the free device memory
the driver reports.
ag_device("cpu")moves buffers to the host first (about twice their device size, as R doubles); the next GPU use uploads them again. - On the CPU the same calls work on an R matrix.
ag_estimate_training_memory(
shapes = list(c(512, 512), c(512, 512), c(512, 512)),
batch_size = 64, optimizer = "adam")Adds up weights, gradients, optimizer state and activations from parameter shapes alone — no model needed, so the question "will this fit" can be asked before building one.
It counts 8 bytes per scalar, not 4, and that is the most useful thing it
says. ag_dtype() controls the precision of forward tensors uploaded into a
ggml buffer; it does not reach backward(), which accumulates $grad in R
closures, nor the optimizer's moments, which are R matrices. So an f32 estimate
of an ag_* training run is half the real figure. Two things follow that are
usually paid for elsewhere: gradient underflow does not occur, so there is no
loss scaling to configure, and optimizer state is already in higher precision
than the forward pass.
For a model that exists, ag_tape_memory() measures rather than estimates.
ggmlR ships with mlr3 learners for tabular classification and regression. After library(ggmlR) (with mlr3 installed), sequential and functional ggmlR networks become available as first-class learners:
library(mlr3)
library(ggmlR)
# Classification on iris (GPU auto-detected via backend = "auto")
task <- tsk("iris")
learner <- lrn("classif.ggml",
epochs = 50L,
batch_size = 16L,
backend = "auto") # "auto" | "cpu" | "gpu"
learner$predict_type <- "prob"
learner$train(task)
pred <- learner$predict(task)
pred$score(msr("classif.logloss"))- Both ggmlR APIs —
model_fncan return aggml_sequential_modelorggml_functional_model. The default builder isggml_default_mlp(), an exported MLP builder you can also use directly. - Vulkan GPU — set
backend = "gpu"(or leave"auto") and the learner trains and predicts on GPU. - Parallel tuning — the learners declare
properties = "marshal"and implement in-memory marshalling (SHA-256-checksummed container), so trained models can be shipped tofuture/callrworkers without file-system round-trips. - Weighted training —
classif.ggmlhonourstask$weights_learner, mapping them tosample_weightinggml_fit(). - Callbacks for tuning — pass
ggml_callback_early_stopping()etc. via thecallbacksparameter to drive early stopping insidemlr3tuning runs. - Custom architectures — set
learner$model_fn <- function(task, n_features, n_out, pars) { ... }to build any ggmlR network you like; the learner handles task → matrix conversion, compilation, training, and prediction.
# Regression with a custom architecture
library(mlr3)
library(ggmlR)
learner <- lrn("regr.ggml", epochs = 100L)
learner$model_fn <- function(task, n_features, n_out, pars) {
ggml_model_sequential() |>
ggml_layer_dense(256L, activation = "gelu", input_shape = n_features) |>
ggml_layer_dropout(0.2) |>
ggml_layer_dense( 64L, activation = "gelu") |>
ggml_layer_dense(n_out, activation = "linear")
}
learner$train(tsk("mtcars"))Only numeric features are supported: wrap the learner in a pipeline
(po("encode") %>>% po("scale") %>>% lrn("classif.ggml")) when the task has
factor columns. mlr3, paradox, R6, and mlr3pipelines are Suggests;
ggmlR only wires up the integration when they are installed.
ggmlR registers a "ggml" engine for parsnip::mlp(), letting you use ggmlR networks inside the tidymodels ecosystem.
library(ggmlR)
library(parsnip)
spec <- mlp(
hidden_units = 64,
epochs = 100,
dropout = 0.2,
learn_rate = 0.001
) |>
set_engine("ggml",
batch_size = 32,
backend = "auto", # "auto" | "cpu" | "gpu"
verbose = 0) |>
set_mode("classification")
fit_obj <- fit(spec, Species ~ ., data = iris)
predict(fit_obj, new_data = iris)
predict(fit_obj, new_data = iris, type = "prob")Regression works the same way:
spec_reg <- mlp(hidden_units = 32, epochs = 200) |>
set_engine("ggml", batch_size = 8, backend = "gpu") |>
set_mode("regression")
fit_reg <- fit(spec_reg, mpg ~ ., data = mtcars)
predict(fit_reg, new_data = mtcars)parsnip, tibble, rlang, and dials are in Suggests — ggmlR only wires up the engine when they are installed.
Run GPU-accelerated operations directly on Seurat objects — no conversion on your side, and no hard dependency: Seurat/SeuratObject stay in Suggests, so ggmlR installs fine without them and the adapter activates only when they are present.
Setup. The adapter needs two packages installed alongside ggmlR — Seurat (the object model and its pipeline) and a couple of light helpers it leans on for speed (Matrix for the sparse graphs, FNN for the kd-tree kNN). They are all Suggests, so install them yourself:
install.packages(c("Seurat", "Matrix", "FNN")) # ggmlR is already installed
library(ggmlR)
library(Seurat) # ggmlR's S3 methods (RunGGML, ggml_extract, ...) activate
# automatically once SeuratObject is on the search pathIf FNN is absent the kNN falls back to the GPU/CPU distance matrix; if Matrix is absent the "neighbors" op is unavailable (the graphs are sparse). The rest works with Seurat alone.
RunGGML() is the one-call, Seurat-style entry point (object in, object out — pipe-friendly, mirrors RunPCA()). Vulkan is used automatically when a GPU is present, with a transparent CPU fallback. The supported operations are the heavy matrix steps of a standard pipeline:
op |
Replaces | What runs on the GPU |
|---|---|---|
"normalize" |
NormalizeData() (LogNormalize) |
sparse log1p(x / colSum × sf) over the stored non-zeros (sparse_lognorm.comp) — the counts stay a dgCMatrix, never densified |
"scale" |
ScaleData() |
per-gene z-score (x − mean) / sd + clamp over the full dense matrix |
"embed" |
RunPCA() |
gene-by-gene covariance multiply (eigendecomposition stays on the CPU — ggml has no eigensolver) |
"umap" |
RunUMAP() |
two custom compute shaders: pairwise_dist.comp (kNN distances, honest f32) and umap_sgd.comp (SGD layout) |
"neighbors" |
FindNeighbors() |
kNN distances feeding the SNN/Jaccard graph; the FNN kd-tree by default, or a fused GPU kNN (knn_tiled.comp) with knn_backend = "vulkan". Exact kNN either way |
"largest_gene" |
percent.Largest.Gene |
per-cell highest-expressed gene + its share of the cell total, over the sparse dgCMatrix CSC slots (no densify); bit-exact with qlcMatrix::colMax |
"normalize" and "scale" write the transformed matrix back into the assay (the data and scale.data layers), so the rest of the Seurat pipeline picks them up unchanged. "embed" and "umap" add a dimensionality reduction; "neighbors" writes the <assay>_nn and <assay>_snn graphs into @graphs, exactly where FindClusters() looks; "largest_gene" adds largest_gene and percent.Largest.Gene columns to meta.data / colData.
A standard Seurat workflow maps onto the adapter like this — five of the heavy steps move to the GPU, and only the final community detection stays on the CPU:
| Standard step | ggmlR | Runs on |
|---|---|---|
NormalizeData() |
RunGGML(op = "normalize") |
GPU |
ScaleData() |
RunGGML(op = "scale") |
GPU |
RunPCA() |
RunGGML(op = "embed") |
GPU matrix multiply (eigensolve on CPU — ggml has none) |
RunUMAP() |
RunGGML(op = "umap") |
GPU (distance + SGD shaders) |
FindNeighbors() |
RunGGML(op = "neighbors") |
kNN on CPU (FNN kd-tree) or GPU (knn_backend = "vulkan", knn_tiled.comp) → sparse SNN |
FindClusters() |
— (use Seurat's) | CPU — iterative graph Louvain/Leiden, a poor GPU fit and already fast on the CPU |
On a Vulkan GPU (AMD RADV) the GPU steps are markedly faster than a naive reference. Indicative numbers at 2000 cells:
| Step | What was sped up | Speed-up |
|---|---|---|
neighbors distance kernel |
tiled f32 vs stats::dist |
up to ~4× |
neighbors GPU kNN (knn_backend = "vulkan") |
fused knn_tiled.comp vs FNN kd-tree |
~2–3× at 25–50k cells, and grows — the kNN cost is linear in n rather than the kd-tree's ~O(n²) in ~10-D |
largest_gene |
sparse CSC argmax vs qlcMatrix::colMax |
~30× (≈20 s → 0.6 s at 53k cells) |
umap pipeline |
kd-tree kNN + sparse fuzzy graph + GPU SGD | ~13× (≈1.45 s → 0.11 s) |
umap SGD shader alone |
one GPU dispatch per epoch | ~10⁹ edge-updates/s |
The whole pipeline A/B is reproducible with inst/examples/seurat_op2_gpu.R,
which runs the classic Seurat route twice — stock CPU vs RunGGML() on the GPU —
on the Kaggle Open Problems – Single-Cell Perturbations counts (18 211 genes ×
240 090 PBMCs), then checks the two arms agree. A representative run (11 %
subsample = 23 279 cells, 2000 HVGs, 50 PCs, --gpu-knn, AMD RADV):
| step | cpu (s) | gpu (s) | speedup |
|---|---|---|---|
largest_gene |
8.92 | 0.28 | 32.2× |
normalize |
2.84 | 1.61 | 1.8× |
scale |
1.84 | 1.94 | 0.9× |
embed (PCA) |
15.92 | 3.03 | 5.3× |
neighbors |
4.43 | 1.45 | 3.1× |
umap |
15.32 | 5.34 | 2.9× |
| TOTAL (GPU ops) | 49.27 | 13.66 | 3.6× |
Every accelerated step matches Seurat to float noise (normalize/scale max abs err
~1e-6, PCA |cor| = 1.0000 over PC1–10, clusters ARI 0.94, largest_gene
top-gene agreement 1.0000). scale comes out ~1× — it is memory-bound with
nothing to accelerate, so it defaults to the CPU even under Vulkan; PCA's
covariance multiply is the biggest matrix-multiply win. See the
single-cell-seurat vignette for the full breakdown.
(Numbers are hardware-dependent; reproduce the UMAP shaders separately with
inst/examples/umap_shaders_bench.R.)
library(Seurat)
# A whole standard pipeline, GPU-accelerated end to end — every heavy step is a
# RunGGML() call; only FindClusters (graph Louvain) stays on Seurat's side:
pbmc <- RunGGML(pbmc, op = "normalize") # -> "data" layer
pbmc <- RunGGML(pbmc, op = "scale") # -> "scale.data" layer
pbmc <- RunGGML(pbmc, op = "embed", n_components = 30,
reduction_name = "pca") # -> reduction "pca"
# UMAP on the PC coordinates — both phases on the GPU (distance + SGD shaders):
pbmc <- RunGGML(pbmc, op = "umap", reduction = "pca", dims = 1:30,
reduction_name = "umap") # -> reduction "umap"
# Neighbour graphs on the PC coordinates -> @graphs, then cluster as usual:
pbmc <- RunGGML(pbmc, op = "neighbors", reduction = "pca", dims = 1:30)
pbmc <- FindClusters(pbmc, graph.name = paste0(DefaultAssay(pbmc), "_snn"))
DimPlot(pbmc, reduction = "umap", group.by = "seurat_clusters")The "umap" op has two phases. The graph/distance phase uses a tiled f32
pairwise-distance kernel that sidesteps mul_mat (whose f16 accumulation reorders
nearest neighbours and corrupts the graph), with a kd-tree kNN and a sparse fuzzy
graph. The SGD layout phase defaults to the single-threaded reference, which
runs in compiled C (R_umap_sgd_cpu) — the earlier interpreted-R loop took
~700 s on ~10k cells; the C version is bit-for-bit identical but finishes in
seconds, so it stays the default for best embedding quality. Passing
sgd_backend = "vulkan" opts into the umap_sgd.comp shader (one dispatch per
epoch, Hogwild updates — faster still, but the lock-free races smear dense-graph
clusters, a known hard problem for async UMAP-SGD). Layout numerics match the CPU
reference to float32 precision; the SNN graph matches FindNeighbors()
bit-for-bit on identical exact kNN.
The adapter is layered, and each layer is a public generic you can call on its own:
| Function | Layer | Responsibility |
|---|---|---|
ggml_extract() |
Extraction | Pull a feature × cell matrix out of a Seurat/dgCMatrix/matrix; handles Seurat v4 (GetAssayData) vs v5 (LayerData) and sparse → dense |
ggml_run() |
Dispatch | Validate against ggml_ops_registry(), route to Vulkan GPU or CPU (auto, with fallback) |
ggml_inject() |
Injection | Write the result back into the object — a reduction (CreateDimReducObject()) for "embed"/"umap", an assay layer for the "normalize"/"scale" transforms, or <assay>_nn/<assay>_snn Graph objects in @graphs for "neighbors" |
# Compose the layers manually (e.g. on a bare matrix, no Seurat needed):
mat <- ggml_extract(expr_matrix) # genes × cells, dense
task <- ggml_task("embed", mat, params = list(n_components = 30))
res <- ggml_run(task) # ggml_result: cells × components
res$embedding
# Check capabilities before dispatch:
ggml_ops_registry() # all supported operations
ggml_ops_registry("embed") # required params + descriptionA SingleCellExperiment (Bioconductor) path with an S4 runGGML() generic is planned next.
Load pre-trained ONNX models from PyTorch, TensorFlow, or other frameworks and run inference on Vulkan GPU or CPU. No Python or external libraries required — ggmlR includes a built-in zero-dependency protobuf parser.
library(ggmlR)
# 1. Load the model
model <- onnx_load("squeezenet.onnx")
model
#> ONNX Model: torch_jit
#> Producer: pytorch 2.0.1
#> IR version: 8 / Opset: 18
#> Nodes: 66 / Weights: 26
# 2. Check expected inputs
onnx_inputs(model)
#> $x.1
#> [1] 1 3 224 224
# 3. Prepare input data (flat numeric vector, row-major NCHW order)
img <- runif(1 * 3 * 224 * 224)
# 4. Run inference — pass a named list matching input names
result <- onnx_run(model, list(x.1 = img))
# 5. Get predictions
scores <- result[[1]]
cat("Predicted class:", which.max(scores) - 1L, "\n")onnx_load() parses the .onnx file, builds a ggml computation graph, and allocates tensors on the specified device. Weights are loaded from the file via memory-mapping (zero-copy).
# Auto-detect device (Vulkan GPU if available, else CPU)
model <- onnx_load("model.onnx")
# Force CPU
model <- onnx_load("model.onnx", device = "cpu")
# Force Vulkan GPU
model <- onnx_load("model.onnx", device = "vulkan")Some models (BERT, RoBERTa, etc.) have dynamic dimensions for batch size or sequence length. Specify fixed shapes at load time:
model <- onnx_load("bert.onnx",
input_shapes = list(
input_ids = c(1L, 128L),
attention_mask = c(1L, 128L)
))If you forget, onnx_load() will tell you which inputs need shapes:
Error: Input 'input_ids' has dynamic shape [?x?].
Specify fixed shape via input_shapes parameter.
# Print overview
model
#> ONNX Model: torch_jit
#> Nodes: 533 / Weights: 199
#> Ops: MatMul, Add, LayerNormalization, Softmax, ...
# Detailed metadata
onnx_summary(model)
# Input names and shapes (what to pass to onnx_run)
onnx_inputs(model)
# Backend placement: GPU vs CPU split, scheduler info
onnx_device_info(model)# Single input model
result <- onnx_run(model, list(x = my_data))
# Multiple inputs
result <- onnx_run(model, list(
input_ids = as.numeric(token_ids),
attention_mask = rep(1, 128)
))
# Result is a named list of output tensors
str(result)
#> List of 1
#> $ output: num [1:1000] 0.00123 0.00045 ...ONNX models expect inputs in row-major order (batch, channels, height, width for images). R matrices are column-major, so you may need to transpose:
# Image classification: model expects [1, 3, 224, 224]
# If you have a 224x224x3 R array:
img_array <- array(runif(224 * 224 * 3), dim = c(224, 224, 3))
# Rearrange to NCHW: [1, 3, 224, 224] — channel first
img_chw <- aperm(img_array, c(3, 1, 2)) # [3, 224, 224]
input <- as.numeric(img_chw) # flat vector, row-major
result <- onnx_run(model, list(x = input))For NLP models, inputs are typically 1D integer sequences:
# BERT-style: token IDs as numeric vector
tokens <- c(101, 2023, 2003, 1037, 3231, 102, rep(0, 122)) # pad to 128
result <- onnx_run(model, list(input_ids = as.numeric(tokens)))# Classification: get top-5 predictions
scores <- result[[1]]
top5 <- order(scores, decreasing = TRUE)[1:5]
cat("Top-5 classes:", top5 - 1L, "\n") # 0-based class indices
cat("Top-5 scores:", scores[top5], "\n")
# Apply softmax if model outputs logits (not probabilities)
probs <- exp(scores) / sum(exp(scores))A loaded model can be run any number of times; weights stay on the GPU and are never re-transferred. Models whose shapes depend on data (NonZero, NonMaxSuppression, TopK — MaskRCNN, RoBERTa) are run in segments, and the segments after the first are rebuilt on every call at the sizes that call measures; the rebuild costs about 4 ms per run on RoBERTa and about 80 ms on MaskRCNN, and memory does not grow from one run to the next:
model <- onnx_load("classifier.onnx")
for (batch in data_batches) {
result <- onnx_run(model, list(x = batch))
# process result...
}All 15 ONNX Model Zoo models below load, run, and are checked value by value
against ONNX Runtime with inst/scripts/ref_check_vs_onnxruntime.sh, which runs
each model through both implementations and diffs the outputs rather than their
length. On the CPU all 15 agree within 1e-3 (worst max|d| about 2e-5), and the
quantised detector MaskRCNN-12-int8 matches bit for bit (max|d| = 0) on both
CPU and Vulkan. On Vulkan five transformer models (XCiT, CaiT, BERT, RoBERTa,
GPT-NeoX) land between 1e-3 and 7e-3: float sums over long reductions taken in a
different order, not a wrong operator.
| Model | Nodes | Key ops |
|---|---|---|
| mnist-8 | 12 | Conv, Relu, MaxPool, Reshape, MatMul |
| squeezenet1.0-8 | 66 | Conv, Relu, MaxPool, Concat, GlobalAveragePool, Softmax |
| adv_inception_v3 (Opset 17/18) | 215 | Conv, BatchNorm, Relu, Concat, AveragePool |
| emotion-ferplus-8 | 52 | Conv, Relu, MaxPool, Gemm, Constant |
| bat_resnext26ts (Opset 18) | 570 | Conv, BatchNorm, SiLU, Concat, Expand, Split |
| bert (Opset 17) | 533 | MatMul, LayerNorm, GELU/Erf, Softmax, Shape, Gather, Where |
| gptneox (Opset 18) | 482 | MatMul, LayerNorm, GELU, Softmax, Shape, Gather |
| MaskRCNN-12-int8 | 3001 | QLinearConv, QLinearMatMul, NonMaxSuppression, RoiAlign, NonZero, TopK |
| roberta-9 | 1180 | MatMul, LayerNorm, Erf, Softmax, Shape, Gather, Cast |
| sageconv (Opset 16) | 24 | MatMul, Add, Mul, Sigmoid, ScatterElements |
| super-resolution-10 | 12 | Conv, Reshape, Transpose |
| botnet26t_256 (Opset 16) | 530 | Conv, BatchNorm, RelPosBias2D (fused custom op), Softmax |
| xcit_tiny | 436 | MatMul, LayerNorm, Softmax, Concat, Transpose |
| cait_xs24_384 (Opset 16) | 1748 | MatMul, LayerNorm, Softmax, Transpose, Gather |
Arithmetic: Add, Sub, Mul, Div, Pow, Mod, Sqrt, Exp, Log, Abs, Neg, Floor, Ceil, Clip, Erf, Sin, Cos, Max, Min.
Comparison / logic: Equal, Greater, GreaterOrEqual, Less, LessOrEqual, And, Or, Not, Xor, Where.
Linear: MatMul (batched), Gemm, Einsum (limited set of patterns).
Convolution: Conv (1D/2D, grouped, depthwise), ConvTranspose (1D/2D), with auto_pad (SAME_UPPER, SAME_LOWER).
Pooling / resize: MaxPool, AveragePool, GlobalAveragePool, Resize/Upsample (nearest, bilinear).
Normalization: BatchNormalization, LayerNormalization, GroupNormalization, InstanceNormalization, RMSNormalization.
Activations: Relu, Sigmoid, Tanh, Gelu, SiLU, Softmax, LeakyRelu, Elu, PRelu, HardSigmoid, HardSwish.
Shape: Reshape, Transpose, Concat, Flatten, Squeeze, Unsqueeze, Expand, Slice, Split, Tile, Pad, Shape, Cast, Identity, EyeLike, Range, CumSum.
Constants: Constant (TensorProto + scalar), ConstantOfShape (INT64/INT32/DOUBLE/FLOAT value).
Indexing: Gather, GatherND, ScatterElements (reduction none/add), ScatterND, TopK (any axis), NonZero.
Reduction: ReduceMean, ReduceSum, ReduceMax, ReduceMin.
Detection: NonMaxSuppression, RoiAlign.
Quantization: DequantizeLinear, QuantizeLinear, QLinearConv, QLinearMatMul, QLinearAdd, QLinearSigmoid, QLinearConcat — QLinearConv and QLinearMatMul with an exact integer accumulator that reproduces ONNX Runtime's int16 pair saturation.
Fused custom ops: RelPosBias2D (BoTNet-style 2D relative position bias).
Pass-through: Dropout.
ggml_custom() adds a graph node computed by a C kernel — for operations ggml
has no graph op for. Kernels are addressed by name, not by a function
pointer passed through R, so a wrong name is a clean error instead of a segfault.
Three kernels ship built in; ggml_custom_ops() lists the registry.
ctx <- ggml_init(16 * 1024 * 1024)
x <- ggml_new_tensor_2d(ctx, GGML_TYPE_F32, 5, 2)
ggml_set_f32(x, c(5, 1, 3, 2, 4,
9, 7, 8, 6, 10))
# One median per row — the output collapses ne[1] to 1
med <- ggml_custom(ctx, "row_median", args = list(x), ne = c(1, 2))
graph <- ggml_build_forward_expand(ctx, med)
ggml_graph_compute(ctx, graph)
ggml_get_f32(med) # 3, 8
ggml_free(ctx)| Kernel | What it does | API path it shows |
|---|---|---|
row_median |
median of each row | one input, caller-chosen output shape |
row_permute |
out[i] = x[perm[i]] within each row (ggml_get_rows() permutes whole rows) |
two inputs, I32 index tensor |
clip_inplace |
clamp into [lo, hi], writing into the input |
ggml_custom_inplace() |
Registering your own kernel. A package that links against ggmlR registers its kernels from C at load time, then references them from R by name:
#include <ggmlR.h> // system.file("include", "ggmlR.h", package = "ggmlR")
static void my_kernel(struct ggml_tensor * dst, int ith, int nth, void * ud) {
const struct ggml_tensor * a = dst->src[0];
/* write only the slice of dst that ith owns */
}
void R_init_mypkg(DllInfo * dll) {
ggmlR_register_custom_op_t reg = (ggmlR_register_custom_op_t)
R_GetCCallable("ggmlR", "ggmlR_register_custom_op");
reg("my_kernel", my_kernel);
}Custom nodes are CPU-only — the kernel is a host function pointer, so no GPU
backend can run it. Computing a graph that would put one on a GPU backend raises
an error rather than failing inside the backend. Kernels may run on any thread
(ith in [0, nth)) and must not call into R; pass n_tasks = 1 for a kernel
that is not thread-safe.
Low-level bindings for the recurrent alternatives to attention — linear in sequence length instead of quadratic.
| Op | Signature |
|---|---|
| Mamba convolution | ggml_ssm_conv(ctx, sx, c) |
| Mamba selective scan | ggml_ssm_scan(ctx, s, x, dt, A, B, C, ids) |
| RWKV-6 | ggml_rwkv_wkv6(ctx, k, v, r, tf, td, state) |
| RWKV-7 | ggml_rwkv_wkv7(ctx, r, w, k, v, a, b, state) |
| Gated linear attention | ggml_gated_linear_attn(ctx, k, v, q, g, state, scale) |
ggml_ssm_scan() and the three RWKV-family ops each return one tensor with
the sequence output and the final recurrent state concatenated. View either half
rather than computing byte offsets by hand:
r <- ggml_ssm_scan(ctx, s, x, dt, A, B, C, ids)
y <- ggml_ssm_scan_output(ctx, r, x) # shaped like x
state <- ggml_ssm_scan_state(ctx, r, s, ids) # [d_state, head_dim, n_head, n_seqs]
w <- ggml_rwkv_wkv6(ctx, k, v, rr, tf, td, st)
out <- ggml_rwkv_output(ctx, w, k) # [S*H, n_tokens]
st2 <- ggml_rwkv_state(ctx, w, k, st) # [S*H, S*n_seqs]Trainable — upstream ggml has no backward pass for any of these, which made
them inference-only; ggmlR adds one for all five, so a Mamba or RWKV block
trains end to end. ggml_ssm_conv() and ggml_ssm_scan() also have backward
Vulkan shaders, so a state-space block trains entirely on the GPU (d_state
128 or 256, Mamba-2 shapes; other shapes compute the backward on the CPU). The
RWKV and GLA backward kernels are CPU-only, as is ggml_gated_linear_attn() in
the forward direction.
The same holds for flash attention: upstream's ggml_flash_attn_back() is a
stub that aborts, so ggml_flash_attn_ext() was inference-only. ggmlR
implements it on CPU and Vulkan and wires it into the autodiff, so attention
trains through the fused op. No ALiBi, logit softcap or attention sinks.
Convolution backward (GGML_OP_IM2COL_BACK) also has a Vulkan shader here,
which upstream lacks. Without it the scheduler moved the backward — and the
nodes feeding it — to the CPU, and a convolution trained on the GPU came out
slower than the same model trained on the CPU outright; with it the graph
stays in one piece on the device. Training a convolution requires an F32
kernel: an F16 one, the usual choice for inference, leaves the backward with no
backend that can run it, and ggmlR says so when the gradient is requested
rather than aborting later inside the scheduler.
Embedding backward (GGML_OP_GET_ROWS_BACK) has a Vulkan shader here as well,
so a trained embedding table stays on the device. On a 30000-word vocabulary
that puts GPU training 1.5–2.5x ahead of the CPU, widening with depth.
See inst/examples/mamba_train_demo.R for a block trained end to end on both
backends, and inst/examples/backward_gpu_demo.R for which backward ops run
where.
Load pre-trained weights from GGUF files (llama.cpp, Hugging Face, etc.) with automatic dequantization. Supports all ggml quantization types (F32, F16, Q4_0, Q8_0, K-quants, IQ, MXFP4, Q1_0, NVFP4).
library(ggmlR)
# Load a GGUF file
g <- gguf_load("model.gguf")
g
#> GGUF file: model.gguf
#> Version: 3
#> Tensors: 291
#> Metadata: 24 key-value pairs
# Inspect metadata (architecture, tokenizer, quant info)
meta <- gguf_metadata(g)
meta[["general.architecture"]]
# List all tensor names
gguf_tensor_names(g)
# Get shape and type for a specific tensor
gguf_tensor_info(g, "blk.0.attn_q.weight")
#> $name: "blk.0.attn_q.weight"
#> $shape: 4096 4096
#> $type: "Q4_0"
# Extract dequantized weights as R numeric array
w <- gguf_tensor_data(g, "blk.0.attn_q.weight")
dim(w)
#> [1] 4096 4096
# Free when done (also freed by GC)
gguf_free(g)Ready-to-run example scripts in inst/examples/:
| Script | Description |
|---|---|
titanic_classification.R |
Binary classification with feature engineering, dropout, stratified split, manual metrics (~82% val accuracy) |
mnist_cnn.R |
CNN image classifier on MNIST |
functional_resnet_cifar.R |
ResNet-style model with skip connections (Functional API) |
functional_text_gru.R |
Text classification with GRU + embedding (Functional API) |
transformer_encoder_demo.R |
Transformer encoder with multi-head attention (autograd) |
dp_train_demo.R |
Data-parallel training across multiple replicas |
benchmark_onnx.R |
GPU vs CPU inference benchmark for ONNX models |
benchmark_ops.R |
Per-op micro-benchmark: every ggml op on CPU and GPU with auto-batching |
profile_onnx_superres_gpu.R |
GPU profiler for SuperResolution ONNX model across input sizes |
mlr3_integration.R |
mlr3 learners: CPU vs GPU comparison, iris + mtcars, 3-fold CV |
tidymodels_integration.R |
parsnip engine: CPU vs GPU comparison, iris + mtcars |
ggml_save_model(model, "my_model.rds")
model <- ggml_load_model("my_model.rds")Measured on AMD Ryzen 5 5600 + AMD RX 9070 (Resizable BAR off), single-image inference, mean of 3 runs:
| Model | CPU (ms) | GPU (ms) | Speedup | CPU FPS | GPU FPS |
|---|---|---|---|---|---|
| Inception V3 | 205.0 | 8.7 | 23.6x | 4.9 | 115.4 |
| MNIST | 0.0 | 0.3 | — | — | 3000.0 |
| SqueezeNet 1.0 | 21.7 | 2.0 | 10.8x | 46.2 | 500.0 |
| SuperResolution | 83.7 | 3.0 | 27.9x | 12.0 | 333.3 |
| EmotionFerPlus | 27.0 | 2.0 | 13.5x | 37.0 | 500.0 |
| Inception V3 Op18 | 200.3 | 11.3 | 17.7x | 5.0 | 88.2 |
| BAT-ResNeXt26ts | 91.7 | 6.3 | 14.5x | 10.9 | 157.9 |
| BERT (Opset17) | 222.0 | 12.0 | 18.5x | 4.5 | 83.3 |
| GPT-NeoX | 1.7 | 3.7 | 0.45x | 600.0 | 272.7 |
| RoBERTa SeqClass | 230.7 | 22.0 | 10.5x | 4.3 | 45.5 |
| CaiT XS24 | 624.0 | 39.7 | 15.7x | 1.6 | 25.2 |
| XCiT Tiny12 P8 | 128.7 | 17.7 | 7.3x | 7.8 | 56.6 |
| BoTNet26t | 122.0 | 4.3 | 28.2x | 8.2 | 230.8 |
| SAGEConv | 110.7 | 12.7 | 8.7x | 9.0 | 78.9 |
| MaskRCNN int8 | 1820.0 | 734.7 | 2.5x | 0.5 | 1.4 |
- MNIST runs below the timer's resolution on the CPU, so no speedup is given.
- GPT-NeoX here is a tiny test model (1.7 MB, ~2 ms on the CPU): the GPU's fixed per-run cost (submission, synchronisation, input upload) exceeds the work itself.
- GPU memory placement follows the device's memory map, not the driver name or version: without Resizable BAR, device buffers are kept in plain VRAM (see Diagnostics for the check and its overrides).
Benchmark scripts: inst/examples/benchmark_onnx.R, inst/examples/profile_onnx_superres_gpu.R
Measured on AMD Ryzen 5 5600 + AMD RX 9070, via llamaR. Model: Ministral-3-3B-Instruct-2512-Q8_0, 50 tokens, avg of 3 runs:
| Backend | Speed (tokens/sec) | Speedup |
|---|---|---|
| CPU (8 threads) | 8.5 | 1.0x |
| GPU (Vulkan) | 108.0 | 12.7x |
Measured on AMD Ryzen 5 5600 + AMD RX 9070, via sd2R:
| Scenario | Resolution | Time (s) |
|---|---|---|
| Direct | 768×768 | 13.94 |
| Tiled VAE | 1024×1024 | 25.32 |
| Highres fix | 2048×1024 | 52.53 |
| img2img | 768×768 | 8.73 |
| Direct | 1024×1024 | 25.40 |
ggmlR is designed GPU-first: Vulkan is auto-detected at build time and, when available, used automatically — on the ONNX benchmark above, 7–28x over the CPU for models large enough to keep the GPU busy (see that section for the exceptions). On machines without a Vulkan-capable GPU the package falls back to CPU transparently — no code changes required.
ggml_vulkan_available() # TRUE if a Vulkan GPU was detected
ggml_vulkan_status() # device name, memory, capabilities
# Dynamic autograd: switch device at runtime
ag_device("gpu") # move subsequent ops to GPU (f16 by default)
ag_device("cpu") # fall back to CPUSupported GPUs: NVIDIA, AMD, Intel, ARM Mali, Qualcomm Adreno.
- Vulkan 1.4 — push constants limit raised to 256 bytes, enabling full 5D tensor parameter blocks in compute shaders without staging buffers.
- Push Descriptors (
VK_KHR_push_descriptor) — when available, descriptors are pushed directly into the command buffer, eliminating descriptor pool allocation overhead. Falls back to descriptor pools on older hardware. - Q4_K flash attention —
GGML_OP_FLASH_ATTN_EXTwith Q4_K key/value tensors now runs fully on GPU (FA_SCALAR and FA_COOPMAT1 paths). Previously Q4_K attention fell back to CPU. Relevant for llamaR with quantized LLMs on AMD/Intel GPU (KHR cooperative matrix). - Subgroup-shuffle mmq (
USE_SUBGROUP_NO_SHMEM) — on wavefront-64 devices (RDNA4, subgroup_size=64) Q4_K / Q5_K / Q6_K weight tiles are loaded directly into registers viasubgroupShuffle, eliminating the shared-memory staging round-trip. ~10-15% token-generation throughput gain on LLaMA 3.x models.
- R ≥ 4.1.0, C++17 compiler
- Optional GPU: Vulkan 1.2+,
libvulkan-dev+glslc(Linux) or Vulkan SDK (Windows) - Platforms: Linux, macOS, Windows (x86-64, ARM64)
MIT