Learning theory: uniform convergence, margins, and why overparameterized models generalize

Why this subject matters now

For two decades the working theory of generalization was a capacity budget. A model class has some complexity measure (parameter count, VC dimension, a norm), the generalization gap scales like the square root of complexity over sample size, and the practitioner's job is to keep capacity matched to data. Deep learning broke this story in public. Zhang, Bengio, Hardt, Recht, and Vinyals (2017) showed that standard architectures can fit ImageNet-scale data with randomized labels, meaning the function class the optimizer searches is expressive enough to memorize noise, so any bound based on the class alone is vacuous for these models. Yet the same networks, trained on real labels, generalize well. Belkin, Hsu, Ma, and Mandal (2019) then showed the test-error curve itself is not the U shape in the textbooks. Past the interpolation threshold it descends a second time. The response of the theory community has been the most productive period in the field since the 1990s. Benign overfitting was characterized exactly for linear regression by Bartlett, Long, Lugosi, and Tsigler at Berkeley, the implicit bias of gradient descent was pinned down by Soudry and collaborators at the Technion with parallel work by Ji and Telgarsky, PAC-Bayes bounds were made numerically nonvacuous for real networks by Dziugaite and Roy, and stability arguments were extended to SGD by Hardt, Recht, and Singer. What a practitioner is expected to know today is therefore double. The classical machinery still matters because it governs test sets, model selection, and every finite-sample guarantee actually deployed, and the modern corrections matter because interview questions and research discussions now routinely turn on why the classical machinery does not explain the models in production. Both halves are derivations, not slogans, and this page works through them.

Concentration of measure from first principles

Everything in this subject rests on one primitive, that an average of independent random quantities is close to its expectation with overwhelming probability. The chain of inequalities that quantifies "overwhelming" starts from almost nothing.

Markov and Chebyshev

Markov's inequality. If \(X \ge 0\) and \(a > 0\), then \(\P(X \ge a) \le \E[X]/a\). The proof is one line of honest accounting. Since \(X \ge 0\),

where the second inequality holds because on the event \(\{X \ge a\}\) the integrand is at least \(a\), and the discarded term is nonnegative. Dividing by \(a\) finishes it. Markov is weak (the tail decays only like \(1/a\)) but it is the engine inside every stronger bound. Each refinement below is Markov applied to a cleverly chosen transformation of \(X\).

Chebyshev's inequality. Apply Markov to the nonnegative variable \((X - \mu)^2\) with \(\mu = \E[X]\) and threshold \(a = \varepsilon^2\).

For a mean of \(n\) i.i.d. variables, \(\Var(\bar X) = \Var(X_1)/n\), so Chebyshev gives a tail decaying like \(1/(n\varepsilon^2)\). This is already enough to prove laws of large numbers, but the polynomial decay is far too slow to survive the union bounds coming later. A union over \(|\mathcal H|\) hypotheses would need the tail to shrink faster than \(1/|\mathcal H|\), and \(1/n\) cannot absorb a large class. Exponential tails can.

The Chernoff method

The trick that upgrades polynomial tails to exponential ones is to apply Markov to \(e^{\lambda X}\) rather than to \(X\) or its square. For any \(\lambda > 0\), the map \(x \mapsto e^{\lambda x}\) is increasing, so the events \(\{X \ge a\}\) and \(\{e^{\lambda X} \ge e^{\lambda a}\}\) are identical, and Markov gives

Since the bound holds for every \(\lambda\), take the best one, \(\P(X \ge a) \le \inf_{\lambda > 0} e^{-\lambda a} \E[e^{\lambda X}]\). Two things make this powerful. First, the moment generating function \(\E[e^{\lambda X}]\) encodes all moments at once, so optimizing \(\lambda\) extracts the strongest polynomial-moment bound simultaneously. Second, for a sum \(S_n = \sum_i X_i\) of independent variables the MGF factorizes, \(\E[e^{\lambda S_n}] = \prod_i \E[e^{\lambda X_i}]\), which turns a bound on one variable into an exponentially compounding bound on the sum. Independence enters the theory exactly here and nowhere else.

Hoeffding's lemma, derived in full

To run the Chernoff method on bounded variables, a bound on the MGF is needed. Hoeffding's lemma. If \(X \in [a, b]\) almost surely and \(\E[X] = 0\), then for every \(\lambda \in \R\),

Define the cumulant generating function \(\psi(\lambda) = \ln \E[e^{\lambda X}]\). Differentiating under the expectation,

These two expressions have a probabilistic reading that does all the work. Define the tilted distribution \(\P_\lambda\) with density proportional to \(e^{\lambda x}\) against the law of \(X\). Then \(\psi'(\lambda) = \E_\lambda[X]\) and \(\psi''(\lambda) = \E_\lambda[X^2] - (\E_\lambda[X])^2 = \Var_\lambda(X)\), so the second derivative of the cumulant function is the variance of \(X\) under the tilted measure. Tilting reweights but cannot move mass outside \([a,b]\), so \(X\) is still supported there under \(\P_\lambda\), and the variance of any variable supported on an interval is at most the squared half-width. Writing \(m = (a+b)/2\) for the midpoint,

since \(|X - m| \le (b-a)/2\) pointwise. So \(\psi''(\lambda) \le (b-a)^2/4\) for all \(\lambda\). Now integrate twice. \(\psi(0) = \ln 1 = 0\), and \(\psi'(0) = \E[X] = 0\) by assumption, so Taylor's theorem with the integral remainder gives, for some intermediate \(\xi\),

Exponentiating both sides is the lemma. Nothing was skipped. The only inputs were boundedness, the zero mean, and the interval variance bound, and the constant \(1/8\) is exactly the \(1/2\) from Taylor times the \(1/4\) from the half-width.

Hoeffding's inequality

Let \(X_1, \dots, X_n\) be independent with \(X_i \in [a_i, b_i]\), and let \(\bar X = \frac1n \sum_i X_i\), \(\mu = \E[\bar X]\). Run the Chernoff method on \(\bar X - \mu = \frac1n \sum_i (X_i - \E X_i)\). Each centered term is mean-zero and supported on an interval of width \(b_i - a_i\), so by independence and the lemma,

Chernoff with threshold \(n\varepsilon\) then gives \(\P(\bar X - \mu \ge \varepsilon) \le \exp(-\lambda n \varepsilon + \frac{\lambda^2}{8}\sum_i (b_i-a_i)^2)\). The exponent is a parabola in \(\lambda\). Setting its derivative to zero, \(-n\varepsilon + \frac{\lambda}{4}\sum_i (b_i - a_i)^2 = 0\), gives \(\lambda^\ast = 4 n \varepsilon / \sum_i (b_i-a_i)^2\), and substituting back,

The same argument applied to \(-X_i\) bounds the lower tail, and a union of the two events gives the two-sided form \(\P(|\bar X - \mu| \ge \varepsilon) \le 2 e^{-2n\varepsilon^2}\) for variables in \([0,1]\). This is the workhorse of the entire subject. The empirical error of one fixed classifier on \(n\) i.i.d. samples is a mean of \([0,1]\) variables, so it sits within \(\varepsilon\) of the true error with probability \(1 - 2e^{-2n\varepsilon^2}\). Note what it does not say. It promises nothing yet about a classifier chosen by looking at the data. That gap is where the rest of the page lives.

Uniform convergence and PAC learning

Why one hypothesis is not enough

The setting for everything that follows is an unknown distribution \(\D\) over examples \(z = (x, y)\), a hypothesis class \(\mathcal H\), a loss \(\ell(h, z) \in [0,1]\) (0-1 loss unless said otherwise), true risk \(L(h) = \E_{z \sim \D}[\ell(h,z)]\), and empirical risk \(\hat L_S(h) = \frac1n \sum_{i=1}^n \ell(h, z_i)\) on a sample \(S\) of \(n\) i.i.d. draws. The learner picks \(\hat h = \argmin_{h \in \mathcal H} \hat L_S(h)\), empirical risk minimization (ERM).

Hoeffding controls \(|\hat L_S(h) - L(h)|\) for one \(h\) fixed before the sample is drawn. But \(\hat h\) is chosen because it looks good on \(S\), which is exactly the conditioning Hoeffding forbids. Among many hypotheses, some will look good by luck, and ERM is a magnet for them. The fix is to demand that the empirical risk be close to the true risk for every hypothesis simultaneously, an event called uniform convergence. If \(\sup_{h \in \mathcal H} |\hat L_S(h) - L(h)| \le \varepsilon\), then ERM is automatically sound. For any competitor \(h^\ast\),

using uniform convergence, then the definition of ERM, then uniform convergence again. So a uniform deviation of \(\varepsilon\) costs at most \(2\varepsilon\) of excess risk, and the whole problem reduces to bounding the probability that the supremum exceeds \(\varepsilon\).

Finite classes and the union bound

For a finite class, the supremum is a maximum over \(|\mathcal H|\) events, and the union bound (\(\P(\cup_j A_j) \le \sum_j \P(A_j)\), immediate by induction from \(\P(A \cup B) = \P(A) + \P(B) - \P(A \cap B) \le \P(A) + \P(B)\)) converts the one-hypothesis Hoeffding bound into a uniform one.

Setting the right side to \(\delta\) and solving for \(\varepsilon\) shows that with probability at least \(1 - \delta\), every \(h \in \mathcal H\) satisfies

The class size enters only through its logarithm, which is why finite-class bounds are surprisingly usable. Concretely, take \(|\mathcal H| = 1000\) models evaluated on \(n = 5000\) held-out examples with tolerance \(\varepsilon = 0.03\). The failure probability is bounded by \(2 \cdot 1000 \cdot e^{-2 \cdot 5000 \cdot 0.0009} = 2000\, e^{-9} = 0.24682\), too weak to be comfortable, but asking instead for the \(\varepsilon\) achievable at \(\delta = 0.01\) gives \(\varepsilon = \sqrt{\ln(200000)/10000} = 0.034937\). A thousand-way model selection on five thousand held-out points costs about 3.5 accuracy points of guarantee, and because the dependence is \(\sqrt{\ln |\mathcal H|}\), evaluating a million models instead of a thousand only moves the guarantee to \(\sqrt{\ln(2 \times 10^8)/10^4} = 0.0437\). These two numbers are the honest theory behind every leaderboard and every hyperparameter sweep.

PAC learning

Valiant (1984) packaged this into a definition that still frames the field. A class \(\mathcal H\) is PAC learnable if there is an algorithm \(A\) and a sample complexity function \(n(\varepsilon, \delta)\) such that for every distribution \(\D\) and every \(\varepsilon, \delta \in (0,1)\), given \(n \ge n(\varepsilon,\delta)\) i.i.d. examples, with probability at least \(1-\delta\) the output satisfies \(L(A(S)) \le \min_{h \in \mathcal H} L(h) + \varepsilon\). "Probably" is the \(\delta\), "approximately correct" is the \(\varepsilon\), and the quantifier order matters. One algorithm must work for all distributions. The realizable case assumes some \(h^\star \in \mathcal H\) has \(L(h^\star) = 0\). The agnostic case (Haussler, and Kearns, Schapire, Sellie) drops that assumption and is the definition given above.

For finite classes the two regimes have different rates, and the realizable one is worth deriving because the argument is different and sharper. Suppose \(L(h^\star) = 0\) and ERM returns some \(\hat h\) with \(\hat L_S(\hat h) = 0\) (it can, since \(h^\star\) achieves zero). Call \(h\) "bad" if \(L(h) > \varepsilon\). A single bad hypothesis survives (labels all \(n\) samples correctly) with probability at most \((1-\varepsilon)^n \le e^{-\varepsilon n}\), using \(1 - u \le e^{-u}\). By the union bound over at most \(|\mathcal H|\) bad hypotheses, the probability that any bad hypothesis survives is at most \(|\mathcal H| e^{-\varepsilon n}\), and setting this to \(\delta\) gives

The realizable rate is \(1/\varepsilon\), the agnostic rate \(1/\varepsilon^2\). Knowing that zero training error is achievable is worth a full factor of \(1/\varepsilon\) in sample complexity, because in the realizable case the bound only needs one-sided protection against impostors rather than two-sided estimation of every risk. This \(1/\varepsilon\) versus \(1/\varepsilon^2\) gap reappears in every framework on this page, including the optimistic rates of modern interpolation theory.

VC dimension and the growth function

Shattering, worked for concrete classes

Infinite classes break the union bound, but they do not break the underlying idea, because what matters is not how many hypotheses exist but how many distinct behaviors they can produce on \(n\) points. For binary classifiers, the restriction of \(\mathcal H\) to a point set \(C = \{x_1, \dots, x_n\}\) is the set of label vectors \(\mathcal H_C = \{(h(x_1), \dots, h(x_n)) : h \in \mathcal H\} \subseteq \{-1,+1\}^n\), the growth function is \(\Pi_{\mathcal H}(n) = \max_{|C| = n} |\mathcal H_C|\), and \(C\) is shattered when \(|\mathcal H_C| = 2^n\), i.e. every labeling of \(C\) is realized by some hypothesis. The VC dimension (Vapnik and Chervonenkis, 1971) is the size of the largest shatterable set, \(d_{\mathrm{VC}}(\mathcal H) = \max\{n : \Pi_{\mathcal H}(n) = 2^n\}\). Establishing \(d_{\mathrm{VC}} = d\) always takes two arguments, exhibiting one set of size \(d\) that is shattered and proving no set of size \(d+1\) can be.

Intervals on the line. \(\mathcal H = \{ h_{[a,b]} : h_{[a,b]}(x) = +1 \text{ iff } x \in [a,b] \}\). Two points \(x_1 < x_2\) are shattered. \([x_1 - 1, x_1]\) gives \((+,-)\), \([x_2, x_2+1]\) gives \((-,+)\), \([x_1, x_2]\) gives \((+,+)\), and any interval to the left of \(x_1\) gives \((-,-)\). Three points \(x_1 < x_2 < x_3\) cannot be. The labeling \((+, -, +)\) requires an interval containing \(x_1\) and \(x_3\) but not \(x_2\), and an interval is convex, so containing the endpoints forces containing the middle. Hence \(d_{\mathrm{VC}} = 2\), matching the parameter count, two endpoints.

Halfspaces in the plane. \(\mathcal H = \{x \mapsto \mathrm{sign}(w\T x + b)\}\) on \(\R^2\). Three points not on a common line are shattered. For each of the 8 labelings, either all labels agree (push the line past the triangle) or the labeling splits the points 2-1, and a line can always separate one vertex of a triangle from the other two. Four points can never be shattered, by a two-case argument. If the four points are in convex position, name them around the hull \(p_1, p_2, p_3, p_4\) and try the alternating labeling \((+,-,+,-)\). A separating halfplane would put \(p_1, p_3\) on one side and \(p_2, p_4\) on the other, but the segments \(p_1 p_3\) and \(p_2 p_4\) are the diagonals of a convex quadrilateral and intersect at some point \(q\). Halfplanes are convex, so \(q\) would lie in both the positive halfplane and its complement, a contradiction. If instead one point lies inside the triangle of the other three, label the triangle \(+\) and the inner point \(-\). The positive halfplane would contain the three vertices, hence their convex hull, hence the inner point, again a contradiction. So \(d_{\mathrm{VC}} = 3\). The same argument through Radon's theorem (any \(d+2\) points in \(\R^d\) split into two disjoint subsets whose convex hulls intersect) gives \(d_{\mathrm{VC}} = d+1\) for halfspaces in \(\R^d\). Label one side of the Radon partition \(+\) and the other \(-\), and the intersecting hulls produce the same contradiction, while shattering \(d+1\) affinely independent points is a direct construction. Parameter count \(d+1\) again, but the match is a coincidence of nice classes. The single-parameter family \(x \mapsto \mathrm{sign}(\sin(\omega x))\) has infinite VC dimension, because for the points \(x_i = 2^{-i}\) a frequency \(\omega\) can be chosen to realize any labeling, encoding the labels in the binary expansion of \(\omega/\pi\). VC dimension measures combinatorial richness, not parameters.

The Sauer-Shelah lemma

The growth function of any class is either \(2^n\) forever (infinite VC dimension) or, once \(n\) passes \(d = d_{\mathrm{VC}}\), polynomial. This dichotomy is the Sauer-Shelah lemma (Sauer 1972, and independently Shelah 1972 in a model-theoretic setting, while Vapnik and Chervonenkis proved a version in 1971).

Proof. By induction on \(n + d\). Fix a set \(S = \{x_1, \dots, x_n\}\) and write \(H = \mathcal H_S\) for the set of realized label vectors. The claim is \(|H| \le \Phi_d(n)\) whenever no subset of \(S\) of size \(d+1\) is shattered by \(H\). Remove the last point. Let \(S' = S \setminus \{x_n\}\), let \(H'\) be the restriction of \(H\) to \(S'\) (project each vector onto its first \(n-1\) coordinates and deduplicate), and let

the patterns on \(S'\) that appear with both labels on \(x_n\). Every vector of \(H\) maps to its projection in \(H'\). The map is 2-to-1 exactly on \(H''\) and 1-to-1 elsewhere, so \(|H| = |H'| + |H''|\), an exact count, not a bound. Now the two key observations. First, \(H'\) is a class on \(n-1\) points with VC dimension at most \(d\) (a shattered set for \(H'\) is shattered by \(H\)), so by induction \(|H'| \le \Phi_d(n-1)\). Second, and this is the step the whole proof turns on, \(H''\) has VC dimension at most \(d - 1\). Suppose \(H''\) shattered some \(T \subseteq S'\) with \(|T| = d\). Every pattern on \(T\) extends to a member of \(H''\), and every member of \(H''\) carries both labels on \(x_n\). Therefore every pattern on \(T \cup \{x_n\}\), meaning any labeling of \(T\) with either label on \(x_n\), is realized by \(H\), so \(H\) shatters the \((d+1)\)-set \(T \cup \{x_n\}\), contradicting \(d_{\mathrm{VC}}(H) \le d\). By induction \(|H''| \le \Phi_{d-1}(n-1)\), and Pascal's rule finishes it.

The base cases are immediate. \(d = 0\) forces \(|H| = 1\) (any two distinct vectors differ somewhere, shattering a single point), and \(n = 0\) gives \(|H| = 1 = \Phi_d(0)\). \(\blacksquare\)

The polynomial form follows from a clean generating-function trick. For \(n \ge d\), multiply and divide by \((d/n)^d\).

where the first inequality uses \((d/n)^{d-i} \le 1\) for \(i \le d\), the binomial theorem gives the equality, and \(1 + u \le e^u\) gives the last step. Hence \(\Phi_d(n) \le (en/d)^d\). Numerically, for \(d = 3\) and \(n = 10\), \(\Phi_3(10) = 1 + 10 + 45 + 120 = 176\), versus \(2^{10} = 1024\) labelings unrestricted and the polynomial cap \((10e/3)^3 = 743.9\). A class either shatters everything or is exponentially thin on large samples. There is no middle ground.

The VC generalization bound

The growth function replaces class size in the uniform bound, through the double-sample symmetrization argument of Vapnik and Chervonenkis. The result, in one standard form, states that with probability at least \(1 - \delta\), for every \(h \in \mathcal H\),

with constants that vary by derivation. The mechanism deserves to be seen even though the bookkeeping is long, because it is the same mechanism behind Rademacher complexity in the next section. Step one is the ghost sample. Draw a second, independent sample \(S'\) of size \(n\). If some \(h\) has \(L(h) - \hat L_S(h) > \varepsilon\), then since \(\hat L_{S'}(h)\) is an average of \(n\) independent \([0,1]\) variables with mean \(L(h)\), Chebyshev gives \(\P(\hat L_{S'}(h) > L(h) - \varepsilon/2) \ge 1/2\) for \(n \ge 8/\varepsilon^2\), so with probability at least \(1/2\) the ghost sample confirms the deviation, \(\hat L_{S'}(h) - \hat L_S(h) > \varepsilon/2\). Therefore

This is the decisive move. The right-hand event mentions only the \(2n\) points of \(S \cup S'\), so the infinite class collapses to at most \(\Pi_{\mathcal H}(2n)\) distinct behaviors. Step two is random swaps. Conditioned on the multiset \(S \cup S'\), swapping \(z_i \leftrightarrow z_i'\) independently per coordinate leaves the joint distribution invariant, and under random swaps each fixed behavior's deviation is a mean-zero bounded sum, so Hoeffding bounds it by \(e^{-n\varepsilon^2/8}\). Step three union-bounds over the \(\Pi_{\mathcal H}(2n)\) behaviors and applies Sauer-Shelah.

Inverting for \(\varepsilon\) gives the boxed bound. Finite VC dimension therefore implies distribution-free uniform convergence, and the converse holds too (the "fundamental theorem of statistical learning", assembled from Vapnik-Chervonenkis and Blumer, Ehrenfeucht, Haussler, and Warmuth 1989). A binary class is PAC learnable if and only if its VC dimension is finite, with sample complexity \(\Theta((d + \ln(1/\delta))/\varepsilon^2)\) agnostic and \(\Theta((d \ln(1/\varepsilon) + \ln(1/\delta))/\varepsilon)\) realizable.

Rademacher complexity

Definition and what it measures

The VC machinery is distribution-free, charging for the worst point set the class could ever meet. Rademacher complexity (Koltchinskii, then Bartlett and Mendelson 2002) charges only for the sample actually drawn. Let \(\sigma_1, \dots, \sigma_n\) be independent uniform \(\pm 1\) signs ("Rademacher variables"). For a class \(\mathcal F\) of real-valued functions and a sample \(S = (z_1, \dots, z_n)\), the empirical Rademacher complexity is

The interpretation is literal. \(\sigma\) is a vector of pure noise, and \(\hat{\mathfrak R}_S\) asks how well the class can correlate with noise on these points. A class that can chase arbitrary sign patterns has complexity near 1, while a class of nearly-constant functions has complexity near 0. "Can fit random labels" is exactly the property Zhang et al. (2017) demonstrated for deep networks, which is why that experiment is best read as a measurement. The empirical Rademacher complexity of the full network class on CIFAR-10-sized samples is essentially 1, so uniform convergence over the whole class cannot explain generalization there. The quantity is computable, by Monte Carlo over \(\sigma\) and for small cases exactly. Both are done below, and the measured values match the theory to three digits.

The generalization bound via symmetrization

Theorem. For \(\mathcal F\) with values in \([0,1]\), with probability at least \(1-\delta\) over \(S \sim \D^n\), every \(f \in \mathcal F\) satisfies

Applied to \(\mathcal F = \{ z \mapsto \ell(h, z) : h \in \mathcal H \}\) this bounds the generalization gap of every hypothesis by twice the Rademacher complexity of the loss class plus a concentration term. The derivation has two halves.

Half one, symmetrization in expectation. Write \(P f = \E[f(z)]\) and \(\hat P_S f = \frac1n \sum_i f(z_i)\). Introduce a ghost sample \(S' = (z_1', \dots, z_n')\), i.i.d. from \(\D\), independent of \(S\). Since \(P f = \E_{S'}[\hat P_{S'} f]\) for each fixed \(f\),

the inequality because a supremum of expectations is at most the expectation of the supremum (Jensen applied to the convex function \(\sup\)). The bracketed quantity is \(\sup_f \frac1n \sum_i (f(z_i') - f(z_i))\). Now the symmetry step. For each \(i\) independently, swapping \(z_i \leftrightarrow z_i'\) leaves the joint distribution of the \(2n\) variables unchanged (they are i.i.d.), and the swap flips the sign of the \(i\)-th summand. Therefore multiplying the \(i\)-th term by an independent random sign \(\sigma_i\) does not change the distribution of the supremum, and

splitting the supremum of a sum into the sum of suprema and using that \(-\sigma_i\) has the same distribution as \(\sigma_i\). So \(\E_S[\sup_f (Pf - \hat P_S f)] \le 2 \mathfrak R_n(\mathcal F)\), meaning the expected worst-case gap is at most twice the expected ability to fit noise.

Half two, concentration. The random variable \(g(S) = \sup_f (Pf - \hat P_S f)\) satisfies the bounded-differences property. Replacing one sample point moves each \(\hat P_S f\) by at most \(1/n\) (values lie in \([0,1]\)), hence moves the supremum by at most \(1/n\). McDiarmid's inequality states that any function of independent variables with bounded differences \(c_i\) concentrates as \(\P(g - \E g \ge t) \le \exp(-2t^2 / \sum_i c_i^2)\). Its proof is a martingale version of the Chernoff-Hoeffding argument above (Azuma's inequality applied to the Doob martingale of \(g\)) and is developed in full in Boucheron, Lugosi, and Massart's monograph. It is the one component this page cites rather than derives. With \(c_i = 1/n\) the deviation term is \(\sqrt{\ln(1/\delta)/2n}\). Applying McDiarmid once to \(g\) and once to \(\hat{\mathfrak R}_S\) (also \(1/n\)-bounded differences) to replace \(\mathfrak R_n\) by its empirical version, then collecting the two \(\delta/2\) events, produces the stated bound with the constant 3. For binary classes, Rademacher complexity also recovers the VC bound. Massart's lemma below gives \(\hat{\mathfrak R}_S \le \sqrt{2 d \ln(en/d) / n}\) via Sauer-Shelah, so the Rademacher route is never worse than VC up to constants, and on easy data it is much better.

Massart's lemma and exact computation for a finite class

Massart's lemma. Let \(A \subset \R^n\) be finite with \(\lVert a \rVert_2 \le r\) for all \(a \in A\). Then

Proof. The Chernoff device again, now on a maximum. For any \(\lambda > 0\), by Jensen's inequality applied to \(\exp\),

Each term factorizes over coordinates, and each \(\sigma_i a_i\) is mean-zero supported on \([-|a_i|, |a_i|]\), so Hoeffding's lemma gives \(\E e^{\lambda \sigma_i a_i} \le e^{\lambda^2 a_i^2 / 2}\) (width \(2|a_i|\), so \((2|a_i|)^2/8 = a_i^2/2\)), and the product is \(e^{\lambda^2 \lVert a \rVert^2/2} \le e^{\lambda^2 r^2 / 2}\). Taking logs,

minimized at \(\lambda = \sqrt{2 \ln|A|}/r\) where the two terms are equal, each \(r\sqrt{2\ln|A|}/2\). Dividing by \(n\) proves the lemma. For a finite class of binary hypotheses, the behavior vectors have \(\lVert a \rVert = \sqrt n\), giving \(\hat{\mathfrak R}_S \le \sqrt{2 \ln |\mathcal H| / n}\), the finite-class bound again but through a completely different door.

These quantities are small enough to compute exactly, which was done for this page (code in the implementation section, all linear algebra in float64 on the GPU). For \(k = 5\) fixed random sign vectors on \(n = 8\) points, enumerating all \(2^8 = 256\) sign assignments gives \(\hat{\mathfrak R}_S = 0.401367\) exactly. Monte Carlo with 1,000 draws gives \(0.400500\) and with 100,000 draws \(0.401973\). Massart's bound evaluates to \(\sqrt{2 \ln 5 / 8} = 0.634318\), loose by a factor of 1.58, which is typical, since Massart prices the worst geometry of 5 vectors and these are random. For the class of one-sided threshold functions restricted to 16 points (17 distinct behaviors by the growth-function count, since \(d_{\mathrm{VC}} = 1\) and Sauer gives \(n + 1\)), exact enumeration of all \(2^{16} = 65{,}536\) sign vectors gives \(\hat{\mathfrak R}_S = 0.342535\), Monte Carlo with 200,000 draws gives \(0.342498\), and Massart's bound is \(\sqrt{2 \ln 17/16} = 0.595106\). A class with 17 behaviors on 16 points can achieve about a third of perfect noise correlation. That measured number, not the class's cardinality, is what enters the generalization bound.

Two structural tools: linear classes and contraction

Two more facts convert Rademacher complexity from a definition into a calculus. First, the complexity of norm-bounded linear functions has a closed-form bound. For \(\mathcal F_B = \{ x \mapsto w\T x : \lVert w \rVert_2 \le B \}\),

where \(R_x = \max_i \lVert x_i \rVert\). The first equality is Cauchy-Schwarz made tight (choose \(w\) parallel to the sum), the middle inequality is Jensen, and the last equality expands the square, where cross terms \(\E[\sigma_i \sigma_j] = 0\) vanish for \(i \ne j\), leaving \(\sum_i \lVert x_i \rVert^2\). No dimension appears anywhere. Only norms do. This single line is the reason SVMs, boosting, and (through the NTK) wide networks can generalize in feature spaces of unbounded dimension.

Second, Talagrand's contraction lemma states that if \(\phi\) is \(\rho\)-Lipschitz, then \(\hat{\mathfrak R}_S(\phi \circ \mathcal F) \le \rho\, \hat{\mathfrak R}_S(\mathcal F)\). Composing with a gentle loss cannot amplify the class's ability to fit noise by more than the Lipschitz constant. The proof (an exchange argument over one coordinate at a time, in Ledoux and Talagrand's book and reproduced in Shalev-Shwartz and Ben-David and in Mohri et al.) is cited rather than derived here. Together the two facts produce margin bounds, next.

Margins: why boosting and SVMs generalize

The margin bound

For a real-valued score \(f(x)\) predicting via \(\mathrm{sign}(f(x))\), the margin of an example is \(y f(x)\), positive when correct, and its size measures the confidence of the call. Fix a margin scale \(\theta > 0\) and define the ramp function \(\phi_\theta(u) = 1\) for \(u \le 0\), \(\phi_\theta(u) = 1 - u/\theta\) for \(0 < u < \theta\), and \(0\) for \(u \ge \theta\). Three pointwise facts do the work, \(\mathbf 1\{u \le 0\} \le \phi_\theta(u) \le \mathbf 1\{u \le \theta\}\), and \(\phi_\theta\) is \((1/\theta)\)-Lipschitz. Then for any class \(\mathcal F\) of scores, chaining the Rademacher bound through the ramp loss gives

with probability \(1 - \delta\), for all \(f \in \mathcal F\) simultaneously. The first inequality is the pointwise domination in expectation, the second is the Rademacher generalization theorem applied to \(\phi_\theta \circ \mathcal F\) (values in \([0,1]\)) followed by contraction with \(\rho = 1/\theta\), and the empirical ramp is bounded by the fraction of training points with margin at most \(\theta\). Reading it off, test error is at most the fraction of small-margin training points plus a complexity term shrunk by the margin scale. For norm-bounded linear scores in a feature space with \(\lVert x \rVert \le R_x\), plugging in the linear-class bound gives

This is the SVM story stated as a theorem. Maximizing the geometric margin minimizes \(B/\theta\) at fixed data radius, and the guarantee depends on dimension not at all, only on how confidently the data is separated relative to its scale. It also disciplines the kernel view from the statistical-learning page. An infinite-dimensional RKHS costs nothing if margins are large, because \(B R_x / \theta\) is the real capacity.

Boosting's margin distribution, derived

AdaBoost (Freund and Schapire 1997) was the puzzle that forced margins into the theory. The algorithm maintains weights \(D_t(i)\) over training examples (initialized uniform), fits a weak learner \(h_t\) with weighted error \(\varepsilon_t = \P_{i \sim D_t}[h_t(x_i) \ne y_i]\), sets \(\alpha_t = \frac12 \ln \frac{1-\varepsilon_t}{\varepsilon_t}\), and updates \(D_{t+1}(i) \propto D_t(i)\, e^{-\alpha_t y_i h_t(x_i)}\) with normalizer \(Z_t\). The observed phenomenon (Schapire, Freund, Bartlett, and Lee 1998, following experiments by Breiman and by Drucker and Cortes) is that training error hits zero within a few rounds, yet test error keeps falling for hundreds more. If capacity-counting were the whole story, adding terms to an already-interpolating ensemble should hurt. The resolution is that the margin distribution keeps improving, and the margin bound above (applied to the convex hull of the weak learners, whose Rademacher complexity is that of the base class, by convexity of the supremum) prices the ensemble by its margins, not its term count.

The empirical half is exactly derivable. Let \(g(x) = \sum_{t=1}^T \alpha_t h_t(x)\) and normalize \(f = g / \sum_t \alpha_t\) so margins lie in \([-1,1]\). For any \(\theta\), start from the indicator bound \(\mathbf 1\{u \le 0\} \le e^{-u}\), applied at \(u = y_i g(x_i) - \theta \sum_t \alpha_t\).

Now the identity that makes AdaBoost analyzable. Unrolling the weight update gives \(D_{T+1}(i) = \frac{1}{n} \cdot \frac{e^{-y_i g(x_i)}}{\prod_t Z_t}\), and since \(D_{T+1}\) sums to 1, \(\frac1n \sum_i e^{-y_i g(x_i)} = \prod_t Z_t\), an exact equality. The normalizer itself evaluates by splitting the expectation over correct and incorrect examples, \(Z_t = (1-\varepsilon_t) e^{-\alpha_t} + \varepsilon_t e^{\alpha_t}\), and with AdaBoost's choice of \(\alpha_t\), \(e^{\alpha_t} = \sqrt{(1-\varepsilon_t)/\varepsilon_t}\), so \(Z_t = 2\sqrt{\varepsilon_t(1 - \varepsilon_t)}\). Combining, with \(e^{\theta \alpha_t} = ((1-\varepsilon_t)/\varepsilon_t)^{\theta/2}\),

At \(\theta = 0\) this is the classic training-error bound \(\prod_t 2\sqrt{\varepsilon_t (1-\varepsilon_t)}\), which is below 1 whenever \(\varepsilon_t \ne 1/2\). The content for \(\theta > 0\) is that if every weak learner beats chance by \(\gamma\) (\(\varepsilon_t \le 1/2 - \gamma\)), each factor is a constant less than 1 for any \(\theta <\) roughly \(\gamma\), so the fraction of training points with margin below \(\theta\) decays exponentially in the number of rounds, long after the fraction below zero has hit exactly zero. Numerically, \(\varepsilon_t = 0.3\) and \(\theta = 0.1\) give a per-round factor of \(2 \cdot 0.3^{0.45} \cdot 0.7^{0.55} = 0.9559\), so 100 further rounds shrink the sub-0.1-margin fraction by \(0.9559^{100} \approx 0.011\). Boosting is a margin-maximization device that never mentions margins.

This page measured the phenomenon directly (full code in the implementation section), running AdaBoost with depth-1 stumps on the UCI breast cancer data (400 training, 169 test points). With the from-scratch implementation, training error reaches zero by round 100, and from round 100 to 1000 the minimum normalized training margin climbs from \(+0.121\) to \(+0.164\) while test error falls from \(4.73\%\) to \(4.14\%\). The scikit-learn implementation on the same split shows the same signature (train error zero from round 100, minimum margin \(0.118 \to 0.155\), test error \(5.33\% \to 3.55\%\)). Nothing about the training labels changes after round 100. What changes is the confidence structure, exactly as the bound predicts. One caveat is owed for honesty. Breiman (1999) constructed ensembles with better minimum margins and worse test error, so minimum margin alone is not the whole distribution. The modern sharpening through Koltchinskii and Panchenko (2002) and the k-th margin bounds of Gao and Zhou keeps the explanation intact with the full margin distribution in place of the minimum.

Algorithmic stability

Bounding the algorithm instead of the class

Uniform convergence bounds the worst hypothesis the class contains, even ones no sane algorithm would return. Stability theory (Bousquet and Elisseeff 2002, with roots in Rogers-Wagner and Devroye-Wagner in the 1970s) bounds only the function the algorithm actually outputs, by asking how much that output depends on any single training example. Let \(A(S)\) denote the hypothesis learned from \(S\), and let \(S^{(i)}\) be \(S\) with example \(z_i\) replaced by an independent copy \(z_i'\). The algorithm has uniform stability \(\beta\) if

Theorem (stability implies generalization, in expectation). \( \E_S\big[ L(A(S)) - \hat L_S(A(S)) \big] \le \beta. \) The proof is a renaming argument with no analysis in it. Write the expected empirical risk as \(\E_S[\hat L_S(A(S))] = \frac1n \sum_i \E_S[\ell(A(S), z_i)]\). For the true risk, draw the ghost point \(z_i'\) and note \(\E_S[L(A(S))] = \E_{S, z_i'}[\ell(A(S), z_i')]\) for any fixed \(i\), since \(z_i'\) is a fresh draw. Now rename. The pair "train on \(S\), test on \(z_i'\)" has the same joint distribution as "train on \(S^{(i)}\), test on \(z_i\)", because both are functions of the same \(n+1\) i.i.d. variables with the roles of \(z_i\) and \(z_i'\) exchanged. Therefore, for each \(i\),

by uniform stability applied at the point \(z = z_i\). Averaging over \(i\) gives the theorem. A high-probability version follows from McDiarmid, since replacing one point moves the generalization gap by at most \(2\beta + 1/n\) when losses are bounded. The Bousquet-Elisseeff form is \(L \le \hat L + \beta + (2n\beta + M)\sqrt{\ln(1/\delta)/2n}\) for losses in \([0, M]\), useful only when \(\beta = o(1/\sqrt n)\). Feldman and Vondrák (2019) tightened the high-probability rate to nearly \(\beta \ln n + \sqrt{1/n}\), one of the quiet sharp results of the modern period.

Regularization buys stability

The canonical stable algorithm is regularized ERM. Let the loss \(\ell(\cdot, z)\) be convex and \(\rho\)-Lipschitz in the parameter \(w\), and let \(A(S) = \argmin_w F_S(w)\) with \(F_S(w) = \hat L_S(w) + \lambda \lVert w \rVert^2\). The regularizer makes \(F_S\) \(2\lambda\)-strongly convex, so for the minimizer \(w_S\) and any \(w\), \(F_S(w) - F_S(w_S) \ge \lambda \lVert w - w_S \rVert^2\) (expand \(F_S(w) \ge F_S(w_S) + \nabla F_S(w_S)\T (w - w_S) + \lambda \lVert w - w_S\rVert^2\) and use \(\nabla F_S(w_S) = 0\)). Write \(u = w_S\), \(v = w_{S^{(i)}}\), \(\Delta = \lVert u - v\rVert\). The empirical objectives differ only in the \(i\)-th term, so

because \(v\) minimizes \(F_{S^{(i)}}\) (making the first difference \(\le 0\)) and each bracketed pair is bounded by \(\rho\Delta\) by Lipschitzness. Combining with strong convexity, \(\lambda \Delta^2 \le F_S(v) - F_S(u) \le 2\rho\Delta/n\), so \(\Delta \le 2\rho / (\lambda n)\), and one more application of Lipschitzness converts parameter distance to loss distance.

Stability decays as \(1/n\), so regularized ERM generalizes at rate \(O(\rho^2/(\lambda n))\) with no mention of dimension, VC dimension, or any property of the hypothesis class. This is the second dimension-free mechanism on this page after margins, and it is the one that survives non-uniform-convergence regimes. Hardt, Recht, and Singer (2016) extended the idea to stochastic gradient descent itself. For convex smooth losses, \(T\) steps of SGD with step sizes \(\eta_t\) are uniformly stable with \(\beta \le \frac{2\rho^2}{n} \sum_t \eta_t\), and for non-convex smooth losses a weaker but nonvacuous bound holds if training is not too long, "train faster, generalize better". Their analysis is a growth-recursion on how far two SGD trajectories, coupled on all but one example, can drift apart. It is the closest thing deep learning has to a classical generalization guarantee that actually references the optimizer.

The implicit bias of gradient descent

Least squares: gradient descent finds the minimum-norm interpolator

When a model has more parameters than constraints, the training loss has many global minima, and "which one does the optimizer pick" becomes the real regularization question. For underdetermined least squares the answer is exact. Minimize \(J(\theta) = \frac12 \lVert X\theta - y\rVert^2\) with \(X \in \R^{n \times p}\), \(p > n\), \(X\) of full row rank, by gradient descent \(\theta_{k+1} = \theta_k - \eta X\T (X\theta_k - y)\) from \(\theta_0 = 0\).

Every update adds a vector of the form \(X\T(\cdot)\), which lies in the row space of \(X\), so by induction every iterate lies in \(\mathrm{row}(X) = \mathrm{span}\{x_1, \dots, x_n\}\). On that \(n\)-dimensional subspace the quadratic is strictly convex (its Hessian \(X\T X\) restricted to the row space has eigenvalues equal to the nonzero eigenvalues of \(X X\T\), all positive), so for \(\eta < 2/\sigma_{\max}^2(X)\) gradient descent converges to the unique minimizer within the row space, which satisfies \(X\theta = y\) exactly since \(p > n\) makes the system solvable. Now decompose any interpolator as \(\theta = \theta_r + v\) with \(\theta_r \in \mathrm{row}(X)\) and \(v \in \mathrm{null}(X)\). Interpolation fixes \(\theta_r\) (two row-space solutions would differ by a row-space null vector, forcing equality), and Pythagoras gives \(\lVert \theta \rVert^2 = \lVert \theta_r \rVert^2 + \lVert v \rVert^2\). The minimum-norm interpolator is therefore exactly the row-space one, \(\theta^\dagger = X\T (XX\T)^{-1} y = X^+ y\), and gradient descent from zero converges to it. The optimizer is a regularizer. Among the infinitely many zero-loss solutions, plain gradient descent selects the one of smallest Euclidean norm, without any penalty term in sight. Every appearance of "min-norm interpolation" in the double-descent section below is implicitly a statement about what gradient descent trained to convergence would find.

Logistic regression on separable data: convergence to the max-margin direction

For classification the analogous result is deeper. On linearly separable data, logistic loss \(\L(w) = \frac1n \sum_i \ln(1 + e^{-y_i w\T x_i})\) has no minimizer, since scaling any separating \(w\) up decreases the loss toward zero, so gradient descent sends \(\lVert w \rVert \to \infty\). The question is the direction. Soudry, Hoffer, Nacson, Gunasekar, and Srebro (2018) proved that for almost every separable dataset, gradient descent with small constant step size satisfies

where \(\hat w\) is the solution of the hard-margin SVM problem \(\min \lVert w \rVert^2\) subject to \(y_i w\T x_i \ge 1\). The mechanism is visible in a derivation-sketch whose every step can be checked, with the full rigor in the paper. Late in training all margins are large, and \(\ln(1 + e^{-u}) = e^{-u}(1 + O(e^{-u}))\), so the loss behaves as the exponential loss \(\sum_i e^{-y_i w\T x_i}\). Posit \(w(t) \approx u \ln t\) for some direction-scale vector \(u\) and ask what makes this self-consistent under gradient flow \(\dot w = \frac1n\sum_i y_i x_i e^{-y_i w\T x_i}\). Substituting, \(e^{-y_i w\T x_i} = t^{-y_i u\T x_i}\), so each example's pull decays polynomially in \(t\) with exponent equal to its margin under \(u\). As \(t\) grows, only the examples of minimal margin \(\min_i y_i u\T x_i\) contribute at leading order. Call that minimal margin 1 (the scale of \(u\) adjusts to make it so, and if it were larger or smaller than 1, \(\dot w\) would decay faster or slower than \(1/t\), contradicting \(\frac{d}{dt}(u\ln t) = u/t\)). Then \(\dot w \approx \frac1t \sum_{i \in \mathrm{SV}} c_i y_i x_i\) with coefficients \(c_i \ge 0\) supported on the minimal-margin set, and consistency demands

These are precisely the KKT conditions of \(\min \frac12\lVert w\rVert^2\) s.t. \(y_i w\T x_i \ge 1\), namely stationarity (\(w = \sum_i \mu_i y_i x_i\)), dual feasibility, primal feasibility, and complementary slackness (only tight constraints carry weight). Gradient descent on plain logistic loss, with no margin term anywhere, converges in direction to the maximum-margin classifier, at the glacially slow rate \(O(1/\ln t)\). Ji and Telgarsky proved matching results independently with different techniques, and Lyu and Li, with parallel work by Chizat and Bach at INRIA/ENS, extended directional convergence to KKT points of the margin problem for deep homogeneous networks, where the result is local rather than global.

The phenomenon is cheap to measure, and was measured for this page (code in the implementation section). On 120 separable 2-D points, the reference max-margin direction was computed by a hard-margin solver (geometric margin \(0.4052\)), and full-batch gradient descent was run on logistic loss from zero with \(\eta = 0.5\). The cosine between \(w(t)/\lVert w(t)\rVert\) and the SVM direction was \(0.986468\) at \(t = 10^2\), \(0.997077\) at \(10^3\), \(0.999857\) at \(10^4\), and \(0.999996\) at \(10^5\), while \(\lVert w \rVert\) grew \(2.72 \to 4.61 \to 8.09 \to 13.07\), the unbounded logarithmic growth the theorem predicts. The unnormalized minimum margin \(\min_i y_i w\T x_i\) grew from \(0.18\) to \(4.86\) over the same horizon. The classifier's decisions froze early, while its confidence and its direction kept sharpening without bound.

Double descent and benign overfitting

The classical U-curve and what replaced it

In the classical picture, as capacity grows, training error falls monotonically while test error traces a U, falling while bias dominates and rising once variance does, with the practitioner's job being to stop at the bottom. In the modern picture, named by Belkin, Hsu, Ma, and Mandal (2019), the U is only the left half. As capacity passes the interpolation threshold, the point where the model family first achieves zero training error, test error spikes, and then falls again, often below the classical sweet spot, as capacity grows further. Nakkiran, Kaplun, Bansal, Yang, Barak, and Sutskever (2019) showed the same shape in deep networks along three different axes (model width, dataset size, and training epochs) and showed label noise amplifies the peak. Neither observation repeals the bias-variance identity, which is exact algebra. What fails is the folklore that capacity is a single knob monotonically trading bias against variance.

The experiment for this page reproduces the full curve in its cleanest known setting, random Fourier features (Rahimi and Recht 2007) with minimum-norm least squares, the miniature of Belkin et al.'s original demonstration. The setup, exactly as run, uses the scikit-learn digits dataset (1,797 8×8 images, 10 classes) with 1,000 training and 797 test points. Features are \(\phi_j(x) = \sqrt{2/N}\cos(w_j\T x + b_j)\) with \(w_j \sim \N(0, 2\gamma I)\), \(b_j\) uniform, and \(\gamma = 0.1089\) from the median heuristic (median squared pairwise distance 9.184). Targets are one-hot, minimum-norm least squares is solved via QR in float64, predictions are by argmax, and five independent feature draws are averaged per \(N\). The full sweep runs in 5.4 s on the NVIDIA H100 80GB in this machine. The measured results are below.

Every signature is present in the measurements. The classical regime bottoms out at \(N = 600\) with \(1.98\%\) test error. At the interpolation threshold \(N = n = 1000\) test error jumps to \(80.4\%\), worse than random guessing among 10 classes, and the overparameterized regime then descends to \(0.90\%\) at \(N = 12{,}000\), less than half the best underparameterized error. The fourth column explains the peak. The norm of the min-norm solution spikes by a factor of \(\sim\)45 at the threshold and then falls as \(N\) grows. At \(N = n\) the feature matrix is square and its smallest singular values are typically tiny, so interpolation requires enormous coefficients that amplify label noise badly. Past the threshold there are many interpolators, the minimum-norm one can spread the fit across abundant directions, and (this is the implicit-bias section paying rent) gradient descent from zero would find exactly that solution. More features here means less effective capacity in the sense that matters, because the selected solution is smoother.

Benign overfitting: the precise conditions

Bartlett, Long, Lugosi, and Tsigler (2020) made "when does interpolating noise not hurt" a theorem for linear regression. The setting is \(y = x\T \theta^\ast + \varepsilon\), with noise variance \(\sigma^2\), covariates with covariance \(\Sigma\) (eigenvalues \(\lambda_1 \ge \lambda_2 \ge \cdots\)), a well-specified subgaussian model, and the minimum-norm interpolator \(\hat\theta = X^+ y\) in the overparameterized regime. The excess risk is controlled by two effective ranks of the covariance tail,

for a universal constant \(b\). Their theorem says there are constants \(b, c\) such that if \(k^\ast < n/c\), then with probability at least \(1 - \delta\),

together with a matching lower bound \(\E[R(\hat\theta)] - R^\ast \ge \frac{\sigma^2}{c}\big( \frac{k^\ast}{n} + \frac{n}{R_{k^\ast}(\Sigma)} \big)\), so the variance characterization is tight. (The proof, a careful analysis of the eigenstructure of random Gram matrices, is in the cited paper. It is genuinely long and is not reproduced here.) Reading the conditions, the risk of the interpolator converges to the noise floor exactly when \(r_0(\Sigma) = o(n)\) (the bias term, total spectral mass concentrated enough that the signal directions are learnable), \(k^\ast = o(n)\) (only a vanishing fraction of directions are "important"), and \(R_{k^\ast}(\Sigma) = \omega(n)\) (the tail spreads its mass over many more than \(n\) comparable directions). The interpretation is the honest headline of the modern theory. Benign overfitting requires significant overparameterization with a long, flat spectral tail, because the tail directions act as a shock absorber, soaking up the noise-fitting at a cost of \(n / R_{k^\ast}\) that vanishes only when the number of unimportant-but-live directions dwarfs the sample size. Concrete spectra make it vivid. \(\lambda_k = k^{-1}\ln^{-\beta}(k+1)\) in infinite dimension is benign precisely when \(\beta > 1\) (summable, but barely, since slower decay breaks \(r_0 = o(n)\) and faster decay like \(\lambda_k = 2^{-k}\) breaks \(R_{k^\ast} = \omega(n)\) because noise concentrates on a few directions), and fixed-dimension isotropic covariates are never benign, since the tail cannot outgrow \(n\). Tsigler and Bartlett (2023) extended the analysis to ridge regression, showing negative or zero regularization can even be optimal under these spectra, and Hastie, Montanari, Rosset, and Tibshirani (2022) computed the exact asymptotic risk curves under proportional scaling, the theoretical version of the measured table above.

The neural tangent kernel, briefly

The bridge from linear theory to networks is linearization. Around initialization \(\theta_0\), a network \(f(x; \theta)\) expands as \(f(x;\theta) \approx f(x;\theta_0) + \nabla_\theta f(x;\theta_0)\T (\theta - \theta_0)\), a linear model in the fixed feature map \(x \mapsto \nabla_\theta f(x;\theta_0)\), whose kernel \(\Theta(x, x') = \nabla_\theta f(x;\theta_0)\T \nabla_\theta f(x';\theta_0)\) is the neural tangent kernel. Jacot, Gabriel, and Hongler (2018, EPFL) proved that as width goes to infinity (with the right scaling), the NTK converges to a deterministic kernel at initialization and stays constant during training, so gradient flow on square loss obeys the linear ODE \(\frac{d}{dt} f_t = -\Theta\, (f_t - y)\) on the training outputs, i.e. infinitely wide networks train as kernel regression with the NTK. The complementary Google Brain line (Lee, Xiao, Schoenholz, Bahri, Novak, Sohl-Dickstein, Pennington 2019) verified that wide-but-finite networks track their linearization closely. Two consequences matter for this page. Generalization of very wide networks inherits the entire kernel/minimum-norm theory above, including benign overfitting analyses transplanted to NTK spectra, and the regime where deep learning beats kernels must involve feature learning, the movement of the kernel during training, which the infinite-width limit freezes out (the "lazy training" critique of Chizat, Oyallon, and Bach). The google/neural-tangents library computes exact infinite-width NTKs for concrete architectures.

PAC-Bayes bounds

The change-of-measure engine

PAC-Bayes (McAllester 1999, building on Shawe-Taylor and Williamson) bounds the risk of a distribution \(Q\) over hypotheses, a "posterior" chosen after seeing data, against a "prior" \(P\) fixed before. Its engine is the Donsker-Varadhan variational formula, derived in three lines. For any measurable \(\phi : \mathcal H \to \R\) and any \(Q\) absolutely continuous with respect to \(P\),

where the first equality just unpacks \(\KL(Q\|P) = \E_Q[\ln \frac{dQ}{dP}]\), the inequality is Jensen for the concave \(\ln\), and the last step is the change of measure \(\E_Q[g \frac{dP}{dQ}] = \E_P[g]\). Rearranged, for every posterior \(Q\) simultaneously,

This is a uniform statement over all \(Q\) obtained with no union bound at all. The \(\KL\) term is the price of choosing \(Q\) after the fact. To make it a generalization bound, choose \(\phi\) to measure the empirical-true gap and control the exponential moment. Take \(\phi(h) = n \cdot \mathrm{kl}(\hat L_S(h) \,\|\, L(h))\), where \(\mathrm{kl}(p\|q)\) is the KL divergence between Bernoulli distributions. Maurer's lemma (2004, a sharp binomial computation, cited not derived) gives \(\E_S \E_P [e^{\phi}] = \E_P \E_S [e^{\phi}] \le 2\sqrt n\) for losses in \([0,1]\). By Markov, with probability at least \(1 - \delta\) over \(S\), \(\E_P[e^{\phi}] \le 2\sqrt n / \delta\). On that event, applying Donsker-Varadhan and then Jensen once more (\(\mathrm{kl}\) is jointly convex, so \(\E_Q\,\mathrm{kl}(\hat L(h)\|L(h)) \ge \mathrm{kl}(\E_Q \hat L \| \E_Q L)\)) gives

Pinsker's inequality \(\,|p - q| \le \sqrt{\mathrm{kl}(p\|q)/2}\,\) (provable by integrating the second derivative of \(\mathrm{kl}\), which is \(\ge 4\) on \([0,1]\)) converts this to the McAllester-style square-root form \(\E_Q L \le \E_Q \hat L + \sqrt{ (\KL(Q\|P) + \ln(2\sqrt n/\delta)) / 2n }\). Three properties distinguish PAC-Bayes from everything else on this page. The complexity term is an algorithm-and-data-dependent divergence rather than a class-level capacity. The bound holds uniformly over posteriors, so \(Q\) may be optimized against the bound itself. And randomized predictors are first-class, which matches SGD's noisy output better than worst-case analysis does.

That last point became an event in 2017, when Dziugaite and Roy took a stochastic neural network on a binarized MNIST task, set \(Q\) to a Gaussian around the SGD solution, and directly optimized the PAC-Bayes bound over \(Q\)'s mean and per-parameter variances, obtaining a certified test-error bound around \(0.16\) for a network whose test error was a few percent, the first nonvacuous generalization certificate for a modern (if small) network, where VC-style counting gives bounds in the millions. The line continued through Pérez-Ortiz, Rivasplata, Shawe-Taylor, and Szepesvári (tighter training-set-split priors) and Lotfi, Finzi, Kapoor, Potapczynski, Goldblum, and Wilson (2022), whose compression-based PAC-Bayes bounds are among the tightest known for image classifiers. PAC-Bayes also quietly connects back to two earlier sections. Choosing \(Q\) and \(P\) as point masses smoothed by noise recovers margin bounds (the Langford-Shawe-Taylor route to SVM bounds), and the stability of the SGD posterior is what keeps \(\KL(Q\|P)\) small, so the three modern frameworks are different renderings of the same underlying insensitivity.

Worked problems

Implementation

Four experiments, each small enough to rerun in seconds, each measuring a theorem from this page. All quoted numbers were produced by these programs on this machine (float64 throughout, with the linear algebra executed on an NVIDIA H100 80GB, where the full double-descent sweep takes 5.4 s wall clock, but nothing here needs more than a laptop CPU). One engineering note from actually running them. At the interpolation threshold the feature matrix is nearly singular by construction, so the min-norm solve must be done with a rank-revealing or orthogonal factorization (QR or SVD, as below), never by forming \(F\T F\), which squares the condition number exactly where the experiment is most sensitive.

Double descent on random Fourier features

The experiment behind the table above is RFF approximation of the RBF kernel on the digits dataset, with minimum-norm least squares on one-hot targets, sweeping the feature count \(N\) through the interpolation threshold \(N = n = 1000\). The three tabs are equivalent. The PyTorch version is the one whose outputs are quoted.

import numpy as np
import torch
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

dev, dt = "cuda", torch.float64

def minnorm_lstsq(F, Y):
    n, N = F.shape                        # F: (n, N), Y: (n, 10)
    if n >= N:                            # overdetermined: unique LS solution
        Q, R = torch.linalg.qr(F)         # Q: (n, N), R: (N, N)
        return torch.linalg.solve_triangular(R, Q.T @ Y, upper=True)
    Q, R = torch.linalg.qr(F.T)           # wide case: F = R^T Q^T
    return Q @ torch.linalg.solve_triangular(R.T, Y, upper=False)
    # returns the MIN-NORM interpolator, i.e. what GD from 0 converges to

X, y = load_digits(return_X_y=True)
X = X / 16.0                              # (1797, 64) in [0, 1]
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=1000,
                                      random_state=0, stratify=y)
Xtr_t = torch.tensor(Xtr, device=dev, dtype=dt)
Xte_t = torch.tensor(Xte, device=dev, dtype=dt)
Ytr = torch.eye(10, device=dev, dtype=dt)[torch.tensor(ytr, device=dev)]
ytr_t, yte_t = torch.tensor(ytr, device=dev), torch.tensor(yte, device=dev)

# median heuristic for the RBF bandwidth exp(-gamma ||x - x'||^2)
idx = np.random.default_rng(0).choice(1000, 200, replace=False)
D2 = ((Xtr[idx, None, :] - Xtr[None, idx, :]) ** 2).sum(-1)
gamma = 1.0 / np.median(D2[D2 > 0])       # measured: 0.1089

for N in (10, 100, 300, 600, 900, 1000, 1100, 2000, 6000, 12000, 24000):
    te = []
    for r in range(5):                    # average 5 feature draws
        g = torch.Generator(device=dev).manual_seed(7000 + 13 * N + r)
        W = torch.randn(64, N, generator=g, device=dev, dtype=dt) \
            * np.sqrt(2 * gamma)          # (64, N)
        b = torch.rand(N, generator=g, device=dev, dtype=dt) * 2 * np.pi
        Ftr = np.sqrt(2.0 / N) * torch.cos(Xtr_t @ W + b)   # (1000, N)
        Fte = np.sqrt(2.0 / N) * torch.cos(Xte_t @ W + b)   # (797, N)
        Th = minnorm_lstsq(Ftr, Ytr)                        # (N, 10)
        te.append(((Fte @ Th).argmax(1) != yte_t).double().mean().item())
    print(f"N={N:>6}  test01={np.mean(te):.4f}")

import numpy as np
import jax
import jax.numpy as jnp
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

jax.config.update("jax_enable_x64", True)

X, y = load_digits(return_X_y=True)
X = X / 16.0                              # (1797, 64)
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=1000,
                                      random_state=0, stratify=y)
Ytr = np.eye(10)[ytr]                     # (1000, 10) one-hot targets

idx = np.random.default_rng(0).choice(1000, 200, replace=False)
D2 = ((Xtr[idx, None, :] - Xtr[None, idx, :]) ** 2).sum(-1)
gamma = 1.0 / np.median(D2[D2 > 0])

for N in (10, 100, 300, 600, 900, 1000, 1100, 2000, 6000, 12000, 24000):
    te = []
    for r in range(5):
        rl = np.random.default_rng(7000 + 13 * N + r)
        W = rl.normal(0.0, np.sqrt(2 * gamma), size=(64, N))
        b = rl.uniform(0, 2 * np.pi, size=N)
        Ftr = jnp.sqrt(2.0 / N) * jnp.cos(jnp.asarray(Xtr) @ W + b)
        Fte = jnp.sqrt(2.0 / N) * jnp.cos(jnp.asarray(Xte) @ W + b)
        # SVD-based lstsq returns the min-norm solution in the wide case
        Th, *_ = jnp.linalg.lstsq(Ftr, jnp.asarray(Ytr), rcond=None)
        te.append(float(jnp.mean(jnp.argmax(Fte @ Th, 1)
                                 != jnp.asarray(yte))))
    print(f"N={N:>6}  test01={np.mean(te):.4f}")

import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

X, y = load_digits(return_X_y=True)
X = X / 16.0
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=1000,
                                      random_state=0, stratify=y)
Ytr = np.eye(10)[ytr]                     # (1000, 10)

idx = np.random.default_rng(0).choice(1000, 200, replace=False)
D2 = ((Xtr[idx, None, :] - Xtr[None, idx, :]) ** 2).sum(-1)
gamma = 1.0 / np.median(D2[D2 > 0])

for N in (10, 100, 300, 600, 900, 1000, 1100, 2000, 6000, 12000, 24000):
    te = []
    for r in range(5):
        rl = np.random.default_rng(7000 + 13 * N + r)
        W = rl.normal(0.0, np.sqrt(2 * gamma), size=(64, N))
        b = rl.uniform(0, 2 * np.pi, size=N)
        Ftr = np.sqrt(2.0 / N) * np.cos(Xtr @ W + b)   # (1000, N)
        Fte = np.sqrt(2.0 / N) * np.cos(Xte @ W + b)   # (797, N)
        # lstsq (SVD) returns the min-norm solution when N > n
        Th, *_ = np.linalg.lstsq(Ftr, Ytr, rcond=None)
        te.append(np.mean((Fte @ Th).argmax(1) != yte))
    print(f"N={N:>6}  test01={np.mean(te):.4f}")

The measured output of the PyTorch version (5-draw means, quoted in full in the table above) is test error \(0.0198\) at \(N=600\), \(0.8040\) at the threshold \(N=1000\), and \(0.0090\) at \(N=12{,}000\). The JAX version on the same seeds reproduces the shape (\(0.0209 / 0.7771 / 0.0100\) with 3 draws). The single most instructive variable to watch is the solution norm, which peaks at \(1685\) at the threshold and settles to \(36\) in the heavily overparameterized regime, the entire double-descent mechanism in one number.

Monte Carlo Rademacher complexity

Estimating \(\hat{\mathfrak R}_S\) is a two-line Monte Carlo once the class's behavior vectors on the sample are enumerated. For small \(n\) the expectation over \(\sigma \in \{\pm1\}^n\) can be computed exactly, which is how the numbers in the Rademacher section were checked.

import numpy as np
import torch

dev, dt = "cuda", torch.float64
rng = np.random.default_rng(42)

n, k = 8, 5
H = torch.tensor(rng.choice([-1.0, 1.0], size=(k, n)),
                 device=dev, dtype=dt)          # (k, n) behavior vectors

# exact: enumerate all 2^n sign vectors
sig = np.array([[(i >> j) & 1 for j in range(n)] for i in range(2 ** n)])
sig = torch.tensor(2.0 * sig - 1.0, device=dev, dtype=dt)   # (256, 8)
exact = ((sig @ H.T).max(dim=1).values.mean() / n).item()

# Monte Carlo with m draws of sigma
m = 100_000
s = torch.tensor(rng.choice([-1.0, 1.0], size=(m, n)),
                 device=dev, dtype=dt)          # (m, n)
mc = ((s @ H.T).max(dim=1).values.mean() / n).item()

massart = float(np.sqrt(2 * np.log(k) / n))
print(f"exact={exact:.6f}  mc={mc:.6f}  massart={massart:.6f}")
# measured: exact=0.401367  mc=0.401973  massart=0.634318

import numpy as np
import jax
import jax.numpy as jnp

jax.config.update("jax_enable_x64", True)
rng = np.random.default_rng(42)

n, k = 8, 5
H = jnp.asarray(rng.choice([-1.0, 1.0], size=(k, n)))   # (k, n)

sig = np.array([[(i >> j) & 1 for j in range(n)]
                for i in range(2 ** n)]) * 2.0 - 1.0     # (256, 8)
exact = float((jnp.asarray(sig) @ H.T).max(axis=1).mean() / n)

m = 100_000
s = jnp.asarray(rng.choice([-1.0, 1.0], size=(m, n)))   # (m, n)
mc = float((s @ H.T).max(axis=1).mean() / n)

massart = float(np.sqrt(2 * np.log(k) / n))
print(f"exact={exact:.6f}  mc={mc:.6f}  massart={massart:.6f}")

import numpy as np

rng = np.random.default_rng(42)
n, k = 8, 5
H = rng.choice([-1.0, 1.0], size=(k, n))                 # (k, n)

sig = np.array([[(i >> j) & 1 for j in range(n)]
                for i in range(2 ** n)]) * 2.0 - 1.0     # (256, 8)
exact = (sig @ H.T).max(axis=1).mean() / n

s = rng.choice([-1.0, 1.0], size=(100_000, n))           # (m, n)
mc = (s @ H.T).max(axis=1).mean() / n

print(f"exact={exact:.6f}  mc={mc:.6f}  "
      f"massart={np.sqrt(2 * np.log(k) / n):.6f}")

The measured values are exact \(0.401367\) over all 256 sign vectors, Monte Carlo \(0.400500\) with 1,000 draws and \(0.401973\) with 100,000, against Massart's \(0.634318\). The same harness on the 17 threshold behaviors over 16 points gives exact \(0.342535\) (all 65,536 sign vectors) versus Monte Carlo \(0.342498\) (200,000 draws) and Massart \(0.595106\). Swapping in the behavior vectors of any classifier restricted to a validation batch turns this into a practical complexity probe. Run on a deep network's top layer it is a direct measurement of the "fits random labels" capacity from Zhang et al.

Margin distribution of a boosted classifier

AdaBoost with depth-1 stumps on the UCI breast cancer dataset (400 train / 169 test), tracking the normalized margin distribution \(y_i g(x_i) / \sum_t \alpha_t\) as rounds accumulate. The Python tab uses scikit-learn (the quoted \(5.33\% \to 3.55\%\) series), while the PyTorch and JAX tabs implement the algorithm from scratch, vectorized over all candidate stumps at once, and reproduce the same signature (identical algorithm with greedier tie-breaking, train error 0 by round 100, minimum margin \(+0.121 \to +0.164\), test \(4.73\% \to 4.14\%\)).

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import AdaBoostClassifier
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=400,
                                      random_state=0, stratify=y)
ypm = 2.0 * ytr - 1.0                     # labels in {-1, +1}

clf = AdaBoostClassifier(
    DecisionTreeClassifier(max_depth=1),  # stumps as weak learners
    n_estimators=1000, algorithm="SAMME", random_state=0)
clf.fit(Xtr, ytr)

alphas = clf.estimator_weights_           # (T,)
Hp = np.array([2.0 * e.predict(Xtr) - 1.0 for e in clf.estimators_])
Ht = np.array([2.0 * e.predict(Xte) - 1.0 for e in clf.estimators_])

for T in (5, 20, 100, 300, 1000):
    a = alphas[:T]
    marg = ypm * (a[:, None] * Hp[:T]).sum(0) / a.sum()   # (400,)
    fte = (a[:, None] * Ht[:T]).sum(0)
    print(f"T={T:>4} train={np.mean(marg <= 0):.4f} "
          f"test={np.mean(np.sign(fte) != 2.0 * yte - 1.0):.4f} "
          f"min_margin={marg.min():+.4f}")

import numpy as np
import torch
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split

dev, dt = "cuda", torch.float64
X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=400,
                                      random_state=0, stratify=y)
Xt = torch.tensor(Xtr, device=dev, dtype=dt)      # (400, 30)
yt = torch.tensor(2.0 * ytr - 1.0, device=dev, dtype=dt)

# enumerate every stump sign(x_j - thr), thr at midpoints per feature
preds = []
for j in range(30):
    v = torch.unique(Xt[:, j])
    thr = (v[:-1] + v[1:]) / 2                    # (n_thr,)
    preds.append((Xt[None, :, j] > thr[:, None]).to(dt) * 2 - 1)
P = torch.cat(preds)                              # (S, 400), S = 11055

w = torch.full((400,), 1 / 400, device=dev, dtype=dt)
f, alpha_sum = torch.zeros(400, device=dev, dtype=dt), 0.0
for t in range(1, 1001):
    corr = P @ (w * yt)                           # (S,) = 1 - 2*w_err
    err = 0.5 - 0.5 * corr
    best = int(torch.argmin(torch.minimum(err, 1 - err)))
    sgn = 1.0 if err[best] <= 1 - err[best] else -1.0
    e = float(min(err[best], 1 - err[best]).clamp(1e-12, 1 - 1e-12))
    a = 0.5 * np.log((1 - e) / e)                 # alpha_t
    h = sgn * P[best]                             # (400,)
    w = w * torch.exp(-a * yt * h); w = w / w.sum()
    f = f + a * h; alpha_sum += a
    if t in (100, 300, 1000):
        marg = yt * f / alpha_sum
        print(f"T={t} train={(marg <= 0).double().mean():.4f} "
              f"min_margin={marg.min():+.4f}")

import numpy as np
import jax
import jax.numpy as jnp
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split

jax.config.update("jax_enable_x64", True)
X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, train_size=400,
                                      random_state=0, stratify=y)
ypm = jnp.asarray(2.0 * ytr - 1.0)                # (400,)

preds = []
for j in range(30):
    v = np.unique(Xtr[:, j])
    thr = (v[:-1] + v[1:]) / 2.0
    preds.append(np.where(Xtr[None, :, j] > thr[:, None], 1.0, -1.0))
P = jnp.asarray(np.concatenate(preds))            # (S, 400)

@jax.jit
def boost_round(w):
    corr = P @ (w * ypm)                          # (S,)
    err = 0.5 - 0.5 * corr
    best = jnp.argmin(jnp.minimum(err, 1 - err))
    sgn = jnp.where(err[best] <= 1 - err[best], 1.0, -1.0)
    e = jnp.clip(jnp.minimum(err[best], 1 - err[best]), 1e-12, 1 - 1e-12)
    a = 0.5 * jnp.log((1 - e) / e)
    h = sgn * P[best]
    w = w * jnp.exp(-a * ypm * h)
    return w / w.sum(), a, h

w = jnp.full(400, 1 / 400)
f, alpha_sum = jnp.zeros(400), 0.0
for t in range(1, 1001):
    w, a, h = boost_round(w)
    f, alpha_sum = f + a * h, alpha_sum + float(a)
    if t in (100, 300, 1000):
        marg = ypm * f / alpha_sum
        print(f"T={t} train={float((marg <= 0).mean()):.4f} "
              f"min_margin={float(marg.min()):+.4f}")

The scikit-learn tab measures, at \(T = 5\), train error \(2.75\%\), test \(7.10\%\), and minimum margin \(-0.366\). At \(T = 100\), train is \(0\), test \(5.33\%\), minimum margin \(+0.118\). At \(T = 300\), test is \(4.73\%\), minimum margin \(+0.143\). At \(T = 1000\), test is \(3.55\%\), minimum margin \(+0.155\) with the 5th-percentile margin at \(0.160\). Nine hundred rounds after the training labels were exhausted as a signal, the test error improved by a third, and the margin distribution is the only training-set quantity that moved.

Implicit bias of gradient descent, watched live

Full-batch gradient descent on logistic loss over separable 2-D data, measuring the cosine between the iterate direction and the hard-margin SVM direction computed independently. The PyTorch tab is the run quoted in the implicit-bias section.

import numpy as np
import torch

dev, dt = "cuda", torch.float64
g = torch.Generator(device=dev).manual_seed(0)
Xp = torch.randn(60, 2, generator=g, device=dev, dtype=dt) * 0.7 \
    + torch.tensor([2.0, 1.0], device=dev, dtype=dt)
Xm = torch.randn(60, 2, generator=g, device=dev, dtype=dt) * 0.7 \
    + torch.tensor([-2.0, -1.0], device=dev, dtype=dt)
X = torch.cat([Xp, Xm])                            # (120, 2), separable
y = torch.cat([torch.ones(60, device=dev, dtype=dt),
               -torch.ones(60, device=dev, dtype=dt)])

# reference max-margin direction: penalized hard-margin solve via LBFGS
w_s = torch.zeros(2, device=dev, dtype=dt, requires_grad=True)
b_s = torch.zeros(1, device=dev, dtype=dt, requires_grad=True)
opt = torch.optim.LBFGS([w_s, b_s], max_iter=500)
def closure():
    opt.zero_grad()
    slack = torch.clamp(1 - y * (X @ w_s + b_s), min=0)
    loss = 0.5 * (w_s ** 2).sum() + 1e6 * (slack ** 2).sum()
    loss.backward(); return loss
for _ in range(20): opt.step(closure)
w_svm = (w_s / torch.linalg.norm(w_s)).detach()    # measured margin 0.4052

Xb = torch.cat([X, torch.ones(120, 1, device=dev, dtype=dt)], 1)
w = torch.zeros(3, device=dev, dtype=dt)           # last coord = bias
for t in range(1, 100001):
    s = torch.sigmoid(-y * (Xb @ w))               # (120,)
    w = w - 0.5 * (-(y * s)[:, None] * Xb).mean(0) # GD on logistic loss
    if t in (100, 1000, 10000, 100000):
        wd = w[:2] / torch.linalg.norm(w[:2])
        print(f"t={t:>6} |w|={torch.linalg.norm(w):.3f} "
              f"cos={float(wd @ w_svm):.6f}")

import numpy as np
import jax
import jax.numpy as jnp

jax.config.update("jax_enable_x64", True)
rng = np.random.default_rng(0)
Xp = rng.standard_normal((60, 2)) * 0.7 + np.array([2.0, 1.0])
Xm = rng.standard_normal((60, 2)) * 0.7 + np.array([-2.0, -1.0])
X = jnp.asarray(np.vstack([Xp, Xm]))               # (120, 2)
y = jnp.asarray(np.r_[np.ones(60), -np.ones(60)])
Xb = jnp.concatenate([X, jnp.ones((120, 1))], 1)   # bias folded in

def loss(w):                                       # mean logistic loss
    return jnp.mean(jnp.log1p(jnp.exp(-y * (Xb @ w))))

step = jax.jit(lambda w: w - 0.5 * jax.grad(loss)(w))

w = jnp.zeros(3)
for t in range(1, 100001):
    w = step(w)
    if t in (100, 1000, 10000, 100000):
        # ||w|| grows ~ log t; direction w[:2]/||w[:2]|| converges to
        # the max-margin separator of THIS sample (compute as at left)
        print(t, float(jnp.linalg.norm(w)))

The PyTorch tab measures a cosine to the max-margin direction of \(0.986468\), \(0.997077\), \(0.999857\), and \(0.999996\) at \(t = 10^2, 10^3, 10^4, 10^5\), with \(\lVert w \rVert\) growing \(2.72, 4.61, 8.09, 13.07\), logarithmic growth with a converging direction, exactly the \(w(t) \approx \hat w \ln t\) shape of the theorem. The convergence rate in angle is \(O(1/\ln t)\), which the numbers also show. Each decade of steps buys roughly a constant factor in \(1 - \cos\).

How it is done in practice

The theory on this page is deployed every day, mostly without being named. The workhorse is the held-out estimate plus Hoeffding. A test set of \(m\) examples estimates any fixed model's error to within \(\sqrt{\ln(2/\delta)/2m}\), which at \(m = 10{,}000\) and \(\delta = 0.05\) is \(\pm 1.36\%\). Production teams that track model quality to a tenth of a point on a ten-thousand-example eval set are, by this arithmetic, reading noise, and the fix (bigger eval sets or paired comparisons, which cancel per-example variance) is finite-class reasoning applied honestly. Model selection over \(k\) candidates on a shared validation set is the union bound in costume. The winner's reported score is biased upward by up to \(\sqrt{\ln(2k/\delta)/2m}\), a few tenths of a point for typical sweeps, which is why serious comparisons re-measure the single chosen model on a fresh split. Repeated leaderboard probing is the same phenomenon compounding. Blum and Hardt's ladder mechanism (2015) formalized how many adaptive queries a test set survives, and the ImageNet-v2 replication study of Recht, Roelofs, Schmidt, and Shankar (2019) measured the ecosystem outcome, absolute accuracies dropped several points on a freshly collected test set while model rankings were preserved, adaptive overfitting to the benchmark distribution being smaller than the community feared but the distribution-shift sensitivity larger.

Where explicit certificates matter, the frameworks split by regime. Finite-class and Hoeffding bounds run inside A/B testing platforms and safety evaluations, where the hypothesis is fixed and the arithmetic is exactly Problem 1. Stability is a design principle more than a computed number. Weight decay, early stopping, dropout, and data augmentation are all interventions whose first-order effect is reducing the output's sensitivity to individual examples, and the Hardt-Recht-Singer analysis is the standard citation for why bounded training time itself regularizes. PAC-Bayes is the only framework currently producing nonvacuous numeric certificates for neural networks, and it remains a research-grade activity (Dziugaite-Roy and successors) rather than a deployment default. The VC and Rademacher machinery earns its keep upstream, in how practitioners reason. The margin bound is why calibrated confidence and margin distributions are monitored in gradient-boosting deployments (XGBoost and LightGBM inherit AdaBoost's margin story with different losses), and the min-norm/implicit-bias results are why "train an overparameterized model to convergence with weak explicit regularization" stopped being malpractice and became the default recipe. The double-descent peak is an operational hazard, not a curiosity. It sits exactly where parameter count matches sample count, and teams that scale models and datasets together can walk into it. Both the Belkin et al. and Nakkiran et al. papers show the practical mitigations, regularize even lightly, or jump capacity well past the threshold rather than stopping on it.

The current research frontier

The sharpest open question remains the one Zhang et al. posed. No known complexity measure simultaneously (a) is computable for real networks, (b) provably bounds their generalization, and (c) correlates with observed generalization across architectures and datasets. Norm-based bounds (Bartlett, Foster, and Telgarsky's spectrally-normalized margin bound, 2017, and Neyshabur, Bhojanapalli, and Srebro's PAC-Bayes analogue, 2018) satisfy (b) but are numerically vacuous at scale and sometimes anti-correlate in the large empirical audits (Jiang, Neyshabur, Mobahi, Krishnan, Bengio 2020, the "fantastic generalization measures" study at Google). Nagarajan and Kolter (2019, CMU) sharpened the negative side. In natural settings, any two-sided uniform-convergence argument, however cleverly restricted to the reachable set, provably cannot explain generalization, because the reachable set is itself too rich, an argument that pushed the field toward stability, PAC-Bayes, and trajectory-dependent analyses. On the constructive side, benign overfitting has moved beyond linear models. Tsigler and Bartlett extended it to ridge, Frei, Chatterji, and Bartlett and separately Cao, Chen, Belkin, and Gu to shallow networks and convolutional cases, and Mallinar, Simon, Abedsoltan, Pandit, Belkin, and Nakkiran (2022) proposed the useful taxonomy of benign, tempered, and catastrophic overfitting, most practical interpolation being tempered, worse than optimal, better than disaster.

The implicit-bias program is filling in the optimizer column, through Lyu and Li's directional convergence for homogeneous networks, Chizat and Bach's mean-field versus lazy dichotomy at INRIA/ENS, Gunasekar, Woodworth, and Srebro on mirror-descent geometries determining which norm gets implicitly minimized, and a growing line on the implicit bias of Adam and of large-step-size dynamics (edge of stability, Cohen et al. 2021) where the gradient-flow picture breaks. PAC-Bayes and compression fused into the current best-certificate technology (Lotfi et al. 2022, NYU, with earlier Zhou, Veeravalli, and colleagues' compression bounds following Arora, Ge, Neyshabur, and Zhang 2018), and the same toolkit was turned on language models to bound memorization versus generalization at token level (Lotfi, Finzi, and collaborators, 2024). Grokking (Power, Burda, Edwards, Babuschkin, Misra 2022, OpenAI), where test accuracy jumps from chance to perfect long after training accuracy saturates, became the fruit-fly of implicit-regularization dynamics, with mechanistic accounts (Nanda et al., Anthropic-adjacent interpretability work, and Liu, Michaud, Tegmark at MIT) connecting the transition to the slow growth of structured representations, a phenomenon the margin-and-norm lens from this page describes surprisingly well. And scaling laws sit at the edge of the subject. Empirical risk curves of the Kaplan et al. and Hoffmann et al. (Chinchilla) form behave like power-law bias-variance tradeoffs in disguise, and putting them on the same footing as the finite-sample theory here, a quantitative theory of when more data beats more parameters, is arguably the field's most consequential open modeling problem.

Open source to read

scikit-learn/scikit-learn, the reference AdaBoost used above. Open sklearn/ensemble/_weight_boosting.py and read _boost_discrete. The weight update, \(\alpha_t\) formula, and early-termination logic map line for line to the derivation in the margins section.

google/neural-tangents, exact infinite-width NTK and NNGP kernels for concrete architectures, built on JAX. Open neural_tangents/_src/stax/linear.py to see how each layer transforms the kernel, then the empirical.py module for finite-width tangent kernels, the fastest way to test whether a network is in the lazy regime.

jax-ml/jax. Beyond the demos above, the autodiff machinery is the right tool for implicit-bias experiments generally. Open jax/_src/api.py for grad/jvp. An NTK is four lines of jax.jvp against parameters.

pytorch/pytorch. torch.func (in torch/_functorch/) provides per-sample gradients and Jacobians, the primitives for measuring stability (how much one example moves the model) and empirical tangent kernels on real architectures.

dmlc/xgboost, the industrial descendant of the boosting story. Open src/objective/regression_obj.cu to see how losses enter only through gradients and Hessians. The margin distribution of any trained model is one predict(output_margin=True) call away.

tomgoldstein/loss-landscape, the standard visualization code for loss-surface geometry (filter-normalized directions), useful for building intuition about flatness claims before believing or disbelieving them. Open plot_surface.py.

probml/pyprobml, executable notebooks for Murphy's books, including double descent, kernel regression, and PAC-Bayes style calculations, and a good source of small, checkable reference implementations to test derivations against.

Common misconceptions

"Deep networks generalize because their VC dimension is small enough." It is not. VC dimension of ReLU networks scales roughly like parameters times depth (times log factors), vastly exceeding dataset sizes in every modern regime, and Zhang et al. verified empirically that the trained classes shatter real-sized samples. Whatever explains generalization must depend on the algorithm and the data, not the architecture's function class alone. That is the precise content of the randomization experiment, and the reason this page spends its second half on stability, implicit bias, and PAC-Bayes.

"Zero training error means overfitting." Interpolation is compatible with optimal rates when the fitting mechanism spreads noise across many low-signal directions. The benign-overfitting theorem gives exact conditions, and the measured double-descent table shows an interpolating model (\(N = 12{,}000\)) beating every non-interpolating one. What matters is never the training error but which interpolator the algorithm selects. "Zero training error" at the interpolation threshold, where the selected solution has an enormous norm, genuinely is catastrophic, and the same words two decades of capacity to the right are benign.

"The bias-variance tradeoff has been disproved." The decomposition is an identity and cannot be disproved. What failed is the assumption that "capacity" moves bias and variance monotonically in opposite directions along a single axis. In the overparameterized regime, growing capacity changes which minimum the optimizer selects, and the min-norm solution's variance decreases with added features, as the falling solution norms in the table show directly. Double descent is bias-variance with the correct accounting of what the estimator actually is.

"AdaBoost should overfit as rounds accumulate, since capacity grows with every term." The ensemble lives in the convex hull of the base class, whose Rademacher complexity does not grow with the number of terms. What grows is the margin distribution, and the derived bound \(\prod_t 2\varepsilon_t^{(1-\theta)/2}(1-\varepsilon_t)^{(1+\theta)/2}\) shows small-margin mass shrinking exponentially in rounds. The measured run (test error falling from 5.3% to 3.6% across 900 post-interpolation rounds) is the phenomenon. Unbounded rounds on noisy data can eventually overfit, but the mechanism is margin decay on mislabeled points, not term-counting.

"Rademacher complexity is a property of the model class, like VC dimension." It is a property of the class and the sample, which is its entire advantage. The same class has different complexity on different data (the measured 0.401 vs 0.343 above), margins and norms enter through the sample, and data-dependent bounds can be tight where distribution-free ones are vacuous. Treating it as a fixed architecture constant erases exactly the term that does the work.

"Uniform stability requires convexity, so it says nothing about deep learning." The clean \(2\rho^2/(\lambda n)\) rate uses strong convexity, but stability itself is a property of the algorithm. Hardt, Recht, and Singer bound the stability of SGD on smooth non-convex objectives as a function of steps and step sizes, and stability-flavored arguments underlie current analyses of early stopping and checkpoint averaging. What is true is that stability bounds degrade with training time in the non-convex case, which is a real limitation, not a disqualification.

"PAC-Bayes is Bayesian inference with a rebranding." The prior in PAC-Bayes needs no relationship to any true data-generating process, the posterior is any distribution the learner likes (typically chosen by optimizing the bound, not by Bayes' rule), and the guarantee is a frequentist tail bound holding uniformly over posteriors. The KL term is a description-length price, not a belief update. The framework is closer to MDL than to Bayesian statistics, which is exactly why it survives misspecification that breaks Bayesian coverage claims.

"Margin maximization fully explains boosting, so bigger minimum margin always means better generalization." Breiman's 1999 arc-gv counterexample achieves a larger minimum margin than AdaBoost with worse test error. The resolution (Reyzin-Schapire, then the k-th margin theory of Gao and Zhou) is that the guarantee depends on the whole margin distribution and the base-class complexity together, and arc-gv buys its minimum margin with more complex base learners. The measured margin experiments above report percentiles, not just the minimum, for exactly this reason.

Self-check

References

1. Shalev-Shwartz, S. and Ben-David, S. (2014). Understanding Machine Learning: From Theory to Algorithms. Cambridge University Press. The cleanest single source for PAC learning, VC theory, Rademacher complexity, and stability. book page
2. Mohri, M., Rostamizadeh, A., and Talwalkar, A. (2018). Foundations of Machine Learning, 2nd ed. MIT Press. The Rademacher-first treatment. The constants used on this page follow it. book page
3. Boucheron, S., Lugosi, G., and Massart, P. (2013). Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press. The full development of the Chernoff method, McDiarmid, and much more. doi:10.1093/acprof:oso/9780199535255.001.0001
4. Vapnik, V. N. and Chervonenkis, A. Ya. (1971). On the uniform convergence of relative frequencies of events to their probabilities. Theory of Probability and its Applications 16(2). The founding paper, with growth functions, symmetrization, and the VC bound. doi:10.1137/1116025
5. Valiant, L. G. (1984). A theory of the learnable. Communications of the ACM 27(11). The PAC framework. doi:10.1145/1968.1972
6. Sauer, N. (1972). On the density of families of sets. Journal of Combinatorial Theory, Series A 13(1). The lemma. Shelah proved a parallel version the same year in model theory. doi:10.1016/0097-3165(72)90019-2
7. Hoeffding, W. (1963). Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association 58(301). doi:10.1080/01621459.1963.10500830
8. Blumer, A., Ehrenfeucht, A., Haussler, D., and Warmuth, M. K. (1989). Learnability and the Vapnik-Chervonenkis dimension. Journal of the ACM 36(4). Finite VC dimension characterizes PAC learnability. doi:10.1145/76359.76371
9. Bartlett, P. L. and Mendelson, S. (2002). Rademacher and Gaussian complexities: risk bounds and structural results. JMLR 3. jmlr.org/papers/v3/bartlett02a
10. Bousquet, O. and Elisseeff, A. (2002). Stability and generalization. JMLR 2. Uniform stability and the theorems of the stability section. jmlr.org/papers/v2/bousquet02a
11. Hardt, M., Recht, B., and Singer, Y. (2016). Train faster, generalize better: stability of stochastic gradient descent. ICML. arXiv:1509.01240
12. Freund, Y. and Schapire, R. E. (1997). A decision-theoretic generalization of on-line learning and an application to boosting. Journal of Computer and System Sciences 55(1). AdaBoost. doi:10.1006/jcss.1997.1504
13. Schapire, R. E., Freund, Y., Bartlett, P., and Lee, W. S. (1998). Boosting the margin: a new explanation for the effectiveness of voting methods. Annals of Statistics 26(5). doi:10.1214/aos/1024691352
14. Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O. (2017). Understanding deep learning requires rethinking generalization. ICLR. arXiv:1611.03530
15. Belkin, M., Hsu, D., Ma, S., and Mandal, S. (2019). Reconciling modern machine-learning practice and the classical bias-variance trade-off. PNAS 116(32). Double descent named and demonstrated. arXiv:1812.11118
16. Nakkiran, P., Kaplun, G., Bansal, Y., Yang, T., Barak, B., and Sutskever, I. (2019). Deep double descent: where bigger models and more data hurt. ICLR 2020. arXiv:1912.02292
17. Bartlett, P. L., Long, P. M., Lugosi, G., and Tsigler, A. (2020). Benign overfitting in linear regression. PNAS 117(48). The effective-rank characterization stated in full above, with proofs. arXiv:1906.11300
18. Soudry, D., Hoffer, E., Nacson, M. S., Gunasekar, S., and Srebro, N. (2018). The implicit bias of gradient descent on separable data. JMLR 19. arXiv:1710.10345
19. Jacot, A., Gabriel, F., and Hongler, C. (2018). Neural tangent kernel: convergence and generalization in neural networks. NeurIPS. arXiv:1806.07572
20. McAllester, D. A. (1999). Some PAC-Bayesian theorems. Machine Learning 37(3). Maurer, A. (2004), A note on the PAC-Bayesian theorem, gives the sharp form derived above. doi:10.1023/A:1007618624809
21. Dziugaite, G. K. and Roy, D. M. (2017). Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data. UAI. arXiv:1703.11008
22. Bartlett, P. L., Foster, D. J., and Telgarsky, M. (2017). Spectrally-normalized margin bounds for neural networks. NeurIPS. arXiv:1706.08498
23. Nagarajan, V. and Kolter, J. Z. (2019). Uniform convergence may be unable to explain generalization in deep learning. NeurIPS. arXiv:1902.04742
24. Rahimi, A. and Recht, B. (2007). Random features for large-scale kernel machines. NeurIPS. The feature map of the double-descent experiment. paper
25. Belkin, M. (2021). Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation. Acta Numerica 30. The survey to read after this page. arXiv:2105.14368