Chapter 18 of 37 · intermediate
Overfitting and the bias-variance tradeoff
What this chapter covers
Every model in this book has a flexibility knob. A decision tree has its depth, a polynomial has its degree, a neural network has its width. Turn the knob up and the model can fit more shapes; the training error drops, and it keeps dropping the harder you turn. The trap is that the number you actually care about — the error on data you haven't seen — does not keep dropping. It falls, bottoms out, and then climbs back up. That U-shaped curve is the single most important picture in supervised learning, and this chapter is about drawing it, understanding it, and naming the two forces that bend it.
The vehicle is polynomial regression, because it makes the knob concrete: fit a line, then a parabola, then a degree-11 curve, and watch each one work on the same data. We build the whole thing from scratch in NumPy — the feature expansion, the least-squares fit, the error metric — then run an empirical bias-variance decomposition, refitting on hundreds of bootstrap resamples to measure how much of the error is the model being too stiff to reach the truth (bias) and how much is it thrashing around chasing noise (variance). The centerpiece is an animation: sweep the degree from a straight line up to a wild wiggle, and watch training error fall toward zero while test error carves its U. This is the chapter the rest of the course keeps pointing back at — regularization, cross-validation, and every hyperparameter you'll ever tune are all attempts to find the bottom of that U.
A bit of history
The tension is old. In the fourteenth century William of Ockham argued that among competing explanations you should prefer the one with the fewest assumptions — Occam's razor, the principle that a simpler hypothesis is likelier to be right. That's the whole intuition behind fighting overfitting, six hundred years early: a model that needs a hundred wiggles to explain your data is probably explaining the noise, not the signal.
The modern statement is younger and sharper. The decomposition of prediction error into bias and variance was folklore in statistical estimation for decades, but the paper that pinned it to machine learning is Stuart Geman, Élie Bienenstock, and René Doursat's 1992 "Neural Networks and the Bias/Variance Dilemma" in Neural Computation. Their argument was that a flexible learner — a big neural net, in their case — can drive its bias to nearly zero, but it pays for that flexibility in variance: the fit swings wildly depending on which training sample it happened to see. You cannot have both at once for free, and choosing where to sit on that tradeoff is, they argued, the central problem of learning from finite data. The name "bias-variance dilemma" comes from that paper, and thirty years later it's still the frame everyone reaches for when a model looks perfect on the training set and falls apart in production.
The intuition
Here is the setup, stripped to one dimension so you can see everything. There's a smooth true function — a sine wave — and I've drawn 40 noisy samples from it. The model never sees the smooth curve; it only sees the dots. Its job is to recover the shape from the dots alone, and how hard it's allowed to try is set by the polynomial degree.
Three fits, same data. A degree-1 line, the degree the test error will pick, and a degree-11 polynomial:
The orange line is hopeless. A straight line cannot bend, so it splits the difference and misses the whole shape — it's too stiff, and it would miss the same way no matter which 40 points you handed it. That stiffness is bias. The purple curve is the opposite failure: it has enough freedom to chase individual points, so it lurches up and down to pass near noise that isn't part of the true signal at all. Hand it a different 40 points and it would lurch somewhere else entirely. That sensitivity is variance. The cyan curve, degree 5, is the one that found the sine without chasing the dots — and the point of this chapter is that "degree 5" isn't a guess, it's the measurable bottom of a curve we're about to draw.
The math
Fix a single test input . The target there is the true function plus noise: , where has mean zero and variance and is the irreducible randomness no model can predict. Our fitted model depends on the particular training set we happened to draw, so it's a random object; imagine drawing many training sets and refitting each time. The quantity that matters is the expected squared error at , averaged over that randomness. It splits into exactly three pieces:
Read the three terms one at a time. Bias is the gap between the average fit and the truth: is the curve you'd get by averaging the fits from all those training sets, and if that average still misses , the model is too rigid to represent the truth no matter how much data you give it. Variance is how much a single fit jitters around that average from one training set to the next — the spread of the cloud of curves. Noise, , is the floor: even a perfect model that nailed exactly would still be off by the noise on each new point, so no amount of cleverness drives the error below it.
The whole tradeoff lives in the first two terms, and they pull in opposite directions as you turn up flexibility. A stiff model (low degree) has high bias and low variance: it can't reach the truth, but it barely moves when the data changes. A flexible model (high degree) has low bias and high variance: its average fit can trace any shape, but any single fit is at the mercy of the particular noise it saw. Total error is their sum plus the floor, so minimizing it means finding the degree where their sum bottoms out — not where either one alone is smallest. That sum is the U.
The framework's reach, and its limits
The bias-variance decomposition is a lens, not an algorithm, and its value is that it turns a vague worry — "is my model overfitting?" — into two questions you can answer separately. If your model does badly on both training and test data, that's bias: it's too simple, and you need a more flexible model or better features, not more data. If it does great on training and poorly on test, that's variance: it's too flexible for the data you have, and the fixes are the opposite — simplify the model, regularize it, or get more data. Knowing which of the two you're fighting is the difference between a fix that helps and one that makes things worse, and I reach for this framing on nearly every model that isn't behaving.
Two honest caveats. First, the clean three-way split is exact only for squared-error regression; for classification and other losses the same story holds in spirit but the algebra isn't as tidy, and "bias" and "variance" become more qualitative. Second — and this is the interesting one — the strict U-curve is a statement about the classical regime where the model has fewer parameters than you have data. Modern deep networks routinely run with far more parameters than data points and, past the point where they interpolate the training set perfectly, test error can fall again — the "double descent" curve Belkin and coauthors documented in 2019. That doesn't repeal the tradeoff; it says the picture past the interpolation threshold is richer than one U. For everything in this book, and most models you'll ship, the single U is the right mental model, and it's the one worth internalizing first.
The data
The dataset is synthetic and seeded, committed as data/synthetic.csv, so the
whole chapter reproduces to the bit. The generative process is deliberately simple:
draw the input uniformly on , set the true target to
— one full period, real curvature — and add Gaussian noise with
standard deviation 0.25. The training set is 40 such points; a separate 4,000-point
draw from the same process is the test set that scores generalization. The noise
variance is 0.0625 in target units squared, and that number is the floor from the
math above: no model, however clever, gets its test error below it.
I chose sine on purpose. A straight line can't follow it, so degree 1 is forced to underfit and the bias is unmistakable. A middle degree traces one period cleanly. A high degree has more freedom than 40 noisy points can constrain, so it wiggles — and because I'm fitting on raw powers, the wiggle gets genuinely wild near the edges of the interval, which is overfitting you can see with your eyes. Here's the data on its own, points and the true curve the model has to recover from them:
Build it, one function at a time
Polynomial regression sounds like a new model, but it's linear regression wearing a costume. The trick is the feature expansion: a degree- polynomial in is a linear model in the columns . Build that matrix — the Vandermonde matrix — and every bit of least-squares machinery from the linear regression chapter applies unchanged. Degree is just how many columns you hand the solver, which is exactly the flexibility knob:
def vandermonde(x, degree, x_lo=0.0, x_hi=1.0):
"""Expand a 1-D input into polynomial features: columns 1, u, u^2, ..., u^p.
This is the whole trick that turns linear regression into curve fitting. A
degree-p polynomial c0 + c1 u + c2 u^2 + ... + cp u^p is a *linear* model in
the columns [1, u, u^2, ..., u^p], so if we build that design matrix we can
fit a curve with the exact machinery we'd use to fit a line. Row i is one
sample; column j is that sample raised to the power j. The result is
(len(x), degree + 1) — degree sets how many columns, i.e. how flexible the
model is allowed to be.
One numerical nicety: we first map the input from its domain [x_lo, x_hi]
onto [-1, 1]. Raw powers of x on [0, 1] get badly correlated at high degree
(x^10 and x^11 are nearly identical there), which makes the least-squares
system ill-conditioned; centering to [-1, 1] keeps the columns well separated
so the fit stays honest up to the degrees we push it to. It changes the basis,
not the model — the fitted curve is the same polynomial either way.
"""
x = np.asarray(x, dtype=float)
u = 2.0 * (x - x_lo) / (x_hi - x_lo) - 1.0
return np.vander(u, N=degree + 1, increasing=True)
The one wrinkle is the rescale to before taking powers. Raw powers of a number in get nearly indistinguishable at high degree — and are almost the same column — and that makes the least-squares system ill-conditioned, which is a numerical problem, not a modeling one. Centering the input keeps the columns separated so the fit stays honest all the way up to degree 11. It changes the basis, not the model.
With the design matrix in hand, the fit is one least-squares solve. We use
np.linalg.lstsq, which minimizes through an SVD — the
numerically stable route, which matters precisely because high-degree Vandermonde
matrices are the ill-conditioned case:
def fit_poly(x, y, degree):
"""Least-squares fit of a degree-p polynomial. Returns the coefficient vector.
Build the Vandermonde design matrix, then solve the linear least-squares
problem min_c ||V c - y||^2 for the coefficients c. We use np.linalg.lstsq,
which solves it through an SVD — the numerically stable way, because at high
degree the Vandermonde columns get badly correlated (x^11 and x^12 look
nearly identical on [0, 1]) and forming the normal equations directly would
lose precision. Same objective as ordinary least squares; only the number of
columns changes.
"""
V = vandermonde(x, degree)
coef, *_ = np.linalg.lstsq(V, y, rcond=None)
return coef
Prediction rebuilds the same feature expansion at the query points and takes the weighted sum. The degree is implied by the length of the coefficient vector, so a caller never passes it twice:
def predict_poly(coef, x):
"""Evaluate the fitted polynomial at new inputs: yhat = V(x) @ coef.
Rebuild the same feature expansion at the query points and take the weighted
sum. degree is implied by len(coef) - 1, so a caller never has to pass it
twice.
"""
degree = len(coef) - 1
return vandermonde(x, degree) @ coef
The metric is mean squared error, and I want to be precise about its unit because the whole chapter is told in it: MSE is the average of the squared misses, so its unit is the target's unit squared — "target units²" everywhere below. Squaring punishes big misses harder than small ones and is what makes the decomposition algebra work out:
def mse(y_true, y_pred):
"""Mean squared error, in target units squared.
The average of the squared misses. Squaring keeps sign out of it and
punishes big misses harder than small ones; the unit is the target's unit
squared, which is what "target units^2" means everywhere in this chapter.
This is the number the whole overfitting story is told in.
"""
y_true = np.asarray(y_true, dtype=float)
y_pred = np.asarray(y_pred, dtype=float)
return float(np.mean((y_true - y_pred) ** 2))
Now the piece that does the real work of this chapter: the empirical bias-variance decomposition. The math above is an expectation over training sets, and here we estimate that expectation the honest way — by actually resampling. Draw many bootstrap resamples of the training set, fit a degree- polynomial on each, and predict at a grid of evaluation points. The spread of those predictions is the variance; the gap between their average and the true function is the bias:
def bias_variance_decompose(x_train, y_train, x_eval, f_true, degree,
n_boot, noise_var, rng):
"""Empirically split test error into bias^2, variance, and noise at one degree.
The recipe follows the decomposition literally. Draw n_boot bootstrap
resamples of the training set (sample the rows with replacement), fit a
degree-p polynomial on each, and predict at every point of a fixed evaluation
grid x_eval. That gives a cloud of predictions at each grid point — one per
resample.
- bias^2 = mean over x of (average prediction - true f(x))^2.
How far the *typical* fit sits from the truth. High for a stiff
model that can't reach the target shape.
- variance = mean over x of the variance of the predictions across resamples.
How much the fit jitters when the training data changes. High
for a flexible model that chases noise.
- noise = the irreducible term, the noise variance baked into the data.
No model can beat it; we pass it in because we know it here.
Returns (bias2, variance, noise, total) where total = bias2 + variance +
noise is the model's expected squared error on a fresh point.
"""
n = len(x_train)
preds = np.empty((n_boot, len(x_eval)))
for b in range(n_boot):
idx = rng.integers(0, n, size=n) # bootstrap resample, with replacement
coef = fit_poly(x_train[idx], y_train[idx], degree)
preds[b] = predict_poly(coef, x_eval)
mean_pred = preds.mean(axis=0) # average fit at each grid point
bias2 = float(np.mean((mean_pred - f_true(x_eval)) ** 2))
variance = float(np.mean(preds.var(axis=0)))
noise = float(noise_var)
total = bias2 + variance + noise
return bias2, variance, noise, total
That function is the experiment the whole tradeoff rests on. Everything else is plotting.
Watch it work
Here is the picture the chapter exists for. Press play and the polynomial degree sweeps from 1 to 11, one real refit per frame. The top panel is the fitted curve over the data, with the true sine dashed behind it. The bottom panel builds the two error curves as the degree climbs — training error in cyan, test error in purple — with a dot marking the current degree. The caption names the degree and both errors in target units².
Watch the two panels together, because the story is in their disagreement. At degree 1 the fit is a straight line that misses the sine badly, and both errors are high — that's pure bias, the model too stiff to bend. As the degree climbs into the middle range the curve snaps onto the sine, and both errors drop together to their lowest point around degree 5: training error 0.041, test error 0.061 in target units². Then the panels split. Training error keeps sliding — 0.038 by degree 11, an ever-closer fit to the 40 points it can see. But test error turns and climbs, gently at first, then steeply: 0.079 at degree 7, 0.154 at degree 8, and 0.817 by degree 11. Up top you can see why — the curve stops tracing the sine and starts lunging toward individual dots, flinging itself off the panel near the edges to thread points that are mostly noise. That widening gap between the cyan and purple lines is overfitting, drawn in real numbers. The model isn't learning anymore; it's memorizing.
The bottom panel is a single realization, though — one training set, one test set — and a single realization is a little jagged. The decomposition smooths it into its two components. Refit on 400 bootstrap resamples at each degree, measure the bias² and variance of the predictions, and add the noise floor, and the U resolves into the two forces underneath it:
This is the bias-variance dilemma as a plot, and the y-axis is log because the two components span orders of magnitude. Bias² (orange) starts high — 0.123 at degree 1, the stiff line's inability to reach the sine — and collapses as soon as the model has enough degree to trace the shape, flattening near 0.003 from degree 3 on. Variance (purple) does the mirror image: it starts near nothing, 0.012 at degree 1, and climbs as flexibility grows, gently through the middle degrees and then explosively — 0.29 at degree 9, 1.6 by degree 11 — as the fit starts swinging wildly between resamples. They cross around degree 4. Total error (cyan) is their sum plus the 0.0625 noise floor, and it's the U: down while bias falls faster than variance rises, up once variance takes over. The bottom of that U, degree 3 here, is the sweet spot the whole apparatus exists to find. Left of it you're bias-limited; right of it, variance-limited.
The full implementation
The whole file, no library, top to bottom — the feature expansion, the fit, the prediction, the metric, and the bootstrap decomposition. This is the code the animation and both charts actually ran:
"""Polynomial regression and the bias-variance decomposition, from scratch.
This chapter is about a single knob — model flexibility — and what it does to
generalization. The vehicle is polynomial regression: fit y with a polynomial of
degree p, turn p up, and watch the fit go from too stiff to draw the data (a
line through a curve) to so loose it threads every noisy point. Somewhere in the
middle is the fit that generalizes, and the whole point is that "somewhere in
the middle" is not a slogan — it's a measurable minimum of the test error.
The mechanism underneath is linear regression. A degree-p polynomial is a linear
model on an expanded feature set: instead of the single column x you fit on the
columns 1, x, x^2, ..., x^p. That expansion is a Vandermonde matrix, and once
you have it the fit is the same least-squares solve as ordinary regression. So
nothing here is new machinery — the flexibility all lives in how many columns
you hand the solver.
The second half is the bias-variance decomposition: refit the model on many
resampled training sets, look at the spread of the predictions, and split the
expected error into bias (how far the average fit sits from the truth), variance
(how much the fit jitters from sample to sample), and the irreducible noise. As
degree climbs, bias falls and variance rises; where they cross is the sweet spot.
Pure NumPy — no ML library in this file. Every function appears in the chapter
one region at a time via the include directives.
"""
import numpy as np
# region: vandermonde
def vandermonde(x, degree, x_lo=0.0, x_hi=1.0):
"""Expand a 1-D input into polynomial features: columns 1, u, u^2, ..., u^p.
This is the whole trick that turns linear regression into curve fitting. A
degree-p polynomial c0 + c1 u + c2 u^2 + ... + cp u^p is a *linear* model in
the columns [1, u, u^2, ..., u^p], so if we build that design matrix we can
fit a curve with the exact machinery we'd use to fit a line. Row i is one
sample; column j is that sample raised to the power j. The result is
(len(x), degree + 1) — degree sets how many columns, i.e. how flexible the
model is allowed to be.
One numerical nicety: we first map the input from its domain [x_lo, x_hi]
onto [-1, 1]. Raw powers of x on [0, 1] get badly correlated at high degree
(x^10 and x^11 are nearly identical there), which makes the least-squares
system ill-conditioned; centering to [-1, 1] keeps the columns well separated
so the fit stays honest up to the degrees we push it to. It changes the basis,
not the model — the fitted curve is the same polynomial either way.
"""
x = np.asarray(x, dtype=float)
u = 2.0 * (x - x_lo) / (x_hi - x_lo) - 1.0
return np.vander(u, N=degree + 1, increasing=True)
# endregion
# region: fit_poly
def fit_poly(x, y, degree):
"""Least-squares fit of a degree-p polynomial. Returns the coefficient vector.
Build the Vandermonde design matrix, then solve the linear least-squares
problem min_c ||V c - y||^2 for the coefficients c. We use np.linalg.lstsq,
which solves it through an SVD — the numerically stable way, because at high
degree the Vandermonde columns get badly correlated (x^11 and x^12 look
nearly identical on [0, 1]) and forming the normal equations directly would
lose precision. Same objective as ordinary least squares; only the number of
columns changes.
"""
V = vandermonde(x, degree)
coef, *_ = np.linalg.lstsq(V, y, rcond=None)
return coef
# endregion
# region: predict_poly
def predict_poly(coef, x):
"""Evaluate the fitted polynomial at new inputs: yhat = V(x) @ coef.
Rebuild the same feature expansion at the query points and take the weighted
sum. degree is implied by len(coef) - 1, so a caller never has to pass it
twice.
"""
degree = len(coef) - 1
return vandermonde(x, degree) @ coef
# endregion
# region: mse
def mse(y_true, y_pred):
"""Mean squared error, in target units squared.
The average of the squared misses. Squaring keeps sign out of it and
punishes big misses harder than small ones; the unit is the target's unit
squared, which is what "target units^2" means everywhere in this chapter.
This is the number the whole overfitting story is told in.
"""
y_true = np.asarray(y_true, dtype=float)
y_pred = np.asarray(y_pred, dtype=float)
return float(np.mean((y_true - y_pred) ** 2))
# endregion
# region: bias_variance
def bias_variance_decompose(x_train, y_train, x_eval, f_true, degree,
n_boot, noise_var, rng):
"""Empirically split test error into bias^2, variance, and noise at one degree.
The recipe follows the decomposition literally. Draw n_boot bootstrap
resamples of the training set (sample the rows with replacement), fit a
degree-p polynomial on each, and predict at every point of a fixed evaluation
grid x_eval. That gives a cloud of predictions at each grid point — one per
resample.
- bias^2 = mean over x of (average prediction - true f(x))^2.
How far the *typical* fit sits from the truth. High for a stiff
model that can't reach the target shape.
- variance = mean over x of the variance of the predictions across resamples.
How much the fit jitters when the training data changes. High
for a flexible model that chases noise.
- noise = the irreducible term, the noise variance baked into the data.
No model can beat it; we pass it in because we know it here.
Returns (bias2, variance, noise, total) where total = bias2 + variance +
noise is the model's expected squared error on a fresh point.
"""
n = len(x_train)
preds = np.empty((n_boot, len(x_eval)))
for b in range(n_boot):
idx = rng.integers(0, n, size=n) # bootstrap resample, with replacement
coef = fit_poly(x_train[idx], y_train[idx], degree)
preds[b] = predict_poly(coef, x_eval)
mean_pred = preds.mean(axis=0) # average fit at each grid point
bias2 = float(np.mean((mean_pred - f_true(x_eval)) ** 2))
variance = float(np.mean(preds.var(axis=0)))
noise = float(noise_var)
total = bias2 + variance + noise
return bias2, variance, noise, total
# endregion
The library version
Nobody hand-rolls a Vandermonde matrix in production. scikit-learn splits the same
job into two composable steps: PolynomialFeatures builds the expanded design
matrix, and LinearRegression solves the least-squares fit on it. Chain them in a
pipeline and you have exactly our fit_poly plus predict_poly:
def sklearn_fit(x, y, degree):
"""Fit a degree-p polynomial with PolynomialFeatures + LinearRegression.
Returns the fitted pipeline. include_bias=True gives the constant column, so
LinearRegression runs with fit_intercept=False — the twin of our fit_poly,
same design matrix, same least-squares solve. We rescale the input to [-1, 1]
first, matching impl.vandermonde, so the coefficients line up one for one.
"""
model = make_pipeline(
PolynomialFeatures(degree=degree, include_bias=True),
LinearRegression(fit_intercept=False),
)
model.fit(_rescale(x), y)
return model
Prediction runs the same pipeline forward:
def sklearn_predict(model, x):
"""Predict with the fitted pipeline at new inputs."""
return model.predict(_rescale(x))
The one subtlety is the intercept. PolynomialFeatures already emits the constant
column, so we set fit_intercept=False to avoid a duplicated bias term, and we
rescale the input to first to match our Vandermonde exactly. With those
two choices the pipeline's coefficient vector lines up one for one with ours, and
since both use the same SVD least-squares solver underneath, they land on the same
curve at every degree — the coefficients agree to about 4e-15 and the predictions to
1e-12. There's nothing for them to disagree about; it's the same linear algebra.
Scratch versus library
Because the two implementations are the same model, the face-off isn't really scratch versus library — it's a check that we built the reference correctly, and we did. At the sweet-spot degree 5 both reach a test MSE of 0.061 in target units², identical to the precision quoted above. The interesting comparison is the one the whole chapter is about: the cost of picking the wrong degree. Here is the held-out test error at the sweet spot next to the same error at the overfit degree, both in target units²:
The overfit degree-11 model scores 0.817 against the sweet spot's 0.061 — 13 times
the error, on the same data, from the same model family, bought purely by turning the
flexibility knob too far. And here's the part that makes overfitting insidious: on the
training data the degree-11 model looks better, 0.038 versus 0.041, so if you judged
it by the numbers it hands you on the data it trained on, you'd ship the worse model
and feel good about it. The only defense is to score on data the model never saw,
which is why every number in this chapter that matters is a test number, and why the
test set is 100 times larger than the training set — so the U it traces is smooth
enough to trust. The exact figures live in results.json, regenerated whenever the
code changes, so the prose and the charts can't drift from what the code produced.
Takeaways
The U-curve is universal. It doesn't care that we used polynomials — swap in tree depth, network width, the number of neighbors in k-NN, the number of training steps before you stop, and you get the same shape: training error sliding toward zero while test error falls, bottoms out, and climbs. Internalize that picture and a huge amount of practical machine learning becomes one question: where is the bottom of my U, and which side of it am I on? Underfitting and overfitting aren't two unrelated problems; they're the two ends of a single dial, and the whole game is landing in the middle.
Once you know which side you're on, the moves are opposite, and this is the payoff of keeping bias and variance separate in your head. If you're bias-limited — bad on train and test both — more data won't help; you need a more flexible model or better features. If you're variance-limited — great on train, poor on test — the fixes all trade a little variance for a little bias: simplify the model, get more data to pin the flexible one down, or regularize it, which is the next thread in this book. Ridge and lasso are exactly this trade made continuous — instead of choosing a degree, you add a penalty that pulls the coefficients toward zero and dial its strength to slide along the U. Cross-validation is how you find the bottom without a held-out set to burn. They're all the same idea wearing different clothes.
One caveat to carry forward so it doesn't surprise you later: the tidy single U is the classical regime, where the model is smaller than the data. Very large models — modern deep networks especially — can push past the point of fitting the training set perfectly and see test error improve again, the double-descent phenomenon. It's real and it's worth knowing about, but it doesn't change the advice for anything you'll build in this course. For polynomials, trees, forests, linear models, and every tabular problem you're likely to meet, the single U is the map, and finding its bottom is the job.