GAM for large datasets using discretized fitting.
BigGAM(
formula,
family=None,
*,
n_discrete=200,
)
BigGAM is a drop-in subclass of ~whittaker.gam.GAM for datasets too large to fit comfortably with the standard dense design matrix (roughly n > 1_000_000). It uses the bam approach of Wood, Li, & Shaddick (2017): each covariate is rounded onto a grid of at most n_discrete representative values, and the smooth basis is evaluated only once per unique (combination of) discretized value(s) rather than once per observation. An index array records, for every observation, which unique bin it fell into (see DiscretizedBlock.indices in build_discretized_model_matrix).
Fitting still runs penalized iteratively reweighted least squares (P-IRLS), exactly as in ~whittaker.gam.GAM.fit, alternating an inner coefficient update with an outer smoothing parameter selection. The difference is purely computational: instead of forming the full n x p design matrix X and computing X'WX and X'Wz directly, bam_fit (in whittaker.fitting.bam) accumulates these quantities per bin — for each smooth’s discretized block, observation weights are aggregated into per-bin totals via numpy.bincount over the bin indices, and the resulting d x d cross-products of unique basis rows (d = number of unique bins) are scattered into the correct p x p block of X'WX (see _compute_XtWX). Because d can be orders of magnitude smaller than n, this reduces the memory needed for the cross-product step from O(n p) to O(d p), and the resulting fit closely approximates the exact (non-discretized) GAM fit on the same data.
BigGAM is a drop-in subclass of GAM: predict(), summary(), plot(), and check() all work the same way as for GAM. The one internal difference is that self._model_matrix.X is an empty array (the dense design matrix is never materialized) so operations that would otherwise reconstruct per-term columns (e.g. smooth_tests()) instead re-expand each smooth’s columns on demand from its DiscretizedBlock via _expand_block_columns.
Parameters
formula: str | Formula-
Model formula as a string (e.g. "y ~ s(x1) + s(x2) + x3"), or an already-parsed Formula object. Same syntax as ~whittaker.gam.GAM.
family: Family | None = None-
Response distribution family, e.g. Gaussian(), Binomial(), Poisson(), Gamma(), or Tweedie(). Defaults to Gaussian().
n_discrete: int = 200-
Maximum number of unique representative values per covariate (or per combination of covariates, for multi-dimensional smooths). Defaults to
200.
Notes
n_discrete controls the accuracy/memory tradeoff directly. Larger values give a discretized grid that more finely resolves each covariate’s range, so the fit approaches the exact (non-discretized) GAM fit at the cost of more unique bins d and therefore more memory and computation in the X'WX accumulation step. Smaller values reduce memory and speed up fitting, but coarsen the covariate resolution: because all observations within a bin share the same basis row, this can slightly bias smooths with high curvature or steep local features, since fine-scale variation within a bin is averaged away. The default of 200 is usually more than enough resolution for typical smooth terms; it rarely needs to be increased unless a covariate has an unusually large number of important local features. The benefit of discretization (versus plain GAM) is only realized once n is much larger than n_discrete, i.e. for large datasets — for small or moderate n, GAM is simpler and just as fast.
Examples
import numpy as np
import whittaker as wt
from whittaker.bam import BigGAM
rng = np.random.default_rng(0)
n = 5_000
x1 = rng.uniform(0, 1, n)
x2 = rng.uniform(0, 1, n)
y = np.sin(2 * np.pi * x1) + x2**2 + rng.normal(scale=0.2, size=n)
model = BigGAM("y ~ s(x1) + s(x2)", n_discrete=100).fit({"x1": x1, "x2": x2, "y": y})
print(model.summary())
GAM fit summary
============================================================
Formula: y ~ s(x1) + s(x2)
Family: Gaussian(link='identity')
Inference: REML
Observations: 0
Coefficients: 19
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.3376 0.0000 0.000 1
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x1) 8.91 9 62254.971 < 1e-16
s(x2) 6.26 7 11059.123 < 1e-16
Total EDF: 16.17
Scale est: 0.040068
Deviance: 199.6941
Null dev: 3148.4427
Dev. expl: 93.7%
GCV score: 0.040198
AIC: -1880.28
BIC: -1774.91
This example uses a modest n for speed, but BigGAM’s memory and speed advantage over plain GAM really shows up once n reaches into the millions, where materializing the full design matrix would be impractical.
Attributes
|
Name
|
Description
|
|
n_discrete
|
Number of discretization grid points per covariate.
|
n_discrete
Number of discretization grid points per covariate.
This is the n_discrete value passed to __init__: the maximum number of unique representative values each covariate (or combination of covariates, for multi-dimensional smooths) is rounded onto before the smooth basis is evaluated. It bounds the number of unique rows d used when accumulating X'WX during fitting (see the class docstring), and therefore controls the tradeoff between fit accuracy and memory/speed.
Methods
|
Name
|
Description
|
|
fit()
|
Fit the BigGAM using discretized P-IRLS.
|
|
smooth_tests()
|
Approximate significance tests for each smooth term.
|
fit()
Fit the BigGAM using discretized P-IRLS.
fit(
data,
*,
smoothing_params=None,
method="fREML",
weights=None,
select=False,
vi_options=None,
mcmc_options=None,
)
Builds a DiscretizedModelMatrix via build_discretized_model_matrix and fits it with bam_fit, which accumulates X'WX and X'Wz from per-bin basis blocks instead of the full design matrix (see the class docstring for how the discretization and cross-product accumulation work).
Parameters
data: dict[str, numpy.ndarray] or InputData-
Column-oriented data as {name: 1-D array} (or any InputData-compatible object). All columns referenced by the formula must be present and of equal length.
smoothing_params: list of float = None-
Fixed smoothing parameters lambda_j, one per smooth term, in formula order. If None (the default), smoothing parameters are selected automatically according to method.
method: str = "fREML"-
Criterion used to select smoothing parameters when smoothing_params is None. One of "fREML" (default, fast discretized REML), "REML", "ML", or "GCV". See ~whittaker.gam.GAM.fit for a description of each criterion.
weights: numpy.ndarray = None-
Observation (prior) weights, shape (n,). Must be strictly positive.
select: bool = False-
If
True, add an extra penalty on each smooth’s null space so the term can be shrunk to exactly zero (double-penalty selection). Defaults to False.
Returns
BigGAM-
Returns
self for method chaining, e.g. model = BigGAM("y ~ s(x)").fit(data).
smooth_tests()
Approximate significance tests for each smooth term.
This is the BigGAM counterpart to ~whittaker.gam.GAM.smooth_tests, adapted to the discretized fitting path: because self._model_matrix.X is never materialized, the per-term design columns needed for the test are reconstructed on demand from each smooth’s DiscretizedBlock via _expand_block_columns (or read directly from parametric_cols for non-smooth columns) rather than sliced out of a dense X. For each smooth term, the coefficient sub-vector, its covariance sub-block, and its expanded design columns are passed to ~whittaker.fitting.inference._smooth_test to obtain an approximate chi-squared test of whether the term is uniformly zero.
Returns
list of SmoothTestResult-
One result per smooth term (or per
by= level), each with term_label, stat, edf, ref_df, and p_value attributes.