Sketched GLM: Temporal Sketches, Pooled AR, and Landmarks

Set up a synthetic dataset

We start with a small 8 x 8 x 3 volume and 120 time-points so the example runs quickly. A two-condition event design is convolved with the canonical HRF to produce the design matrix.

TR <- 2; Tlen <- 120; dim3 <- c(8L, 8L, 3L)
space4d <- NeuroSpace(c(dim3, Tlen))
maskVol <- LogicalNeuroVol(array(TRUE, dim3), NeuroSpace(dim3))

onsets  <- seq(20, 200, by = 30)
conds   <- factor(rep(c("A", "B"), length.out = length(onsets)))
events_df <- data.frame(onset = onsets, condition = conds, run = 1L)

sframe <- sampling_frame(blocklens = Tlen, TR = TR)
em <- event_model(onset ~ hrf(condition), data = events_df,
                  block = ~ run, sampling_frame = sframe)
X <- as.matrix(design_matrix(em))
p <- ncol(X); V <- prod(dim3)

Simulate data with AR(1) noise

The observed data is Y = X B + noise where the noise follows an AR(1) process with rho = 0.3. This gives the sketched solvers declared temporal dependence to estimate and whiten.

task_cols <- which(grepl("condition|hrf", colnames(X), ignore.case = TRUE))
B_true <- matrix(0, p, V)

# Landmark interpolation assumes spatial structure. Construct smooth task
# maps in the same voxel order used by the array-packing loop below.
coords_sim <- expand.grid(
  z = seq_len(dim3[3]), y = seq_len(dim3[2]), x = seq_len(dim3[1])
)[, c("x", "y", "z")]
xyz01 <- sweep(as.matrix(coords_sim), 2, 1, "-")
xyz01 <- sweep(xyz01, 2, pmax(dim3 - 1, 1), "/")
smooth_maps <- rbind(
  sin(pi * xyz01[, 1]) * cos(pi * xyz01[, 2]) * (0.5 + xyz01[, 3]),
  cos(pi * xyz01[, 1]) * sin(pi * xyz01[, 2]) * (1.0 - 0.5 * xyz01[, 3])
)
# Scale each map by its regressor's peak height so that the response peaks
# at 1.75 signal units per unit of map value, whatever the HRF's raw
# amplitude; the noise level below is chosen relative to that peak.
for (j in seq_along(task_cols)) {
  B_true[task_cols[j], ] <- smooth_maps[1L + (j - 1L) %% nrow(smooth_maps), ] *
    (1.75 / max(abs(X[, task_cols[j]])))
}

Y <- X %*% B_true + ar1_noise(Tlen, V, sd = 0.15)

Pack into an fmri dataset

We reshape the simulated matrix into a 4D NeuroVec for fmri_lm().

vec  <- NeuroVec(arr, space4d)
dset <- neurovec_frame(scans = list(vec), mask = maskVol,
                         TR = TR, event_table = events_df)

Choose the approximation engine

fmrireg has two built-in approximation paths with different controls. engine = "latent_sketch" is the sketched GLM engine used in most of this vignette; it is configured by lowrank_control(). engine = "rrr_gls" is the reduced-rank-regression GLS engine; it is configured through engine_args and is the path that exposes block-bootstrap task standard errors. The rrr_gls inference surface is deliberately narrow: standard errors and post-hoc contrasts apply to event/task coefficients, not baseline or nuisance terms.

Build parcels

Parcels are used for parcel-pooled AR whitening. The landmark extension below is coordinate-based and does not enforce parcel boundaries. Here we k-means cluster voxel coordinates solely to demonstrate AR pooling.

coords  <- expand.grid(x = seq_len(dim3[1]),
                       y = seq_len(dim3[2]),
                       z = seq_len(dim3[3]))
parcel_fit <- kmeans(
  coords, centers = 40, iter.max = 1000L, nstart = 10L,
  algorithm = "Lloyd"
)
parcels <- ClusteredNeuroVol(maskVol, parcel_fit$cluster)

Fit with Global AR + SRHT sketch

engine = "latent_sketch" activates the sketched path. lowrank_control() sets the sketch method and size; here we use SRHT with m = 8p rows. pooling = "global" estimates one shared AR(1) coefficient across the mask.

low_srht <- lowrank_control(
  parcels     = parcels,
  time_sketch = list(method = "srht", m = min(8L * p, Tlen))
)

fit_srht_global <- fmri_lm(
  onset ~ hrf(condition), block = ~ run,
  dataset    = dset,
  engine     = "latent_sketch",
  lowrank    = low_srht,
  control    = fmri_lm_control(noise = noise_spec("ar1", pooling = "global"))
)

Fit with Parcel AR + SRHT sketch

Switching to pooling = "parcel" estimates a separate AR coefficient per parcel, then shrinks them toward the global mean (shrink_c0 = 100). This stabilizes whitening in small or noisy parcels.

fit_srht_group <- fmri_lm(
  onset ~ hrf(condition), block = ~ run,
  dataset    = dset,
  engine     = "latent_sketch",
  lowrank    = low_srht,
  control    = fmri_lm_control(
    noise = noise_spec("ar1", pooling = "parcel", parcels = parcels,
                       shrink_c0 = 100L)
  )
)

Choosing a sketch

All three sketch-and-solve methods give calibrated inference; they differ in cost. In this engine the design is narrow (p of about 10-30 columns) and the response is wide (1e4-1e5 voxels or components), so exact least squares costs a single X'Z product. A sketch pays off only when applying it to Z is cheaper than that product. With p = 20, m = 8p and 5e4 response columns, measured on an Apple M3 Max with Accelerate BLAS:

Prefer "countsketch" or "gaussian" for speed. SRHT’s fast Walsh-Hadamard transform costs O(T log T) per column, so it trails both at these sizes; use it when you want its orthonormal rows. The sketches save more time relative to exact OLS as T grows.

The former method = "ihs" (iterative Hessian sketch) is deprecated. Each IHS iteration needs the exact gradient X'(Z - XM), a full X'Z pass, so in this regime it can never beat one exact solve. It now computes the exact OLS fit of the (whitened) data, with T - rank(X) residual degrees of freedom; without AR it reproduces the joint exact fit fmri_lm(..., strategy = "chunkwise"). It ignores iters and tol and says so once per session. For exact results, omit engine = "latent_sketch".

Compare sketching to an exact OLS baseline

The AR examples above illustrate whitening behavior, but they are not a like-for-like comparison to exact OLS because the estimators differ. For a direct sketch-fidelity check, we disable AR, fit the same model with and without sketching, and compare the task-effect coefficients only. In this moderate-SNR setting, correlations above 0.9 indicate that the sketch is preserving the task part of the GLM solution well.

fit_exact_iid <- fmri_lm(
  onset ~ hrf(condition), block = ~ run,
  dataset = dset
)

fit_srht_iid <- fmri_lm(
  onset ~ hrf(condition), block = ~ run,
  dataset = dset,
  engine = "latent_sketch",
  lowrank = low_srht
)

B_exact_iid <- t(fit_exact_iid$result$betas$data[[1]]$estimate[[1]])

cor_srht_iid_task <- cor(
  as.numeric(B_exact_iid[task_cols, , drop = FALSE]),
  as.numeric(fit_srht_iid$betas_fixed[task_cols, , drop = FALSE])
)
cor_srht_global_task <- cor(
  as.numeric(B_exact_iid[task_cols, , drop = FALSE]),
  as.numeric(fit_srht_global$betas_fixed[task_cols, , drop = FALSE])
)
cor_srht_group_task <- cor(
  as.numeric(B_exact_iid[task_cols, , drop = FALSE]),
  as.numeric(fit_srht_group$betas_fixed[task_cols, , drop = FALSE])
)

stopifnot(cor_srht_iid_task > 0.9)

list(
  cor_srht_iid_task = cor_srht_iid_task,
  cor_srht_global_task = cor_srht_global_task,
  cor_srht_group_task = cor_srht_group_task
)
#> $cor_srht_iid_task
#> [1] 0.9687722
#> 
#> $cor_srht_global_task
#> [1] 0.9740778
#> 
#> $cor_srht_group_task
#> [1] 0.9642624

The exact iid-versus-sketch comparison is the meaningful fidelity check here, and it should stay high. The AR-whitened fits usually track the same task effects reasonably well, but they are solving a different estimation problem, so their correlations to exact OLS are not expected to match the iid sketch baseline.

Standard errors follow the estimator

A sketch-and-solve fit uses only m of the information in T scans, so its coefficients vary more than the full-data ones: by about sqrt(T / m) in standard deviation. Its reported standard errors are the conditional-on-sketch standard errors of that estimator, and its residual degrees of freedom are those of the sketched residual (about m - p), so t-statistics and p-values stay calibrated. The same holds for every sketch-and-solve method.

se_task <- function(fit) {
  as.matrix(standard_error(fit, type = "estimates"))[, seq_along(task_cols)]
}
sketch_fits <- list(srht = fit_srht_iid)
for (method in c("gaussian", "countsketch")) {
  sketch_fits[[method]] <- fmri_lm(
    onset ~ hrf(condition), block = ~ run,
    dataset = dset,
    engine = "latent_sketch",
    lowrank = lowrank_control(
      time_sketch = list(method = method, m = min(8L * p, Tlen))
    )
  )
}
se_ratio <- vapply(sketch_fits, function(fit) {
  mean(se_task(fit) / se_task(fit_exact_iid))
}, numeric(1))

stopifnot(all(se_ratio > 1))

data.frame(
  method = names(sketch_fits),
  sketch_rows = vapply(sketch_fits, function(f) f$sketch$m, numeric(1)),
  residual_df = vapply(sketch_fits, function(f) f$result$rdf, numeric(1)),
  se_ratio_to_exact = se_ratio,
  sqrt_T_over_m = sqrt(Tlen / vapply(sketch_fits, function(f) f$sketch$m, numeric(1))),
  row.names = NULL
)
#>        method sketch_rows residual_df se_ratio_to_exact sqrt_T_over_m
#> 1        srht          16    9.956096          3.356971      2.738613
#> 2    gaussian          16    9.358252          3.165295      2.738613
#> 3 countsketch          16    9.356405          3.195592      2.738613

Reduced-rank GLS with bootstrap standard errors

Use rrr_gls when you want a low-rank task fit with task-level standard errors from either the conditional analytic approximation or a residual block bootstrap. The bootstrap is useful when short-TR data leave temporal autocorrelation that ordinary iid standard errors do not represent well.

fit_rrr_boot <- fmri_lm(
  onset ~ hrf(condition), block = ~ run,
  dataset = dset,
  control = fmri_lm_control(noise = noise_spec("ar1")),
  engine = "rrr_gls",
  engine_args = list(
    rank_mode = "energy",
    energy_keep = 0.99,
    se_mode = "bootstrap",
    bootstrap_n = 20L,
    bootstrap_block_size = 8L,
    bootstrap_seed = 1L
  )
)

For a real analysis, increase bootstrap_n substantially; 200 is a reasonable starting point. The vignette uses a smaller value so the example remains quick.

rrr_se <- standard_error(fit_rrr_boot, type = "betas")

data.frame(
  engine = attr(fit_rrr_boot, "engine"),
  se_mode = fit_rrr_boot$rrr$se_mode,
  rank_used = fit_rrr_boot$rrr$rank_used,
  finite_standard_errors = all(is.finite(as.matrix(rrr_se)))
)
#>    engine   se_mode rank_used finite_standard_errors
#> 1 rrr_gls bootstrap         2                   TRUE

Landmark extension (full-voxel, global AR)

Here we isolate the spatial approximation: both reference and landmark fits use the same SRHT plan size and global AR specification. The only difference is that the landmark fit solves at 48 of 192 voxels and interpolates the rest. The simulated task maps are spatially smooth because interpolating independent voxel effects would be an ill-posed target.

# Choose landmarks (L) and neighbors (k)
L <- 48L; k_nn <- 8L
low_lm <- lowrank_control(
  parcels = parcels,
  landmarks = L,
  k_neighbors = k_nn,
  time_sketch = list(method = "srht", m = min(8L*p, Tlen))
)

fit_lm_srht <- fmri_lm(onset ~ hrf(condition), block = ~ run, dataset = dset,
                       engine = "latent_sketch",
                       lowrank = low_lm,
                       control = fmri_lm_control(
                         noise = noise_spec("ar1", pooling = "global")
                       ))

# Compare task coefficients to the like-for-like full-voxel sketch.
task_reference <- as.numeric(
  fit_srht_global$betas_fixed[task_cols, , drop = FALSE]
)
task_landmarks <- as.numeric(
  fit_lm_srht$betas_fixed[task_cols, , drop = FALSE]
)
cor_landmarks <- cor(task_reference, task_landmarks)
rmse_landmarks <- sqrt(mean((task_reference - task_landmarks)^2))
nrmse_landmarks <- rmse_landmarks / stats::sd(task_reference)

stopifnot(cor_landmarks > 0.85, nrmse_landmarks < 0.60)

c(
  correlation = cor_landmarks,
  RMSE = rmse_landmarks,
  NRMSE = nrmse_landmarks
)
#> correlation        RMSE       NRMSE 
#>   0.8899257   2.3383580   0.4839667

The correlation checks spatial ordering, while NRMSE detects scale errors that correlation alone would miss. The thresholds above are an explicit contract for this smooth synthetic field, not a general guarantee for real images. Validate both metrics on representative pilot data before adopting landmarks in an analysis.

# Compare one task-effect slice on a common scale, then show signed error.
task_index <- task_cols[1]
slice_index <- ceiling(dim3[3] / 2)
reference_volume <- array(
  fit_srht_global$betas_fixed[task_index, ], dim = dim3
)
landmark_volume <- array(
  fit_lm_srht$betas_fixed[task_index, ], dim = dim3
)
reference_slice <- reference_volume[, , slice_index]
landmark_slice <- landmark_volume[, , slice_index]
error_slice <- landmark_slice - reference_slice

effect_limit <- max(abs(c(reference_slice, landmark_slice)))
error_limit <- max(abs(error_slice))
effect_palette <- hcl.colors(64, "Blue-Red 3")
error_palette <- hcl.colors(64, "Blue-Red 3")

old_par <- par(mfrow = c(1, 3), mar = c(3, 3, 3, 1))
image(
  reference_slice, zlim = c(-effect_limit, effect_limit),
  col = effect_palette, asp = 1, axes = FALSE,
  main = "Full-voxel reference"
)
image(
  landmark_slice, zlim = c(-effect_limit, effect_limit),
  col = effect_palette, asp = 1, axes = FALSE,
  main = "Landmark extension"
)
image(
  error_slice, zlim = c(-error_limit, error_limit),
  col = error_palette, asp = 1, axes = FALSE,
  main = "Error: landmark - reference"
)

par(old_par)

The common color scale in the first two panels makes attenuation visible; the third panel localizes signed interpolation error. Smooth regions are recovered most closely, while gradients and boundaries show the largest departures. This is why a whole-map correlation alone is insufficient.

Notes:

• Landmarks are selected by k-means over voxel coordinates and mapped to nearest mask voxels.
• The extension uses coordinate k-NN heat-kernel weights with row-normalization. It can smooth across anatomical or parcel boundaries.
• Residual variance is propagated from landmark residuals through squared interpolation weights; this is an approximation, not a calibrated voxelwise-inference guarantee.

What the comparisons establish

• The iid SRHT comparison isolates temporal-sketch error from AR-model changes.
• The landmark comparison isolates spatial interpolation while holding the temporal sketch and AR specification fixed.
• Results on smooth synthetic fields do not establish accuracy at sharp boundaries or for spatially irregular effects.
• No runtime claim is made here. Benchmark the exact workload, including data access and preprocessing, before choosing an approximation for speed.

When to use which knob

• Reduce T with a temporal sketch only after validating task-coefficient error against the unsketched estimator.
• Reduce V with landmarks only when coefficient fields are spatially smooth enough to satisfy a declared map-error tolerance.
• Omit landmarks when voxelwise detail or sharp boundaries are central to the scientific question.

Pros and cons

Pros

• Temporal and spatial approximation controls are explicit in lowrank_control().
• Global and parcel-pooled AR models are expressed through the same public noise_spec() interface.
• Correlation and scale-sensitive error can be validated separately.

Cons

• Landmark extension introduces spatial shrinkage that is data-dependent and can cross boundaries; if you need strict voxelwise inference, omit landmarks.
• Sketch-and-solve (SRHT, Gaussian, CountSketch) trades precision for data touched: standard errors grow by about sqrt(T / m) and residual degrees of freedom shrink to about m - p. Choose m for the precision you need, or fit exactly (omit engine = "latent_sketch") when full-data precision matters.

Practical defaults

• Time sketch: list(method = "countsketch", m = min(8*p, T)) or list(method = "gaussian", m = min(8*p, T)) for speed; SRHT gives the same inference at somewhat higher cost. "ihs" is deprecated and fits exact OLS.
• AR order: 1 is a solid default; use 2 if residuals show lag-2 structure.
• Parcel AR: use noise_spec("ar1", pooling = "parcel", parcels = parcels); tune shrinkage for the target data rather than assuming one universal value.
• Landmarks: choose L and k_neighbors by a declared correlation and NRMSE contract on representative data.

Reproducibility tips

• Set a seed before building SRHT plans or k-means landmarks (set.seed(...)).
• Keep the sketch method and m fixed across runs to compare apples-to-apples.
• Validate on a small ROI first: check correlations vs exact, and histogram of t-maps in null tissue.

Next

• vignette("a_09_linear_model", package = "fmrireg") — fMRI Linear Model (GLM)
• vignette("group_analysis", package = "fmrireg") — Group Analysis