| Component | Version |
|---|---|
| Whittaker | 0.1.dev390+g931406871.d20260812 |
| Python | 3.13.12 |
| NumPy | 2.5.2 |
| SciPy | 1.18.0 |
| OS | Darwin arm64 |
| CPU | arm |
| CPU count | 12 |
Benchmarks
This page reports the performance of Whittaker’s core operations on a single machine. All timings are wall-clock medians over three runs (one run for the largest sizes and MCMC). The benchmarks establish a baseline for the current pure-Python/NumPy/SciPy implementation before any future Rust or JAX acceleration work.
Environment
Sample size scaling
The most fundamental performance question is how fit time grows with the number of observations. This benchmark fits a Gaussian GAM with a single cubic regression spline (y ~ s(x, bs='cr')) using REML across a range of sample sizes. Cubic regression splines use a banded basis matrix, so timing reflects the general P-IRLS cost without the O(n^2) distance-matrix overhead of TPRS.
| n | Fit time (s) |
|---|---|
| 500 | 0.0025 |
| 1,000 | 0.0025 |
| 2,000 | 0.0029 |
| 5,000 | 0.0037 |
| 10,000 | 0.0056 |
Fit time grows roughly as O(n \cdot k^2) where k is the basis dimension (default 10). The dominant cost is the Cholesky factorization of the penalized normal equations at each P-IRLS step.
Smooth type comparison
Different smooth bases have different construction and penalty costs. This benchmark fixes n = 3{,}000 and compares the main smooth types on a Gaussian response with REML.
| Smooth | Fit time (s) |
|---|---|
| TPRS | 3.5078 |
| Cubic RS | 0.0033 |
| P-spline | 0.0029 |
| Cyclic CRS | 0.0031 |
| Tensor (te) | 6.9691 |
Cubic regression splines, P-splines, and cyclic splines are the fastest because they use banded basis matrices. TPRS requires an O(n^2) eigendecomposition of the distance matrix, making it substantially slower at this sample size. Tensor products are the most expensive because they construct a Kronecker product penalty over two marginal bases.
Family comparison
Non-Gaussian families require iterative fitting (P-IRLS) with multiple Newton steps, so they are inherently more expensive than the Gaussian case. This benchmark fixes n = 3{,}000 with a cubic regression spline and compares across response distributions.
| Family | Fit time (s) |
|---|---|
| Gaussian | 0.0035 |
| Poisson | 0.0070 |
| Binomial | 0.0073 |
| Gamma | 0.0080 |
| Neg. Binomial | 0.0087 |
Gaussian fitting converges in a single P-IRLS step (it is a direct solve). Non-Gaussian families typically require 3–8 iterations, with the Negative Binomial being the most expensive due to its additional scale parameter estimation.
Fitting method comparison
Whittaker supports several smoothness selection and inference methods. This benchmark compares them on a Poisson GAM with a cubic regression spline at n = 2{,}000 to show the relative cost of each approach.
| Method | Fit time (s) |
|---|---|
| REML | 0.0064 |
| GCV | 0.0128 |
| ML | 0.0070 |
| VI | 0.0232 |
| MCMC | 2.2017 |
REML, GCV, and ML have similar costs because they all use P-IRLS with outer smoothness selection. VI adds optimization of the variational parameters on top of the P-IRLS initialization. MCMC is the most expensive because it runs hundreds of leapfrog steps per sample, but it provides exact posterior inference.
Prediction throughput
Once a model is fitted, prediction speed matters for serving. This benchmark fits a single Gaussian GAM on n = 5{,}000 observations and then times predict() across a range of prediction-set sizes.
| n_new | Predict time (ms) |
|---|---|
| 100 | 0.39 |
| 500 | 0.47 |
| 1,000 | 0.60 |
| 5,000 | 1.18 |
| 10,000 | 2.03 |
Prediction is dominated by design matrix construction (evaluating the basis functions at the new covariate values) and a single matrix–vector multiply. It scales linearly with the prediction set size.
BigGAM vs. GAM scaling
BigGAM uses discretised covariates and compressed sufficient statistics to avoid building the full n \times p design matrix. This makes it slower than GAM at small sample sizes (due to discretization overhead) but faster at large ones.
| n | GAM (s) | BigGAM (s) | Faster |
|---|---|---|---|
| 1,000 | 0.0027 | 0.0028 | GAM |
| 2,000 | 0.0028 | 0.0027 | BigGAM |
| 5,000 | 0.0039 | 0.0032 | BigGAM |
| 10,000 | 0.0058 | 0.0039 | BigGAM |
| 25,000 | 0.0118 | 0.0059 | BigGAM |
BigGAM’s discretized design matrix trades a small fixed overhead for dramatically better scaling. At small sample sizes the overhead may make BigGAM slightly slower than GAM, but it quickly becomes the faster option as n grows and dominates for large-scale fitting.