Performance Optimizations for sparse long and wide random matrices

[aschneuw]

Performance Evaluation

Compute Node

## System Information

### OS
Description:    Ubuntu 22.04.5 LTS

### Kernel
6.8.0-45-generic

### CPU
Architecture:                         x86_64
CPU(s):                               32
On-line CPU(s) list:                  0-31
Model name:                           AMD Ryzen Threadripper PRO 5955WX 16-Cores
NUMA node0 CPU(s):                    0-31

### Memory
               total
Mem:           503Gi

### R
R version 4.5.1 (2025-06-13) -- "Great Square Root"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu

Comparison

The aim of this performance evaluation is to compare the two following gamm4 versions:

• current (2.6 / 2.7): current vanilla gamm4 - Optimizer is always nloptwrap since the optimizer argument is not propagated properly.
• mod (2.8 ?): modified version of gamm4 with bug fixed, scikit-sparse support for data covariance matrix Cholesky decomposition and performance-optimized version of Fabian Scheipl's trick

Heat Stress Dataset

More info about the data: https://aschneuw.github.io/heatstress/materials/

model <- gamm4(
  formula= milk ~ 1 + s(thi_mean_t0_3d, by=parity, k=10) + s(days_in_milk_t, by=parity, k=15) + parity + year,
  random =~((1 | zip_month) + (1 | farmIdLocationSample) + (1 | animalId)),
  data = data,
  drop.unused.levels = TRUE,
  REML = TRUE,
  control = lmerControl(calc.derivs = FALSE, optCtrl = list(optimizer="nloptwrap", maxfun = 10000)),
  verbose = 10,
  python_cholmod = TRUE
)

Version	Samples	Animals	Farms	ZIPxMonth	Optimizer	scikit-sparse	Execution Time [s]
2.7	500'000	23'675	4'241	16'056	nloptwrap	-	21'536
2.7	734’685	23'675	4'302	16'648	nloptwrap	-	crash
mod	500'000	23'675	4'241	16'056	nloptwrap	False	5'240
mod	500'000	23'675	4'241	16'056	nloptwrap	True	4'316
mod	734’685	23'675	4'302	16'648	nloptwrap	False	crash
mod	734’685	23'675	4'302	16'648	nloptwrap	True	6'089
mod	734’685	23'675	4'302	16'648	BOBYQA	True	6'067

Random Data Example from Documentation

Benchmarking a mixed model example from the documentation (with more samples and factor levels) with microbench.

Current

gamm4_example <- function() {
    set.seed(0) 
    
    n = 20000
    m = 2000
    dat <- gamSim(1, n = n, scale = 2)  # simulate 4-term additive truth
    
    # Add 3,000-level random effect 'fac'
    dat$fac <- fac <- as.factor(sample(1:m, n, replace = TRUE))
    dat$y <- dat$y + model.matrix(~fac - 1) %*% rnorm(m) * 0.5
    
    # Fit gamm4 model
    br <- gamm4(y ~ s(x0) + x1 + s(x2), data = dat, random = ~(1 | fac))
    plot(br$gam, pages = 1)
    
    # Model summaries
    summary(br$gam)  # summary of gam
    summary(br$mer)  # underlying mixed model
}

Unit: seconds
            expr      min       lq     mean   median       uq     max neval
 gamm4_example() 17.02194 17.29564 18.51402 19.16437 19.44374 19.6596   100

Mod

gamm4_example <- function() {
    set.seed(0) 
    
    n = 20000
    m = 2000
    dat <- gamSim(1, n = n, scale = 2)  # simulate 4-term additive truth
    
    # Add 3,000-level random effect 'fac'
    dat$fac <- fac <- as.factor(sample(1:m, n, replace = TRUE))
    dat$y <- dat$y + model.matrix(~fac - 1) %*% rnorm(m) * 0.5
    
    # Fit gamm4 model
    br <- gamm4(y ~ s(x0) + x1 + s(x2), data = dat, random = ~(1 | fac), python_cholmod=TRUE)
    plot(br$gam, pages = 1)
    
    # Model summaries
    summary(br$gam)  # summary of gam
    summary(br$mer)  # underlying mixed model
}

Unit: seconds
            expr      min      lq     mean   median       uq      max neval
 gamm4_example() 8.155614 8.32772 9.079993 9.511128 9.655152 9.924096   100