ggmlR — Neural Networks for R

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:

Also serves as the backend engine for llamaR (LLM inference) and sd2R (Stable Diffusion).

Installation

install.packages("ggmlR")

GPU (Vulkan) support is auto-detected at build time.

Ubuntu / Debian — to enable GPU:

sudo apt install libvulkan-dev glslc

Windows — install Rtools and optionally the Vulkan SDK for GPU support.

Build options

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")

Linux (detailed)

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-dev

2. 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-sdk

3. Verify the GPU is visible to Vulkan:

vulkaninfo --summary

4. 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()'

Windows build options

Important: on Windows, R ignores configure.args / --configure-args (they are only honoured by the Unix ./configure path). 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.

Diagnostics

GGMLR_LOG_DEVICE — when set to a non-empty value other than 0, enables Vulkan telemetry logging. Two messages are emitted:

1. 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

2. 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.R checks the backward ops one at a time.

3. 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 overrides

   The 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=1 keeps the window (the previous behaviour), GGML_VK_DISABLE_HOST_VISIBLE_VIDMEM=1 never 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")

CPU performance: multithreaded BLAS/LAPACK

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 with OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 for a clean valgrind report; CRAN’s own machines use single-threaded reference BLAS, so it does not appear there.

GPU linear algebra (drop-in for %*%)

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 with device = "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-3 relative 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.

Sequential API

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)

Available layers (Sequential)

Available losses

"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.

CNN example (MNIST)

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")

Training from a generator

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.

Your own training loop

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$snapshot()                 # independent copy: the target network
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) tr$update_target(target)       # hard copy (DQN)
  # or every step: tr$update_target(target, tau = 0.005)   # Polyak (DDPG, SAC)
}

ggml_snapshot_free(target)
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.

Do not emulate this by calling ggml_fit() once per batch: every call builds a fresh optimizer context, so the Adam moments restart and each step is the optimizer’s first.

Sequential models and functional models with one input and one output are supported.

Masked losses. With loss_mask, the loss becomes a weighted MSE and $step() takes a sample_weight. loss_mask = 1 gives one weight per sample; loss_mask = <outputs> gives one per output element, so a 0/1 matrix trains only the outputs whose label is known — one action of a Q-network, say. The masked terms and their gradients are removed exactly; the loss is divided by the total element count, not by the number of active ones.

tr <- ggml_trainer(q_net, batch_size = 512L, loss_mask = 18L)  # 18 actions
tr$step(x, y, sample_weight = mask)     # mask: [512, 18], 1 where y is known

Predicting during training. $predict(x) runs the current weights in eval mode (no dropout) on a scheduler of its own, so it never disturbs the training graph or the optimizer state, and x may have any number of rows. Use it for validation inside the loop rather than ggml_predict(tr$model(), x).

Target networks. $snapshot() returns an independent model — its own weights in its own buffer, its own scheduler — that later steps do not move and that outlives $free(). $update_target(target, tau) moves it towards the current weights, tau * theta + (1 - tau) * theta_target: tau = 1 is the hard copy, a small tau the Polyak average. Batch-norm running statistics are copied whole, never averaged, and a target with different weights is refused. Release a snapshot with ggml_snapshot_free().

Functional API

Wire layers into arbitrary graphs — residual connections, multi-input/output, shared weights.

Residual (skip) connection

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)

Embedding + GRU + skip connection (NLP)

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].

Multi-input model

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))

Multi-output model

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]

Training several heads at once

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_loss

loss 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 / transformer encoder block

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 2

Note 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.

ResNet-like image classifier

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)

Shared layers (Siamese / weight sharing)

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)

Differences from Keras

evaluate() masks the keras3 generic

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)

Dynamic Autograd Engine (PyTorch-style)

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()

Transformer encoder block

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()
}

Data-parallel training

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 0

Multi-GPU: tensor & pipeline parallelism

For 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.

Benchmark: which split strategy to use

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

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:

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 plain ggml_vulkan_shutdown() (no hard) 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 build hard = TRUE does 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 with ggml_vulkan_hard_exit_available().

When to use it — and when not. hard = TRUE only 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 need hard = 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 as DeadKernelError: Kernel died even though all results were already computed and written. In a notebook, either drop the call or use plain ggml_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)

Autograd op reference

Fused attention

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):

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.

Gradient checkpointing

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:

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 as a single graph

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 breakdown

GGMLR_AG_BWD_PROF=1 turns the same per-stage timings on inside any script; read them with ggmlR:::ag_backward_profile_report().

When the GPU is worth it on the ag_* path

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").

⚠️ The table below predates device residency and has not been re-measured. It was taken when every 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):

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. Set ag_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 goes

This 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.

Keeping tensors on the device

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 reset with_grad_tape() performs at every step;
• optimizer moments — Adam’s m and v live 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=0

Numerically 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 now

ag_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 shape

Device buffers

A 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 R rows x capacity view; the name follows ggml_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 to 0..capacity-1.
• Inside ag_capture() pass the indices as an argument and list the buffer in params: 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.

Budgeting a training run

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.

mlr3 Integration

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"))

Features

• Both ggmlR APIs — model_fn can return a ggml_sequential_model or ggml_functional_model. The default builder is ggml_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 to future / callr workers without file-system round-trips.
• Weighted training — classif.ggml honours task$weights_learner, mapping them to sample_weight in ggml_fit().
• Callbacks for tuning — pass ggml_callback_early_stopping() etc. via the callbacks parameter to drive early stopping inside mlr3 tuning 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.

tidymodels / parsnip Integration

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.

Single-cell GPU Acceleration (Seurat)

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 path

If 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:

"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.

What runs where

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:

On a Vulkan GPU (AMD RADV) the GPU steps are markedly faster than a naive reference. Indicative numbers at 2000 cells:

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):

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:

# 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 + description

A SingleCellExperiment (Bioconductor) path with an S4 runGGML() generic is planned next.

ONNX Model Import

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.

Quick start

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")

Loading models

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")

Dynamic shapes

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.

Inspecting the model

# 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)

Running inference

# 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 ...

Preparing input data

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)))

Interpreting outputs

# 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))

Repeated inference

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...
}

Tested models

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.

Supported ONNX ops (90+)

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.

Custom Operations

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)

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.

State-Space (Mamba) and RWKV Ops

Low-level bindings for the recurrent alternatives to attention — linear in sequence length instead of quadratic.

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.

GGUF Pre-trained Weights

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)

Examples

Ready-to-run example scripts in inst/examples/:

Save / Load

ggml_save_model(model, "my_model.rds")
model <- ggml_load_model("my_model.rds")

File format. ggml_save_model() writes format 3, which adds to format 2 a functional model’s ggml_freeze_weights() marks, the compiled weight_decay and predict_precision, and fixed batch-norm statistics. ggml_load_model() reads formats 2 and 3 and refuses anything newer with a message to update ggmlR. Releases before 0.8.7 accept any format from 2 up, so they load a format-3 file and ignore its frozen marks and fixed batch-norm statistics without a word: load format-3 files with 0.8.7 or later.

Checks on load. A truncated or non-RDS file stops at ggml_load_model(). Every saved tensor is compared by name and size with the one it is about to fill, before anything is written. A file from another architecture, or a damaged one, then stops with the layer and tensor named, instead of being copied in partly or falling back to random init. This check runs when the graph is first built: the first ggml_predict(), ggml_fit() or ggml_trainer() call.

ggml_freeze_weights() and ggml_unfreeze_weights() on anything other than a ggml_* model (an ag_* module, say) stop with a message naming what is supported.

Fine-tuning a saved model

ggml_load_model() returns a compiled model with the saved weights, so training goes on from where it stopped with either ggml_fit() or ggml_trainer():

model <- ggml_load_model("my_model.rds", backend = "auto")

# a few more epochs on new data
model <- ggml_fit(model, x_new, y_new, epochs = 5L, batch_size = 32L)

# or your own loop
tr <- ggml_trainer(model, batch_size = 512L)
tr$set_lr(5e-4)
for (i in seq_len(steps)) tr$step(x_batch(i), y_batch(i))
model <- tr$model()
tr$free()

ggml_save_model(model, "my_model_ft.rds")

• Only the weights are saved, not the optimizer state. The Adam moments start from zero, so the first steps after loading are Adam’s first steps. A lower rate or a short warm-up keeps them from knocking the weights about.
• The architecture comes from the file, along with the optimizer, loss and metrics given to ggml_compile(). Inputs and outputs must keep the shapes the model was saved with.
• Files written by ggml_save_weights() are not accepted by ggml_load_model(). Rebuild and compile the same architecture, then call ggml_load_weights(model, path).

Freezing layers. ggml_freeze_weights() marks layers as not trainable, and ggml_unfreeze_weights() undoes it. The flag is read when the training graph is built, so it can be set on a loaded model before ggml_fit() or ggml_trainer():

model <- ggml_load_model("my_model.rds")
model <- ggml_freeze_weights(model, from = 1L, to = 3L)  # sequential: by index
# model <- ggml_freeze_weights(model, layer_names = c("enc1", "enc2"))  # by name
model <- ggml_fit(model, x_new, y_new, epochs = 5L)

prefix selects every layer whose name starts with it – a whole block, a head, the embeddings – and combines with layer_names; a name or prefix that matches nothing is an error. For functional models only names and prefixes select layers; without them every layer is frozen. Weight decay (weight_decay in ggml_compile() or ggml_trainer()) reaches trainable weights only, so frozen layers are not decayed either. ggml_save_model() keeps the frozen marks for both kinds of model, so a reloaded model trains the same layers as before. Functional files saved before this was added load with nothing frozen.

What freezing reaches:

Layers without weights (activations, dropout, pooling, add, flatten) have nothing to freeze.

Batch-norm statistics. Freezing a sequential batch-norm stops gamma and beta only: in training it still normalizes each batch by its own mean and variance, and ggml_fit() still recalibrates its running statistics. ggml_freeze_bn_stats() fixes them – the layer then uses the stored statistics in training as at prediction, and they are no longer recalibrated – while gamma and beta keep training unless frozen too. That is the usual setting for fine-tuning on small batches; ggml_unfreeze_bn_stats() undoes it. The state is saved with the model. Fix them after the statistics have been calibrated by a fit. A functional batch_norm is per-sample RMS normalization with no statistics, and the call says so.

model <- ggml_freeze_bn_stats(model)                     # all batch-norm layers
model <- ggml_freeze_weights(model, prefix = "block1_")  # and, if wanted, weights

ONNX Benchmark: GPU (Vulkan) vs CPU

Measured on AMD Ryzen 5 5600 + AMD RX 9070 (Resizable BAR off), single-image inference, mean of 3 runs:

• 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

LLM Inference Benchmark

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:

Flux Image Generation Benchmark (10 steps)

Measured on AMD Ryzen 5 5600 + AMD RX 9070, via sd2R:

GPU Acceleration

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 CPU

Supported GPUs: NVIDIA, AMD, Intel, ARM Mali, Qualcomm Adreno.

Vulkan optimizations

• 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_EXT with 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 via subgroupShuffle, eliminating the shared-memory staging round-trip. ~10-15% token-generation throughput gain on LLaMA 3.x models.

Precision of matrix multiplies on Vulkan

On a Vulkan device with fp16 support, MUL_MAT keeps its tiles in f16 by default. A single row goes through an exact f32 matrix-vector kernel, larger batches through the tiled one, so predictions carry a relative error around 1e-3 that depends on the batch size: the same input can come out slightly differently alone and in a batch of 64. Training is not hurt by noise of that size; comparisons against a reference, tests and reproducible evaluation are.

model <- ggml_compile(model, optimizer = "adam", loss = "mse",
                      predict_precision = "f32")      # every prediction
ggml_predict(model, x, matmul_precision = "f32")      # one call
tr$predict(x, matmul_precision = "f32")              # a trainer's predict

"f32" marks every multiply GGML_PREC_F32; the result then matches the CPU to about 1e-6. The argument of the call wins over the model’s predict_precision, which wins over "default". Training graphs are not affected. The throughput cost depends on the device and driver: on an RX 9070 under RADV the f32 tiles ran at about a seventh of the f16 ones for large multiplies (inst/examples/benchmark_coopmat.R measures it on yours). The CPU backend is exact either way.

System Requirements

• 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)

See Also

• llamaR — LLM inference in R
• sd2R — Stable Diffusion in R
• ggml — underlying C library

License

MIT