Empirical Risk Minimization and Generalization

Training reduces a loss measured on the available sample; prediction requires that performance carry over to new observations. Generalization theory studies how the sample loss relates to the population risk. A small gap alone is not enough: both losses could be large. We will compare model families, regularization, and repeated selection using NumPy and scikit-learn. Run the code blocks in order; the examples assume the means, sampling variation, and model-fitting vocabulary introduced in the preceding articles.

Risk and empirical risk

For a loss \(\ell\) and a distribution \(\mathcal{D}\) over pairs \((x, y)\), the risk of a model \(f\) is \(R(f) = \mathbb{E}_{(x,y) \sim \mathcal{D}}[\ell(f(x), y)]\). In a real problem it cannot be computed directly, because \(\mathcal{D}\) is unknown. The empirical risk is its sample average \(\hat{R}(f) = \frac{1}{n}\sum_i \ell(f(x_i), y_i)\), and empirical risk minimization is the strategy of choosing the \(f\) that minimizes it. If a model’s squared losses on three rows are 1, 0, and 4, its empirical risk is \((1+0+4)/3 = 5/3\); for a fixed model with finite expected loss, the sample average converges to the risk as the number of iid observations grows. The set \(\mathcal{F}\) written below is simply the collection of models the search is allowed to choose from — every polynomial of degree at most 9, say, or every random forest with the given hyperparameters.

Take a fixed \(f\), chosen before seeing the data, and evaluate it on independent, identically distributed draws from \(\mathcal{D}\). If the expected loss is finite, \(\hat{R}(f)\) is an unbiased estimate of \(R(f)\); if the loss also has finite variance \(\sigma_\ell^2\), the standard error of that estimate is \(\sigma_\ell/\sqrt{n}\). A stronger claim — that the probability of \(|\hat{R}(f) – R(f)| > \epsilon\) decreases exponentially with \(n\) — requires an additional condition, such as a bounded loss. The unbiasedness is no longer guaranteed once \(f\) is selected using the same data, because then \(f\) depends on the sample and the two are not independent. That does not mean such an estimate must be biased; it means nothing here rules it out.

Selecting a model because it scored well on a sample can select favorable sampling noise along with genuine performance. A separate evaluation sample avoids reusing that same noise for selection and reporting, provided it is independent and represents the target distribution. Training loss remains useful as an optimization diagnostic; independence is needed for the simple unbiased evaluation argument above, not for the loss to have any meaning.

Capacity and the generalization gap

A uniform bound applies to every candidate in a class at once. Suppose \(\mathcal F\) is fixed before the sample is drawn. Under suitable sampling and loss assumptions, a bound of the form \(R(f)\leq\hat R(f)+\text{complexity}(\mathcal F,n,\delta)\) holds with probability at least \(1-\delta\) over the sample, simultaneously for every \(f\in\mathcal F\). Here \(\delta\) is the allowed probability that the guarantee fails. Covering every candidate also covers one selected from the class using that sample. Choosing the class itself from the data needs additional justification.

For a fixed finite class of \(M\) models, iid samples, and losses in [0, 1], one such complexity term is \(\sqrt{\log(M/\delta)/(2n)}\). With \(M=1000\), \(n=1000\), and \(\delta=0.05\), it is about 0.0704. Thus, with probability at least 95%, every candidate’s risk is at most its sample loss plus 0.0704. This is an upper bound, not a prediction of the actual gap. Infinite classes need other measures, such as VC dimension or Rademacher complexity; their derivations are beyond this introduction.

The following polynomial experiments use squared loss with Gaussian noise, which is unbounded. They therefore do not verify this bounded-loss formula. They illustrate how the fitted model and its generalization gap change as the available curves or fitting rule change.

The experiment below changes the degree of a polynomial model. A degree-2 model predicts \(c_0+c_1x+c_2x^2\): each input becomes a row \((1,x,x^2)\), and fitting chooses the coefficients. This is the matrix-and-weights setup from the linear algebra article. The curve can bend with \(x\) while remaining linear in its coefficients. Increasing the degree adds powers of \(x\), allowing more possible curves.

import numpy as np
import warnings

rng = np.random.default_rng(2)
x_tr = np.sort(rng.uniform(-1, 1, 60))
y_tr = np.sin(3 * x_tr) + 0.2 * rng.normal(size=60)
x_te = np.sort(rng.uniform(-1, 1, 2000))
y_te = np.sin(3 * x_te) + 0.2 * rng.normal(size=2000)

for deg in (1, 3, 9, 20, 40):
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        c = np.polyfit(x_tr, y_tr, deg)
    tr = np.mean((np.polyval(c, x_tr) - y_tr) ** 2)
    te = np.mean((np.polyval(c, x_te) - y_te) ** 2)
    print(f"degree {deg:2d}  train {tr:.4f}  test {te:10.3e}  gap {te - tr:10.3e}")
    for warning in caught:
        print(f"  {warning.category.__name__} at degree {deg}: {warning.message}")
# degree  1  train 0.1638  test  2.148e-01  gap  5.097e-02
# degree  3  train 0.0383  test  4.448e-02  gap  6.136e-03
# degree  9  train 0.0349  test  6.463e-02  gap  2.968e-02
# degree 20  train 0.0268  test  9.683e+01  gap  9.681e+01
# degree 40  train 0.0138  test  6.102e+12  gap  6.102e+12
# RankWarning at degree 40: Polyfit may be poorly conditioned

Training error falls with degree here. The minimum achievable training loss cannot increase as the class grows, since a higher-degree polynomial can represent everything a lower-degree one can — but that is a statement about exact minimization, and what polyfit returns at high degree is a numerical approximation that may not attain it. Among the degrees tested, test error is lowest at degree 3 and becomes extremely large at degrees 20 and 40; degree 40 fits 60 points with 41 parameters and produces test error of order \(10^{12}\).

The code captures warnings around each fit and prints them beside its result. In the recorded run, degree 40 reports RankWarning: some singular values of the scaled polynomial design matrix fall below the fitting cutoff. Its columns contain powers of the inputs, which can become nearly linearly dependent on this interval. This numerical rank warning does not quantify how much of the test error is statistical overfitting versus numerical sensitivity. Scientific notation reduces the printed precision but does not guarantee reproduction across NumPy versions or linear algebra backends; the high-degree results remain environment-sensitive.

Capacity is not the number of parameters

When a hypothesis class is very rich relative to the sample size, a worst-case bound can become too loose to guarantee useful predictive performance. That is not a prediction that such a model will fail to generalize — it says the bound alone cannot certify that it will succeed. The next example fits random forests to a structured target and to independent random labels. A forest fits many trees on resampled data with randomized feature choices, then combines their predicted class probabilities.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

rng = np.random.default_rng(0)
X = rng.normal(size=(500, 10))
structured = (X[:, 0] + X[:, 1] > 0).astype(int)
noise = rng.integers(0, 2, 500)

for name, y in (("real structure", structured), ("random labels", noise)):
    X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.4, random_state=0)
    m = RandomForestClassifier(n_estimators=200, random_state=0).fit(X_tr, y_tr)
    print(f"{name:15s} train {accuracy_score(y_tr, m.predict(X_tr)):.3f}"
          f"  test {accuracy_score(y_te, m.predict(X_te)):.3f}")
# real structure  train 1.000  test 0.960
# random labels   train 1.000  test 0.515

The same estimator with the same hyperparameters achieves perfect training accuracy on both. Fitting one random labeling does not prove the class can fit every labeling, but the fitted forests have enough flexibility to interpolate these particular training samples. The two forests are not matched on size: the same hyperparameters produce different numbers of nodes on different labels. The structured problem gives test accuracy 0.960; random labels give 0.515, close to their expected accuracy of 0.5.

The experiment shows that fitting one random labeling and generalizing on structured data are compatible; worst-case capacity alone does not distinguish the two outcomes, because it is the same in both runs while the results are not. What differs is the interaction between the data and the algorithm’s inductive bias — the preferences and assumptions it uses to predict cases it has not seen. That bias is unchanged across the two runs — trees split on axis-aligned feature thresholds, and the forest averages their class probabilities — but the random labels contain no input-label relationship for it to exploit.

This observation, made forcefully for deep networks by Zhang and colleagues, is a limit on what worst-case capacity measures explain rather than a refutation of the bounds themselves. Capacity, the fitting algorithm, and the structure of the data all matter. The experiment changes the labels while holding the model settings fixed; it does not compare classes with equal parameter counts.

No free lunch

Fix a finite input space and suppose every binary labeling of it is equally likely. Then under zero–one loss, every learning algorithm has the same expected error on the inputs it did not see in training. The reason is short: knowing the training labels says nothing about the rest, so each unseen label remains a fair coin no matter what the learner does with what it saw. The result is provable by symmetry and can be checked exhaustively on a small enough problem.

import itertools

points = list(itertools.product([0, 1], repeat=3))    # 8 inputs
train, test = [0, 1, 2, 3], [4, 5, 6, 7]

def majority(labels):    return [int(sum(labels) * 2 > len(labels))] * 4
def always_zero(labels): return [0] * 4
def parity(labels):      return [sum(points[i]) % 2 for i in test]

scores = {"majority": 0, "always 0": 0, "parity": 0}
for bits in itertools.product([0, 1], repeat=8):      # all 256 target functions
    tr_labels = [bits[i] for i in train]
    te_labels = [bits[i] for i in test]
    for name, learner in (("majority", majority), ("always 0", always_zero),
                          ("parity", parity)):
        scores[name] += sum(a == b for a, b in zip(learner(tr_labels), te_labels))

for name, correct in scores.items():
    print(f"{name:9s} {correct}/1024 = {correct / 1024:.4f}")
# majority  512/1024 = 0.5000
# always 0  512/1024 = 0.5000
# parity    512/1024 = 0.5000

All three score exactly 0.5. Majority uses the four observed labels; always-zero ignores them; parity predicts whether an input has an odd number of 1s, also ignoring the observed labels. Each method is evaluated on four held-out inputs for each of 256 possible labelings, giving 1,024 predictions. This is an exhaustive enumeration, not a simulation — the equality is exact.

The theorem is often misread as saying no algorithm is better than another. What it says is narrower and more useful: superiority cannot come from the algorithm alone, it has to come from a match between the algorithm’s assumptions and the problems you actually face. In a particular application, assumptions such as smoothness, sparsity, or local structure may be useful. The uniform-labeling average does not tell us which assumptions hold there, or which algorithm will be best under the actual distribution and constraints.

The practical form of this is that “which algorithm is best” has no general answer, and benchmark leaderboards measure a match to one collection of problems. Prior knowledge can narrow the candidate families; validation can then compare their performance on the problem at hand.

Exercises

1. The learning curve. Hold the model capacity fixed at degree 9 and vary the training set size from 15 to 1000, over 200 draws at each size. Report the mean training error, the median test error, and the median of the per-draw gap.

You should get: a gap that becomes much smaller over the sample sizes tested, while the training error rises to meet the test error.

import numpy as np

rng = np.random.default_rng(0)
x_te = np.sort(rng.uniform(-1, 1, 5000))
y_te = np.sin(3 * x_te) + 0.2 * rng.normal(size=5000)

for n in (15, 30, 60, 250, 1000):
    trs, tes = [], []
    for _ in range(200):
        xs = rng.uniform(-1, 1, n)
        ys = np.sin(3 * xs) + 0.2 * rng.normal(size=n)
        c = np.polyfit(xs, ys, 9)
        trs.append(np.mean((np.polyval(c, xs) - ys) ** 2))
        tes.append(np.mean((np.polyval(c, x_te) - y_te) ** 2))
    trs, tes = np.asarray(trs), np.asarray(tes)
    print(f"n={n:5d}  mean train {trs.mean():.4f}  median test {np.median(tes):.4f}"
          f"  median gap {np.median(tes - trs):.4f}")
# n=   15  mean train 0.0131  median test 9.6369  median gap 9.6149
# n=   30  mean train 0.0260  median test 0.0879  median gap 0.0682
# n=   60  mean train 0.0341  median test 0.0506  median gap 0.0180
# n=  250  mean train 0.0386  median test 0.0420  median gap 0.0036
# n= 1000  mean train 0.0397  median test 0.0407  median gap 0.0011

2. Regularization with a fixed parameterization. Fit a degree-15 polynomial to 30 points with ridge penalties from 0 to 10, and report the training error, test error, and coefficient norm at each.

You should get: a large change in test error while the polynomial feature set and number of fitted parameters stay fixed.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures

rng = np.random.default_rng(1)
x_tr = np.sort(rng.uniform(-1, 1, 30))
y_tr = np.sin(3 * x_tr) + 0.2 * rng.normal(size=30)
x_te = np.sort(rng.uniform(-1, 1, 5000))
y_te = np.sin(3 * x_te) + 0.2 * rng.normal(size=5000)

P = PolynomialFeatures(15, include_bias=False)
A = P.fit_transform(x_tr.reshape(-1, 1))
B = P.transform(x_te.reshape(-1, 1))

def report(label, m):
    tr = np.mean((m.predict(A) - y_tr) ** 2)
    te = np.mean((m.predict(B) - y_te) ** 2)
    print(f"{label:16s}  train {tr:.4f}  test {te:9.4f}"
          f"  |w| {np.linalg.norm(m.coef_):10.1f}")

report("no penalty", Ridge(alpha=0.0, solver="svd").fit(A, y_tr))
for lam in (1e-12, 1e-6, 1e-3, 1e-1, 10.0):
    report(f"lambda {lam:8.1e}", Ridge(alpha=lam).fit(A, y_tr))
# no penalty        train 0.0110  test 1369.5112  |w|    83279.1
# lambda  1.0e-12   train 0.0111  test 1213.0929  |w|    74373.7
# lambda  1.0e-06   train 0.0211  test    1.2931  |w|      167.7
# lambda  1.0e-03   train 0.0253  test    0.0484  |w|        8.1
# lambda  1.0e-01   train 0.0428  test    0.0605  |w|        2.6
# lambda  1.0e+01   train 0.2822  test    0.3208  |w|        0.5

3. Validation-set overfitting. Score \(m\) useless models on 200 validation rows and report the best score observed, for \(m = 1, 10, 100, 1000\). All models have true accuracy 0.5, and their validation correctness indicators are independent in this toy setup. The code draws their counts directly from a binomial distribution instead of fitting models.

You should get: a best-of-many validation score that is well above the truth, with no model having learned anything.

import numpy as np

rng = np.random.default_rng(2)
for m in (1, 10, 100, 1000):
    best = [ (rng.binomial(200, 0.5, size=m) / 200).max() for _ in range(2000) ]
    print(f"{m:5d} candidates: best validation accuracy {np.mean(best):.4f}"
          f"  (true accuracy 0.5000)")
#     1 candidates: best validation accuracy 0.4998  (true accuracy 0.5000)
#    10 candidates: best validation accuracy 0.5548  (true accuracy 0.5000)
#   100 candidates: best validation accuracy 0.5885  (true accuracy 0.5000)
#  1000 candidates: best validation accuracy 0.6143  (true accuracy 0.5000)

References

• Zhang, C., Bengio, S., Hardt, M., Recht, B., & Vinyals, O. (2017). Understanding deep learning requires rethinking generalization. ICLR.
• Shalev-Shwartz, S., & Ben-David, S. (2014). Understanding Machine Learning: From Theory to Algorithms. Chapter 4 discusses uniform convergence for finite hypothesis classes.
• Wolpert, D. H. (1996). The Lack of A Priori Distinctions Between Learning Algorithms. Neural Computation.