Generalized Additive Models (GAM)
Sums of smooth functions, fit by backfitting — coordinate descent for curves.
On this page
Beginner: a generalized additive model (GAM) predicts y as a plain SUM of separate curves, one per feature — y = β₀ + f₁(x₁) + f₂(x₂) + ⋯ + fₚ(xₚ) — instead of a single straight-line combination of features (plain OLS) or one giant multivariate curved surface fit over all features at once. Each fⱼ is free to bend and wiggle however the data for that one feature demands; the model as a whole is just those bent curves stacked on top of each other.
Intermediate: why not just fit one big smooth surface over all features jointly, using the same spline machinery from the previous lesson but in higher dimensions? Two practical reasons. First, a fully general smooth surface over p features needs enough data to cover a p-dimensional neighborhood around every point you'd like to predict at, and the amount of data needed for that explodes with p — the same curse of dimensionality that shows up everywhere in ML. Past about 2–3 features, this becomes both statistically intractable to fit well and impossible to even draw. Second, and just as important: even if you somehow fit that surface perfectly, you couldn't read it. A GAM sidesteps both problems by restricting the surface to a very particular, much smaller family — sums of independent one-dimensional curves — each of which is exactly as easy to fit and to plot as the single-feature smoothing spline from the previous lesson.
Advanced: the additive restriction is precisely what buys back interpretability. Because the model is a sum, holding every other feature fixed and varying just xⱼ changes the prediction by exactly fⱼ(xⱼ) — no other term moves. That means you can plot fⱼ on its own, on a single 2D axis, and read off precisely how that one feature affects the prediction, independent of whatever the other features happen to be doing. A single p-dimensional smooth surface offers no such decomposition — you cannot, in general, ask "what does feature 3 alone contribute?" of an arbitrary multivariate function, because its effect on the output can depend on every other feature's value too. This is also exactly the price a GAM pays, revisited in the Pitfalls below: it can only represent effects that decompose this cleanly, by construction.
Each fⱼ is typically represented the same way the previous lesson represented a single smooth function — a spline or basis expansion — and each carries its own smoothing penalty λⱼ controlling how wiggly that particular feature's curve is allowed to be. Setting every fⱼ to a straight line recovers plain linear regression exactly; a GAM is a strict generalization of OLS, not a different model family bolted on beside it.
Fitting a GAM means minimizing a penalized sum of squared errors over all p smooth functions jointly:
That is a single optimization problem over p unknown functions at once — in general no closed form exists for all of them simultaneously. But now fix every function fₖ for k ≠ j at whatever their current values are, and minimize over fⱼ alone. Every penalty term except λⱼ∫fⱼ''² is now a constant (it doesn't involve fⱼ), so it drops out of the minimization, and the data-fit term can be regrouped by pulling every already-fixed fₖ(xᵢₖ) to the other side:
Read that carefully: it is exactly the single-feature smoothing-spline problem from the previous lesson — minimize squared error plus a roughness penalty over one function of one variable — with the ordinary target y replaced by the partial residual r⁽ʲ⁾ = y - β₀ - Σₖ≠ⱼ fₖ(xₖ): whatever is left of y after every other feature's current contribution has been subtracted out. Nothing new needs to be derived to solve this — it's the previous lesson's machinery, called on a different target vector.
That single fact is the whole algorithm. Backfitting cycles through the features one at a time: compute the partial residual for feature j (removing every other feature's current fit from y), refit fⱼ as a 1D smoothing spline against that partial residual, move to feature j+1, and repeat the whole cycle until none of the fⱼ change appreciably between passes:
This is precisely coordinate descent (Module 2, 2.2.5) — optimize one coordinate with everything else held fixed, cycle through all coordinates, repeat — except each "coordinate" here is not a single scalar but an entire one-dimensional function. Every other structural fact from that lesson carries over unchanged: each individual update is an exact minimizer of the joint objective over its one coordinate (function), so the joint penalized objective can only decrease or stay exactly the same at every single half-step, for the same reason a scalar coordinate-descent step can never make the overall objective worse — you are, by definition, replacing one term with whatever minimizes the full objective over that term alone, holding everything else fixed. Under mild conditions (in particular, the featurewise smoothers behaving like well-defined linear projections), this monotone non-increasing sequence of objective values is bounded below by zero and is therefore guaranteed to converge, by the same "objective can't decrease forever" argument Module 2 used for coordinate descent and for EM — a different algorithm, the same proof skeleton, because both are instances of "optimize a piece at a time, holding the rest fixed."
Where this is used: this is not a simplified textbook sketch — it is exactly how production GAM-fitting software works internally. pyGAM and R's mgcv package both fit their smooth terms by cycling through features and refitting each one's smoother against the current partial residual, repeating until the fitted values stop moving — backfitting is the actual solver, not a theoretical warm-up for it.
Each half-iteration refits exactly one curve against the partial residual left after removing the other curve's current contribution (dots), while the other curve holds still — the highlighted panel shows which one is currently updating. Training MSE is recomputed live from whatever is on screen, and it only ever drops or holds flat, exactly as the derivation above claims.
The synthetic surface here has a real interaction: x₁'s effect on y grows in amplitude as x₂ increases. Drag the x₂ slider — the solid green curve (the true effect of x₁ at that x₂) visibly changes shape, but the dashed red curve — the best possible additive fit β₀+f₁(x₁)+f₂(x₂), found the same alternating-refit way backfitting works above — can only slide vertically, because f₁'s shape is frozen once fit. The RMSE readout tracks how far apart the two curves are at the chosen x₂.
#include <cmath>
#include <cstdio>
#include <vector>
#include <random>
// Same idea as the Python tab: average the target over nearby points along x.
std::vector<double> movingAverageSmoother(
const std::vector<double>& x,
const std::vector<double>& target,
double bandwidth
) {
int n = static_cast<int>(x.size());
std::vector<double> fitted(n, 0.0);
for (int i = 0; i < n; ++i) {
double sumW = 0.0, sumWT = 0.0;
for (int k = 0; k < n; ++k) {
if (std::fabs(x[k] - x[i]) <= bandwidth) {
sumW += 1.0;
sumWT += target[k];
}
}
fitted[i] = sumWT / sumW;
}
return fitted;
}
double meanOf(const std::vector<double>& v) {
double s = 0.0;
for (double val : v) s += val;
return s / v.size();
}
void backfit(
const std::vector<double>& x1,
const std::vector<double>& x2,
const std::vector<double>& y,
double bandwidth,
int maxIters,
double tol,
double& beta0,
std::vector<double>& f1,
std::vector<double>& f2
) {
int n = static_cast<int>(y.size());
beta0 = meanOf(y);
f1.assign(n, 0.0);
f2.assign(n, 0.0);
for (int iter = 0; iter < maxIters; ++iter) {
std::vector<double> totalBefore(n);
for (int i = 0; i < n; ++i) totalBefore[i] = beta0 + f1[i] + f2[i];
std::vector<double> r1(n);
for (int i = 0; i < n; ++i) r1[i] = y[i] - beta0 - f2[i];
f1 = movingAverageSmoother(x1, r1, bandwidth);
double shift1 = meanOf(f1);
for (double& v : f1) v -= shift1;
beta0 += shift1;
std::vector<double> r2(n);
for (int i = 0; i < n; ++i) r2[i] = y[i] - beta0 - f1[i];
f2 = movingAverageSmoother(x2, r2, bandwidth);
double shift2 = meanOf(f2);
for (double& v : f2) v -= shift2;
beta0 += shift2;
double maxChange = 0.0;
for (int i = 0; i < n; ++i) {
double totalAfter = beta0 + f1[i] + f2[i];
maxChange = std::max(maxChange, std::fabs(totalAfter - totalBefore[i]));
}
if (maxChange < tol) {
std::printf("Converged after %d iterations (max change=%.2e)\n", iter + 1, maxChange);
break;
}
}
}
int main() {
std::mt19937 rng(0);
std::uniform_real_distribution<double> unif(-3.0, 3.0);
std::normal_distribution<double> noise(0.0, 0.3);
int n = 300;
std::vector<double> x1(n), x2(n), y(n);
for (int i = 0; i < n; ++i) {
x1[i] = unif(rng);
x2[i] = unif(rng);
double trueF1 = 1.5 * std::sin(1.1 * x1[i]);
double trueF2 = 0.4 * x2[i] * x2[i] - 1.0;
y[i] = trueF1 + trueF2 + noise(rng);
}
double beta0;
std::vector<double> f1, f2;
backfit(x1, x2, y, 1.0, 25, 1e-6, beta0, f1, f2);
std::printf("beta0 = %.3f\n", beta0);
std::printf("fitted f1 sample: %.3f %.3f %.3f\n", f1[0], f1[1], f1[2]);
std::printf("fitted f2 sample: %.3f %.3f %.3f\n", f2[0], f2[1], f2[2]);
return 0;
}- Epidemiology and public health — modeling how disease risk depends on age, exposure level, and other continuous covariates, each allowed its own possibly nonlinear shape (risk rising sharply past a threshold age, say), while the model stays simple enough that a policy audience can be shown "here is exactly how risk changes with age, holding exposure fixed" as a single plotted curve.
- Credit risk scoring and insurance pricing — regulators in these industries typically require that a pricing or risk model's per-feature effects be auditable and explainable to a human reviewer, not just accurate. A GAM's separate, plottable
fⱼcurves are often exactly what satisfies that requirement in a way a black-box model with entangled feature interactions cannot. - Ecological and environmental modeling — species abundance or pollutant concentration as a sum of smooth effects of several separate covariates (temperature, rainfall, elevation), where ecologists specifically want to inspect each covariate's individual response curve rather than a single opaque prediction.
- Medical dose-response curves — modeling a patient outcome as a smooth, possibly non-monotonic function of drug dosage, additively combined with smooth effects of age or other continuous clinical measurements, while still being able to hand a clinician one readable dosage-response curve.
- Energy demand forecasting — daily electricity demand as a sum of smooth effects of temperature and of time-of-year, each fit independently and each interpretable on its own, which is exactly the situation the previous lesson's single-feature spline can't scale to without becoming either a black box or intractable to fit.
- Any setting that needs the previous lesson's per-feature flexibility at scale — whenever a plain linear model underfits a genuinely curved relationship for several features at once, but a fully general multivariate smooth surface would both overfit and become impossible to explain, a GAM is the standard middle path.
- Assuming a GAM's purely additive structure captures everything relevant. It fundamentally cannot represent an interaction — an effect of
x₁that itself depends on the value ofx₂— because the model is, by construction, a sum of functions that each see only one feature. Adding an explicit interaction smooth term,f₁₂(x₁, x₂), fixes this for the pairs of features that need it, but that term is itself a 2D (or higher) smooth surface — it reintroduces exactly the curse-of-dimensionality trouble a GAM was chosen to avoid, just scoped down to however many features are actually interacting instead of all of them at once. - Forgetting that a GAM has more hyperparameters to tune than the previous lesson's single smoothing spline, not fewer — every
fⱼgets its own smoothing penaltyλⱼ, so a model with 10 features has 10 separate wiggliness knobs to choose (usually via cross-validation or the automatic selection built intomgcv/pyGAM), not one. - Running backfitting without an identifiability constraint and being confused by component functions that "drift" — since only the SUM
f₁(x₁) + f₂(x₂)is pinned down by the data, a constant could in principle be added tof₁and subtracted fromf₂without changing any prediction at all, leaving each individual curve's vertical position ambiguous. The standard fix, used explicitly in the derivation and diagram above, is to constrain eachfⱼto average to zero over the training data and let the interceptβ₀absorb every constant term instead — that single constraint per feature is what makes each plotted curve's shape meaningful and reproducible rather than an arbitrary vertical shift.
Going deeper
The backfitting algorithm's convergence rate depends heavily on how correlated the features are — and this is a precise, not a vague, analogy to plain coordinate descent's behavior on ill-conditioned problems (Module 2). When two features x₁ and x₂ are strongly correlated, a large share of the variation in y that either f₁ or f₂ alone could explain is ambiguous between them — much like two nearly-collinear columns in a design matrix make a coordinate-descent step on one coefficient barely move the loss, because the other coefficient can absorb almost the same effect. Backfitting responds the same way: it takes many more cycles to settle when refitting f₁ on the partial residual only reveals a small amount that f₂ couldn't already explain, and vice versa. Just as ill-conditioning slows scalar coordinate descent without breaking its convergence guarantee, correlated features slow backfitting's convergence without breaking the monotone-non-increasing-objective argument from the derivation above — it still converges, just more slowly, and diagnosing "why is this GAM fit so slow" often comes down to checking feature correlations first.
A colleague fits a GAM with y = β₀ + f₁(x₁) + f₂(x₂) and finds it fits training data noticeably worse than a black-box model, even after tuning both smoothing penalties. Age (x₁) and income (x₂) are the two features. What is the single most likely structural explanation, and how would you fix it without abandoning the GAM's interpretability?
The most likely explanation is a real interaction: perhaps income's effect on the outcome genuinely differs by age bracket (e.g. extra income matters much more for young adults than for retirees). A purely additive model cannot represent that by construction, no matter how the two 1D smoothing penalties are tuned -- tuning wiggliness only changes the shape of each independent curve, not whether the model can represent effects that depend on both features jointly. The fix that preserves most of the interpretability is to add an explicit interaction smooth term f12(age, income) for just this one suspected pair, rather than reverting to a fully unrestricted multivariate model -- you pay the curse-of-dimensionality cost only for the one interaction that actually matters, and every other feature stays in its own clean, individually plottable additive term.
A GAM keeps the previous lesson's per-feature spline flexibility while staying interpretable, by restricting the model to a SUM of one-dimensional smooth functions instead of one multivariate surface. That restriction is what makes backfitting possible in the first place: holding every other feature's function fixed collapses the fitting problem for one feature down to exactly the single-feature smoothing-spline problem already solved, applied to a partial residual — which is why backfitting is coordinate descent over functions, inherits its convergence guarantee, and is genuinely what production software runs. The price for all of this is the additive assumption itself: a GAM cannot see interactions between features unless you explicitly add a term for them.