# Large datasets Standard GAM fitting builds the full n \times p design matrix, which becomes memory-prohibitive when n is in the millions. Whittaker provides three scalable backends that use discretised covariates and streaming aggregation to fit GAMs on large datasets without materializing the full matrix: | Class | Data source | Key feature | |----|----|----| | [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) | In-memory dict | Discretised P-IRLS on NumPy arrays | | [PolarsGAM](../reference/PolarsGAM.md#whittaker.PolarsGAM) | Polars DataFrames, LazyFrames, or files | Lazy streaming via Polars | | [DuckDBGAM](../reference/DuckDBGAM.md#whittaker.DuckDBGAM) | DuckDB tables or SQL queries | SQL-native streaming via Arrow | All three extend [GAM](../reference/GAM.md#whittaker.GAM), so prediction, summary, diagnostics, and serialization work the same way as a regular GAM once fitted. # BigGAM: discretised fitting [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) implements the discretised P-IRLS algorithm of Wood, Li & Shaddick (2017). Instead of storing one row per observation, each covariate is discretised into a grid of [n_discrete](../reference/BigGAM.md#whittaker.BigGAM.n_discrete) unique values (default 200), and the sufficient statistics are accumulated over the grid. This reduces memory from O(n \cdot p) to O(d \cdot p) where d \ll n. ``` python import numpy as np import whittaker as wk # Simulate a dataset large enough to demonstrate discretised fitting rng = np.random.default_rng(23) n = 2_000 x1 = rng.uniform(0, 2 * np.pi, n) x2 = rng.uniform(0, 1, n) y = np.sin(x1) + 2 * x2 + rng.normal(0, 0.3, n) data = {"x1": x1, "x2": x2, "y": y} # BigGAM uses the same formula syntax model = wk.BigGAM("y ~ s(x1) + s(x2)", n_discrete=200) model.fit(data, method="fREML") 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.9970 0.0000 0.000 1 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x1) 8.60 9 11090.453 < 1e-16 s(x2) 1.00 2 7182.779 < 1e-16 Total EDF: 10.60 Scale est: 0.093291 Deviance: 185.5926 Null dev: 1886.2454 Dev. expl: 90.2% GCV score: 0.093787 AIC: 942.28 BIC: 1001.63 The `method="fREML"` (fast REML) is the default and recommended smoothing parameter selection method for [BigGAM](../reference/BigGAM.md#whittaker.BigGAM). It exploits the discretised structure for efficient computation of the REML criterion. ## Controlling the discretisation grid The [n_discrete](../reference/BigGAM.md#whittaker.BigGAM.n_discrete) parameter controls the resolution of the covariate grid. Larger values give a more faithful approximation but use more memory and time. ``` python # Compare discretisation resolutions for nd in [50, 200]: m = wk.BigGAM("y ~ s(x1) + s(x2)", n_discrete=nd) m.fit(data, method="fREML") print(f"n_discrete={nd:3d}: EDF = {m.edf_total:.1f}, scale = {m.scale:.4f}") ``` n_discrete= 50: EDF = 10.6, scale = 0.0939 n_discrete=200: EDF = 10.6, scale = 0.0933 > **Tip: Choosing n_discrete** > > For most datasets, the default of 200 is a good balance. You rarely need more than 500. If the true function is very smooth, even 50 can work well. Increase [n_discrete](../reference/BigGAM.md#whittaker.BigGAM.n_discrete) if the model diagnostics suggest the discretisation is too coarse, but check `model.check()` first since basis dimension inadequacy is more common than discretisation error. ## Predictions Predictions work exactly like a standard GAM: ``` python x_new = np.linspace(0, 2 * np.pi, 100) pred = model.predict({"x1": x_new, "x2": np.full(100, 0.5)}) print(f"Prediction shape: {pred.values.shape}") ``` Prediction shape: (100,) ``` python import altair as alt # Compare BigGAM fit with the true function true_vals = np.sin(x_new) + 2 * 0.5 plot_data = [ {"x1": float(x_new[i]), "y": float(pred.values[i]), "source": "BigGAM"} for i in range(len(x_new)) ] + [ {"x1": float(x_new[i]), "y": float(true_vals[i]), "source": "Truth"} for i in range(len(x_new)) ] alt.Chart({"values": plot_data}).mark_line(strokeWidth=2).encode( x=alt.X("x1:Q", title="x1"), y=alt.Y("y:Q", title="f(x1) at x2=0.5"), color=alt.Color("source:N"), strokeDash=alt.condition( alt.datum.source == "Truth", alt.value([4, 4]), alt.value([0]) ), ).properties(width="container", height=300, title="BigGAM fit vs. truth (2,000 observations)") ``` # PolarsGAM: fitting from Polars and files [PolarsGAM](../reference/PolarsGAM.md#whittaker.PolarsGAM) extends [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) to accept Polars DataFrames, LazyFrames, or file paths. Data is streamed in chunks, so the full dataset never needs to be in memory at once. ## From a Polars DataFrame ``` python import polars as pl # Create a Polars DataFrame df = pl.DataFrame({"x1": x1, "x2": x2, "y": y}) model_pl = wk.PolarsGAM("y ~ s(x1) + s(x2)", chunk_size=2_000) model_pl.fit(df, method="fREML") print(f"Rows processed: {model_pl.n_rows}") print(f"Chunk size: {model_pl.chunk_size}") print(f"EDF: {model_pl.edf_total:.1f}") ``` Rows processed: 2000 Chunk size: 2000 EDF: 10.6 ## From a LazyFrame or file [PolarsGAM](../reference/PolarsGAM.md#whittaker.PolarsGAM) can scan Parquet, CSV, IPC (Arrow), and NDJSON files directly. The file is read lazily via Polars' streaming engine, so only one chunk is in memory at a time: ``` python # Save to Parquet for demonstration import tempfile, pathlib tmpdir = pathlib.Path(tempfile.mkdtemp()) parquet_path = tmpdir / "data.parquet" df.write_parquet(parquet_path) # Fit directly from the file model_file = wk.PolarsGAM("y ~ s(x1) + s(x2)", chunk_size=2_000) model_file.fit(str(parquet_path), method="fREML") print(f"EDF from Parquet: {model_file.edf_total:.1f}") ``` EDF from Parquet: 10.6 > **Note: Supported file formats** > > | Extension | Reader | > |------------------|---------------------| > | `.parquet` | `pl.scan_parquet()` | > | `.csv` | `pl.scan_csv()` | > | `.ipc`, `.arrow` | `pl.scan_ipc()` | > | `.ndjson` | `pl.scan_ndjson()` | > > The format is detected from the file extension. # DuckDBGAM: fitting from SQL [DuckDBGAM](../reference/DuckDBGAM.md#whittaker.DuckDBGAM) streams data from DuckDB tables or SQL queries via DuckDB's Arrow interface. This is ideal when the data lives in a database or when you want to filter and transform data with SQL before fitting. ``` python import duckdb # Create a DuckDB connection and load data conn = duckdb.connect() conn.execute("CREATE TABLE obs AS SELECT * FROM read_parquet(?)", [str(parquet_path)]) # Fit from a table name model_duck = wk.DuckDBGAM("y ~ s(x1) + s(x2)", chunk_size=2_000) model_duck.fit("obs", conn, method="fREML") print(f"EDF from DuckDB: {model_duck.edf_total:.1f}") del model_pl, model_file, df ``` EDF from DuckDB: 10.6 ## SQL queries Use [fit_query()](../reference/DuckDBGAM.md#whittaker.DuckDBGAM.fit_query) to fit from an arbitrary SQL query: ``` python # Fit on a filtered subset via SQL model_filtered = wk.DuckDBGAM("y ~ s(x1)", chunk_size=2_000) model_filtered.fit_query( "SELECT x1, y FROM obs WHERE x2 > 0.5", conn, method="fREML", ) print(f"EDF (filtered): {model_filtered.edf_total:.1f}") ``` EDF (filtered): 9.0 > **Tip: When to use DuckDBGAM** > > Use [DuckDBGAM](../reference/DuckDBGAM.md#whittaker.DuckDBGAM) when: > > - Your data is already in DuckDB (a common analytics stack) > - You want to pre-filter or transform data with SQL before fitting > - The dataset is too large for a Polars scan (DuckDB's out-of-core engine handles very large files) # Smoothing parameter selection All three scalable backends support the same selection methods: | Method | Description | |-----------|-------------------------------------------------| | `"fREML"` | Fast REML (default, recommended for large data) | | `"REML"` | Restricted maximum likelihood | | `"ML"` | Maximum likelihood | | `"GCV"` | Generalised cross-validation | `fREML` is specifically optimised for discretised fitting and is significantly faster than standard REML on large datasets. # Smooth selection The `select=True` option enables double-penalty smooth selection, which can shrink entire terms to zero: ``` python # Add a noise variable noise = rng.uniform(0, 1, n) data_noise = {"x1": x1, "x2": x2, "noise": noise, "y": y} model_select = wk.BigGAM("y ~ s(x1) + s(x2) + s(noise)", n_discrete=200) model_select.fit(data_noise, method="fREML", select=True) print(model_select.summary()) ``` GAM fit summary ============================================================ Formula: y ~ s(x1) + s(x2) + s(noise) Family: Gaussian(link='identity') Inference: REML Observations: 0 Coefficients: 28 Parametric coefficients: Term Estimate Std.Err t value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) 0.9970 0.0000 0.000 1 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x1) 8.59 9 11095.808 < 1e-16 s(x2) 1.00 2 7188.896 < 1e-16 s(noise) 0.65 1 1.781 0.182 Total EDF: 11.24 Scale est: 0.093209 Deviance: 185.3698 Null dev: 1886.2454 Dev. expl: 90.2% GCV score: 0.093735 AIC: 941.16 BIC: 1004.10 The noise term should receive a very low EDF, effectively being selected out. # Non-Gaussian families Scalable backends work with all response families: ``` python # Poisson counts rng = np.random.default_rng(23) n_pois = 2_000 x = rng.uniform(0, 2 * np.pi, n_pois) mu = np.exp(0.5 + 0.8 * np.sin(x)) y_pois = rng.poisson(mu).astype(float) model_pois = wk.BigGAM("y ~ s(x)", family=wk.Poisson(), n_discrete=200) model_pois.fit({"x": x, "y": y_pois}, method="fREML") print(f"Poisson BigGAM EDF: {model_pois.edf_total:.1f}") ``` Poisson BigGAM EDF: 8.2 # Which backend to choose | Scenario | Backend | Why | |----|----|----| | Data fits in memory | [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) | Simplest API, dict input | | Data in Polars or files | [PolarsGAM](../reference/PolarsGAM.md#whittaker.PolarsGAM) | Lazy streaming, zero-copy | | Data in DuckDB / SQL | [DuckDBGAM](../reference/DuckDBGAM.md#whittaker.DuckDBGAM) | SQL filtering, out-of-core | | Incremental batches | [StreamingGAM](../reference/StreamingGAM.md#whittaker.StreamingGAM) | Online updates, decay | For datasets under ~100K rows, the standard [GAM](../reference/GAM.md#whittaker.GAM) is usually fast enough. Above that, [BigGAM](../reference/BigGAM.md#whittaker.BigGAM) or its streaming variants keep fitting times manageable without sacrificing accuracy. ``` python # Clean up conn.close() import shutil shutil.rmtree(tmpdir) ``` You can now choose the right scalable backend for your data source and fit GAMs on datasets that would not fit in memory with the standard [GAM](../reference/GAM.md#whittaker.GAM) class. # Where to go next - **[Streaming and online fitting](streaming.md)**: incremental fitting when data arrives in batches over time. - **[Model fitting](fitting.md)**: how P-IRLS and smoothness selection work in the standard pipeline. - **[Saving and loading models](serialization.md)**: persist a fitted large-data model for deployment.