Chapter 4 — Generalization: Overfitting, Bias–Variance and Regularization

Chapter 3 ended with a classifier that fits ten loan files and a 20,000-applicant book, and a question: will it work on next year's applicants? That is the question of generalization, and it applies to every model in this book. By the end of this chapter you can say exactly why a model that wins on its training data can lose on new data; split its error into three parts and compute each by hand; measure how flexible a model really is; read a learning curve; derive what ridge, lasso, elastic net, early stopping, noise injection and dropout do to a fit; tune any knob with cross-validation and report an honest number; explain why billion-parameter models that fit their training data perfectly still generalize; and spot the shifts and leaks that make a test score lie. One idea carries all of it: a model is judged on days it has never seen, and its error there splits into noise you can't beat, bias you chose and variance you paid for — every cure for overfitting trades a little bias for a lot less variance.

We study it on a regression curve, because with squared error the answer splits exactly into three parts we can compute and plot; Section 3 says what changes for classifiers. Three earlier pieces are used, not re-taught: train/test optimism, polynomial features and the ridge/lasso geometry (Ch 1, Sections 9–10); gradient descent, implicit regularization and AdamW (Ch 2, Sections 4, 6, 8); and penalized logistic regression with C = 1/\lambda, where shrinkage pulls probabilities toward the base rate (Ch 3, Sections 9 and 12).

How to read it. The chapter is a chain of 14 questions in three acts: The problem (Sections 1–4: what generalization is, the U-curve, the decomposition, capacity), The cures (Sections 5–9: more data, ridge, lasso and elastic net, cross-validation, the other regularizers) and The modern view (Sections 10–14: double descent, shift and leakage, then building, using and scaling it). Each section opens with its question and a 30-second version, and ends with the question it leaves open.

Two passes. On the first pass, read only the question, the 30-second version and the figures of each section. You'll follow the whole arc in under an hour. On the second pass, read everything, doing the Remember box at the end of each section from memory before you look. Practice lives in its own tabs: Derivations (what to prove on a whiteboard) and Interview Bank (every question by level, with the intuition behind each answer).

1. What Is Generalization? Training Error vs the Error That Matters

The question: four curves fit last season's six days of bike rides. Which one will be closest on July 20 next year — and how can we know before next year arrives?

The 30-second version: what we care about is risk: the average loss on new days drawn from the same process. What we can compute is empirical risk: the average loss on the days we trained on. The difference is the generalization gap, and on average it is never negative, because the fit chose its curve to look good on exactly those days. On our six days the degree-5 curve has training error exactly 0, yet on next year's six days it is six times worse than the cubic (0.258 vs 0.041). The fix for measuring it: hold out data the model never saw. Even a perfect model can't beat the noise floor \sigma^2, and what's left above it splits three ways — Family, Sample, Search: the family can't reach the truth, six days can't pin down the best member, and the optimizer may not find it. First read: the toy table, risk vs empirical risk, Family-Sample-Search.

The problem

A city bike-share operator plans its rebalancing trucks and staff from a forecast of daily rides. The target is y = daily rides minus the annual average, in thousands (0.5 means 500 rides above an average day). The only feature is time of year, measured from the spring equinox: x \in[0, 1), with x = 0 on March 20, 0.25 near the summer peak (June 20), 0.5 on September 20 and 0.75 near the winter trough (December 20). For hand work we use m = 12x, months since March 20.

Last season's pilot logged one day every two months. Next year's planning days sit between them:

The planning question for the whole chapter: how many rides above average should we staff for on July 20 next year (T3)? We'll use these twelve days throughout, so every formula can be checked by hand.

Four curves through six days

Fit polynomials of degree 0 (a flat line), 1 (a line), 3 (a cubic) and 5 by least squares. Six points and six coefficients means degree 5 passes through every training day.

Training error falls at every step and hits zero at degree 5. Test error falls to the cubic and then rises again: the degree-5 curve's test error is 6.27 times the cubic's. Look at Figure 4.1(a): to pass through all six points the degree-5 curve bends hard between them and shoots to 1.83 at the left edge of the plot, half a month before the season starts. On July 20 the cubic and degree 5 look equally good (0.791 vs 0.801), so one day can't tell them apart; next year's six days can. (Degree 4 also nearly interpolates, train MSE 0.0005, by a coincidence of these six values.)

Vocabulary you'll be asked about

Two traps in that table. First, "generalization error" means the risk in some texts and the gap in others; say which. Second, a "test error" computed on one test set is itself an estimate of the risk, with a standard error, as the next subsection shows.

Risk vs empirical risk

R

risk: a curve's average loss on new days drawn from the same process

f

a candidate curve: a rule turning time of year into predicted rides

\mathbb{E}_{(X,Y)}

average over a fresh random day: its time of year X and its rides Y

X

a new day's time of year, random (capital = random variable)

Y

that new day's rides above average, in thousands; random too

\ell

the loss for one prediction; here squared error (y - \hat y)^2

\hat R_n

empirical risk: the curve's average loss on the n training days (hat = estimated from data)

n

the number of training days (6 in the toy)

\sum_{i=1}^{n}

add up over the training days i = 1 to n

y_i

the rides logged on training day i

x_i

the time of year of training day i

For any fixed curve f chosen before seeing the data, the training error is an unbiased estimate of its risk: \mathbb{E}[\hat{R}_n(f)] = R(f). The trouble is that we don't evaluate a fixed curve; we evaluate \hat{f}, the curve the data chose. ERM picks \hat{f} because it has the lowest training error in the family, so its training error is the minimum of many noisy scores, and a minimum of noisy scores is biased low:

\mathbb{E}

average over many possible training seasons, each with fresh noise

\hat f

the curve ERM picked: lowest training error in the family, chosen after seeing the data

f_{\mathcal F}

the best member of the family \mathcal F: the one with the lowest risk on new days

Here f_{F} is the best member of the family. The first inequality holds because \hat{f} beats f_{F} on the training data by definition; the second because f_{F} is the best on new data. So the expected gap is never negative (Card 1). The more curves the family offers, the deeper that minimum dips, and the larger the gap.

The test set is an estimate too

A test MSE is an average of squared errors, so it has a standard error: the standard deviation of the squared errors divided by \sqrt{n_{\text{test}}}.

\widehat{\text{MSE}}_{\text{test}}

the test MSE we measured on next year's days: itself an estimate of the risk

y_j

the rides on test day j

x_j

the time of year of test day j

n_{\text{test}}

the number of test days (6 next year)

The cubic's 0.0411 on six test days has an SE of 0.0078 — about 19% of the estimate. Later we simulate 10,000 test days, where the same cubic-type fit scores 0.1137 with an SE of 0.0016. Report test errors with their SE, especially when two models are close (Q2).

Family, Sample, Search

Above the noise floor, the excess risk of what we actually ship splits into three parts (Bottou and Bousquet, 2008). With f^* the Bayes predictor, f_{F} the best member of the family, \hat{f}_n the exact ERM fit on n rows, and \tilde{f} the fit the optimizer actually returned:

\tilde f

the fit the optimizer actually returned: the curve we ship (tilde = what we really got)

f^*

the Bayes predictor: the best possible curve, \mathbb{E}[Y \mid X = x] (star = the optimum)

\hat f_n

the exact ERM fit on n training days: what a perfect optimizer would return

A richer family shrinks the first term and grows the second; that tension is the whole chapter (Q1, Q3). Search error is zero for closed-form fits and becomes a regularizer in its own right when we stop early (Section 9).

Same story for classifiers. Nothing above depends on squared error. For the loan classifier of Ch 3, \ell is 0-1 loss (accuracy) or log loss, risk is the error rate or cross-entropy on next year's applicants, and Ch 3's held-out split was an estimate of it. What is special about squared error is the exact three-way split of Section 3.

Remember: we're graded on new days, not old ones. Risk vs empirical risk; the gap is the price of choosing your curve after seeing the data, so it's never negative on average. A test score has an SE. Above the noise floor: Family, Sample, Search.

Next question → why does a more flexible curve keep winning on the training days but eventually lose on new ones? (Section 2)

2. Underfitting and Overfitting: the U-Curve

The question: as we raise the polynomial degree from 0 to 20 on 30 days of rides, training error keeps falling. Why does test error fall, then rise?

The 30-second version: a flexible family contains every curve a simpler one does, so its training error can only go down — that's a guarantee, not good news. Test error falls while extra flexibility captures real shape (degree 0 → 3: 0.582 → 0.106), then rises once the extra flexibility fits the noise (degree 15: 0.379; degree 18: 8.54). Too stiff = underfitting: training and test error are both high. Too flexible = overfitting: training error is low and the gap is large. Hook: train tells bias, gap tells variance. Ch 1's optimism formula predicts the training side once the curve is flexible enough: \mathbb{E}[\text{train MSE}] = \sigma^2(1 - p/n) (degree 9: 0.0612 simulated vs 0.0600). First read: Figure 4.1(b), the diagnosis table.

The season data

Six days are good for hand checks and bad for statistics. From here on we also simulate a full season: 30 days, one every 12 days or so (x_i = (i - 0.5)/30), with the same truth and noise as the toy (revealed in Section 3: f(x) = \sin(2\pi x), \sigma = 0.3). Because it's simulated, we can do what no real team can: rerun the season 500 times with fresh noise, and score every fit on 10,000 new days drawn uniformly over the year. Replicate 0 is "this season", the one dataset a real team would have. All numbers come from sims/sim_season.py (seed 404).

We fit degrees 0–20 by least squares on Legendre polynomial features instead of raw powers x, x^2, \ldots. In exact arithmetic the fits are identical (both span the same polynomials), but raw powers of x \in[0, 1] are nearly collinear — x^{14} and x^{15} are almost the same column — so the solve loses digits. Legendre polynomials are close to orthogonal on the interval, which keeps the condition number small (Ch 1, Section 3).

Why training error must fall. The degree-d family is nested inside the degree-(d + 1) family: set the new coefficient to zero and you're back. Least squares finds the best member, so it can always do at least as well as before. Training error is non-increasing in flexibility for any nested family fitted by ERM — trees with more leaves, networks with more units, kNN with smaller k. It says nothing about new data.

The U

(Means over 500 seasons. Degrees 16 and up have heavy-tailed test errors — a few seasons explode at the edges — so their means are dominated by those seasons; Figure 4.1(b) uses a log scale.)

Underfit vs overfit: the diagnosis

The hook works because training error mostly measures how far the family is from the truth (it can't fit the shape), while the gap mostly measures how much the fit chased this sample's noise. Section 3 makes both halves exact (Q5).

Optimism, recapped and checked

Ch 1 showed that 10 parameters on 50 rows give expected training MSE 0.80\sigma^2 and test MSE 1.20\sigma^2 (Ch 1, Drill 30). The general statement, for a correct linear model with p coefficients, fixed design and new noise at the same x's (Card 3):

\mathbb{E}

average over many reruns of the season with fresh noise

\sigma^2

noise variance: day-to-day luck no curve can predict (0.09 here)

p

the number of fitted coefficients (10 for degree 9)

n

the number of training days (30 in the season)

x

the training days' times of year; the test here reuses the same ones

Check it on the season. Degree 9 has p = 10 coefficients and its bias is negligible, so theory says train MSE = 0.09 \times(1 - 10/30) = 0.0600; the simulation gives 0.0612. Theory for test is 0.09 \times(1 + 10/30) = 0.120; on random new days the simulation gives 0.124, close because the 30 design days fill the year evenly. At degree 3 the theory says 0.078 for training but the simulation gives 0.083: the formula assumes the model is correct, and the cubic still carries a little bias (Section 3).

For random x (each new day at a new, random time), the test error is larger. With Gaussian features and an intercept-free model it is \sigma^2(1 + p/(n - p - 1)), which blows up as p approaches n. Section 6 meets this on the station data and Section 10 explains the blow-up (Card 3, Q15, Drill 2).

Wild edges

Why do high degrees fail so badly? A polynomial that must pass near 30 noisy points uses its spare coefficients to bend, and the bending is cheapest where data are thinnest: between points and beyond the outermost ones. Over 500 seasons, the degree-15 fit's variance is 0.037 on June 20 (x_0 = 0.25, mid-data) and 1.61 on March 16 (x_0 = 0.99, past the last design day) — 44 times larger. The same effect drives extrapolation failures under shift (Section 11; Q4).

Remember: training error always falls with flexibility; test error is a U. Train tells bias, gap tells variance. Optimism for least squares: \sigma^2(1 - p/n) on train, \sigma^2(1 + p/n) on test, gap 2\sigma^2p/n. Polynomials go wild at the edges.

Next question → "too stiff" and "too wiggly" are words. Can we split the test error into exactly those two pieces? (Section 3)

3. The Bias–Variance Decomposition

The question: if we could rerun the season 500 times, how much of each curve's error on July 20 is aim, how much is wobble, and how much is luck we can never remove?

The 30-second version: treat the fitted curve as random — a different season gives a different curve. At any day x_0, expected squared error = \sigma^2 (noise) + \text{bias}^2 (how far the average curve is from the truth) + variance (how much the curves scatter around their average). Hook: Noise you can't beat, Bias you chose, Variance you paid for. The proof is one add-and-subtract and two cross terms that vanish. On the toy, by hand, the line at July 20: 0.09 + 0.274 + 0.020 = 0.384; the cubic: 0.09 + 0.001 + 0.052 = 0.143. Over 500 seasons, degree 1 is almost all bias (0.199 vs 0.006), degree 15 almost all variance (\approx 0 vs 0.288). First read: Figure 4.2(a–b), the hand computation, the hook.

Now we reveal the truth: f(x) = \sin(2\pi x), so July 20 (x = 1/3) has f = 0.866, and the noise is \varepsilon \sim N(0, 0.3^2), \sigma^2 = 0.09. Imagine many seasons, each with the same 30 days and fresh noise. Each season gives a different fitted curve \hat{f}, so \hat{f}(x_0) is a random variable. Its average over seasons, \mathbb{E}[\hat{f}(x_0)], is where the procedure aims; its spread is how steady it is. This is Ch 1's estimator picture (bias = aim, variance = steady hand; Ch 1, Section 5), applied to a prediction instead of a coefficient, plus one new term for the new day's own noise.

The derivation

A new day at x_0 has y = f(x_0) + \varepsilon with \mathbb{E}[\varepsilon] = 0, \operatorname{Var} (\varepsilon) = \sigma^2, and \varepsilon independent of the training data. Write f, \hat{f} and \bar{f} = \mathbb{E}[\hat{f}(x_0)] for short. Add and subtract f and \bar{f} (Card 2):

y

the rides on a new day at x_0 (e.g. July 20)

\hat f

the fitted curve's prediction at x_0: random, because it depends on which season we got

f

the true curve's value at x_0, \sin(2\pi x_0)

\bar f

the average prediction over all possible seasons, \mathbb{E}[\hat f(x_0)] (bar = average): where the procedure aims

\varepsilon

the new day's noise: luck with mean 0 and variance \sigma^2

Square and take expectations over the training data and the new noise. There are three squares and three cross terms. The cross terms with \varepsilon vanish because the new day's noise is independent of the training data and has mean zero. The cross term between bias and wobble vanishes because the bias is a constant and \mathbb{E}[\bar{f} - \hat{f}] = 0 (Q13):

\mathbb{E}

average over training seasons and the new day's noise

x_0

the day we forecast, e.g. July 20 (x_0 = 1/3)

\sigma^2

noise variance, 0.09: the floor no curve can beat

\text{Var}

variance: how much the prediction scatters across seasons

Noise you can't beat, Bias you chose, Variance you paid for. The noise is set by the world. The bias is set by your choice of family and penalty. The variance is the price of flexibility, paid in data.

By hand on the toy, at July 20

Keep the six training days fixed and imagine fresh noise. For least squares, \hat{y} is linear in y, so \mathbb{E}[\hat{y}] is the fit to the true values f, and the variance follows from Ch 1's prediction-variance formula (Ch 1, Section 8; Q14).

The line. Fitting a line to the true values (0.5, 1, 0.5, −0.5, −1, −0.5) gives slope -12/70 = -0.1714 per month and intercept 1.0286, so at m = 4 the line aims at \mathbb{E}[\hat{y}] = 0.343. The truth is 0.866, so \text{bias} = -0.523, \text{bias}^2 = 0.274. With centered months (m - 6 ranges over −5, −3, −1, 1, 3, 5; S_{\text{mm}} = 70):

\hat y(4)

the fitted line's predicted rides at month m = 4 (July 20)

n

the number of training days, 6

S_{mm}

the spread of the training months, \sum (m_i - 6)^2 = 70

\hat y

the line's prediction for July 20

The same computation for the other degrees (toy_* keys; each confirmed by 200,000 simulated noise draws):

Degree 5 aims almost perfectly (\text{bias}^2 \approx 10^{-6}) and still loses to the cubic, because it pays 0.083 in variance against the cubic's 0.052. In Section 1 the cubic and degree 5 gave nearly the same July 20 prediction on this season (0.791 vs 0.801); the decomposition says that over many seasons degree 5 lands farther from the truth on average. One lucky season can't show that.

The 500-season simulation

On the season, average the three terms over the 10,000 test days:

The sum matches the simulated test MSE within 0.00015 for degrees 1–15. At degree 0 the match is only within 0.006, because the 10,000 test days' own noise happens to correlate slightly with degree 0's huge bias. For degree 15 the simulation reports \text{bias}^2 = 0.001, but that is mostly variance leaking in: an average of 500 noisy curves is itself noisy by about \operatorname{Var} /500, and the exact \text{bias}^2 is \approx 0. Degree 1 is almost all bias; degree 15 is almost all variance; degree 3 sits at the bottom of the U, where \text{bias}^2 has collapsed and variance is still small.

Bias and variance depend on where you look

The table averages over the year. Point by point:

Variance is largest at the edges for every degree, and enormously so for degree 15. The cubic's bias is also concentrated at the edge, where a cubic can't quite match the sine's turn back to zero.

Three traps

1. Bias–variance describes a procedure, not a curve. "This fitted model has high variance" is loose talk. Variance belongs to the recipe — family, n, penalty, optimizer — run on many datasets. One fitted curve has an error, not a variance.
2. You can't compute it on real data. You never see f or 500 seasons. The proxies are bootstrap refits (to see variance: how much do predictions move when you resample rows?) and learning curves (Section 5: the plateau shows bias, the gap shows variance).
3. This "bias" is neither the intercept nor data bias. ML calls the intercept a bias term (Ch 1, Section 2), and Section 11's biased samples are a third meaning. Here bias is \mathbb{E}[\hat{f}(x_0)] - f(x_0).

A second knob: k-nearest neighbors

k-nearest neighbors (kNN) regression predicts y at x_0 by averaging the y's of the k training days closest to x_0. It has no coefficients at all, yet it has a complexity knob. Averaging k independent noisy values gives variance \sigma^2/k; reaching further out to collect them gives bias, because days farther from x_0 have different true values (Card 5):

\hat f(x_0)

the kNN forecast for day x_0

k

the number of neighboring days averaged

\sum_{i \in N_k(x_0)}

add up over the k training days closest to x_0

y_i

the rides on training day i

x_i

the time of year of training day i

f

the true curve, \sin(2\pi x), evaluated at the day in brackets

On the season (sims/sim_knn.py), k = 1 has variance 0.091 and \text{bias}^2 0.002 (test 0.183); k = 30 predicts the season mean everywhere (\text{bias}^2 0.495, test 0.582); k = 5 is best (test 0.124, \text{bias}^2 0.015, variance 0.018). Small k is flexible, large k is stiff: the knob runs the opposite way to degree (Drill 1). Ch 8 owns kNN; here it's a second knob for Section 4.

Beyond squared error

For 0-1 loss there is no additive decomposition. Worse, variance can reduce error: if the average prediction at x_0 is on the wrong side of the decision boundary, a fit that wobbles sometimes lands on the right side (Friedman, 1997; Domingos, 2000). For log loss an analogous decomposition exists (Heskes, 1998). The trade-off intuition survives — flexibility buys aim and costs steadiness — but the arithmetic is special to squared error (Q38).

Remember: expected test error = \sigma^2 + \text{bias}^2 + \text{variance}. Noise you can't beat, Bias you chose, Variance you paid for. The cross terms die because the new noise is independent and the wobble has mean zero. It describes a procedure, not a curve. Variance peaks at the edges.

Next question → "flexibility" moved bias and variance. How do we measure flexibility when parameter count gives the wrong answer? (Section 4)

4. Capacity: Degrees of Freedom, VC Dimension, and Why Counting Parameters Misleads

The question: a degree-15 polynomial with a penalty, kNN with k = 5, and a one-parameter sine wave all break the "count the parameters" rule. How do we measure how flexible a model really is?

The 30-second version: measure flexibility by how hard the fit chases the noise: degrees of freedom, \text{df} = \sum_i \operatorname{Cov} (\hat{y}_i, y_i)/\sigma^2. For any linear smoother \hat{y} = Sy this is \operatorname{tr} (S): p for least squares (Ch 1's \operatorname{tr} (H) = p), \sum d^2/(d^2 + \lambda) for ridge, n/k for kNN — kNN has no parameters but k = 1 uses all 30 df. df is exactly what the training error hides: optimism = 2\sigma^2\cdot \text{df}/n, which is why Mallows' Cp and AIC add it back. Plot test error against df and the degree, \lambda and k knobs trace the same U. For classifiers, VC dimension counts the largest set of points a family can label every possible way: 3 for lines in the plane; infinite for \operatorname{sign} (\sin(\omega x)), which has one parameter. Hook: Count it, Correlate it, Shatter it. First read: the df table, Figure 4.2(c), the sine example.

Why counting fails. Three examples. kNN with k = 1 has no fitted parameters, yet it reproduces every training point exactly — as flexible as a model can be. A degree-15 polynomial with a strong ridge penalty has 16 coefficients but behaves like a much smaller model, because the penalty won't let most of them move. And \operatorname{sign} (\sin(\omega x)), one parameter, can classify any number of suitably placed points in every possible way. Counting coefficients measures storage, not flexibility.

df as covariance

The right measure asks: if I nudge y_i, how much does the fit's \hat{y}_i follow it? A fit that copies its training targets has a lot of capacity; one that ignores them has none (Efron, 2004; Card 4):

\frac{1}{\sigma^2}

divide by the noise variance \sigma^2 (0.09), so df is measured in units of noise

\sum_{i=1}^{n}

add up over all n training days

n

the number of training days (30)

\text{Cov}

covariance over noise redraws: how much two quantities move together

\hat y_i

the fit's prediction for training day i

y_i

the rides observed on day i

It rebuts "kNN has no parameters, so it can't overfit" (Q16) and explains a model card's "15 features, 7.8 df" (Q18). You can estimate it by simulation: refit on many noise draws and compute each covariance. On the season, Monte Carlo gives 3.98 for degree 3 (theory 4), 9.99 for degree 9 (theory 10), 6.02 for kNN with k = 5 (theory 6), 30.05 for k = 1 (theory 30) and 12.45 for ridge on degree 15 at \lambda = 5.62 (theory 12.44).

Linear smoothers and tr(S)

Most fits in this chapter are linear smoothers: \hat{y} = Sy for a matrix S that depends on the x's and the knob but not on y. Then \operatorname{Cov} (\hat{y}, y) = S\sigma^2, so \text{df} = \operatorname{tr} (S) exactly.

Optimism = rent, and the information criteria

The expected gap between in-sample test error (same x's, new noise) and training error is exactly 2\sigma^2\cdot \text{df}/n for any fit, linear or not (Card 4):

\mathbb{E}

average over noise redraws at the same 30 days

\frac{2\sigma^2}{n}

the rent per df: twice the noise variance \sigma^2, spread over the n days

C_p

Mallows' Cp: training error plus rent, an estimate of test error

\hat\sigma^2

the noise variance estimated from the data (hat = estimated)

For kNN with k = 5 on the season, the formula gives 2 \times 0.09 \times 6/30 = 0.036; simulation gives 0.0359 (Drill 12). Mallows' Cp adds this rent back to the training error to estimate test error. It's Ch 1's "adjusted R^2 charges rent", made exact. The information criteria charge rent on the log-likelihood scale (Gaussian errors, \sigma unknown):

\text{AIC}_c

AIC corrected for small samples: extra rent when n/\text{df} is small

BIC's rent per df is \ln 30 = 3.40 at n = 30 against AIC's 2. On this season, Cp (with \sigma^2 = 0.09 known), AICc and BIC all pick degree 3. Plain AIC picks degree 15: with \hat{\sigma}^2 estimated from the shrinking RSS, n \ln(\text{RSS}/n) keeps falling faster than 2 per df, and at n = 30 AIC's rent is too cheap. AICc's small-sample correction fixes it. Trap: AIC is an asymptotic approximation; with n/\text{df} small, use AICc or cross-validation.

Three knobs, one axis

Put degree (\text{df} = d + 1), ridge \lambda on degree 15 (df from 16 down to 1.4) and kNN k (\text{df} = 30/k) on one df axis, and the three U's line up: test error is low between about 4 and 12 df and rises on both sides. The knob is just a dial; df is the quantity that sets bias and variance. (The curves don't coincide — the families differ — but their shapes and their best region do.)

VC dimension, light

For classifiers, capacity is often measured combinatorially. A family shatters a set of points if, for every one of the 2ᵐ ways of labeling them \pm, some member of the family gets all labels right. The VC dimension is the size of the largest set it can shatter (Vapnik and Chervonenkis, 1971; Card 19).

• Lines in the plane: any 3 points not on a line can be shattered. No 4 points can: of the 16 labelings, a line realizes only 14 (the two "diagonal" labelings of four points in convex position, or "inner point vs the rest" for a triangle plus an interior point, fail). \text{VC} = 3. In d dimensions, linear classifiers have \text{VC} = d + 1 (Drill 18).
• \operatorname{sign} (\sin(\omega x)): one parameter, \text{VC} = \infty. Place points at x_i = 10^{-i}; for any labeling there is an \omega that gets every one right (checked for 5 points in sims/exact.py).

The classic bound has the shape "\text{risk} \le \text{training risk} + \text{roughly } \sqrt{\text{VC} \cdot \log n / n}". It's loose — often vacuous for real models — but it says the right thing: what matters is capacity per data point. A related measure, Rademacher complexity, asks how well the family can fit random \pm 1 signs on your actual data; Section 10's random-label experiment is this measure in action.

The trap the next act exploits. Capacity is a property of the fitting procedure, not only the family. The same degree-15 family has 16 df under least squares and 12.4 under ridge at \lambda = 5.62; early stopping and the minimum-norm rule (Sections 9–10) shrink it further without changing a single parameter. Hook: Count it, Correlate it, Shatter it — parameter count (crude), covariance df (right for smoothers), VC (worst case for classifiers).

Remember: capacity = how hard the fit chases the noise. \text{df} = \sum \operatorname{Cov} (\hat{y}_i, y_i)/\sigma^2 = \operatorname{tr} (S): p for OLS, \sum d^2/(d^2 + \lambda) for ridge, n/k for kNN. Optimism = 2\sigma^2\cdot \text{df}/n, so \text{Cp} = \text{train} + \text{rent}. Plain AIC under-charges at small n. Count it, Correlate it, Shatter it.

Next question → we can see overfitting and measure capacity. What's the first cure — and how do we know whether to spend on more data or a better model? (Section 5)

5. More Data: Reading Learning Curves

The question: a vendor offers ten times more days of ride data for $200k. Will it help the degree-1 curve? The degree-12 curve? How do we decide before paying?

The 30-second version: a learning curve plots training and validation error against the number of training rows. With more data, variance shrinks (roughly \sigma^2\cdot \text{df}/n) but bias doesn't move. So two curves that meet high above the noise floor = high bias: more data won't help; get a more flexible model or better features (degree 1 stalls at 0.286 forever). A wide gap that is still closing = high variance: more data helps (degree 12 goes from 0.146 at n = 30 to 0.091 at n = 1{,}000, overtaking the cubic from n = 200). Hook: more data cures variance, never bias. Learning curves often follow a power law, \text{error} \approx a + b\cdot n^{-\alpha}, which lets you price data before buying it — the same curve at enormous scale is an LLM scaling law (Section 14). First read: Figure 4.3, the decision table.

What a learning curve is

Pick a model. For a grid of training sizes n, fit on n rows and record two numbers: the training error and the error on held-out data (a validation set, or cross-validation at each n). Plot both against n, usually on a log axis. Training error rises with n (more points are harder to fit exactly); validation error falls (more points pin the fit down). They converge to the same level: the model's best achievable risk, \sigma^2 + \text{bias}^2.

Here we simulate it (sims/sim_learning_curves.py): n from 8 to 1,000 days spread evenly over the year, 200 replicates at each n, and the expected test error computed exactly over new noise (\sigma^2 plus the mean of (f - \hat{f})^2 over 10,000 random days).

The three shapes

• Degree 1 (high bias). The curves meet by n \approx 100 at 0.286, far above the noise floor of 0.09. That level is exact: the best possible line through a full sine wave leaves \text{bias}^2 = 0.5 - 3/\pi^2 = 0.196, plus \sigma^2 = 0.09. Ten times more data buys nothing.
• Degree 3 (about right). Small gap at n = 30, converging to 0.094 — just above the floor, because a cubic leaves a little bias.
• Degree 12 (high variance). A large gap at n = 30 that keeps closing; by n = 200 degree 12 matches the cubic and then beats it, converging to the floor 0.09 itself, because a degree-12 polynomial can follow the sine almost perfectly. Notice also n = 15: test error spikes to 16.3 when the 13-coefficient fit has barely more points than coefficients. That's double descent showing up inside a learning curve (Section 10).

The decision table: more data or more model?

So for the vendor: if the current model is the line, decline the data and buy a better feature (a seasonal term). If the current model is degree 12, the data would help — but a degree-3 model on today's data is already at 0.107. Q6 and Q36 play this out.

How much data?

For a least-squares smoother with a fixed design, the average variance is \sigma^2\cdot \text{df}/n (Card 18). So the excess error over the model's floor shrinks like 1/n, and you can solve for n. Degree 12 has \text{df} = 13: to get its variance down to 0.01 you need n \approx 0.09 \times 13/0.01 = 117 days; to 0.001, about 1,170. The simulation agrees: excess × n stays between 1.14 and 1.27 from n = 50 to 1,000, close to \sigma^2\cdot \text{df} = 0.09 \times 13.

Power-law extrapolation

Real learning curves rarely come with a known df, but they often follow

n

the number of training days (rows)

a

the asymptote: noise floor plus \text{bias}^2, where the curve flattens out

b

the scale of the extra error a small dataset carries

\alpha

how fast extra data pays off; about 1 for the 1/n variance rate

where a is the asymptote (noise floor plus \text{bias}^2), and \alpha says how fast extra data pays (Hestness et al., 2017). Fitting degree 12's test curve on n \ge 50 gives a = 0.0903, b = 1.82, \alpha = 1.09 — the asymptote is the noise floor, and \alpha \approx 1 is the 1/n variance rate. With three runs at n, 4n and 16n you can fit it by hand: the successive differences shrink by 4^\alpha, so \alpha = \log_4 of their ratio (Drill 16, Q39). Fit the curve before you buy data, then ask what the predicted error at 10n is worth.

Two traps

1. A learning curve on a random split can't show you shift. If next year's days, or the vendor's days, come from a different distribution, the curve built on this year's split says nothing about them. More data from the wrong distribution doesn't help (Section 11).
2. Near the interpolation threshold, more data can briefly hurt. Degree 12's spike at n = 15 is the sample-wise version of double descent (Section 10). Don't read one point of a learning curve in isolation.

Remember: more data cures variance, never bias. Curves meet high → change the model; curves gapped → get data. Excess error \approx \sigma^2\cdot \text{df}/n. Fit a + b\cdot n^{-\alpha} before you buy.

Next question → data is expensive. Can we cut variance without it, by restraining the model? (Section 6)

6. Ridge in Depth: Shrinkage, Effective Degrees of Freedom and the Bayesian View

The question: the degree-15 curve has test error 0.379. How exactly does the ridge penalty tame it, how much bias does it cost, and what \lambda is best?

The 30-second version: recap Ch 1: ridge minimizes \lVert y - Xw\rVert^2 + \lambda \lVert w\rVert^2, with solution (X^\top X + \lambda I)^{-1}X^\top y. Rotate into the SVD basis and ridge is simple: along each direction of the data it multiplies the least-squares answer by d^2/(d^2 + \lambda), where d is that direction's singular value. Directions the data pins down firmly (big d) are barely touched; weakly determined directions — where all the variance lives — are shrunk hardest. Hook: shrink the weak directions most. Effective df = \sum d^2/(d^2 + \lambda) slides smoothly from p to 0. Bias grows with \lambda, variance falls, and some \lambda > 0 always beats least squares (Hoerl and Kennard, 1970). Bayesian view: with a N(0, \tau^2) prior on each weight, the posterior mean is ridge with \lambda = \sigma^2/\tau^2 — noise over prior — and that \lambda is the best possible on average. On the season, ridge brings degree 15 from 0.379 to 0.137 (df 12.4), but the cubic still wins at 0.106: a penalty is not a substitute for the right features. On 100 stations × 80 features, ridge cuts test MSE from 22.2 (worse than predicting the mean, 14.5) to 5.05. First read: the one-weight example, Figure 4.4(a–b), the Bayesian line.

One weight, by hand

On the toy, center the months (m - 6) and keep the intercept unpenalized (it stays at \bar{y} = 0). Ridge on the slope minimizes \sum(y_i - w(m_i - 6))^2 + \lambda w^2, so

\hat w_\lambda

ridge's slope for the toy line at penalty \lambda, in thousand rides per month

S_{my}

how months and rides move together: \sum (m_i - 6)\,y_i = -13.6

S_{mm}

the spread of the months: \sum (m_i - 6)^2 = 70

\lambda

the ridge penalty strength; 0 gives least squares

\lambda = 70 = S_{\text{mm}} halves the slope. The shrink factor 70/(70 + \lambda) is exactly d^2/(d^2 + \lambda) with d^2 = S_{\text{mm}} = 70, the one squared singular value of the centered month column. And look at July 20: shrinking the line makes its prediction worse (further from 0.866), because the line's problem is bias, not variance. Ridge cures variance.

How much shrinkage is best for the slope itself? Its MSE as an estimate of the best-line slope \beta = -12/70 is (\lambda \beta/(70 + \lambda))^2 + 0.09\cdot 70/(70 + \lambda)^2, minimized at \lambda^* = \sigma^2/\beta^2 = 3.06, where the MSE falls only from 0.001286 to 0.001232 — 4%. A strong, well-measured signal needs little shrinkage (Drill 4).

The SVD picture

With X = UDV^\top (Ch 1, Section 3), rotate the weights into the directions V. In those directions X^\top X is diagonal, so ridge decouples into one-weight problems like the one above (Card 6):

\hat{\mathbf w}_\lambda

ridge's vector of weights at penalty \lambda (bold = vector of all of them)

X

the design matrix: one row per training day, one column per feature

(X^\top X + \lambda I)^{-1}

the inverse (it undoes the matrix) of the penalized cross-product matrix

X^\top X

the features' cross-products; diagonal in the SVD basis

X^\top

the transpose of X: rows become columns

I

the identity matrix: 1s on the diagonal, 0s elsewhere

\lambda

the ridge penalty strength

\mathbf y

the vector of all training days' rides

\sum_{j}

add up over the data's directions j, one per singular value

\mathbf v_j

direction j in weight space (column j of V)

\mathbf u_j^\top\mathbf y

how much of the rides vector points along direction j

\mathbf u_j

direction j across the training days (column j of U)

d_j

singular value j: how much the data spread along direction j

Least squares is the same sum with every factor equal to 1. Why shrink the weak directions most? Along direction j, the least-squares coefficient has variance \sigma^2/d_j^2: a direction with little spread in the data is barely measured, so its coefficient is mostly noise. Ridge's factor d_j^2/(d_j^2 + \lambda) is near 1 when d_j^2 \gg \lambda and near 0 when d_j^2 \ll \lambda — it keeps what the data measured and discards what it didn't (Q17).

The season's degree-15 design (15 standardized Legendre columns, intercept separate) has 13 strong directions with d^2 between 26.7 and 47.6, and two weak ones with d^2 = 1.47 and 0.59. At the best \lambda = 5.62, the strong directions keep 83–89% of their least-squares value; the two weak ones keep 21% and 10%. Those two directions are where degree 15's wiggles live.

Effective degrees of freedom

Ridge is a linear smoother, so its df is the trace of its smoother matrix (Section 4):

\text{tr}

trace: the sum of a matrix's diagonal entries

It slides from p at \lambda = 0 to 0 as \lambda \to \infty; add 1 for the unpenalized intercept. On the season it runs from 16 down to 1.4, and at the best \lambda it is 12.4 — the two weak directions are nearly switched off and the strong ones are trimmed a little. For Drill 5's design (d^2 = 100, 10, 1, 0.1, \lambda = 1): 0.990 + 0.909 + 0.5 + 0.091 = 2.49.

Bias and variance in closed form

Substitute y = X\beta + \varepsilon (Card 7):

\mathbb{E}

average over fresh noise on the same design

\boldsymbol\beta

the true weight vector the data were generated from

\text{Var}

covariance matrix of the weights across noise redraws

\sigma^2

noise variance

Ridge is biased toward zero for every \lambda > 0. Direction by direction, with \theta_j = v_j^\top \beta the true coefficient in that direction:

\text{MSE}_j(\lambda)

the squared error of ridge's coefficient along direction j, at penalty \lambda

\theta_j

the true coefficient along direction j, \mathbf v_j^\top\boldsymbol\beta

\lambda_j^*

the best penalty for direction j alone (star = the optimum)

\text{Bias}^2 rises from 0, variance falls from \sigma^2/d_j^2: that is the "regularization slider" (Q8). Each direction's best \lambda is noise over \text{signal}^2, which is Drill 3 (\sigma^2 = 0.25, \theta = 1, d^2 = 1: \lambda^* = 0.25, MSE 0.20 vs OLS 0.25).

On the season (500 seasons, degree 15):

Some λ > 0 always beats least squares

Sum the per-direction MSE and differentiate at \lambda = 0. The \text{bias}^2 term starts with zero slope (it's quadratic in \lambda) while the variance term starts falling:

\frac{d}{d\lambda}

how fast the total coefficient error changes as you nudge \lambda: a slope, not a fraction to divide

\Big|_{\lambda = 0}

evaluated at no penalty, i.e. at least squares

So a small enough penalty always lowers the total MSE of the coefficients, whatever \beta is (Hoerl and Kennard, 1970; Card 7). On the toy slope (d^2 = 70) the slope is -2 \times 0.09/70^2 = -3.67 \times 10^{-5}, matching a numerical derivative. This is the regression cousin of James–Stein (Ch 1, Drill 33): unbiasedness is not optimal.

The Bayesian view: λ = noise over prior

Put a prior on the weights: each w_j \sim N(0, \tau^2) independently, and y \mid w \sim N(Xw, \sigma^2I). The log-posterior is, up to a constant,

\lVert \mathbf y - X\mathbf w\rVert^2

squared length of the miss vector: the sum of squared misses

\mathbf w

the weight vector, now treated as random with a prior

\tau^2

prior variance: how big we expect each weight to be before seeing data

\lVert \mathbf w\rVert^2

squared length of the weight vector: the sum of squared weights

\mathbb{E}[\mathbf w \mid \mathbf y]

the posterior mean: the average weights given the observed rides

The posterior is Gaussian, so its mode (MAP, Ch 1's statement) and its mean coincide: ridge with \lambda = \sigma^2/\tau^2. Read it as noise over prior: loud noise or a tight prior → shrink more. And because the posterior mean minimizes expected squared error averaged over the prior, no estimator — linear or not — has lower average risk when the prior is right (Card 8). Ch 6 formalizes MAP and Bayes estimators.

The dense station example checks it (Worked Problem 1 in the Interview Bank): 80 descriptors, each with a small true effect drawn from N(0, 0.05), noise \sigma^2 = 4, 500 random draws of the whole problem. Average test risk is minimized at exactly \lambda = 80 = 4/0.05, where it is 5.20; least squares averages 20.7, \lambda ten times too big 6.21 and ten times too small 6.30. The curve is flat near the optimum: getting \lambda within a factor of 3 is most of the win.

Practical notes, one line each

• Standardize the features first and leave the intercept unpenalized (Ch 1, Section 10); sklearn's Ridge(alpha) uses this chapter's sum-of-squares scaling, so \text{alpha} = \lambda.
• X^\top X + \lambda I is always invertible for \lambda > 0, so p > n is fine.
• Compute via SVD (one decomposition gives the whole \lambda path) or the augmented-data trick: append \sqrt{\lambda}\cdot I as extra rows of X and zeros to y, then run ordinary least squares — ridge is least squares on fake data that says "every weight is zero" (on the season the two agree exactly).
• Principal component regression (PCR) keeps the top k directions and drops the rest, a factor of 1 or 0; ridge is its soft version.
• The dual form \hat{w} = X^\top(XX^\top + \lambda I)^{-1}y needs only inner products between rows, which is the door to kernels (Ch 11).
• Shrinkage in a classifier pulls predicted probabilities toward the base rate (Ch 3, Section 9, "Shrink"), so a ridge-penalized logistic model needs a calibration check.

Ridge on the 80 station features

The second running dataset: the operator is opening in a new city and models each existing station's average daily rides (standardized) from 80 candidate descriptors — population within 250 m, 500 m and 1 km, transit stops, offices, elevation, bike lanes, and so on — measured on n = 100 stations (sims/sim_wide.py). The descriptors come in 16 correlated blocks of 5 (within-block correlation 0.8). Only the first two blocks matter (10 features, each with true weight 0.5), and the noise has \sigma = 2, so the irreducible error is 4 and predicting the mean scores 14.5.

Least squares is worse than predicting the mean. With random new stations and an intercept, its exact expected test error is \sigma^2(1 + 1/n)(1 + p/(n - p - 2)) = 22.0 (checked at 21.96 over 3,000 simulated designs); at the same 100 stations it would be only \sigma^2(1 + p/n) = 7.24. Eighty coefficients from 100 rows leaves tiny singular values, and new stations land in exactly those weak directions. Ridge (about 27 effective df) spreads weight across each correlated block instead of betting on single members (Q30).

Remember: ridge = least squares with each direction multiplied by d^2/(d^2 + \lambda). Shrink the weak directions most. \text{df} = \sum d^2/(d^2 + \lambda). \text{Bias}^2 rises from 0, variance falls, and some \lambda > 0 always wins. \lambda = noise over prior, \sigma^2/\tau^2. A penalty doesn't fix the wrong family (0.137 vs the cubic's 0.106).

Next question → ridge never sets a weight to zero. When we want a short list of features, what does the L1 penalty do in depth — and what does it cost? (Section 7)

7. Lasso and Elastic Net in Depth: Sparsity, Paths and Instability

The question: the expansion team wants the handful of station features that drive demand. Lasso gives a short list — but a different list every time we resample. Why, and what do we do?

The 30-second version: recap Ch 1: L1 kills, L2 shrinks. In depth: lasso is the MAP estimate under a Laplace prior (a sharp peak at zero). Its optimality condition says a weight stays exactly zero unless its feature's correlation with the current residual exceeds \lambda/2, so above \lambda_{\max} = 2\cdot \max_j\lvert x_j^\top y\rvert every weight is zero, and as \lambda falls features enter one at a time along a piecewise-linear path. With orthonormal features lasso soft-thresholds: subtract a fixed amount and stop at zero. Two costs: survivors are shrunk by the full threshold (biased even when large), and among correlated features lasso picks one almost arbitrarily — 100 bootstrap refits gave 100 different lists. Ridge spreads weight across the group; elastic net groups them. Hook: lasso picks, ridge spreads, elastic net groups. On the stations lasso keeps 14.6 features with 6.2 false positives; elastic net does slightly better on test error (4.67 vs 4.78, better in 82% of replicates). And never report naive p-values for the survivors: under a pure-noise truth, lasso-then-OLS finds a "significant driver" in 39% of runs. First read: the path figure, picks/spreads/groups, the post-selection warning.

The Laplace prior and MAP

Put an independent Laplace prior on each weight, p(w_j) \propto \exp(-\lvert w_j\rvert/b). Its log is -\lvert w_j\rvert/b, so the MAP estimate under Gaussian noise minimizes \lVert y - Xw\rVert^2 + \lambda \lVert w\rVert_1 with \lambda = 2\sigma^2/b (Card 9). The Laplace density has a sharp peak at zero, and that kink is what lets the mode sit exactly at zero.

Trap: the Bayesian lasso's posterior mean is never exactly zero. Sparsity comes from taking the mode of a kinked posterior, not from the prior "believing" weights are zero.

Optimality conditions, λ_max and soft-thresholding

\lvert w\rvert has no derivative at 0, so the optimality condition uses a subgradient: any value in [-1, 1] at w_j = 0. Setting the subgradient of \lVert y - Xw\rVert^2 + \lambda \lVert w\rVert_1 to contain zero gives the KKT conditions (Card 9):

\mathbf x_j^\top(\mathbf y - X\hat{\mathbf w})

how strongly station feature j lines up with the current residual (actual minus fitted rides)

\lambda

the lasso penalty strength

\text{sign}

the sign: +1 if positive, -1 if negative

\hat w_j

lasso's weight on station feature j

In words: an active feature's correlation with the current residual is pinned at exactly \pm \lambda/2; an inactive feature's correlation is below \lambda/2 (Q19). Start at w = 0, where the residual is y: every weight stays at zero as long as \lambda \ge \lambda_{\max} = 2\cdot \max_j\lvert x_j^\top y\rvert. On station replicate 0, \lambda_{\max} = 512 (sklearn alpha 2.56), and at that \lambda the fit is all zeros.

With orthonormal features (X^\top X = I), the problem splits by coordinate. With z_j = x_j^\top y the least-squares coefficient:

z_j

the least-squares coefficient for feature j (with orthonormal features, \mathbf x_j^\top\mathbf y)

\big(|z_j| - \tfrac{\lambda}{2}\big)_+

shrink \lvert z_j\rvert by \lambda/2 but stop at 0 (the _+ keeps only positive values)

That is soft-thresholding: move every coefficient \lambda/2 toward zero and stop at zero. Drill 10: z = (3, 0.8, -1.5) with \lambda = 2 gives (2, 0, -0.5), and \lambda_{\max} = 6. Compare ridge, which divides by (1 + \lambda) and never reaches zero, and best-subset selection, which hard-thresholds (keep or kill, no shrink). For general X, coordinate descent applies soft-thresholding one coordinate at a time, and proximal gradient descent applies it after each gradient step (Ch 2, Section 12).

The regularization path

As \lambda falls from \lambda_{\max}, features enter one at a time; between entry (and exit) events every coefficient moves linearly in \lambda. LARS (least-angle regression) computes the whole piecewise-linear path for about the cost of one least-squares fit. When p > n, lasso can have at most n non-zeros.

On station replicate 0 the first features to enter are 10, 9, 2, 4, 5, 8, 7, 3, 6 — nine of the ten true features — then noise features 71, 73 and 19. True feature 1 never enters at the CV penalty: its block-mates already carry its signal. Lasso's degrees of freedom are the number of non-zeros (Zou, Hastie and Tibshirani, 2007), so its rent in Cp is easy (Section 4).

Bias of the survivors. Soft-thresholding shrinks every survivor by \lambda/2, even a huge, obviously real one. Lasso is therefore biased toward zero on exactly the coefficients you most want to read. Two fixes, one line each: the relaxed lasso refits the selected set with less or no penalty; the adaptive lasso reweights the penalty by 1/\lvert\text{initial estimate}\rvert, shrinking big coefficients less.

Instability under correlation

Inside a block of five features with correlation 0.8, all five carry nearly the same information. Lasso's penalty prefers one large weight to five small ones of the same total, so it picks one or two and zeroes the rest — and which ones depends on the noise. On replicate 0, 100 bootstrap refits produced 100 different selected sets. The selection frequency of true feature 1 was 11% under lasso and 53% under elastic net; averaged over the ten true features, 76% vs 95%; over the 70 noise features, 9% vs 14%.

Stability selection (Meinshausen and Bühlmann, 2010) turns this into a method: refit on many half-samples and keep features selected in more than a threshold fraction of them. Report selection frequencies, not one list (Q20).

Elastic net and the grouping effect

Elastic net adds both penalties (Zou and Hastie, 2005; Card 10):

\mathbf w

the vector of weights on the 80 station descriptors

\lVert \mathbf y - X\mathbf w\rVert^2

the sum of squared misses over the 100 stations

\lambda_1

strength of the L1 (lasso) penalty

\lVert\mathbf w\rVert_1

L1 norm: the sum of the absolute weights

\lambda_2

strength of the L2 (ridge) penalty

\lVert\mathbf w\rVert_2^2

squared L2 length: the sum of the squared weights

The L2 part makes the objective strictly convex, so the solution is unique, and it rewards splitting weight evenly across correlated features. Two identical columns get identical weights. In Drill 11 a feature is duplicated (\lVert x\rVert^2 = 100, OLS on one copy = 1, each \lambda = 20): ridge gives (0.455, 0.455), elastic net (0.409, 0.409), while lasso's total is 0.9 split in any way with the same sign — sklearn happens to return (0.9, 0). For two standardized features with correlation \rho, the elastic-net weights differ by at most (\lVert y\rVert_1/\lambda_2)\cdot \sqrt{2(1 - \rho)}, so highly correlated features move together. In sklearn the mix is l1_ratio (1 = lasso, 0 = ridge).

So the three penalties treat a correlated block differently. Lasso picks, ridge spreads, elastic net groups.

Which wins when?

Elastic net beat lasso in 82% of sparse replicates. Lasso's 1-SE choice (Section 8) keeps 8.6 features at a small cost in error (5.19). The pattern to say out loud: lasso and elastic net win when few features matter; ridge wins when many small effects add up; elastic net is the safe default when features are correlated (Q7).

Post-selection inference

After lasso picks features, it's tempting to refit least squares on the survivors and report their p-values. Those p-values are invalid: the same data chose the features because they looked strong. A pure-noise check (sims/sim_postsel.py, 400 replicates; true \beta = 0, so every "driver" is false): lasso-CV selected at least one feature in 44% of runs, and in 39% of all runs at least one survivor had a naive p < 0.05 — against the 5% a valid test allows. Among selected features, 41% looked significant. Sample splitting fixes it: select on 50 stations, test on the other 50; the share of runs with a "significant" feature drops to 7.5% (2.0% with a Bonferroni correction for the number tested). Other fixes: selective inference, which conditions on the selection event (Lee, Sun, Sun and Taylor, 2016); simultaneous intervals over all possible selections (Berk et al., 2013); stability selection. Ch 1, Q66 raised the warning; this is its size (Q31).

What to tell the expansion team

"Population density near the station and transit access matter; we can't say which radius, because the radii carry the same information. Here are the groups that matter, and how often each feature is selected when we resample. For prediction we'll use elastic net. If you need effect sizes with p-values, we'll estimate them on held-out stations."

Remember: Laplace prior → sparse mode. KKT: zero unless \lvert x_j^\top \text{residual}\rvert reaches \lambda/2, so \lambda_{\max} = 2\cdot \max \lvert x_j^\top y\rvert. Orthonormal → soft-threshold by \lambda/2. Lasso picks, ridge spreads, elastic net groups. Survivors are shrunk too. Never trust naive p-values after selection.

Next question → every cure has a knob — degree, k, \lambda. How do we set it using only the training data, and how far can we trust the number we get? (Section 8)

8. Choosing λ: Cross-Validation Done Right

The question: we have only this season's 30 days. How do we pick the degree (or \lambda) without a test set, and is the CV error we report an honest estimate of next year's error?

The 30-second version: K-fold cross-validation: split the rows into K folds; for each fold, fit on the other K − 1 and score on it; average. Pick the knob with the lowest average — or, to prefer simplicity, the simplest model within one standard error of the best (1-SE rule). The CV number has three flaws, the 3 Ls: Less data — each fit sees only (K − 1)/K of the rows, so CV is a little pessimistic (on our six days, 3-fold CV prefers the line, 0.300 vs the cubic's 0.559, though the cubic is far better on next year's days); Lots of noise — folds overlap, so the estimate is noisy and its naive standard error too small; Lucky minimum — the minimum over many settings is optimistic, more so the more settings you try. Report an honest number with nested CV: tune inside, test outside. Everything that learns from data — scaling, feature selection, imputation — goes inside the fold. First read: the toy 3-fold table, Figure 4.6, the 3 Ls, tune inside / test outside.

K-fold by hand

Split the six toy days into three folds, each holding one warm and one cold day: {A, D}, {B, E}, {C, F}. Fit on four days, score on the other two, rotate:

CV picks the line, which would staff July 20 for 0.389 when the truth is 0.866. Why? A cubic fitted on four points interpolates them exactly and swings wildly at the held-out day (fold {A, D}: 0.981). That's the first L, Less data, at its most extreme: CV scores each model at n = 4, and at n = 4 the line is better (Drill 15 sizes this pessimism for a known learning curve). With six days, no validation method can recover the cubic.

LOOCV and its shortcut

Leave-one-out CV (LOOCV) is K = n. For any linear smoother at a fixed knob, it needs no refitting (Ch 1, Q69; Card 12):

y_i

the rides on day i

\hat y_{(-i)}

the prediction for day i from a fit that left day i out

\hat y_i

the prediction for day i from the fit on all days

S_{ii}

day i's leverage: the diagonal entry of S, its pull on its own fit

n

the number of days

\sum_{i=1}^{n}

add up over every day i

\sum_i

add up over every day i

\frac1n

divide by n: with the sum, the average over the days

S

the smoother matrix that maps observed rides to fitted rides, \hat{\mathbf y} = S\mathbf y

\text{tr}

trace: the sum of the diagonal, here the fit's df

For the toy line, the leverages S_{ii} are 0.524, 0.295, 0.181 (symmetric): the edge days A and F have the most pull, so their LOO residuals are inflated most. LOOCV by degree is 0.787, 0.258, 0.977, 0.674 for degrees 0–3 — again the line. For the mean (degree 0), every S_{ii} = 1/n, so \text{LOOCV} = (n/(n - 1))^2 \times \text{train MSE} = 1.44 \times 0.547 = 0.787 (Drill 9). The same shortcut makes ridge's LOOCV free at every \lambda (it matches brute-force refitting to 3 \times 10^{-11} on the season).

GCV (Craven and Wahba, 1979) replaces each S_{ii} by the average \operatorname{tr} (S)/n. That's cheap and rotation-invariant, but it fails when leverages are very unequal. On this season GCV picks \lambda = 0.32 for the degree-15 ridge, with test MSE 0.52, while LOOCV and 10-fold both pick 5.62 (test 0.137): the edge days of a degree-15 design have far more leverage than average, and GCV can't see them.

Choosing K

At n = 30 million rows, a single holdout of a few hundred thousand rows already has a tiny SE; K-fold is for when data is scarce (Q24).

The season CV curves

On this season (replicate 0), 10-fold CV at degree 3 is 0.120 with SE 0.041; 5-fold gives 0.117, LOOCV 0.120. All three put the minimum at degree 3, and the 1-SE rule also picks 3. (The true test error on this season is lowest at degree 5, 0.110 vs 0.114 for degree 3 — CV can't resolve a difference that small.) Reshuffling the 10 folds 20 times picks degree 3 eighteen times and degree 5 twice. A single 80/20 split is much worse: over 20 random seeds it picked 10 different degrees, from 3 to 15 (Q9).

The 1-SE rule

Compute the SE of the CV estimate at each setting as the standard deviation of the K fold errors divided by \sqrt{K}. Find the minimum. Then choose the simplest model whose CV error is within one SE of the minimum (Hastie, Tibshirani and Friedman). "Simplest" must be defined: lowest df, lowest degree, or largest \lambda. The logic: settings within one SE of the best are statistically indistinguishable, and the simpler one has lower variance and is easier to defend. Over 500 seasons, plain 10-fold picks degree 3 in 267 seasons and the 1-SE rule in 436, lowering the mean test error of the chosen model from 0.112 to 0.109 (the oracle, always degree 3, gets 0.106). For the degree-15 ridge on this season, 1-SE picks \lambda = 10 (df 11.2, test 0.130) over the minimum's 5.62 (Q21).

The 3 Ls in numbers

Over 500 seasons (sims/sim_cv.py):

Two findings deserve a sentence each. First, the naive SE is too small because the K fold errors are correlated — their training sets overlap — so they are not K independent measurements (Bengio and Grandvalet, 2004; Drill 17: correlation 0.3 between 10 folds understates the SE by 1.92× if the fold SD is known, 2.30× if estimated from the folds). Second, a season's CV estimate barely tracks that season's true test error: the correlation at degree 3 is 0.03. CV estimates the average error of the procedure over datasets better than the error of the model in hand (Bates, Hastie and Tibshirani, 2023).

The Less-data and Lucky-minimum Ls pull in opposite directions, and which wins depends on how many settings you search (Drill 14: the best of ten equal models, each with validation noise SD 0.01, reports 0.185 against a true 0.200). On the stations (sims/sim_nested.py), a 19-value \lambda search reports best_score_ 0.09 ± 0.06 above the true error (Less data wins, 300 replicates); a 760-configuration search over \lambda and feature subsets reports it 0.83 ± 0.15 below (Lucky minimum wins, 40 replicates).

Nested CV: tune inside, test outside

To get an honest number for a procedure that includes tuning, put the whole tuning loop inside an outer CV:

1. Split the rows into outer folds.
2. For each outer fold, run the full tuning procedure (inner K-fold over every setting) on the other folds only, refit the winner, and score it on the held-out outer fold.
3. Average the outer scores. That is the estimate for "tune by CV, then fit" — the procedure you will ship.

The outer fold never influenced any choice, so its score is free of the Lucky minimum. In the two searches above, nested CV is +0.19 ± 0.06 and +0.25 ± 0.19 above the truth: honest to slightly pessimistic both times, because each outer fit sees less data (Varma and Simon, 2006; Cawley and Talbot, 2010). Hook: tune inside, test outside (Q10, Q23). With hundreds of configurations scored on one validation set, the winner's curse is large (Q35).

Leaks inside CV

CV is only honest if every step that learns from the data happens inside the training folds. Fitting a feature selector, a target encoder, an imputer or a scaler on all rows and then cross-validating lets the held-out fold influence the model. The damage depends on how supervised the step is. Supervised screening is catastrophic: choosing genes by their correlation with the label on all samples made random labels look predictive (Ambroise and McLachlan, 2002), and Section 12 reproduces it. An unsupervised scaler barely matters (and for least squares with an intercept, not at all). The fix is mechanical: a pipeline (Section 13; Q22).

Analytic alternatives. When n is tiny or refits are expensive, Cp, AICc, BIC and GCV (Section 4) estimate test error from one fit. They assume a correct df and, for GCV, similar leverages — the season showed both failing (plain AIC chose degree 15, GCV chose \lambda = 0.32). CV assumes only that the rows are exchangeable. For grouped or time-ordered rows, ordinary K-fold breaks that assumption; Ch 5 owns grouped and time-series validation design.

Remember: fit on K − 1 folds, score on the last, rotate, average. 1-SE: the simplest model within one standard error of the best. The 3 Ls: Less data (pessimistic), Lots of noise (naive SE too small), Lucky minimum (optimistic, grows with the number of settings). Tune inside, test outside. Every learned step goes inside the fold.

Next question → penalties are one way to restrain a model. What else works — and why do so many of them turn out to be ridge in disguise? (Section 9)

9. Other Regularizers: Stopping Early, Adding Noise, Dropping Out, Averaging

The question: an engineer stops gradient descent early, another adds noise to the inputs, a third averages 50 models. Are these hacks, or the same idea?

The 30-second version: every regularizer cuts variance by restraining what the fit can do with the data. Hook: PANDAS — Penalty, Augmentation, Noise injection, Dropout, Averaging, Stopping early. Three are ridge in disguise for a linear model: adding Gaussian noise of variance s^2 to every input is ridge with \lambda = n\cdot s^2 (Bishop, 1995); dropout is ridge with a per-feature penalty (Wager, Wang and Liang, 2013); and stopping gradient descent from zero after t steps multiplies each SVD direction by 1 - (1 - \eta d^2)^t, a filter shaped like ridge's d^2/(d^2 + \lambda) with \lambda roughly 1/(\eta t). Averaging B models with correlation \rho cuts variance to \rho \sigma^2 + (1 - \rho)\sigma^2/B — only the uncorrelated part averages away. First read: PANDAS, Figure 4.4(c), the three equivalences in words.

PANDAS

P was Sections 6–7. One reminder: for plain SGD, adding \lambda \lVert w\rVert^2 to the loss and "weight decay" are the same; for Adam they are not, which is why AdamW decouples the decay from the adaptive step (Loshchilov and Hutter, 2019; Ch 2, Section 8; Q37).

S: early stopping as spectral filtering

Run gradient descent on \tfrac 12 \lVert y - Xw\rVert^2 from w_0 = 0 with step \eta. In the SVD basis each direction evolves on its own, and after t steps the fit is (Card 13):

X\mathbf w_t

the fitted values after t gradient-descent steps from \mathbf w_0 = 0

t

the number of gradient-descent steps taken

\sum_j

add up over the data's directions j

\mathbf u_j^\top\mathbf y

how much of the rides vector lies along direction j

\mathbf u_j

direction j across the training days

\eta

the step size (learning rate)

d_j

singular value j: the data's spread along direction j

\lambda

the ridge penalty strength

Strong directions (large \eta d^2) are fitted in a few steps; weak directions creep up slowly and are still near zero when you stop. That's "shrink the weak directions most" again, done by the clock instead of by a penalty. Both filters run from 0 to 1 as d^2 grows, and they match at the half-way point when \lambda = \ln 2/(\eta t), which is why \lambda \approx 1/(\eta t) is the usual rough map (Ali, Kolter and Tibshirani, 2019). Drill 7 (\eta = 0.01, t = 100, d^2 = 10, 1, 0.1): gradient descent keeps 1.000, 0.634, 0.095; ridge with \lambda = 1 keeps 0.909, 0.5, 0.091.

On the season's degree-15 design with \eta = 1/\lambda_{\max} = 0.021, test error falls from 0.176 after one step to its minimum 0.130 at t = 3, then climbs back to least squares' 0.379 as the weak directions fill in. Ridge's minimum is 0.137. The two curves have the same U and nearly the same floor, but the 1/(\eta t) map is rough: at matched t they differ by up to 0.10. With a ten-times smaller step the gap shrinks to 0.077, and the best constant is \lambda \approx 0.6/(\eta t) (gap 0.047). Treat "early stopping ≈ ridge" as a statement about shape, not a formula.

In practice you don't know the best t: hold out a validation split, track its loss, stop after it has failed to improve for a patience of several checks, and restore the best weights. On a 20/10 split of this season, validation loss bottomed at step 4, patience stopped the run at step 42, and the restored weights scored 0.159 on new days — against 2.07 had it run to convergence on 20 days (Q11).

N: noise injection is ridge

Replace each training input x_i by \tilde{x}_i = x_i + \delta_i with \delta_i \sim N(0, s^2I), fresh each time. In expectation over the noise (Card 14):

\mathbb{E}

average over the injected input noise

\mathbf y

the vector of training rides

\tilde X

the inputs with fresh Gaussian noise added (tilde = jittered copy)

X

the clean inputs (the design matrix)

\mathbf w

the weight vector

n

the number of training rows

s^2

the variance of the injected input noise

\lVert\mathbf w\rVert^2

squared length of the weights: the sum of squared weights

\lVert \mathbf y - X\mathbf w\rVert^2

the ordinary sum of squared misses on clean inputs

The cross term dies because \delta has mean zero; the \delta^2 term leaves n\cdot s^2\lVert w\rVert^2. So training on jittered inputs is ridge with \lambda = n\cdot s^2 for the sum loss (\lambda = s^2 for the mean loss — say which). On the season's standardized degree-15 design with s = 0.1 and 200 noisy copies of each row, least squares on the copies and ridge with \lambda = 30 \times 0.01 = 0.3 give coefficients within 0.003 of each other (the largest is 0.68), and test errors of 0.534 and 0.533 (Drill 6, Q25). With finitely many copies the match is approximate, since the copies' cross terms don't cancel exactly.

D: dropout is per-feature ridge

Inverted dropout multiplies each input by \delta/q with \delta \sim \text{Bernoulli}(q): kept with probability q and scaled up so the expected input is unchanged. In expectation (Card 15):

\mathbb{E}

average over the random dropout masks

\tilde X

the inputs after dropout: each kept with probability q and scaled by 1/q

q

keep probability: the chance each input survives

\sum_j

add up over features j

\sum_i

add up over training rows i

x_{ij}

feature j's value on row i

w_j

the weight on feature j

Each weight gets its own ridge penalty, proportional to its feature's energy \sum_ix_{ij}^2. On the season's design with q = 0.8, the Monte Carlo dropout objective (8.769 ± 0.006) matches the formula (8.758). With standardized features on 100 rows, \sum_ix_{ij}^2 = 100 and the equivalent ridge \lambda is (0.2/0.8) \times 100 = 25 (Drill 19; Q26). Dropout must be switched off at prediction time (model.eval() in PyTorch); Ch 17 owns dropout in networks.

A: augmentation. Augmentation adds transformed copies of training rows whose label you are sure doesn't change: flips and crops of images, small time shifts of audio, synonym swaps in text. It encodes an invariance — "a flipped cat is still a cat". For ride demand, the defensible augmentations are few: jittering a date by a day, perhaps. Flipping summer into winter is not one. A wrong invariance is bias you added on purpose. Mixup blends pairs of rows and their labels, which acts like a smoothness penalty between training points.

A: averaging

Average B fits, each with variance \sigma_{f}^2 and pairwise correlation \rho:

\text{Var}

variance across possible training seasons

B

the number of fits averaged

\frac1B

divide by B: with the sum, the average of the B fits

\sum_{b=1}^{B}

add up over the B fits

\hat f_b

fit number b (e.g. 1-NN on bootstrap sample b)

\rho

the correlation between any two fits' predictions

\sigma_f^2

the variance of a single fit's prediction

Only the uncorrelated part averages away; \rho \sigma_{f}^2 is a floor no number of models removes. Bagging (Breiman, 1996) gets the B fits by refitting on bootstrap samples. On the season, 1-NN (the most variable model we have) goes from test 0.184 to 0.136 with B = 50 bootstrap fits; its variance falls from 0.091 to 0.044, bias barely moves, and bagging wins in every one of 200 seasons. The bootstrap fits have correlation \rho = 0.42 and variance 0.102 each, so the formula predicts 0.0437 against 0.0439 observed.

Trap: bagging is not automatically safe. Bagging a degree-12 polynomial explodes (median test MSE 341 against 0.151 for a single fit), because many bootstrap samples miss the edge days and a degree-12 polynomial extrapolates wildly past them. Average models whose instability is benign. Ch 10 owns ensembles.

Implicit regularization and label smoothing. Some regularization comes free with the optimizer. SGD's gradient noise prefers flat regions of the loss (Ch 2, Section 6). Gradient descent started at zero on an underdetermined least-squares problem converges to the minimum-norm solution among all perfect fits — a hidden ridge with \lambda \to 0 that Section 10 depends on. For classifiers, label smoothing trains on targets like 0.9 / 0.1 instead of 1 / 0, which stops logits growing without bound (Ch 3; transformers, Ch 21).

Remember: PANDAS: Penalty, Augmentation, Noise, Dropout, Averaging, Stopping early. For a linear model: noise = ridge with \lambda = n\cdot s^2; dropout = per-feature ridge with ((1 - q)/q)\sum_ix_{ij}^2; early stopping filters each direction by 1 - (1 - \eta d^2)^t, roughly ridge with \lambda \approx 1/(\eta t). Averaging leaves \rho \sigma^2.

Next question → all of this says: never let the model fit the training data perfectly. Yet the largest modern models fit their training data perfectly and still generalize. How? (Section 10)

10. Double Descent and Benign Overfitting

The question: on the station data, as we add features past the number of stations, test error explodes — and then comes back down. Why does the U-curve have a second act?

The 30-second version: when features p exceed rows n, infinitely many weight vectors fit the training data exactly (interpolation). Least squares via the pseudo-inverse — and gradient descent started from zero — picks the one with the smallest norm. Right at p = n there is exactly one way to pass through every point, and it is the wildest: the smallest singular value of X is near zero and the fit's variance blows up (pure noise: \sigma^2(1 + p/(n - p - 1)) below n, \sigma^2(1 + n/(p - n - 1)) above). Past the threshold, extra features give the minimum-norm rule room to find a calm interpolator, and test error falls again: double descent. Hook: at the threshold there's one way to fit everything, and it's the wildest; past it, minimum norm picks the calmest. Interpolating noise is benign when the signal lives in a few strong directions and many weak directions soak up the noise. A tuned ridge penalty removes the spike. Bias–variance still holds at every p; what failed was measuring capacity by parameter count. First read: Figure 4.7, the hook, benign overfitting in words.

Interpolation and the minimum-norm solution

With p > n and full row rank, Xw = y has infinitely many solutions. The pseudo-inverse picks the one with the smallest \lVert w\rVert (Card 16):

\hat{\mathbf w}

the minimum-norm weights: the interpolating fit least squares returns when p > n

X

the design: 100 stations by p features

X^{+}

the pseudo-inverse of X: the best "undo" when X has no true inverse

\mathbf y

the vector of station rides

X^\top

the transpose of X: rows become columns

(XX^\top)^{-1}

the inverse of the n \times n matrix of inner products between stations

\arg\min

the \mathbf w that makes what follows smallest

\lVert\mathbf w\rVert

the length (norm) of a weight vector

\mathbf w

any candidate weight vector

Gradient descent from w_0 = 0 lands on the same point, because every gradient X^\top(Xw - y) is a combination of the rows of X: the iterates never leave the row space, and the only interpolator in the row space is the minimum-norm one (Ch 1, Q67). On the double-descent design below with p = 300, gradient descent's distance to the pseudo-inverse solution is 0.48 after 100 steps, 0.0007 after 1,000, and 10^{-12} at convergence. Started from a random point instead, it converges to a different interpolator, offset by exactly the starting point's component outside the row space. The minimum-norm solution is also the \lambda \to 0 limit of ridge.

The spike at p = n

The double-descent design (sims/sim_double_descent.py, 150 replicates): n = 100 stations. Features are added in order: 20 strong "core descriptors" (variance 1, true weights decreasing like 1/j, \lVert \beta \rVert = 1), then 980 weak "map-embedding" features (variance 0.02, true weight 0). Noise \sigma = 0.5, so predicting the mean scores 1.25.

(Medians, because at p = n the mean is dominated by a few nearly singular draws.) The classical U ends at p = 10. As p approaches n, the smallest singular value of X approaches zero: one direction of the data has almost no spread, the coefficient along it is \sigma/d-sized noise, and \lVert \hat{w}\rVert explodes. At p = n = 100 the median test error is about 155 times its classical best. Then, as weak features pile up, the smallest singular value grows again, the norm falls, and test error descends to 0.351 at p = 1{,}000 — close to the classical best, with zero training error.

For pure noise (\beta = 0, isotropic Gaussian features, \sigma^2 = 1) the risk is exact (Card 16; Worked Problem 2; Drill 8):

\mathbb{E}

average over random designs and noise draws

\sigma^2

noise variance (1 in the pure-noise case)

p

the number of features

n

the number of training rows (100 stations)

At n = 100: 2.02 at p = 50, 11.0 at 90, 24.75 at 95, 26.0 at 105, 12.1 at 110, 2.01 at 200, 1.11 at 1,000; simulation agrees (2.02, 11.05, 24.7, 25.5, 12.0, 2.01, 1.11). With signal and p > n the minimum-norm fit also has bias \lVert \beta \rVert^2(1 - n/p): it can only see the part of \beta in the n-dimensional row space.

The second descent, and when the biggest model wins

In the core + embedding design the second descent returns near, but not below, the classical best. In a flat design (400 isotropic features, signal spread evenly, \lVert \beta \rVert = 1, \sigma = 0.5) every model with p < 100 is worse than predicting the mean (1.25) — each sees only a sliver of the signal — and the largest model is the best of all: 1.09 at p = 400 (sims/sim_dd_flat.py).

Ridge removes the spike

The last row of the table: with \lambda chosen by CV, ridge's test error stays between 0.309 and 0.354 for every p, with no spike (Nakkiran, Venkat, Kakade and Ma, 2021). The spike is a variance explosion along near-zero singular directions, and ridge's factor d^2/(d^2 + \lambda) kills exactly those directions. Weight decay and early stopping do the same in networks.

Benign overfitting in words

Why can a model that fits noise exactly predict well? The minimum-norm interpolator has to put the noise somewhere. If there are many weak directions — 980 low-variance embedding features — it spreads the noise thinly across all of them, and each new station's prediction picks up only a little of it. Meanwhile the few strong directions carry the signal and are estimated well. Benign overfitting needs both: signal concentrated in a few strong directions, plus a long tail of many weak directions that absorb noise harmlessly (Bartlett, Long, Lugosi and Tsigler, 2020; Hastie, Montanari, Rosset and Tibshirani, 2022). With too few weak directions, the noise concentrates and overfitting is harmful — the spike (Q28).

Three faces of double descent

The spike appears along three knobs: model-wise (parameters \approx n; Belkin, Hsu, Ma and Mandal, 2019), sample-wise (rows \approx parameters, model fixed: more data can hurt) and epoch-wise (training time, when a network first fits its training set; Nakkiran et al., 2019). The sample-wise face (sims/sim_dd_null.py, sample-wise arm: p = 100 fixed, flat signal): min-norm test error rises from 1.01 at n = 50 to 2.91 at n = 90 and 6.69 at n = 105, then falls to 0.756 at n = 150 and 0.500 at n = 200. Adding rows toward n = p made things worse. Degree 12's spike at n = 15 in Section 5 is the same effect (Q29).

Random labels. Zhang, Bengio, Hardt, Recht and Vinyals (2017) trained standard image networks on labels assigned at random. They reached zero training error anyway, so their capacity by any count is enough to memorize anything — yet on real labels the same networks generalize. Capacity of the family can't explain it; the procedure (gradient descent from small initial weights, the architecture, the data) decides which of the many perfect fits you get. This is Section 4's trap at full scale.

How this squares with Act I. The decomposition is never violated: at every p, test error = \sigma^2 + \text{bias}^2 + \text{variance}. What broke was using the parameter count as the capacity measure. The minimum-norm interpolator past the threshold has far fewer effective df than p, because minimum norm is a constraint. The classical advice still holds where most practical models live — near p \approx n, or with too few weak directions to absorb noise — so regularize there.

Remember: p > n → infinitely many perfect fits; pseudo-inverse and GD from zero pick the minimum-norm one. At the threshold, one way to fit everything — the wildest; past it, minimum norm picks the calmest. Spike at p = n, second descent after; more data can hurt near n = p. Tuned ridge removes the spike. Benign = few strong directions + many weak ones to absorb the noise.

Next question → all of this assumed next year's days come from the same process as this year's. What if they don't? (Section 11)

11. When Train and Test Differ: the Assumptions Behind Generalization, Shift and Leakage

The question: the pilot season logged mostly spring and summer days; the model will run all year — and e-bikes launch next spring. Which assumption breaks, how do we see it, and how do we correct it?

The 30-second version: a test error predicts future error only under three promises — hook Same, Separate, Single-use: future data come from the same distribution; test rows are separate from training rows (no duplicates, no shared groups, no information from the future); and the test set is used once, not tuned on. Break "Same" and you have distribution shift, in three kinds — X moves, Y moves, the rule moves: covariate shift (which days you see changes), label shift (how common each outcome is changes), concept drift (the relationship itself changes). For covariate shift, reweight training rows by w(x) = p_{\text{test}}(x)/p_{\text{train}}(x) — it helps when the model is misspecified (the line: 0.314 → 0.293 at n = 200), barely matters when the model is right (the cubic), and costs effective sample size: with 80% of days in summer, \text{ESS} = 64\% of n. Break "Separate" and you get fake generalization: four docks report the same day, a random split puts copies on both sides, and 1-NN's test error looks 98 times better than on truly new days. First read: the S-S-S table, X / Y / rule, the weighting result and ESS.

The three promises: remember, see, test, fix

Every generalization claim so far assumed i.i.d. rows. That single assumption is three promises, and each has a fingerprint, a test and a fix:

These map onto Ch 1's Who, What, How, When (Ch 1, Section 7): WHO is in the data drives covariate shift; WHEN drives drift and leakage. Q32 walks a production failure through all three promises.

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_predict
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(0)
n_tr, n_te = 300, 300
summer = rng.random(n_tr) < 0.8                              # pilot: 80% of days in the summer half
x_tr = np.where(summer, rng.uniform(0, 0.5, n_tr), rng.uniform(0.5, 1, n_tr))
x_te = rng.uniform(0, 1, n_te)                               # deployment: the whole year

# SEE + TEST (Same): adversarial validation - can a classifier tell train days from deployment days?
X = np.polynomial.legendre.legvander(2 * np.r_[x_tr, x_te] - 1, 3)[:, 1:]
z = np.r_[np.zeros(n_tr), np.ones(n_te)]                     # 1 = deployment
p = cross_val_predict(LogisticRegression(), X, z, cv=5, method="predict_proba")[:, 1]
print("adversarial AUC:", round(roc_auc_score(z, p), 3))     # ~0.5 means no detectable shift

# FIX: importance weights w(x) = p_test(x) / p_train(x) from the domain classifier, then check the ESS
w = (p[:n_tr] / (1 - p[:n_tr])) * (n_tr / n_te)
print("mean weight, summer vs winter:", w[summer].mean().round(2), w[~summer].mean().round(2))
print("ESS fraction:", round(w.sum() ** 2 / (w ** 2).sum() / n_tr, 2))

# TEST (Separate): no group may sit on both sides of a split
day_train, day_test = {1, 3, 5, 7}, {7, 9, 11}
print("days on both sides:", day_train & day_test)

The AUC comes out well above 0.5, so the shift is detectable. The estimated weights are smoother than the true 0.625 and 2.5, because a cubic logistic model can't draw a sharp step at September 20; for the same reason the printed ESS fraction (0.8) is more flattering than the true 0.64 derived below. (The full simulation, sims/sim_shift.py, gets an AUC of 0.649 against a theoretical 0.65, estimated weights of 0.67 and 2.27, and an AUC of 0.49 when there is no shift.)

X moves, Y moves, the rule moves

Covariate shift on the running example

This is data bias in the sense of Ch 1, Section 7: the rows are an unrepresentative sample of the days the model will serve. Training days come from the mixture 0.8 \times \text{Uniform}(0, 0.5) + 0.2 \times \text{Uniform}(0.5, 1) (80% spring–summer); deployment is uniform over the year. The densities are 1.6 and 0.4 against 1, so the true weights are 0.625 for summer days and 2.5 for winter days. Weighted least squares with those weights minimizes the training loss as if the days had come from the deployment distribution (Card 17):

\mathbb{E}_{\text{train}}

average over days drawn like the pilot's (80% spring and summer)

\mathbb{E}_{\text{test}}

average over deployment days, spread evenly over the year

w

a day's importance weight: how much more common such days are at deployment than in the pilot

X

a random day's time of year

Y

that day's rides

\ell

the loss for one prediction (squared error)

f

the fitted curve

x

a particular time of year

p_{\text{test}}

how common day x is at deployment (density 1 everywhere)

p_{\text{train}}

how common day x is in the pilot: 1.6 in spring and summer, 0.4 otherwise

Weighting matters when the model is misspecified: a line can't fit the whole sine, so which days it fits best decides where its errors land, and the unweighted line spends its effort on summer (Shimodaira, 2000; Q12). Weighting wins in 90% of replicates for the line at n = 200. When the model can represent the truth, both fits converge to it and weighting only adds variance: for degree 5 it wins in only 39% of replicates.

Effective sample size and overlap

Weights cost data. The effective sample size is the number of equally weighted rows that would give the same variance:

\sum_i

add up over training days i

w_i

training day i's importance weight

n

the number of training days

So 30 pilot days are worth 19.2 days of deployment-like data (Drill 20). For a Gaussian mean shift of \mu standard deviations, \text{ESS}/n = e^{-\mu^2}: 37% at \mu = 1, under 2% at \mu = 2 (Drill 13), and at that point the weights are so heavy-tailed that even a simulation struggles to estimate the ESS. Two practical rules: clip or self-normalize extreme weights, and check overlap. If the pilot had no winter days at all, no weight can fix it: the line's winter error is 1.94 and degree 5's median is 8,349, pure extrapolation. Hook: weight by test over train; check the ESS (Q33).

Estimating the weights: train a domain classifier on train vs target rows (as in the code), then w(x) = [P(\text{target} \mid x)/P(\text{train} \mid x)] \times(n_{\text{train}}/n_{\text{target}}). Cross-validate with the same weights — importance-weighted CV (Sugiyama, Krauledat and Müller, 2007).

Label shift and concept drift

Label shift is a classifier's problem: if the base rate of the outcome changes but p(x \mid y) stays, adjust predicted probabilities as in Ch 3's prior-shift correction (Ch 3, Card 13), estimating the new base rate from the confusion matrix on unlabeled target data (black-box shift estimation; Lipton, Wang and Smola, 2018).

Concept drift: next spring e-bikes make every seasonal swing 1.3 times larger. The cubic trained on this season scored 0.106; on the new truth it scores 0.153, and the extra 0.045 is exactly the new bias (0.3 \times \sin)^2, averaged over the year. No reweighting fixes it, because p(y \mid x) itself moved. Retrain on recent windows, weight recent days more, and monitor residuals over time (Ch 28).

Leakage: fake generalization

Each season day is reported by four docks, each adding a tiny reporting error (SD 0.02): 120 rows that are really 30 days. A random 80/20 row split puts twins of almost every test row (99.7%) into training (sims/sim_leak.py, 200 replicates):

The random split says 1-NN is 98 times better than it is: the nearest neighbor of each test row is its own twin. Degree 15 looks 6.5 times better. Grouping by day removes the leak; for degree 15 the grouped error is heavy-tailed (holding out edge days forces extrapolation from 24 days), so quote its median. Ch 5 owns the full leakage taxonomy (Q34).

Remember: a test score is a promise: Same, Separate, Single-use. X moves, Y moves, the rule moves. Weight by p_test/p_train when the model is misspecified; check ESS and overlap. Duplicates across the split = fake generalization.

Next question → we have the whole theory. Can we build the U-curve, the decomposition, ridge and CV ourselves in a page of NumPy? (Section 12)

12. Build It From Scratch

The question: can we reproduce the U-curve, the bias–variance split, ridge by SVD, and K-fold CV in about 50 lines of NumPy?

The 30-second version: build Legendre polynomial features; solve ridge by SVD with an unpenalized intercept; wrap any fit in a shuffled K-fold loop; rerun the season many times to measure \text{bias}^2 and variance. The numbers match Sections 2, 3, 6 and 8. First read: run the code once.

Four functions, NumPy only, each mapping to a section: poly (Section 2), fit_ridge (Section 6, Card 6), kfold_cv (Section 8, Card 11) and bias_variance (Section 3, Card 2). Least squares is fit_ridge with \lambda = 0. Q27 asks you to write this in 15 minutes.

import numpy as np

rng = np.random.default_rng(404)
f = lambda x: np.sin(2 * np.pi * x)                        # the true seasonal curve
sigma, n = 0.3, 30
x_test = rng.uniform(0, 1, 10_000)                         # next year's days
y_test = f(x_test) + rng.normal(0, sigma, 10_000)
x = (np.arange(1, n + 1) - 0.5) / n                         # one day every ~12 days
Y = f(x) + rng.normal(0, sigma, (500, n))                   # 500 possible seasons; Y[0] is this season

def poly(x, d):                                             # Legendre features on [0, 1]
    return np.polynomial.legendre.legvander(2 * x - 1, d)

def fit_ridge(X, y, lam):
    """Ridge on the non-constant columns (standardized), intercept unpenalized; returns predict(Xnew)."""
    Z = X[:, 1:]
    mu, sd = Z.mean(0), (Z.std(0) if Z.shape[1] else np.ones(0))
    U, s, Vt = np.linalg.svd((Z - mu) / sd, full_matrices=False)
    w = Vt.T @ (s / (s ** 2 + lam) * (U.T @ (y - y.mean())))   # shrink each direction by s^2/(s^2 + lam)
    return lambda Xn: y.mean() + ((Xn[:, 1:] - mu) / sd) @ w

def kfold_cv(x, y, d, lam, K=10, seed=0):
    folds = np.array_split(np.random.default_rng(seed).permutation(len(x)), K)
    err = 0.0
    for te in folds:
        tr = np.setdiff1d(np.arange(len(x)), te)
        err += np.sum((y[te] - fit_ridge(poly(x[tr], d), y[tr], lam)(poly(x[te], d))) ** 2)
    return err / len(x)

def bias_variance(d, lam, R=200):
    grid = np.linspace(0, 1, 201)
    P = np.array([fit_ridge(poly(x, d), Y[r], lam)(poly(grid, d)) for r in range(R)])
    bias2 = np.mean((P.mean(0) - f(grid)) ** 2)
    var = np.mean(P.var(0))
    y_new = f(grid) + rng.normal(0, sigma, (R, len(grid)))  # fresh days to score on
    return bias2, var, sigma ** 2, bias2 + var + sigma ** 2, np.mean((y_new - P) ** 2)

print("degree  train   test   (this season)")
for d in [0, 1, 3, 9, 15]:
    pred = fit_ridge(poly(x, d), Y[0], 0.0)
    print(f"{d:6d} {np.mean((Y[0] - pred(poly(x, d))) ** 2):7.4f} {np.mean((y_test - pred(poly(x_test, d))) ** 2):7.4f}")
cv = [kfold_cv(x, Y[0], d, 0.0) for d in range(16)]
print("10-fold CV picks degree", int(np.argmin(cv)))
for d in [1, 3, 15]:
    print(f"degree {d}: bias2 %.4f  var %.4f  noise %.2f  sum %.4f  simulated test %.4f" % bias_variance(d, 0.0))
lams = np.logspace(-4, 3, 29)
lam = lams[np.argmin([kfold_cv(x, Y[0], 15, l) for l in lams])]
pred = fit_ridge(poly(x, 15), Y[0], lam)
print(f"ridge degree 15: CV lambda {lam:.3g}, test MSE {np.mean((y_test - pred(poly(x_test, 15))) ** 2):.4f}")

It runs in under a second, and every printed value is recorded in numbers.json (scratch_*). What to notice:

• The U on one season. Train / test MSE for degrees 0, 1, 3, 9, 15: 0.754 / 0.583, 0.387 / 0.296, 0.092 / 0.114, 0.070 / 0.136, 0.031 / 0.888. This season is one of the bad ones for degree 15 (its 500-season mean is 0.379). 10-fold CV picks degree 3.
• The decomposition adds up. Degree 1: 0.200 + 0.005 + 0.09 = 0.295 vs simulated 0.296. Degree 3: 0.005 + 0.011 + 0.09 = 0.106 vs 0.106. Degree 15: 0.001 + 0.474 + 0.09 = 0.565 vs 0.565. Degree 1 is bias, degree 15 is variance.
• Degree 15's variance is 0.474 here, not Section 3's 0.288, because this code scores on a grid that includes the endpoints 0 and 1, where degree 15 is wildest (Section 3: variance peaks at the edges). Same procedure, different place to look.
• Ridge by SVD shrinks each direction by s^2/(s^2 + \lambda) — the s / (s ** 2 + lam) line — and at the CV-chosen \lambda = 5.62 brings degree 15 from 0.888 to 0.137 on this season.

Common mistake: selecting features on all rows, then cross-validating. Take 50 rows, 1,000 pure-noise features and a pure-noise y (no fit can beat an MSE of 1.0). Keep the 10 features most correlated with y, then 5-fold CV least squares on them: the CV MSE comes out 0.54 — "R^2 = 0.46" from noise — and below 1.0 in 99% of 200 replicates. Screening inside each training fold gives 1.43, the honest answer (sims/sim_leak_cv.py). By contrast, fitting the scaler on all 30 days before CV changes the least-squares CV curve by exactly 0 (an affine rescaling of columns doesn't change a least-squares fit with an intercept) and ridge's by about 0.003. Supervised steps are the ones that leak (Section 8).

Try it: set sigma = 0.6: every variance grows fourfold (it is proportional to \sigma^2) while \text{bias}^2 stays put. Then set n = 300 and add 12 to the bias–variance loop: degree 12's variance collapses below the cubic's \text{bias}^2, and it overtakes the cubic. More data cures variance.

Remember: features → fit (SVD, shrink by s^2/(s^2 + \lambda)) → K-fold → rerun the world to see bias and variance.

Next question → how does this look in production code, and which library defaults bite? (Section 13)

13. Use It Like a Pro

The question: how do practitioners tune, validate and regularize in scikit-learn and PyTorch, and which defaults cause silent mistakes?

The 30-second version: put every learned step in a Pipeline and tune with GridSearchCV or the built-in RidgeCV/LassoCV/ElasticNetCV; report a nested-CV or held-out score, never best_score_. Three gotchas, hook the 3 Ss: Scaler inside the pipeline; Shuffle your folds (KFold does not shuffle by default); the Same \lambda is not the same alpha — Ridge uses \lVert y - Xw\rVert^2 + \alpha \lVert w\rVert^2, while Lasso/ElasticNet use (1/2n)\lVert y - Xw\rVert^2 + \alpha \lVert w\rVert_1, so a lasso alpha is \lambda/(2n) in this chapter's convention. In PyTorch, weight_decay is L2 in SGD and decoupled in AdamW; model.eval() turns dropout off; early stopping restores the best checkpoint. First read: the gotcha table.

The tuning pipeline

import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge, RidgeCV
from sklearn.model_selection import GridSearchCV, KFold, cross_val_score

rng = np.random.default_rng(404)                             # same draws as Sections 2-12
x_test = rng.uniform(0, 1, 10_000)
y_test = np.sin(2 * np.pi * x_test) + rng.normal(0, 0.3, 10_000)
x = (np.arange(1, 31) - 0.5) / 30
y = np.sin(2 * np.pi * x) + rng.normal(0, 0.3, (500, 30))[0]   # this season
X, Xt = x.reshape(-1, 1), x_test.reshape(-1, 1)

pipe = make_pipeline(PolynomialFeatures(include_bias=False), StandardScaler(), Ridge())
grid = {"polynomialfeatures__degree": range(1, 13), "ridge__alpha": np.logspace(-4, 3, 8)}
inner = KFold(10, shuffle=True, random_state=0)               # Shuffle explicitly: KFold() does not
gs = GridSearchCV(pipe, grid, cv=inner, scoring="neg_mean_squared_error").fit(X, y)
print("chosen:", gs.best_params_)
print("best_score_ (chosen on itself):", round(-gs.best_score_, 4))
print("test MSE on 10,000 new days:", round(np.mean((y_test - gs.predict(Xt)) ** 2), 4))

outer = KFold(5, shuffle=True, random_state=1)                # tune inside, test outside
nested = cross_val_score(GridSearchCV(pipe, grid, cv=inner, scoring="neg_mean_squared_error"),
                         X, y, cv=outer, scoring="neg_mean_squared_error")
print("nested CV MSE:", round(-nested.mean(), 4))

rcv = make_pipeline(PolynomialFeatures(15, include_bias=False), StandardScaler(),
                    RidgeCV(alphas=np.logspace(-4, 3, 29))).fit(X, y)   # efficient LOO when cv=None
print("RidgeCV alpha on raw powers:", rcv[-1].alpha_)

It picks degree 3 with alpha 10^{-4}; best_score_ is 0.122, the test MSE on new days 0.113, and the nested estimate 0.160 (its five outer folds range from 0.010 to 0.385 — with 30 rows, any single number is rough). The nested number estimates the procedure; best_score_ was chosen on the same folds it reports. Two lessons hide in the last line. RidgeCV with cv=None uses the exact LOO shortcut of Section 8 at no extra cost. And it picks alpha = 0.0018 on standardized raw powers, where Section 6's Legendre features picked 5.62 for the same data: the scale of \lambda depends on the basis, so never carry a \lambda across feature definitions. LassoCV and ElasticNetCV(l1_ratio=[.1, .5, .9, 1]) work the same way; they fit the whole path with warm starts, and a ConvergenceWarning means raise max_iter or standardize. validation_curve and learning_curve in sklearn.model_selection produce Figures 4.6(a) and 4.3 for any estimator.

Gotchas: defaults that bite

PyTorch: weight decay, dropout, early stopping

try:
    import torch                                              # illustrative: runs only where PyTorch is installed
except ImportError:
    torch = None
if torch is not None:
    torch.manual_seed(0)
    X, y, Xv, yv = torch.randn(200, 15), torch.randn(200, 1), torch.randn(100, 15), torch.randn(100, 1)
    model = torch.nn.Sequential(torch.nn.Linear(15, 64), torch.nn.ReLU(), torch.nn.Dropout(0.1), torch.nn.Linear(64, 1))
    opt = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.1)       # decoupled weight decay
    loss = torch.nn.functional.mse_loss
    best, best_state, bad = float("inf"), None, 0
    for epoch in range(500):
        model.train(); opt.zero_grad(); loss(model(X), y).backward(); opt.step()   # dropout on while training
        model.eval()                                                               # dropout off to evaluate
        with torch.no_grad():
            v = loss(model(Xv), yv).item()
        if v < best:
            best, bad, best_state = v, 0, {k: t.clone() for k, t in model.state_dict().items()}
        elif (bad := bad + 1) >= 10:                                               # patience: 10 checks
            break
    model.load_state_dict(best_state)                                              # restore the best checkpoint

torch.optim.SGD(weight_decay=...) adds \lambda w to the gradient, which is L2. Adam(weight_decay=...) does the same, but Adam then rescales it per weight, so it is no longer ridge; AdamW shrinks the weights directly (Q37).

Before you ship

• [ ] Looked at the learning curve: is the problem bias or variance? (Section 5)
• [ ] Reported a nested-CV or untouched-holdout score with its SE, not best_score_ (Section 8)
• [ ] Considered the 1-SE choice (Section 8)
• [ ] Every learned step inside the pipeline (Sections 8, 12)
• [ ] Ran adversarial validation against recent or target data (Section 11)
• [ ] If features will be interpreted: reported selection frequencies and held-out inference (Section 7)

Remember: pipeline everything, tune inside, report outside. The 3 Ss: Scaler inside, Shuffle explicitly, Same \lambda \ne same alpha. \lambda's scale depends on the basis.

Next question → where does generalization control earn its keep, and what does it look like in a model with a trillion training tokens? (Section 14)

14. In the Real World, and How It Connects to LLMs

The question: where do these ideas decide real outcomes, and what do overfitting and regularization look like in large language models?

The 30-second version: every production model is tuned by the ideas here: ridge and lasso in demand forecasting, credit scorecards and genomics; heavy L2 in ads; shrinkage in experiment readouts. LLMs flipped the regime: they have far more parameters than the classical threshold allows, yet see most tokens about once, so pretraining behaves like the "more data" cure at enormous scale, and scaling laws are learning curves — loss falls as a power law in parameters, data and compute. Weight decay (AdamW) remains standard; dropout is often off in large pretraining; fine-tuning on small datasets brings classical overfitting back, handled by early stopping, low-rank adapters (LoRA) as capacity control, and RLHF's KL penalty, which pulls the policy toward a reference model the way ridge pulls weights toward zero. First read: the two tables.

Where it earns its keep

How this chapter becomes an LLM

Pretraining is the "more data" regime, where variance is small and scaling laws describe the bias-like floor; fine-tuning on thousands of examples is the "30 days" regime, where every PANDAS letter earns its keep again.

Remember: the cures scale. Scaling laws are learning curves, weight decay is ridge, the RLHF KL penalty is shrinkage toward a prior, LoRA is capacity control. Fine-tuning brings classical overfitting back.

We've controlled overfitting. How do we measure how good a model is — honestly enough to convince someone else? That is Chapter 5.

Next → The whole chapter on one page (Section 15). Then practice: the Derivations tab, then the Interview Bank.

15. The Chapter on One Page

One map of the whole chapter and every memory hook in one table. Use it before an interview, or to check a second read: if any row feels unfamiliar, go back to its section.

The memory map

Every hook in one place

Remember: if you can redraw the memory map on a blank page and say the hook for every row of the table, you have the chapter. Now practice.

16. Practice

The chapter ends here. Practice lives in two tabs, so this page stays a story you can read straight through.

• Ch 4 · Derivations — 19 derivation cards built on one toolkit, S-E-A-T (rotate with the SVD, take the Expectation over noise or prior, Add and subtract the mean, use the Trace trick), about 24 formulas to memorize, and a mastery tracker.
• Ch 4 · Interview Bank — 39 questions in four levels, each answered in four beats; two worked research-level problems (Bayes-optimal ridge on the dense stations, Figure 4.10; the spike at p = n, Figure 4.11) solved with one recipe, Rotate, filter, split, sum; 20 drills; and about 38 rapid-fire rows for the day before.

Suggested order: finish the chapter's second read → Derivations Cards 1–5 → Interview Bank Levels 1–2 → Cards 6–19 → Levels 3–4.

We've controlled overfitting. Next question: how do we measure how good a model is — honestly enough to convince someone else? That is Chapter 5, Evaluation.