Keras-like API in ggmlR

ggmlR provides a high-level Keras-like interface for building and training neural networks on top of the ggml C tensor library. This vignette walks through the main APIs using the built-in iris and mtcars datasets.

Installation

# install.packages("ggmlR")   # once on CRAN
library(ggmlR)

1. Sequential API

The sequential API mirrors keras.Sequential: stack layers with |>, compile, fit, evaluate, predict.

1.1 Prepare data — iris (3-class classification)

data(iris)
set.seed(42)

x <- as.matrix(iris[, 1:4])
x <- scale(x)                        # standardise

# one-hot encode species
y <- model.matrix(~ Species - 1, iris)  # [150 x 3]

idx   <- sample(nrow(x))
n_tr  <- 120L
x_tr  <- x[idx[1:n_tr],  ]
y_tr  <- y[idx[1:n_tr],  ]
x_val <- x[idx[(n_tr+1):150], ]
y_val <- y[idx[(n_tr+1):150], ]

1.2 Build and compile

model <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_batch_norm() |>
  ggml_layer_dropout(0.3, stochastic = TRUE) |>
  ggml_layer_dense(16L, activation = "relu") |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(
    optimizer = "adam",
    loss      = "categorical_crossentropy",
    metrics   = c("accuracy")
  )

print(model)

Available losses:

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.

MAE and Huber matter when the targets contain outliers: MSE squares the error, so one wild value dominates the gradient, while MAE weights every point equally and Huber is quadratic only near zero. Binary cross-entropy treats each output as its own Bernoulli, which is what makes it the loss for multi-label targets rather than for one-of-N classification.

1.3 Train

model <- ggml_fit(
  model,
  x_tr, y_tr,
  epochs           = 100L,
  batch_size       = 32L,
  validation_split = 0.0,   # we supply our own val set below
  verbose          = 0L
)

1.4 Evaluate and predict

score <- ggml_evaluate(model, x_val, y_val, batch_size = 16L)
cat(sprintf("Val loss: %.4f   Val accuracy: %.4f\n", score$loss, score$accuracy))

probs   <- ggml_predict(model, x_val, batch_size = 16L)
classes <- apply(probs, 1, which.max)
true    <- apply(y_val, 1, which.max)
cat("Confusion matrix:\n")
print(table(true = true, predicted = classes))

1.5 Save and reload

path <- tempfile(fileext = ".rds")
ggml_save_model(model, path)
cat(sprintf("Saved: %s (%.1f KB)\n", path, file.size(path) / 1024))

model2 <- ggml_load_model(path)
score2 <- ggml_evaluate(model2, x_val, y_val, batch_size = 16L)
cat(sprintf("Reloaded model accuracy: %.4f\n", score2$accuracy))

Only the weights are saved, not the optimizer state. On load, every saved tensor is checked by name and size against the layer it fills, so a file from another architecture stops with the layer named. Files are written in format 3, and ggml_load_model() refuses a format newer than it knows.

1.6 Embedding as the first layer

ggml_layer_embedding() can open a sequential model: give it input_shape = the sequence length, and the input is a matrix of token indices, one row per sample. Indices are 0-based, 0 .. vocab_size - 1; NA, fractions and out-of-vocabulary indices are errors, not silently clamped.

set.seed(3)
vocab <- 20L; seq_len <- 6L
tok <- matrix(sample(0:(vocab - 1L), 200L * seq_len, replace = TRUE),
              200L, seq_len)
lab   <- as.integer(rowMeans(tok < 10L) > 0.5)   # mostly low tokens?
y_tok <- cbind(1 - lab, lab) * 1.0

model_emb <- ggml_model_sequential() |>
  ggml_layer_embedding(vocab_size = vocab, dim = 8L, input_shape = seq_len) |>
  ggml_layer_flatten() |>
  ggml_layer_dense(2L, activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

model_emb <- ggml_fit(model_emb, tok, y_tok,
                      epochs = 30L, batch_size = 20L, verbose = 0L)
cat(sprintf("embedding model accuracy: %.4f\n",
            tail(model_emb$history$train_accuracy, 1)))

1.7 Training from a generator

With generator =, ggml_fit() pulls batches from a function instead of slicing a matrix held in memory, so only the current batch is resident. The function takes no arguments and returns list(x, y) with exactly batch_size rows, or NULL to end the epoch. A finite generator has to rewind itself before returning NULL, or the next epoch gets nothing:

make_gen <- function(x, y, bs) {
  pos <- 0L
  function() {
    if (pos + bs > nrow(x)) { pos <<- 0L; return(NULL) }   # rewind, end epoch
    i <- pos + seq_len(bs)
    pos <<- pos + bs
    list(x[i, , drop = FALSE], y[i, , drop = FALSE])
  }
}

model_gen <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

model_gen <- ggml_fit(model_gen, generator = make_gen(x_tr, y_tr, 24L),
                      epochs = 50L, batch_size = 24L, verbose = 0L)

A generator that never returns NULL — augmentation, a simulator — needs steps_per_epoch to say where an epoch ends:

sampler <- function() {
  i <- sample(nrow(x_tr), 24L)
  list(x_tr[i, , drop = FALSE], y_tr[i, , drop = FALSE])
}

model_gen <- ggml_fit(model_gen, generator = sampler, steps_per_epoch = 5L,
                      epochs = 50L, batch_size = 24L, verbose = 0L)
cat(sprintf("generator model, val accuracy: %.4f\n",
            ggml_evaluate(model_gen, x_val, y_val, batch_size = 10L)$accuracy))

validation_generator and validation_steps validate the same way. The arguments that need the whole data set — validation_split, validation_data, class_weight, sample_weight — are rejected with a generator rather than ignored.

1.8 Your own training loop: ggml_trainer()

When progress is not counted in passes over a data set — reinforcement learning, curricula — ggml_trainer() returns a context that keeps the weights, the Adam moments and the graph alive, and you step it yourself. $step(x, y) runs one forward and backward pass and returns the loss; $eval(x, y) measures without training. Both take batches of exactly batch_size rows; $predict(x) takes any number.

A step normalizes each batch by its own statistics and does not update the running ones a sequential batch norm predicts with, so a model with batch norm warns when the trainer is created. $calibrate_bn(x) sets them from representative inputs; call it before predicting, snapshotting or saving.

set.seed(5)
model_tr <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_batch_norm() |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

tr <- ggml_trainer(model_tr, batch_size = 30L)

for (i in seq_len(300L)) {
  b    <- sample(nrow(x_tr), 30L)
  loss <- tr$step(x_tr[b, , drop = FALSE], y_tr[b, , drop = FALSE])
}
cat(sprintf("last step loss: %.4f   val loss: %.4f\n",
            loss, tr$eval(x_val, y_val)))

tr$calibrate_bn(x_tr)
p_tr <- tr$predict(x_val)
cat(sprintf("trainer val accuracy: %.4f\n",
            mean(apply(p_tr, 1, which.max) == apply(y_val, 1, which.max))))

For target networks, $snapshot() returns an independent model that later steps do not move, and $update_target(target, tau) moves it towards the current weights (tau = 1 is a hard copy, a small tau the Polyak average):

target <- tr$snapshot()
b <- sample(nrow(x_tr), 30L)
tr$step(x_tr[b, , drop = FALSE], y_tr[b, , drop = FALSE])
tr$update_target(target, tau = 0.1)
ggml_snapshot_free(target)

model_tr <- tr$model()   # weights back into a model object
tr$free()

With ggml_trainer(loss_mask = ) the loss becomes a weighted MSE and $step() takes a sample_weight, so only the outputs whose label is known are trained — one action of a Q-network, say.

1.9 Fine-tuning: freezing layers and batch-norm statistics

ggml_freeze_weights() marks layers as not trainable — by index (from, to), by layer_names, or by prefix, which selects every layer whose name starts with it (a whole block, say prefix = "block1_"). ggml_unfreeze_weights() undoes it. A name or prefix that matches no layer stops with an error rather than being ignored.

Freezing a batch-norm layer stops gamma and beta only. ggml_freeze_bn_stats() also fixes its running statistics: training then normalizes with the stored statistics, as prediction does, and they are no longer recalibrated — the usual setting for fine-tuning on small batches.

# `model` from 1.2: dense, batch norm, dropout, dense, dense
model_ft <- ggml_freeze_weights(model, from = 1L, to = 2L)
model_ft <- ggml_freeze_bn_stats(model_ft)
model_ft <- ggml_fit(model_ft, x_tr, y_tr,
                     epochs = 20L, batch_size = 32L, verbose = 0L)

tryCatch(ggml_freeze_weights(model, prefix = "encoder_"),
         error = function(e) cat("Error:", conditionMessage(e), "\n"))

1.10 Weight decay and prediction precision

weight_decay in ggml_compile() (and ggml_trainer()) is decoupled weight decay, AdamW-style, for "adam" and "sgd"; frozen weights are not decayed.

On a Vulkan device with fp16 support, matrix multiplies keep their tiles in f16 by default, so predictions carry a relative error around 1e-3 that depends on the batch size. predict_precision = "f32" makes every prediction match the CPU to about 1e-6; matmul_precision does the same for one ggml_predict() or $predict() call. Training is unaffected.

model_wd <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy",
               weight_decay = 0.01, predict_precision = "f32")

model_wd <- ggml_fit(model_wd, x_tr, y_tr,
                     epochs = 100L, batch_size = 32L, verbose = 0L)
p_wd <- ggml_predict(model_wd, x_val, batch_size = 10L)
p_16 <- ggml_predict(model_wd, x_val, batch_size = 10L,
                     matmul_precision = "default")   # this call only
cat(sprintf("max |f32 - default|: %.2e (exactly 0 on a CPU backend)\n", max(abs(p_wd - p_16))))

2. Functional API

The functional API lets you build DAGs with multiple inputs/branches/outputs. Use ggml_input() to declare inputs, pipe through layers, then call ggml_model().

2.1 Single-input DAG — residual-style network on iris

inp  <- ggml_input(shape = 4L)

h    <- inp |> ggml_layer_dense(32L, activation = "relu") |> ggml_layer_rms_norm()
h    <- h   |> ggml_layer_dense(16L, activation = "relu")
out  <- h   |> ggml_layer_dense(3L,  activation = "softmax")

model_fn <- ggml_model(inputs = inp, outputs = out) |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy",
               metrics = c("accuracy"))

model_fn <- ggml_fit(model_fn, x_tr, y_tr,
                     epochs = 100L, batch_size = 32L, verbose = 0L)

score_fn <- ggml_evaluate(model_fn, x_val, y_val, batch_size = 16L)
cat(sprintf("Functional model accuracy: %.4f\n", score_fn$accuracy))

Normalization in the functional API is per-sample RMS normalization (ggml_layer_rms_norm()), with no batch statistics. A functional ggml_layer_batch_norm() is the same layer under an older name, not the mean/variance batch norm of the sequential model in 1.2.

2.2 Multi-input model — mtcars regression

Predict mpg from two input groups: engine features and transmission features.

data(mtcars)
set.seed(7)

# Input group 1: engine (disp, hp, wt)
# Input group 2: transmission / gearbox (cyl, gear, carb, am)
engine <- as.matrix(scale(mtcars[, c("disp","hp","wt")]))
trans  <- as.matrix(scale(mtcars[, c("cyl","gear","carb","am")]))
y_mpg  <- matrix(scale(mtcars$mpg), ncol = 1L)          # [32 x 1]

# small dataset — use all for training, evaluate on same data for demo
x1 <- engine;  x2 <- trans

inp1 <- ggml_input(shape = 3L, name = "engine")
inp2 <- ggml_input(shape = 4L, name = "transmission")

branch1 <- inp1 |> ggml_layer_dense(16L, activation = "relu")
branch2 <- inp2 |> ggml_layer_dense(16L, activation = "relu")

merged  <- ggml_layer_add(list(branch1, branch2))        # element-wise add
out_reg <- merged |>
  ggml_layer_dense(8L, activation = "relu") |>
  ggml_layer_dense(1L)

model_reg <- ggml_model(inputs = list(inp1, inp2), outputs = out_reg) |>
  ggml_compile(optimizer = "adam", loss = "mse")

model_reg <- ggml_fit(model_reg,
                      x = list(x1, x2), y = y_mpg,
                      epochs = 200L, batch_size = 16L, verbose = 0L)

preds <- ggml_predict(model_reg, x = list(x1, x2), batch_size = 16L)
cat(sprintf("Pearson r (scaled mpg): %.4f\n", cor(preds, y_mpg)))

ggml_layer_add() requires the branches to have the same width. ggml_layer_concatenate() does not — it stacks them along an axis, keeping each branch’s features distinct instead of summing them. Both are trainable: gradients flow back into every branch.

br1 <- inp1 |> ggml_layer_dense(8L,  activation = "relu")
br2 <- inp2 |> ggml_layer_dense(12L, activation = "relu")   # widths differ

out_cat <- ggml_layer_concatenate(list(br1, br2), axis = 0L) |>
  ggml_layer_dense(1L)

model_cat <- ggml_model(inputs = list(inp1, inp2), outputs = out_cat) |>
  ggml_compile(optimizer = "adam", loss = "mse")
model_cat <- ggml_fit(model_cat, x = list(x1, x2), y = y_mpg,
                      epochs = 200L, batch_size = 16L, verbose = 0L)

p_cat <- ggml_predict(model_cat, x = list(x1, x2), batch_size = 16L)
cat(sprintf("Concatenate variant, Pearson r: %.4f\n", cor(p_cat, y_mpg)))

2.3 Multi-output model — one trunk, two heads

A model with several outputs trains all of them. Each head gets its own loss and weight, and the optimizer minimizes sum(loss_weights[i] * loss_i). Losses may be mixed: below a classification head and a regression head share one trunk.

Pass y as a list, one matrix per head, matched by output name.

set.seed(11)

# Two targets from the same cars: efficient or not (class), and mpg (value).
y_cls <- as.integer(mtcars$mpg > median(mtcars$mpg))
y_cls <- cbind(1 - y_cls, y_cls) * 1.0         # one-hot [32 x 2]
x_all <- as.matrix(scale(mtcars[, c("disp","hp","wt","cyl")]))

inp_mo <- ggml_input(shape = 4L, name = "features")
trunk  <- inp_mo |> ggml_layer_dense(16L, activation = "relu", name = "trunk")
head_c <- trunk  |> ggml_layer_dense(2L, activation = "softmax", name = "class")
head_v <- trunk  |> ggml_layer_dense(1L, name = "value")

model_mo <- ggml_model(inputs = inp_mo, outputs = list(head_c, head_v)) |>
  ggml_compile(optimizer = "adam",
               loss         = list(class = "categorical_crossentropy",
                                   value = "mse"),
               loss_weights = c(class = 1.0, value = 0.5))

model_mo <- ggml_fit(model_mo, x_all,
                     y = list(class = y_cls, value = y_mpg),
                     epochs = 150L, batch_size = 16L, verbose = 0L)

# Per-head losses, so a head that stops learning is visible on its own
# rather than hidden inside the total.
h <- model_mo$history
cat(sprintf("class loss: %.4f -> %.4f\n",
            h$train_class_loss[1], tail(h$train_class_loss, 1)))
cat(sprintf("value loss: %.4f -> %.4f\n",
            h$train_value_loss[1], tail(h$train_value_loss, 1)))

ggml_predict() returns a list, one matrix per head, and ggml_evaluate() reports a <output>_loss for each.

2.4 Attention and a transformer block

ggml_layer_attention() is multi-head scaled dot-product attention. It takes and returns a sequence node of shape c(seq_len, d_model), so blocks stack and a residual merge needs no reshaping.

time_distributed = TRUE applies one dense kernel at every position, sharing weights across the sequence — the position-wise feed-forward sublayer that completes a transformer block. (Without it a dense layer would flatten the sequence into a single vector.)

set.seed(13)
n_seq <- 128L; S <- 6L; D <- 8L

# A task that needs mixing across positions: is the sequence sum positive?
xs   <- array(runif(n_seq * S * D, -1, 1), dim = c(n_seq, S, D))
lab  <- as.integer(apply(xs, 1, mean) > 0)
y_sq <- cbind(1 - lab, lab) * 1.0

inp_t <- ggml_input(shape = c(S, D), name = "seq")

# Attention sublayer + residual
att <- inp_t |> ggml_layer_attention(d_model = D, n_heads = 2L, name = "mha")
h1  <- ggml_layer_add(list(inp_t, att))

# Position-wise feed-forward sublayer + residual
ff  <- h1 |> ggml_layer_dense(D * 2L, activation = "relu",
                              time_distributed = TRUE, name = "ff1")
ff  <- ff |> ggml_layer_dense(D, time_distributed = TRUE, name = "ff2")
h2  <- ggml_layer_add(list(h1, ff))

out_t <- h2 |> ggml_layer_flatten() |>
  ggml_layer_dense(2L, activation = "softmax", name = "cls")

model_t <- ggml_model(inputs = inp_t, outputs = out_t) |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

model_t <- ggml_fit(model_t, xs, y_sq, epochs = 30L, batch_size = 32L,
                    verbose = 0L)
cat(sprintf("transformer block accuracy: %.4f\n",
            tail(model_t$history$train_accuracy, 1)))

causal = TRUE masks keys after the query, making the layer autoregressive (GPT-style). Cross-attention takes queries from one node and keys/values from another, which may have a different length:

dec <- inp_t |> ggml_layer_attention(D, n_heads = 2L, causal = TRUE)

ctx_seq <- ggml_input(shape = c(10L, D))
x_cross <- ggml_apply(list(dec, ctx_seq), ggml_attention(D, n_heads = 2L))

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.

3. Callbacks

Callbacks plug into ggml_fit() via the callbacks argument.

3.1 Early stopping

cb_stop <- ggml_callback_early_stopping(
  monitor   = "val_loss",
  patience  = 15L,
  min_delta = 1e-4
)

model_cb <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

model_cb <- ggml_fit(model_cb, x_tr, y_tr,
                     epochs           = 300L,
                     batch_size       = 32L,
                     validation_split = 0.1,
                     callbacks        = list(cb_stop),
                     verbose          = 0L)

3.2 Learning-rate schedulers

# Cosine annealing
cb_cosine <- ggml_schedule_cosine_decay(T_max = 100L, eta_min = 1e-5)

# Step decay: halve LR every 30 epochs
cb_step <- ggml_schedule_step_decay(step_size = 30L, gamma = 0.5)

# Reduce on plateau
cb_plateau <- ggml_schedule_reduce_on_plateau(
  monitor  = "val_loss",
  factor   = 0.5,
  patience = 10L,
  min_lr   = 1e-6
)

model_lr <- ggml_model_sequential() |>
  ggml_layer_dense(32L, activation = "relu", input_shape = 4L) |>
  ggml_layer_dense(3L,  activation = "softmax") |>
  ggml_compile(optimizer = "adam", loss = "categorical_crossentropy")

model_lr <- ggml_fit(model_lr, x_tr, y_tr,
                     epochs           = 150L,
                     batch_size       = 32L,
                     validation_split = 0.1,
                     callbacks        = list(cb_cosine),
                     verbose          = 0L)

4. Autograd Engine (ag_*)

The autograd engine provides PyTorch-style dynamic computation graphs. Tensors are column-major: shape is [features, batch].

4.1 Prepare data (iris, col-major)

# ag_device()/ag_dtype() are package-level state, so a section that relies on a
# particular device has to say so rather than inherit whatever ran before it.
ag_device("cpu")
ag_dtype("f32")

# transpose: rows = features, cols = samples
x_tr_ag  <- t(x_tr)    # [4, 120]
y_tr_ag  <- t(y_tr)    # [3, 120]
x_val_ag <- t(x_val)   # [4, 30]
y_val_ag <- t(y_val)

4.2 Build a model with ag_sequential

ag_mod <- ag_sequential(
  ag_linear(4L,  32L, activation = "relu"),
  ag_batch_norm(32L),
  ag_dropout(0.3),
  ag_linear(32L, 16L, activation = "relu"),
  ag_linear(16L,  3L)
)

params <- ag_mod$parameters()
opt    <- optimizer_adam(params, lr = 1e-3)

4.3 Training loop

BS <- 32L
n  <- ncol(x_tr_ag)

ag_train(ag_mod)
set.seed(42)

for (ep in seq_len(150L)) {
  perm <- sample(n)
  for (b in seq_len(ceiling(n / BS))) {
    idx <- perm[seq((b-1L)*BS + 1L, min(b*BS, n))]
    xb  <- ag_tensor(x_tr_ag[, idx, drop = FALSE])
    yb  <- y_tr_ag[, idx, drop = FALSE]

    with_grad_tape({
      loss <- ag_softmax_cross_entropy_loss(ag_mod$forward(xb), yb)
    })
    grads <- backward(loss)
    opt$step(grads)
    opt$zero_grad()
  }

  if (ep %% 50L == 0L)
    cat(sprintf("epoch %d  loss %.4f\n", ep, as.matrix(loss)[1]))
}

4.4 Inference and accuracy

ag_eval(ag_mod)

# forward in chunks, apply softmax manually
ag_predict_cm <- function(model, x_cm, chunk = 64L) {
  n   <- ncol(x_cm)
  out <- matrix(0.0, nrow(as.matrix(model$forward(ag_tensor(x_cm[,1,drop=FALSE])))), n)
  for (s in seq(1L, n, by = chunk)) {
    e  <- min(s + chunk - 1L, n)
    lg <- as.matrix(model$forward(ag_tensor(x_cm[, s:e, drop = FALSE])))
    ev <- exp(lg - apply(lg, 2, max))
    out[, s:e] <- ev / colSums(ev)
  }
  out
}

probs_ag <- ag_predict_cm(ag_mod, x_val_ag)          # [3, 30]
pred_ag  <- apply(probs_ag, 2, which.max)
true_ag  <- apply(y_val_ag, 2, which.max)            # col-major: rows = classes
acc_ag   <- mean(pred_ag == true_ag)
cat(sprintf("Autograd val accuracy: %.4f\n", acc_ag))

4.5 LR scheduler + gradient clipping

ag_mod2 <- ag_sequential(
  ag_linear(4L,  64L, activation = "relu"),
  ag_batch_norm(64L),
  ag_dropout(0.3),
  ag_linear(64L, 32L, activation = "relu"),
  ag_linear(32L,  3L)
)
params2 <- ag_mod2$parameters()
opt2    <- optimizer_adam(params2, lr = 1e-3)
sch2    <- lr_scheduler_cosine(opt2, T_max = 150L, lr_min = 1e-5)

ag_train(ag_mod2)
set.seed(42)

for (ep in seq_len(150L)) {
  perm <- sample(n)
  for (b in seq_len(ceiling(n / BS))) {
    idx <- perm[seq((b-1L)*BS + 1L, min(b*BS, n))]
    xb  <- ag_tensor(x_tr_ag[, idx, drop = FALSE])
    yb  <- y_tr_ag[, idx, drop = FALSE]

    with_grad_tape({
      loss2 <- ag_softmax_cross_entropy_loss(ag_mod2$forward(xb), yb)
    })
    grads2 <- backward(loss2)
    clip_grad_norm(params2, grads2, max_norm = 5.0)
    opt2$step(grads2)
    opt2$zero_grad()
  }
  sch2$step()
}

ag_eval(ag_mod2)
probs2 <- ag_predict_cm(ag_mod2, x_val_ag)
acc2   <- mean(apply(probs2, 2, which.max) == true_ag)
cat(sprintf("ag + cosine + clip  val accuracy: %.4f\n", acc2))

4.6 Manual model from raw ag_param

For full control, build a network directly from parameters:

make_net <- function(n_in, n_hidden, n_out) {
  W1 <- ag_param(matrix(rnorm(n_hidden * n_in) * sqrt(2/n_in),  n_hidden, n_in))
  b1 <- ag_param(matrix(0.0, n_hidden, 1L))
  W2 <- ag_param(matrix(rnorm(n_out * n_hidden) * sqrt(2/n_hidden), n_out, n_hidden))
  b2 <- ag_param(matrix(0.0, n_out, 1L))

  list(
    forward    = function(x)
      ag_add(ag_matmul(W2, ag_relu(ag_add(ag_matmul(W1, x), b1))), b2),
    parameters = function() list(W1=W1, b1=b1, W2=W2, b2=b2)
  )
}

set.seed(1)
net    <- make_net(4L, 32L, 3L)
opt_r  <- optimizer_adam(net$parameters(), lr = 1e-3)

for (ep in seq_len(200L)) {
  perm <- sample(n)
  for (b in seq_len(ceiling(n / BS))) {
    idx <- perm[seq((b-1L)*BS+1L, min(b*BS, n))]
    xb  <- ag_tensor(x_tr_ag[, idx, drop = FALSE])
    yb  <- y_tr_ag[, idx, drop = FALSE]
    with_grad_tape({ loss_r <- ag_softmax_cross_entropy_loss(net$forward(xb), yb) })
    gr <- backward(loss_r)
    opt_r$step(gr);  opt_r$zero_grad()
  }
}

probs_r <- ag_predict_cm(net, x_val_ag)
acc_r   <- mean(apply(probs_r, 2, which.max) == true_ag)
cat(sprintf("Raw ag_param val accuracy: %.4f\n", acc_r))

5. GPU / device selection

# Use Vulkan GPU if available, fall back to CPU
device <- tryCatch({
  ag_device("gpu")
  "gpu"
}, error = function(e) "cpu")

cat("Running on:", device, "\n")

# Mixed precision (f16 on GPU, f32 on CPU)
ag_dtype(if (device == "gpu") "f16" else "f32")

After calling ag_device() and ag_dtype(), all subsequent ag_param and ag_tensor calls use the selected device and dtype. The sequential API (ggml_fit) picks up the backend automatically.

Both settings are package-level state that outlives this document, so put them back when done:

ag_device("cpu")
ag_dtype("f32")

API summary

Task Sequential / Functional Autograd
Build model ggml_model_sequential() + ggml_layer_* ag_sequential() or raw ag_param
Functional DAG ggml_input() → ggml_model() —
Multi-output outputs = list(...), loss = per head manual
Attention ggml_layer_attention(), ggml_attention() ag_multihead_attention()
Compile ggml_compile() manual optimizer
Train ggml_fit() (data or generator =), ggml_trainer() for your own loop with_grad_tape + backward + opt$step
Fine-tune ggml_freeze_weights(), ggml_freeze_bn_stats() —
Evaluate ggml_evaluate() manual loop
Predict ggml_predict() manual forward
Save / load ggml_save_model() / ggml_load_model() ag_save_model() / ag_load_model(), ag_state_dict() / ag_load_state_dict()
Callbacks ggml_callback_early_stopping(), ggml_schedule_* lr_scheduler_*, clip_grad_norm
Device automatic ag_device() / ag_dtype()