# 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


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


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


<style>
  #altair-viz-10d29884ff5d4466add02083c2911697.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-10d29884ff5d4466add02083c2911697.vega-embed details,
  #altair-viz-10d29884ff5d4466add02083c2911697.vega-embed details summary {
    position: relative;
  }
</style>


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


<style>
  #altair-viz-bdfce2d19da64a2484ae5081a807eac6.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-bdfce2d19da64a2484ae5081a807eac6.vega-embed details,
  #altair-viz-bdfce2d19da64a2484ae5081a807eac6.vega-embed details summary {
    position: relative;
  }
</style>


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


<style>
  #altair-viz-703af072c1d14966ba77a71fd8275c60.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-703af072c1d14966ba77a71fd8275c60.vega-embed details,
  #altair-viz-703af072c1d14966ba77a71fd8275c60.vega-embed details summary {
    position: relative;
  }
</style>


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


<style>
  #altair-viz-bbd758fc0dd14f0cb8016223f363b2ae.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-bbd758fc0dd14f0cb8016223f363b2ae.vega-embed details,
  #altair-viz-bbd758fc0dd14f0cb8016223f363b2ae.vega-embed details summary {
    position: relative;
  }
</style>


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


<style>
  #altair-viz-2ed124a7886248f3be609feb15276ed1.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-2ed124a7886248f3be609feb15276ed1.vega-embed details,
  #altair-viz-2ed124a7886248f3be609feb15276ed1.vega-embed details summary {
    position: relative;
  }
</style>


|  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](../reference/BigGAM.md#whittaker.BigGAM) uses discretised covariates and compressed sufficient statistics to avoid building the full n \times p design matrix. This makes it slower than [GAM](../reference/GAM.md#whittaker.GAM) at small sample sizes (due to discretization overhead) but faster at large ones.


<style>
  #altair-viz-b02bdf6717db4e238d51e72bc7e65d0e.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-b02bdf6717db4e238d51e72bc7e65d0e.vega-embed details,
  #altair-viz-b02bdf6717db4e238d51e72bc7e65d0e.vega-embed details summary {
    position: relative;
  }
</style>


|      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](../reference/BigGAM.md#whittaker.BigGAM)'s discretized design matrix trades a small fixed overhead for dramatically better scaling. At small sample sizes the overhead may make [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) slightly slower than [GAM](../reference/GAM.md#whittaker.GAM), but it quickly becomes the faster option as n grows and dominates for large-scale fitting.