---
cluster: "B"
day: "B.6"
contributors:
- Aleksander Czejdo
- Tudor Dimofte
- Brianna Grado-White
- Charles Renshaw-Whitman
summary: "An introduction to physics-inspired analysis and inquiry in deep learning, particularly focusing on techniques from statistical mechanics, statistical field theory, and (a bit of) dynamical systems."
title: "Physics of Deep Learning"
---
$\gdef\R{\mathbb{R}}\gdef\E{\mathbb{E}}\gdef\bv#1{\boldsymbol{#1}}\gdef\ntk{K_{\mathrm{NTK}}}\gdef\fixme#1{\textcolor{blue}{[#1]}}\gdef\CN{{\mathcal N}}\gdef\CL{{\mathcal L}}\gdef\CO{{\mathcal O}}\gdef\pd{\partial}\gdef\expval#1{\left\langle #1 \right\rangle}\gdef\dv#1#2{\frac{d #1}{d #2}}\gdef\abs#1{\left| #1 \right|}\gdef\arraystretch{1.25}$
**Contents**
- [-1. Prerequisites](#1-prerequisites)
- [0. A one-day teaching plan](#0-a-one-day-teaching-plan)
- [0.1 Schedule](#01-schedule)
- [0.2 Main references](#02-main-references)
- [1. Intro reading: learning mechanics](#1-intro-reading-learning-mechanics)
- [2. Statistical mechanics of learning](#2-statistical-mechanics-of-learning)
- [2.1 The Bayesian model of learning](#21-the-bayesian-model-of-learning)
- [2.2 Why the exponential of a loss is a likelihood](#22-why-the-exponential-of-a-loss-is-a-likelihood)
- [2.3 Partition functions, free energy, and cumulants](#23-partition-functions-free-energy-and-cumulants)
- [2.4 A physics example: the Ising model](#24-a-physics-example-the-ising-model)
- [2.5 Phases and phase transitions](#25-phases-and-phase-transitions)
- [2.6 The Ising perceptron: a first-order transition in learning](#26-the-ising-perceptron-a-first-order-transition-in-learning)
- [2.7 Comparison with QFT](#27-comparison-with-qft)
- [3. Large width at initialization: the Gaussian process limit](#3-large-width-at-initialization-the-gaussian-process-limit)
- [3.1 Notation](#31-notation)
- [3.2 Normalizing the network](#32-normalizing-the-network)
- [3.3 Gaussian processes](#33-gaussian-processes)
- [3.4 The NNGP correspondence](#34-the-nngp-correspondence)
- [3.5 NNGP via functional integrals](#35-nngp-via-functional-integrals)
- [3.6 L/N](#36-ln)
- [3.7 Edge of chaos](#37-edge-of-chaos)
- [4. The NTK and large width during training](#4-the-ntk-and-large-width-during-training)
- [4.1 Gradient flow and the NTK](#41-gradient-flow-and-the-ntk)
- [4.1.1 Geometry of gradient flows](#411-geometry-of-gradient-flows)
- [4.2 Features](#42-features)
- [4.3 Constant NTK at large width](#43-constant-ntk-at-large-width)
- [4.4 Lazy learning](#44-lazy-learning)
- [4.5 Summary](#45-summary)
- [4.6 Bayesian learning at large width: kernel learning](#46-bayesian-learning-at-large-width-kernel-learning)
- [5. Mean-field scaling](#5-mean-field-scaling)
- [5.1 Mean-field theories in physics](#51-mean-field-theories-in-physics)
- [5.1.1 Warm-up: the Curie–Weiss magnet](#511-warm-up-the-curieweiss-magnet)
- [5.1.2 The general shape of a mean-field model](#512-the-general-shape-of-a-mean-field-model)
- [5.2 Mean-field scaling for one-hidden-layer networks](#52-mean-field-scaling-for-one-hidden-layer-networks)
- [5.2.1 Initialization](#521-initialization)
- [5.2.2 The loss is a mean-field energy](#522-the-loss-is-a-mean-field-energy)
- [5.2.3 Equations of motion and the mean field](#523-equations-of-motion-and-the-mean-field)
- [5.2.4 From particles to measures: the continuity equation](#524-from-particles-to-measures-the-continuity-equation)
- [5.2.5 The McKean–Vlasov equation and self-consistency](#525-the-mckeanvlasov-equation-and-self-consistency)
- [5.3 The NTK is no longer constant: feature learning](#53-the-ntk-is-no-longer-constant-feature-learning)
- [5.4 The scaling dial and Maximal Update Parameterization](#54-the-scaling-dial-and-maximal-update-parameterization)
- [5.5 What is this framework good for?](#55-what-is-this-framework-good-for)
- [5.6 Grokking as a lazy-to-rich transition: a reading exercise](#56-grokking-as-a-lazy-to-rich-transition-a-reading-exercise)
- [5.7 Bayesian grokking: a first-order transition (alternative reading)](#57-bayesian-grokking-a-first-order-transition-alternative-reading)
- [6. Non-equilibrium physics, MSRJD, and DMFT](#6-non-equilibrium-physics-msrjd-and-dmft)
- [6.1 From stochastic dynamics to path integrals: the MSRJD formalism](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism)
- [6.2 Double descent as a dynamical phase transition](#62-double-descent-as-a-dynamical-phase-transition)
- [6.2.1 Model definition and formalism](#621-model-definition-and-formalism)
- [6.2.2 Overview of results via training and generalization error](#622-overview-of-results-via-training-and-generalization-error)
- [6.3 Dynamical mean-field theory (DMFT)](#63-dynamical-mean-field-theory-dmft)
- [6.3.1 Self-consistency](#631-self-consistency)
- [6.3.2 Recent applications in deep learning](#632-recent-applications-in-deep-learning)
- [7. Scaling laws](#7-scaling-laws)
- [7.1 The empirical laws](#71-the-empirical-laws)
- [7.2 Joint laws, bottlenecks, and the Chinchilla revision](#72-joint-laws-bottlenecks-and-the-chinchilla-revision)
- [7.3 Where do the laws come from?](#73-where-do-the-laws-come-from)
- [7.4 A survey of theoretical proposals](#74-a-survey-of-theoretical-proposals)
- [7.5 The quantization model: exponents from the statistics of skills](#75-the-quantization-model-exponents-from-the-statistics-of-skills)
- [7.6 Emergence: smooth laws from staggered cliffs](#76-emergence-smooth-laws-from-staggered-cliffs)
- [7.7 Exponents from geometry and from spectra](#77-exponents-from-geometry-and-from-spectra)
- [7.7.1 Geometry: the data-manifold argument](#771-geometry-the-data-manifold-argument)
- [7.7.2 Spectra: scaling laws from kernel eigenvalues](#772-spectra-scaling-laws-from-kernel-eigenvalues)
- [7.8 A solvable large-N field theory](#78-a-solvable-large-n-field-theory)
- [7.9 All the knobs at once: a dynamical theory](#79-all-the-knobs-at-once-a-dynamical-theory)
- [7.10 Beyond pretraining: post-training and test-time scaling](#710-beyond-pretraining-post-training-and-test-time-scaling)
- [A. The replica method](#a-the-replica-method)
- [A.1 What the annealed approximation gets wrong](#a1-what-the-annealed-approximation-gets-wrong)
- [A.2 Replicas](#a2-replicas)
- [A.3 A variant: the spherical perceptron](#a3-a-variant-the-spherical-perceptron)
- [B. From linear regression to kernel regression](#b-from-linear-regression-to-kernel-regression)
- [B.1 Linear regression](#b1-linear-regression)
- [B.2 Linear regression with nonlinear features](#b2-linear-regression-with-nonlinear-features)
- [B.3 Kernel regression](#b3-kernel-regression)
- [B.4 Kernel ridge regression](#b4-kernel-ridge-regression)
- [B.5 Linearized networks are kernel ridge regressors](#b5-linearized-networks-are-kernel-ridge-regressors)
- [C. Solutions to exercises](#c-solutions-to-exercises)
## -1. Prerequisites
We don't assume any physics background, though reading some of the suggested reference material may be helpful. We do assume:
- a basic understanding of probability theory and Bayesian inference;
- familiarity with simple neural networks (MLP's);
- undergraduate-level calculus and linear algebra;
- and, in a few places, differential geometry.
By the end of these notes, you should be able to:
- Use the statistical-mechanics dictionary for Bayesian learning — energy, entropy, partition function, free energy, cumulants, order parameters, quenched averages — and explain why the exponential of a loss function is a likelihood.
- Solve the Ising model in the mean-field approximation; explain what a phase transition is, why finite systems show smooth sigmoids and growing susceptibility peaks instead of true kinks, and how the Ising perceptron realizes a first-order transition to perfect generalization, metastable lag included.
- Derive the Gaussian-process description of a wide network at initialization, and explain the roles of the fan-in normalization, the depth-to-width ratio, and the edge of chaos.
- Derive the NTK, and explain why it freezes at large width — so that training reduces to kernel regression, with no feature learning.
- Derive, heuristically, the mean-field limit of one-hidden-layer networks — neurons as interacting particles, self-consistency, the McKean–Vlasov equation — and explain how its *moving* kernel encodes feature learning; place the NTK and mean-field limits on the scaling dial, and connect it to $\mu$P and hyperparameter transfer.
- Compare the two theoretical accounts of grokking — an equilibrium first-order transition vs. delayed lazy-to-rich dynamics — and propose an experiment that would distinguish them.
- Recognize the MSRJD/DMFT path-integral formulation of training dynamics, the role of correlation and response functions, and double descent as a dynamical phase transition.
- Derive the three scaling laws of the quantization model, explain "emergence" as the aggregation of staggered cliffs, and survey where scaling exponents come from: skill statistics, data geometry, and kernel spectra.
## 0. A one-day teaching plan
These notes contain far more than can be taught in a day. This section records one selection of material that was taught in a single day at the Iliad Intensive (September 2026, day B.6), as a worked example of how the notes can be used. The day alternates three lectures, given from the slide deck that accompanies these notes, with exercise sessions in which participants pick one exercise from the notes and work on it, alone or in small groups. It closes with a short reading and a discussion.
The deck is the file `slides.tex` in the same folder as these notes. The website build compiles it and offers it in the "Slides" row at the top of this chapter's web page; it can also be opened directly at [iliad-team.github.io/iliad-intensive/downloads/qft/qft-slides.pdf](https://iliad-team.github.io/iliad-intensive/downloads/qft/qft-slides.pdf). Its three parts are the three lectures below: frames 1–20 (Part I), 21–41 (Part II), and 42–52 (Part III). The lectures deliberately leave derivations to the exercise sessions, so the deck is sparse by design; the corresponding sections of these notes fill in what the slides state.
### 0.1 Schedule
| **time** | **session** | **where in these notes** |
| --- | --- | --- |
| 10:00–11:00 | Lecture, Part I: Bayesian learning and phase transitions | [Sections 2.1–2.5](#21-the-bayesian-model-of-learning) |
| 11:00–11:15 | Break | |
| 11:15–12:00 | Exercises, Part I (choose one) | [Exercises 2.1](#ex-cumulants), [2.2](#ex-anatomy) and [2.3](#ex-ising) |
| 12:00–12:30 | Lecture, Part IIa: large width at initialization; the NTK | [Sections 3.1–3.4](#31-notation), [Section 3.6](#36-ln), [Sections 4.1–4.2](#41-gradient-flow-and-the-ntk) |
| 12:30–1:30 | Lunch | |
| 1:30–2:00 | Lecture, Part IIb: large width during training; mean-field scaling | [Sections 4.3–4.4](#43-constant-ntk-at-large-width), [Sections 5.2–5.4](#52-mean-field-scaling-for-one-hidden-layer-networks) |
| 2:00–3:00 | Exercises or lab, Part II (choose one) | [Exercises 3.3](#ex-ln-llms), [4.5](#ex-kernel-ridge), [4.4](#ex-ntk-lab) and [5.3](#ex-mf-lab) |
| 3:00–3:45 | Lecture, Part III: neural scaling laws, in theory and practice | [Sections 7.1–7.5](#71-the-empirical-laws) |
| 3:45–4:15 | Break | |
| 4:15–5:00 | Exercise, Part III | [Exercise 7.1](#ex-quanta) |
| 5:00–5:30 | Reading on "learning mechanics" | [Section 1](#1-intro-reading-learning-mechanics) |
| 5:30–5:50 | Discussion in small groups | the prompts of [Section 1](#1-intro-reading-learning-mechanics) |
| 5:50–6:00 | Feedback | |
### 0.2 Main references
The day is based on the following references, all of which appear elsewhere in these notes. For Part I: Tong's lectures on statistical field theory and on kinetic theory ([Tong 2012](#bib-tongsft); [Tong 2012](#bib-tongkinetictheory)), the physics behind the lecture; Engel and Van den Broeck ([Engel & Van den Broeck 2001](#bib-engel2001statistical)) for the perceptron tradition; and Gyorgyi ([Gy"orgyi 1990](#bib-gyorgyi1990first)) for the Ising perceptron transition. For Part II: Hanin's lectures ([Hanin 2026](#bib-hanin2026mdl)) and Roberts, Yaida, and Hanin ([Roberts et al. 2022](#bib-roberts2022principles)) on large-width limits and the $L/N$ expansion; Ringel et al. ([Ringel et al. 2025](#bib-ringel2025applications)) on statistical field theory in deep learning; Spiliopoulos, Sowers, and Sirignano ([Spiliopoulos et al. 2025](#bib-spiliopoulos2025foundations)) for rigorous NTK and mean-field limits; Lin ([Lin 2025](#bib-lin2025feature)) for experiments identifying features through the empirical NTK; and the two accounts of grokking ([Kumar et al. 2024](#bib-kumar2024grokking); [Rubin et al. 2024](#bib-rubin2024grokking)), which were not lectured but are the natural follow-up reading ([Sections 5.6](#56-grokking-as-a-lazy-to-rich-transition-a-reading-exercise) and [5.7](#57-bayesian-grokking-a-first-order-transition-alternative-reading)). For Part III: Kaplan et al. ([Kaplan et al. 2020](#bib-kaplan2020scaling)) and Henighan et al. ([Henighan et al. 2020](#bib-henighan2020scaling)) for the empirical laws; Hoffmann et al. ([Hoffmann et al. 2022](#bib-hoffmann2022training)) and Porian et al. ([Porian et al. 2024](#bib-porian2024resolving)) for Chinchilla and the reconciliation with Kaplan; Michaud et al. ([Michaud et al. 2023](#bib-michaud2023quantization)) for the quantization model; and Maloney, Roberts, and Sully ([Maloney et al. 2022](#bib-maloney2022solvable)) and Bordelon, Atanasov, and Pehlevan ([Bordelon et al. 2024](#bib-bordelon2024dynamical)) for the solvable field theories.
A complementary introduction to the physics of learning, taught at the August 2026 Intensive by Renshaw-Whitman, can be found in ([Renshaw-Whitman 2026](#bib-renshaw2026iliad)); [Appendix B](#b-from-linear-regression-to-kernel-regression) is adapted from it.
## 1. Intro reading: learning mechanics
The review ([Simon et al. 2026](#bib-simon2026scientific)) argues that a scientific theory of deep learning — a "mechanics of learning" — is already emerging, and it organizes the evidence into five converging strands: analytically solvable settings, insightful limits, simple empirical laws, disentangled hyperparameters, and universal phenomena. The body of the paper develops the first two strands in detail and samples the remaining three. These notes will do much the same, so it is worth keeping the five strands in mind as a map of everything that follows.
The paper's starting observation deserves to be savored, because it sets the tone for the whole subject. A deep learning system is completely specified by four explicit ingredients: an architecture (a composition of simple linear and nonlinear maps), a dataset, an objective, and a learning rule (gradient-based updates from a random initialization). Nothing is hidden: every parameter, every gradient, every intermediate activation can be measured, at every step, to whatever precision we like. In the language of the reading, deep learning directly exposes its own equations of motion — *the central challenge is not opacity but complexity*. This places the field in very different territory from the genuinely opaque sciences: in neuroscience or cell biology the microscopic equations of motion are unknown and the state is largely unmeasurable, while in high-energy physics we must infer the microscopic laws indirectly, from scattering amplitudes. The better analogy is fluid dynamics: there, too, we possess the microscopic law — for deep learning, the update equation is five lines of code — and what is missing is the *effective* description, the analogue of Navier–Stokes.
Read **Section 2** (only) of ([Simon et al. 2026](#bib-simon2026scientific)), and think about the following:
1. *Rank the physics analogies in Table 1 of the reading* (solvable models $\leftrightarrow$ hydrogen atom; limits $\leftrightarrow$ thermodynamic limit; empirical laws $\leftrightarrow$ Kepler and Planck; hyperparameters $\leftrightarrow$ Reynolds number; universality $\leftrightarrow$ critical phenomena). Which are the most compelling? Which do you think are most useful in understanding deep learning?
2. *Discretization and other approximations.* The reading's Discretization Hypothesis proposes that practical networks are best understood as noisy, finite approximations to models of infinite size — like a finite-difference discretization of a PDE — so that finite-size effects typically *cost* performance rather than enable it. Given what you've seen so far, does this seem accurate? The hypothesis is falsifiable: it would be refuted by a finite-size effect that delivers a general benefit unachievable in any limit. What are some possible experiments/setups that could falsify it?
3. *What actually matters for alignment?* The reading wants falsifiable, quantitative, *average-case* predictions rather than worst-case bounds — physics rather than classical learning theory. What is lost? E.g. for alignment, does one want worst-case guarantees or forecasting of typical capabilities? Which of the five lines of evidence is most relevant for alignment?
You may want to revisit these questions at later stages as you go through these notes, and learn more about physics-inspired approaches in deep learning.
## 2. Statistical mechanics of learning
There is a large and well-developed body of work applying physics to learning theory. Much of it consists of *statistical mechanics* applied to *Bayesian learning*. The classic textbook is Engel and Van den Broeck ([Engel & Van den Broeck 2001](#bib-engel2001statistical)); modern reviews include ([Zdeborov'a & Krzakala 2016](#bib-zdeborova2016statistical); [Bahri et al. 2020](#bib-bahri2020statistical)); and an excellent reference for the physics itself is David Tong's lecture course on statistical field theory ([Tong 2012](#bib-tongsft)), which we follow closely in [Sections 2.4–2.5](#24-a-physics-example-the-ising-model).
In this section we focus on the Bayesian model of learning. We first recall what it is, why it is a reasonable idealization of how networks are actually trained, and why the exponential of a loss function deserves to be called a likelihood ([Sections 2.1–2.2](#21-the-bayesian-model-of-learning)). We then switch to physics. We introduce partition functions, free energies, and cumulants ([Section 2.3](#23-partition-functions-free-energy-and-cumulants)), and work through the Ising model in the mean-field approximation, building up to first- and second-order phase transitions ([Sections 2.4–2.5](#24-a-physics-example-the-ising-model)). The payoff is a learning problem that is the Ising magnet's twin, the Ising perceptron, which shows a first-order transition to perfect generalization; it is left as a self-contained walk-through exercise ([Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning)). [Section 2.7](#27-comparison-with-qft) is a translation guide for readers arriving from quantum field theory. The ideas introduced here — order parameters, free energies, saddle points, quenched averages, phase transitions — will keep returning throughout the notes.
### 2.1 The Bayesian model of learning
Consider a network $f_{\theta}$ with parameters $\theta\in\R^{P}$, trained on $n$ samples $(x^{\mu},y^{\mu})_{\mu=1}^{n}$ to minimize a loss
a sum over samples of a per-sample loss $\ell$. In practice, the network is initialized at random, $\theta_{0}\sim\pi(\theta)$, and trained with some version of gradient descent,
$$
\Delta\theta \approx -\eta\,\nabla_{\theta}\CL(\theta)\,.
$$
A common idealization of this process, and the one this section is about, is *Bayesian learning*. It says that after training, the parameters should be regarded as drawn from the posterior distribution
$$
p_{\beta}(\theta) = \frac{1}{Z}\,e^{-\beta\,\CL(\theta)}\,\pi(\theta)\,,\qquad Z = \int d\theta\, e^{-\beta\,\CL(\theta)}\,\pi(\theta)\,.
$$
Here the initialization distribution $\pi(\theta)$ plays the role of the Bayesian *prior*, and $e^{-\beta\CL(\theta)}$ is interpreted as the *likelihood* of the training data given the parameters $\theta$ ([Section 2.2](#22-why-the-exponential-of-a-loss-is-a-likelihood) explains why this is reasonable). The positive constant $\beta$, the *inverse temperature*, controls how strongly the data constrain the parameters. We will see in a moment that it can be read as a rough measure of how much training has been done.
*An interpretation of $\beta$.* It is convenient to write the loss as $n$ times an average over samples,
$$
\CL(\theta) = n\,L_{n}(\theta)\,,\qquad L_{n}(\theta) := \frac{1}{n}\sum_{\mu=1}^{n}\ell\big(f_{\theta}(x^{\mu}),y^{\mu}\big)\,.
$$
The *empirical loss* $L_{n}$ stays finite as $n\to\infty$, converging to the *population loss* $L(\theta) := \E_{(x,y)}\,\ell(f_{\theta}(x),y)$, the expectation over the true data distribution. Thus for large $n$ the posterior is
and the sample size $n$ always accompanies the inverse temperature, in the combination $n\beta$. Increasing the number of samples is the same as lowering the temperature. If, moreover, samples are drawn fresh at each training step, then $n$ is also a rough measure of training time. Altogether: *more data $=$ more training $=$ larger $n\beta$ $=$ colder*. We will often write $e^{-\beta L}$ for short, but it is worth remembering that the effective $\beta$ grows with $n$.
*Bayes vs. gradient descent.* We do not train networks by sampling from [(2.3)](#bayes-p); we train them by gradient descent. Why, then, is the Bayesian model a useful idealization? There are two connections.
The first is heuristic: the analogy just described, $n\to\infty$ (more samples, longer training) $\leftrightarrow$ $\beta\to\infty$ (lower temperature). As we will see, sweeping a parameter like $\beta$ can take a system through phase transitions, and the analogy predicts that training can do the same. It is only an analogy. In particular, when a system is driven *dynamically* through a first-order transition, the transition is delayed by *hysteresis* ([Section 2.5](#25-phases-and-phase-transitions)), so timing predictions from the Bayesian model should be treated with care. A genuinely dynamical formalism for training, the MSRJD path integral, appears in [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft).
The second connection is exact, and is what makes the Bayesian model testable in practice. Replace the batch noise of SGD by constant, isotropic Gaussian noise; the result is called *Langevin dynamics*:
Here $\xi$ is Gaussian noise with mean $0$ and covariance $I$, $\eta$ is the learning rate, and $T$ is a parameter called the *temperature* (or diffusion constant). Each discrete algorithm on the left has a continuous-time limit on the right, obtained by shrinking the step size and scaling up the number of steps: gradient descent becomes *gradient flow*, and Langevin dynamics becomes the *Langevin equation*, a stochastic differential equation. The key fact is that if one starts from any reasonable distribution of parameters and evolves with the Langevin equation (for a sufficiently nice $\CL$), the distribution converges at late times to the Gibbs distribution
(See ([Tong 2012](#bib-tongkinetictheory), Chapter 3) for a physics-focused discussion of this convergence.) A prior can be included by training on the modified loss $\CL(\theta) - \beta^{-1}\log\pi(\theta)$, which reproduces [(2.3)](#bayes-p). Thus one can *engineer* the Bayesian posterior by Langevin training, and this is how theoretical predictions derived from [(2.3)](#bayes-p) are usually tested. Note that the Langevin route gives a different Bayes-to-gradient connection from the first: here training time $t\to\infty$ at fixed $T$ gives the posterior at $\beta = 1/T$, while $T\to 0$ at fixed $t$ gives gradient flow.
**Remark.** SGD itself is more complicated than Langevin dynamics. The covariance $\Sigma(\theta)$ of its effective noise depends on the position in the loss landscape, and at late times SGD need not converge to any stationary distribution. Despite many qualitative and empirical similarities between Bayesian learning and SGD/Adam training, the precise relation between them remains an open problem. We will compare the two repeatedly: for wide networks in [Section 4](#4-the-ntk-and-large-width-during-training) (gradient flow) versus [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning) (Bayesian); for toy models of grokking in [Section 5](#5-mean-field-scaling) ([Rubin et al. 2024](#bib-rubin2024grokking)); and, for toy models of superposition ([Elhage et al. 2022](#bib-elhage2022superposition)) through the lens of singular learning theory, in ([Chen et al. 2023](#bib-chen2023dynamical)).
### 2.2 Why the exponential of a loss is a likelihood
The Bayesian model treats $e^{-\beta\CL(\theta)}$ as the probability of the training data given the parameters. For the loss functions actually used in deep learning, this is not a metaphor: the losses were *designed* as negative log-likelihoods, so that minimizing the loss is maximum-likelihood estimation. Two examples cover most of practice.
*Mean squared error.* Suppose the labels are the network's output plus Gaussian noise, $y = f_{\theta}(x) + \epsilon$ with $\epsilon\sim\CN(0,\sigma^{2})$. Then the probability density of a label given the input and the parameters is
$$
p(y\,|\,x,\theta) = \frac{1}{\sqrt{2\pi\sigma^{2}}}\,\exp\Big[-\frac{1}{2\sigma^{2}}\big(y-f_{\theta}(x)\big)^{2}\Big]\,,
$$
and for $n$ independent samples,
The MSE loss is the negative log-likelihood of a Gaussian noise model, and the inverse temperature is the inverse noise variance. Noisier labels mean a hotter posterior, which constrains the parameters less.
*Cross-entropy.* For classification, the network outputs a probability distribution $p_{\theta}(y\,|\,x)$ over $K$ classes, typically the softmax of its $K$ output values. The cross-entropy loss is $\ell = -\log p_{\theta}(y\,|\,x)$, so that
*exactly*, with $\beta = 1$: the exponential of the cross-entropy loss *is* the likelihood, with no noise model to choose. The same is true of the next-token loss of a language model, which is exactly the log-probability the model assigns to the training corpus.
*The general pattern.* Whenever a loss is a negative log-likelihood, $e^{-\CL}$ is the likelihood and [(2.3)](#bayes-p) at $\beta=1$ is the textbook Bayesian posterior. A value $\beta\neq 1$ gives a *tempered* posterior. It arises naturally when the noise level is not known (as in [(2.9)](#mse-likelihood)), and it is also a useful theoretical knob: as the interpretation of $n\beta$ above shows, temperature and sample size are interchangeable in the Bayesian model. There is also a useful limiting case. Some losses are hard constraints, $\ell = 0$ if the prediction is right and $\ell=\infty$ if it is wrong. Then $e^{-\ell}$ is a step function, the posterior is uniform on the parameters that fit every training sample, and no $\beta$ appears at all. This is the *zero-temperature* limit, and it is the setting of the Ising perceptron in [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning).
### 2.3 Partition functions, free energy, and cumulants
The Bayesian posterior [(2.3)](#bayes-p) looks exactly like a *thermal* system in equilibrium with a heat bath, the basic object of statistical mechanics. The dictionary is:
| **Bayesian learning** | **statistical mechanics** |
| --- | --- |
| parameters $\theta\in\R^{P}$ | microscopic state of a system |
| | (*e.g.* the spins of $P$ electrons) |
| loss $\CL(\theta)$ | energy $E(\theta)$, or Euclidean action $S(\theta)$ |
| posterior $p(\theta)\propto e^{-\beta\CL(\theta)}$ | Boltzmann distribution at temperature $T=1/\beta$ |
| number of samples $n$, or training time | inverse temperature $\beta$ |
| #parameters $P$ | #degrees of freedom $N$ |
To include the prior $\pi(\theta)$, simply redefine the energy as $E(\theta) = \CL(\theta) - \beta^{-1}\log\pi(\theta)$, so that $e^{-\beta E}= e^{-\beta\CL}\pi$. In what follows we write $E(\theta)$ for whatever sits in the exponent, and $N$ for the number of degrees of freedom, as physicists do; for a network, $N=P$.
*The thermodynamic limit.* Statistical mechanics is about systems with very many degrees of freedom, and one studies their behavior as $N\to\infty$: the *thermodynamic limit*. In this limit one distinguishes *extensive* quantities, which grow in proportion to $N$ (the total energy of a gas), from *intensive* ones, which stay finite (the temperature, or the energy per particle). Both notions transfer directly to learning: with $P$ parameters and $n$ samples, the loss [(2.1)](#sec2-loss) is extensive in $n$ and the empirical loss $L_{n}$ is intensive. Many of the approximations in these notes are thermodynamic limits of some kind, sending the number of parameters, the width, or the number of samples to infinity.
*Partition function and free energy.* The normalization constant in [(2.3)](#bayes-p),
is called the *partition function*. In Bayesian terms it is the marginal likelihood of the data. It depends on $\beta$ and on any other hyperparameters of the problem. Its logarithm defines the *free energy*[^1] and the free energy *density*,
There are three reasons to care about $F$ rather than $Z$. First, $F$ is extensive and $\varphi$ intensive: the free energy density has a limit as $N\to\infty$, while $Z$ itself is exponentially large or small. Second, "$F$ is minimized": at large $N$ the integral [(2.11)](#z-def) is dominated by the regions of parameter space with the smallest local free energy, so the system's macroscopic state is found by minimizing $F$. [Section 2.4](#24-a-physics-example-the-ising-model) makes this precise. Third, derivatives of $F$ compute expectation values, as follows.
*Sources and cumulants.* Suppose we care about the expectation values of some *observables* $\CO_{i}(\theta)$, functions of the parameters. Add *source terms* $J^{i}$ to the exponent,
$$
Z(\beta,J) = \int d\theta\; e^{-\beta E(\theta) + \sum_i J^i\,\CO_i(\theta)}\,,\qquad F(\beta,J) = -\log Z(\beta,J)\,.
$$
Differentiating under the integral,
and in general the $p$-th derivative is the $p$-th *cumulant*,
$$
\kappa_{p}(\CO_{i_1},\ldots,\CO_{i_p}) = -\frac{\pd^{p} F}{\pd J^{i_1}\cdots\pd J^{i_p}}\bigg|_{J=0}\,:
$$
the mean, the covariance, the third central moment, and from $p=4$ on, combinations of moments that vanish for Gaussian variables. The sourced free energy is the generating function of cumulants. Cumulants will return in field-theory clothing in [Section 3.3](#33-gaussian-processes), where they measure how far a network at initialization is from a Gaussian process.
If the only observable we care about is the energy itself, no source is needed: $\beta$ already couples to $E$, and
We will meet these derivatives again shortly: discontinuities in them are how phase transitions are defined.
**Exercise 2.1 (★) (Partition functions, sources, and cumulants).** Take a Boltzmann distribution $p_{\beta}(\theta)\propto e^{-\beta E(\theta)}$, an observable $\CO(\theta)$, and the sourced partition function $Z(\beta,J) = \int d\theta\, e^{-\beta E(\theta)+J\CO(\theta)}$, with $F = -\log Z$. All expectations are in $p_{\beta}$, that is, at $J=0$.
**(a)** Show that $\E(\CO) = -\pd_{J}F|_{J=0}$ and $\mathrm{Var}(\CO) = -\pd_{J}^{2}F|_{J=0}$.
**(b)** The third cumulant is defined by continuing the pattern, $\kappa_{3}(\CO) := -\pd_{J}^{3}F|_{J=0}$. Compute it in terms of the moments $\E(\CO^{3})$, $\E(\CO^{2})$, $\E(\CO)$, and show that it equals the third *centered* moment $\E[(\CO-\E\CO)^{3}]$.
**(c)** Show that for a Gaussian $\CO$, $\log Z(J)$ is exactly quadratic in $J$, so that all cumulants with $p\geq 3$ vanish. Compute $\kappa_{4}$ in terms of centered moments and check that it is *not* the fourth centered moment. (Moral: cumulants measure the failure of higher moments to be determined by lower ones in the Gaussian pattern.)
**(d)** Verify [(2.16)](#e-cumulants).
*Quenched averages.* In a learning problem there is a second source of randomness besides the posterior: the training set itself. Quantities like the free energy depend on which $n$ samples were drawn, and at large $n$ one expects them to concentrate around their average over draws. Physicists call randomness that is fixed before the system equilibrates *quenched disorder*, by analogy with impurities frozen into an alloy when it is quenched in water; the training data are quenched disorder for the posterior. The *quenched free energy* is
and it is the free energy, not the partition function, that should be averaged, because it is the free energy that is extensive and whose derivatives give expectation values. Averaging $\log Z$ is usually hard, and the *annealed approximation* swaps the order,
$$
\E_{\rm data}\log Z \;\approx\; \log\E_{\rm data}Z\,,
$$
which is an upper bound by Jensen's inequality and often a good approximation. We will use it in [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning). The exact treatment goes through the *replica method*, sketched in [Appendix A](#a-the-replica-method); the classic reference is ([Engel & Van den Broeck 2001](#bib-engel2001statistical)).
### 2.4 A physics example: the Ising model
The analogy with thermal physics opens up a bag of well-developed approximation methods. We illustrate two of them, the *mean-field approximation* and the *saddle-point approximation at large $N$*, on the most famous model in statistical mechanics. For the physics of this example, well beyond what we cover here, see ([Tong 2012](#bib-tongsft), Chapter 1).
*The model.* Let $\Lambda$ be a graph with $N$ vertices and undirected edges $\langle ij\rangle$; typically $\Lambda$ is a regular lattice in $D$ dimensions, and every vertex has the same number $d$ of neighbors ($d=2D$ for a hypercubic lattice). At each vertex $i$ sits a *spin* $s_{i}\in\{\pm 1\}$. A state of the system is a choice of all $N$ spins, so there are $2^{N}$ states, and the energy of a state is
The *coupling* $J>0$ rewards neighboring spins for aligning with each other, and the *external magnetic field* $B$ rewards spins for aligning with it. At inverse temperature $\beta$ the Boltzmann distribution is $p_{\beta}(s) = e^{-\beta E(s)}/Z$ with $Z(\beta,J,B) = \sum_{s\in\{\pm1\}^N}e^{-\beta E(s)}$. Dimensionally, only the combinations $\beta J$ and $\beta B$ matter.
The simplest macroscopic quantity describing a state is its *magnetization*, the average spin,
It is our first example of an *order parameter*: a single number that summarizes an $N$-dimensional microscopic state for the purpose of the questions we want to ask. The question here is whether the spins align ($|m|$ close to $1$, an *ordered* phase) or not ($m$ near $0$, a *disordered* phase), and how this depends on $\beta$, $J$, and $B$.
*Mean-field approximation.* We would like to push the distribution $p_{\beta}(s)$ over $2^{N}$ states forward to a distribution over the single variable $m$. The interaction term stands in the way, since it couples each spin to its particular neighbors. The mean-field approximation replaces the neighbors of each spin by their average: suppose that each spin feels only the *average* field of all the others, so that in the neighbor sum we may replace each $s_{j}$ by $m$,
$$
\sum_{\langle ij\rangle}s_{i} s_{j} = \frac{1}{2}\sum_{i} s_{i} \underbrace{\sum_{j\sim i} s_j}_{\approx\, d\,m}\;\approx\; \frac{d\,m}{2}\sum_{i} s_{i} = \frac{N d}{2}\,m^{2}\,.
$$
(The factor $\frac{1}{2}$ compensates for counting each edge twice.) The energy becomes a function of $m$ alone,
This is an uncontrolled approximation for a lattice in low dimension, where a spin's few neighbors fluctuate along with it, but it turns out to be qualitatively right for $D\geq 2$ and quantitatively right for $D\geq 4$ ([Tong 2012](#bib-tongsft), Chapter 1). It becomes exact for a model in which every spin couples weakly to every other; that model, the Curie–Weiss magnet, returns in [Section 5.1](#51-mean-field-theories-in-physics) as the prototype of mean-field limits of neural networks.
*Entropy, and the reduction to one variable.* With the energy depending only on $m$, the sum over states can be organized by the value of $m$:
$$
Z \;\approx\; \sum_{m}\Omega(m)\; e^{N\beta\left(\frac{dJ}{2}m^2 + Bm\right)}\,,\qquad \Omega(m) := \#\{s : m(s) = m\} = \binom{N}{\frac{1+m}{2}N}\,,
$$
since a state with magnetization $m$ has exactly $\frac{1+m}{2}N$ spins up. Stirling's formula gives
the *entropy* per site at magnetization $m$: the log of the number of microscopic states compatible with the macroscopic value $m$. It is largest at $m=0$, where $h = \log 2$ and almost all of the $2^{N}$ states live, and it vanishes at $m=\pm1$, where exactly one state lives. Replacing the sum over $m$ by an integral (the spacing is $2/N$), the partition function becomes a one-dimensional integral,
The function $s(m)$ is called the effective *action*, or the free-energy density at fixed $m$: energy minus temperature times entropy, per site. *We have reduced the problem from $N$ variables to one.* This is the general shape of every calculation in this section: push the microscopic distribution forward to an order parameter, at the price of an entropy term that counts how many microstates sit at each value.
*Saddle points.* At large $N$ the integrand $e^{-N\beta s(m)}$ is exponentially peaked at the *minima* of $s$, and the integral is dominated by a window of width $\sim 1/\sqrt{N}$ around the lowest one. This is the saddle-point approximation (also called Laplace's method), and it gives
$$
\varphi(\beta,J,B) = -\frac{1}{N}\log Z \;\xrightarrow{N\to\infty}\; \beta\,\min_{m} s(m)\,.
$$
The free energy is minimized, in exactly the sense promised in [Section 2.3](#23-partition-functions-free-energy-and-cumulants). The minimizer $\bar m$ is the magnetization the system actually displays. Setting $\pd s/\pd m = 0$ and using $h'(m) = -\tanh^{-1}(m)$,
$$
dJ\,\bar m + B = \beta^{-1}\tanh^{-1}(\bar m)\,,\qquad\text{that is,}\qquad \boxed{\;\bar m = \tanh\big(\beta\,(dJ\,\bar m + B)\big)\;}\,.
$$
This transcendental equation is best solved graphically ([Figure 1](#fig-ising-tanh)). Depending on $\beta$, $J$, and $B$ it has one solution, or three; when there are three, the outer two are minima of $s$ and the middle one is a maximum.
Graphical solution of the mean-field equation [(2.27)](#ising-sc): intersections of $\tanh\beta(dJ\bar m+B)$ with the diagonal. Left: high temperature, one solution at $\bar m=0$. Middle: low temperature at $B=0$, three solutions, of which $\bar m=0$ is a maximum of $s$ and $\pm\bar m$ are degenerate minima. Right: low temperature with $B>0$, three solutions, with the positive one the global minimum.
*Self-consistency (an aside).* There is another route to [(2.27)](#ising-sc), which is the form in which mean-field theory usually appears and which we will use in [Section 5](#5-mean-field-scaling). Look at the effective field felt by a single spin $s_{i}$:
$$
B^{\rm eff}_{i} = -\frac{\pd E}{\pd s_{i}}= J\sum_{j\sim i}s_{j} + B \;\approx\; dJ\,m + B\,.
$$
A single spin in this field has partition function $Z_{i} = \sum_{s_i=\pm1}e^{\beta B^{\rm eff}s_i}= 2\cosh\beta(dJm+B)$ and mean $\E(s_{i}) = \beta^{-1}\pd_{B}\log Z_{i} = \tanh\beta(dJm+B)$. If the approximation is to be consistent, the average that comes out must equal the average we put in, $m = \E(s_{i})$, and that is [(2.27)](#ising-sc). Solve the one-body problem in the average environment, then demand that the average be reproduced: this is the prototype of every mean-field theory.
### 2.5 Phases and phase transitions
Now we analyze the structure of the minima of $s(m)$, which control where the probability concentrates, as $\beta$, $J$, and $B$ are varied. Taylor-expanding the entropy, $h(m) = \log 2 - \frac{m^{2}}{2}- \frac{m^{4}}{12}- \ldots$, the action is
This expansion of a free energy in powers of an order parameter is called a *Landau theory*, and the qualitative analysis below depends only on its form, not on the microscopic model behind it.
*The ordered phase and a first-order transition.* Take $dJ > \beta^{-1}$: the coupling beats the temperature, the quadratic coefficient in [(2.29)](#landau) is negative, and $s(m)$ is a double well ([Figure 2](#fig-ising-wells)). The system wants to be ordered, with neighboring spins aligned; but aligned which way? The field decides. For $B<0$ the global minimum is near $m=-1$, for $B>0$ near $m=+1$, and at $B=0$ the two minima are exactly degenerate. As $B$ is swept through $B_{c} = 0$, the magnetization the system displays jumps from $\bar m\approx -1$ to $\bar m\approx +1$. This is a *first-order phase transition*.
The action $s(m)$ of [(2.25)](#ising-action) in the ordered phase ($\beta dJ = 1.5$), for a negative, zero, and positive external field. The two minima tilt with $B$; the global minimum (orange) switches sides at $B=0$. The suppressed minimum (grey) survives as a metastable state.
The structure of this transition is universal, in the sense that it follows from two isolated minima crossing and not from any detail of the model. Consider the expected magnetization, which by [(2.14)](#r-cumulants) is a derivative of the free energy with respect to the field:
$$
\E(m) = \frac{1}{N}\sum_{i}\E(s_{i}) = -\frac{1}{\beta}\frac{\pd\varphi}{\pd B}\,.
$$
At finite $N$, the two minima of $s$ are two competing *phases*, and $\E(m)$ is a probability-weighted average of their magnetizations. As $N\to\infty$ the saddle-point approximation puts all the weight on the lower minimum, so $\E(m)$ develops a discontinuity at $B_{c}$, and $\varphi$ itself develops a kink ([Figure 3](#fig-ising-firstorder), left). *A discontinuity in a first derivative of the free-energy density, in the thermodynamic limit, is the definition of a first-order transition.* Its second derivative, the *susceptibility*
$$
\chi := \frac{\pd\,\E(m)}{\pd B}= -\frac{1}{\beta}\frac{\pd^{2}\varphi}{\pd B^{2}}= N\beta\,\mathrm{Var}(m)\,,
$$
is a peak whose height grows like $N$ and whose width shrinks like $1/N$: a nascent delta function ([Figure 3](#fig-ising-firstorder), right). Genuine non-analyticity can occur only at $N=\infty$: at finite $N$ the partition function is a finite sum of analytic functions of $B$, so everything is smooth, and what one sees is a rapid *crossover* that sharpens as $N$ grows. Watching the peak in $\chi$ grow with $N$, rather than saturate, is the standard numerical diagnostic that separates a true transition from a crossover. [Exercise 2.2](#ex-anatomy) works out the shapes of these curves for a general first-order transition; they are worth remembering whenever "emergent capabilities" in models of growing size are under discussion ([Section 7.6](#76-emergence-smooth-laws-from-staggered-cliffs)).
Fingerprints of a first-order transition, in the mean-field Ising model at $\beta dJ = 1.5$, computed from the one-dimensional integral [(2.25)](#ising-action) at finite $N$. Left: the expected magnetization as a function of the field is a sigmoid of width $\propto 1/N$, sharpening into a step. Right: the susceptibility $\chi = \pd_{B}\E(m)$ is a peak of height $\propto N$ and width $\propto 1/N$.
*Metastability and hysteresis.* The suppressed minimum in [Figure 2](#fig-ising-wells) does not disappear when it stops being the global one. A system that starts there and evolves by local dynamics (flipping one spin at a time) must climb over the barrier to reach the true minimum, which takes a time exponential in $N$. In practice the system stays in the *metastable* phase until the field is pushed well past $B_{c}$, and the transition observed dynamically *lags* the thermodynamic one. This is *hysteresis*, familiar from real magnets. It is the reason the "$\beta\leftrightarrow$ training time" analogy of [Section 2.1](#21-the-bayesian-model-of-learning) must be handled with care at first-order transitions, and we will meet it again in the guise of delayed generalization (grokking) in [Section 5.7](#57-bayesian-grokking-a-first-order-transition-alternative-reading).
*A second-order transition.* An *$n$-th order* phase transition is a discontinuity in the $n$-th derivative of the free-energy density as $N\to\infty$. The Ising model has a second-order transition too, at fixed $B=0$, as the temperature $T=\beta^{-1}$ crosses $T_{c} = dJ$ ([Figure 4](#fig-ising-secondorder)). Here nothing crosses anything: the single basin at $m=0$ flattens and splits into two symmetric ones. The quadratic coefficient in [(2.29)](#landau) changes sign at $T_{c}$, and the minima sit at
so the order parameter $\bar m \propto (T_{c} - T)^{1/2}$ moves *continuously* but not smoothly; $\varphi$ and $\pd_{T}\varphi$ are continuous, and the discontinuity sits in $\pd_{T}^{2}\varphi$ (the heat capacity). The point $(T,B) = (T_{c},0)$ is a *critical point*. Below it the two ordered states $\pm\bar m$ are exact mirror images, and the system must pick one: this is *spontaneous symmetry breaking*, the $\mathbb{Z}_{2}$ symmetry $s\to -s$ of the energy at $B=0$ being broken by the state. The exponent $\frac{1}{2}$ in [(2.32)](#ising-exponent) is an example of a *critical exponent*. Mean-field theory predicts it in every dimension, and it is wrong in low dimension: the exact value is $\frac{1}{8}$ in $D=2$ (Onsager), and there is no transition at all in $D=1$. Mean-field exponents become exact for $D\geq 4$. Understanding why is the subject of the renormalization group, for which see ([Tong 2012](#bib-tongsft)).
The second-order transition of the mean-field Ising model at $B=0$. Left: the action $s(m)$ for temperatures above, at, and below $T_{c}= dJ$; a single basin flattens and splits into two symmetric ones. Right: the resulting magnetization $\bar m(T)$, continuous at $T_{c}$ but with the square-root singularity [(2.32)](#ising-exponent) (dashed).
*Back to learning.* In the Bayesian model of [Section 2.1](#21-the-bayesian-model-of-learning), the natural control parameter is $n\beta$, the amount of data or training. The Ising analogue is sweeping $\beta$ at fixed $B=0$, which sees the second-order transition. First-order transitions in learning arise when two qualitatively different families of solutions compete, one entropically favored and one favored by the data; the Ising perceptron of [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning) is the cleanest example. Before turning to it, the following exercise establishes what *every* first-order transition looks like at finite $N$.
**Exercise 2.2 (★) (Anatomy of a first-order transition).** Consider any system whose partition function has the form $Z(B) = \int d\theta\, e^{-N E(\theta,B)}$, depending on *some* hyperparameter $B$. Here $B$ could be the inverse temperature, an external field, or the load $\alpha$ of [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning); we have absorbed $\beta$ into $E$ to keep the setup general. Suppose that $E(\theta,B)$ has exactly two isolated minima, at $\theta_{1}$ and $\theta_{2}$. Then at large $N$ the integral is dominated by two competing phases. Letting $U_{1},U_{2}$ be small neighborhoods of $\theta_{1},\theta_{2}$,
with $\varphi_{i}(B) = -\frac{1}{N}\log\int_{U_i}d\theta\,e^{-NE(\theta,B)}= E(\theta_{i},B) + o(1)$ the free-energy densities of the two basins. Suppose further that as $B$ crosses a critical value $B_{c}$ the two minimum values cross: $\varphi_{1}(B_{c}) = \varphi_{2}(B_{c})$ with $\varphi_{1}'(B_{c})\neq\varphi_{2}'(B_{c})$, where a prime denotes $\pd_{B}$.
**(a)** *Occupation.* Show that the probability of finding the system in basin $i$ is
a logistic function of $N\,\Delta\varphi'\,\Delta B$. Over what window of $B$ does the system switch allegiance?
**(b)** *First derivative.* Show from [(2.33)](#two-branch-z) that
the probability-weighted average of the branch slopes. Sketch it at finite $N$ and at $N=\infty$: a sigmoid of height $\Delta\varphi'$ sharpening into a jump. The kink in $\varphi$ itself is the defining discontinuity of a first-order transition. Show also, directly from the integral, that $\pd_{B}\varphi = \E\big[\pd_{B} E(\theta,B)\big]$: the observable that jumps is the one conjugate to the control parameter. When $B=\beta$ it is the energy, and the jump is the *latent heat*; for the Ising model with $B$ the field, it is the magnetization.
**(c)** *Second derivative.* Using $\pd_{B} p_{1} = -\pd_{B} p_{2} = -N p_{1}p_{2}\,\Delta\varphi'$ (verify this), show that
Sketch this curve. Near $B_{c}$, where $p_{1}p_{2}\approx\frac{1}{4}$, the first term is a peak of height $\sim N(\Delta\varphi')^{2}/4$ and width $\sim 1/(N\Delta\varphi')$, which converges to a delta function as $N\to\infty$. Show also that when $\CO(\theta) = \pd_{B}E(\theta,B)$ does not depend on $B$, $-\pd_{B}^{2}\varphi = N\,\mathrm{Var}(\CO)\geq 0$. (For $B=\beta$ this is the heat capacity; for the Ising field it is the susceptibility.)
**(d)** *Plot it.* Take the minimal model $\varphi_{1}(B) = a\,(B-B_{c})$, $\varphi_{2} = 0$ and plot $\varphi$, $\pd_{B}\varphi$, and $-\pd_{B}^{2}\varphi$ for a few values of $N$. (Use a numerically stable $\log(e^{x}+e^{y})$.) Watching the peak grow with $N$ rather than saturate is the standard diagnostic separating a true transition from a smooth crossover.
### 2.6 The Ising perceptron: a first-order transition in learning
We close the physics part of this section with a complete statistical-mechanics analysis of a learning problem, built to run in parallel with the Ising magnet. It is the simplest model we know in which *more data* produces a genuine phase transition to *perfect generalization*, and it goes back to Gyorgyi ([Gy"orgyi 1990](#bib-gyorgyi1990first)); see ([Engel & Van den Broeck 2001](#bib-engel2001statistical), Chapter 7). It is presented as a walk-through exercise; solutions are in [Appendix C](#c-solutions-to-exercises), but try each part before reading them.
*Setup.* A *teacher* is a fixed vector of binary weights $\theta^{*}\in\{\pm1\}^{N}$. Inputs are Gaussian, $x\sim\CN(0,I_{N})$, and the teacher labels them by
$$
y(x) = \mathrm{sign}(\theta^{*}\cdot x)\,.
$$
A *student* is another binary vector $\theta\in\{\pm1\}^{N}$, guessing labels the same way, $f_{\theta}(x) = \mathrm{sign}(\theta\cdot x)$. This model is called the *Ising perceptron*, because its weights take the values of Ising spins. We observe $n$ examples $(x^{\mu},y^{\mu})$ with $y^{\mu} = y(x^{\mu})$, and do Bayesian learning at zero temperature with a uniform prior, in the sense of [Section 2.2](#22-why-the-exponential-of-a-loss-is-a-likelihood): the per-sample loss is $0$ if the student classifies the example correctly and $\infty$ otherwise, so the posterior is uniform on the students that classify *every* example correctly. Since teacher and student agree on $x$ exactly when $y(x)\,(\theta\cdot x)>0$,
$$
p(\theta) \;\propto\; \prod_{\mu=1}^{n}\Theta\big(y^{\mu}\,(\theta\cdot x^{\mu})\big)\,,\qquad Z = \frac{1}{2^{N}}\sum_{\theta\in\{\pm1\}^N}\;\prod_{\mu=1}^{n}\Theta\big(y^{\mu}\,(\theta\cdot x^{\mu})\big)\,,
$$
with $\Theta(u) = 1$ for $u>0$ and $0$ otherwise. With the factor $1/2^{N}$, $Z$ is the fraction of all students consistent with the data. Both $N$ and $n$ will be taken large, at fixed ratio: the *load*
$$
\alpha := \frac{n}{N}\,,
$$
the number of examples per parameter, is the control parameter of the problem.
**Exercise 2.3 (★) (The Ising perceptron).** **(a)** *The order parameter is a magnetization.* Define the *overlap* $R(\theta) = \frac{1}{N}\theta\cdot\theta^{*}\in[-1,1]$ and the alignment variables $s_{i} := \theta_{i}\theta_{i}^{*}\in\{\pm1\}$. Show that $R = \frac{1}{N}\sum_{i} s_{i}$: the overlap *is* the magnetization of the spins $s_{i}$. Learning the teacher means magnetizing them to $R=1$.
**(b)** *Generalization error.* Show that a student at overlap $R$ misclassifies a fresh Gaussian input with probability
$$
\varepsilon(R) = \frac{\arccos R}{\pi}\,.
$$
[*Hint:* project $x$ onto the plane spanned by $\theta$ and $\theta^{*}$; the projection is an isotropic two-dimensional Gaussian, and the two decision boundaries are diameters meeting at angle $\arccos R$. When do the signs disagree?]
**(c)** *Entropy: the same $h$ as the magnet.* Show that the number of students at overlap $R$ is $\binom{N}{\frac{1+R}{2}N}$, so that by Stirling's formula
$$
\begin{aligned}\frac{1}{N}\log\#\{\theta : R(\theta) = R\} \;\to\; h(R) &= \log 2 - \frac{1-R}{2}\log(1-R) \\ &\quad - \frac{1+R}{2}\log(1+R)\,,\end{aligned}
$$
*identical* to the magnet's entropy [(2.24)](#ising-entropy) with $m\to R$.
**(d)** *Push forward to one variable.* Average $Z$ over datasets. The examples are independent, and each is classified correctly by a student at overlap $R$ with probability $1-\varepsilon(R)$, so
$$
\E_{x} Z = \frac{1}{2^{N}}\sum_{\theta}\big(1-\varepsilon(R(\theta))\big)^{n}\,.
$$
Trade the sum over $\theta$ for an integral over $R$ using (c), and show that at fixed load $\alpha$ and large $N$,
This is the exact analogue of the magnet's [(2.25)](#ising-action): an entropy term counting states, and a "field" term rewarding alignment. The aligning field is now the *data*, and its strength is the load $\alpha$. (Averaging $Z$ rather than $\log Z$ is the annealed approximation [(2.18)](#annealed-swap). It preserves all the qualitative structure here and shifts the numbers only slightly; [Appendix A](#a-the-replica-method) explains how to do the exact calculation.)
**(e)** Plot $s(R)$ for several values of $\alpha$ and locate its critical points. Do you expect a phase transition? At which $\alpha$?
**(f)** *The teacher basin.* $R=1$ is the perfectly aligned state, where the student is the teacher. Show that $s(1) = \log 2$ for *every* $\alpha$ (zero entropy, but zero data cost), and that $R=1$ is always a local minimum: writing $\delta = 1-R$, deviating from the teacher gains entropy $\sim\frac{\delta}{2}\log(1/\delta)$ but costs data-fit $\alpha\varepsilon\sim\alpha\sqrt{2\delta}/\pi$ (expand $\arccos$), and $\sqrt{\delta}\gg\delta\log(1/\delta)$. This differs from the magnet, where the aligned state sat at the bottom of a smooth well and was not always a minimum.
**(g)** *The interior basin, the crossing, and the transition.* At $\alpha=0$, $s$ is minimized at $R=0$ with $s(0) = 0 < s(1)$: the disordered cloud of typical students dominates. As $\alpha$ grows, this interior minimum $R^{*}(\alpha)$ drifts toward $1$ and its value rises. Numerically, find the load $\alpha_{c}$ at which the two basins exchange dominance, $s(R^{*}(\alpha_{c})) = s(1)$. You should find $\alpha_{c}\approx 1.45$, with the interior minimum still at $R^{*}\approx 0.70$, that is, generalization error $\varepsilon^{*}\approx 0.25$. Conclude that at $\alpha_{c}$ the typical student jumps from $25\%$ error to *zero* error: a first-order transition to perfect generalization. Then connect to [Exercise 2.2](#ex-anatomy) with $\alpha$ as the control parameter: what do $\pd_{\alpha}\varphi$ and $-\pd_{\alpha}^{2}\varphi$ look like at finite $N$, and what plays the role of the latent heat? Finally, find the load at which the interior minimum disappears altogether, and explain what happens in between in the light of the metastability discussion of [Section 2.5](#25-phases-and-phase-transitions).
**(h)** *The dictionary.* Complete the table:
| **Ising magnet** | **Ising perceptron** |
| --- | --- |
| spins $s_{i}\in\{\pm1\}$ | alignments $s_{i}= \theta_{i}\theta_{i}^{*}\in\{\pm1\}$ |
| magnetization $m$ | overlap $R$ |
| entropy $h(m)$ | entropy $h(R)$ *(the same function)* |
| aligning field $B$ | ? |
| knob swept across the transition | ? |
| ordered phase $\E(m)\approx\pm1$ | ? |
| jump in $\E(m)$ as a function of $B$ | ? |
Also articulate one structural *difference*: in the magnet the two competing minima are symmetry partners; what are they here, and which ingredient of (c) made a sharp transition possible at all? (What would change for spherical weights $\theta\in\R^{N}$, $|\theta|^{2} = N$, where the slice at overlap $R$ has volume $\propto(1-R^{2})^{N/2}$? [Appendix A.3](#a3-a-variant-the-spherical-perceptron) works this case out.)
### 2.7 Comparison with QFT
This short subsection is a translation guide, aimed mainly at readers arriving from high-energy or condensed-matter physics — if the words "action" and "path integral" mean more to you than "posterior" and "marginal likelihood," read on; everyone else may safely skip ahead.
Nearly everything in these notes lives, in the first instance, in *Euclidean equilibrium statistical mechanics*: probability distributions of Boltzmann–Gibbs form $P \propto e^{-\beta E}$, partition functions, free energies. For a field theorist, a partition function with sources,
$$
Z[J] = \int \mathcal{D}\phi\; e^{-S[\phi] \,+ \int J\,\phi}\,,
$$
is a moment generating functional; its logarithm generates connected correlation functions (cumulants); the action $S$ encodes the microscopic rules of the system, and a quadratic action means a free (Gaussian) theory, with everything else treated by perturbing around it. In these notes the role of the field $\phi$ is played variously by the weights $\theta$ of a network, by an order parameter such as the magnetization $m$ or the overlap $R$ of [Sections 2.4](#24-a-physics-example-the-ising-model) and [2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning), or — most usefully — by the network function $f(x)$ itself, with the input $x$ playing the role of the spacetime coordinate. Functional integrals over the network function are developed carefully in [Sections 3.3](#33-gaussian-processes) and [3.5](#35-nngp-via-functional-integrals), and the reader coming from QFT will find many familiar ingredients there dressed up as learning theory: Wick's theorem, sources, cumulants, effective actions.
If your instincts were trained on Lorentzian signature, picture everything here as already Wick-rotated. There is no $i/\hbar$: the weights $e^{-S}$ are real and positive, as after the standard continuation
$$
Z_{\rm QFT}= \int \mathcal{D}\phi\; e^{\frac{i}{\hbar} S[\phi]}\;\;\xrightarrow{\;\tau\,=\,it\;}\;\; Z_{\rm E}= \int \mathcal{D}\phi\; e^{- S_E [\phi]/\hbar}\,.
$$
The role of $\hbar$ — the parameter controlling the size of fluctuations around saddle points — is played here by temperature, or by $1/N$ at large width, or by $1/n$ at large sample size.
When we eventually study the *dynamics* of training as a field theory in [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft), a genuinely real-time formalism does appear — a path integral over trajectories with a causal structure and a response field. That's the closest object in these notes to a Lorentzian field theory, and is precisely the Martin–Siggia–Rose construction of non-equilibrium statistical mechanics.
For readers who would like physics references alongside the learning-theory ones: David Tong's lecture notes ([Tong 2012](#bib-tongstatphysics); [Tong 2012](#bib-tongsft); [Tong 2012](#bib-tongkinetictheory); [Tong 2006](#bib-tongqft)) are a friendly introduction to statistical physics, statistical field theory, kinetic theory, and QFT respectively; Lancaster and Blundell ([Lancaster & Blundell 2014](#bib-lancaster2014quantumgifted)) is a gentle conceptual route into QFT; Altland and Simons ([Altland & Simons 2010](#bib-altlandcondensedmatterfield2010)) bridges toward condensed-matter field theory (and contains the Langevin/MSRJD material used in [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft)); Pathria ([Pathria & Beale 2011](#bib-pathriastatisticalmechanics2011)) is a standard graduate statistical mechanics text; Goldenfeld ([Goldenfeld 2018](#bib-goldenfeldlecturesphasetransitions2018a)) is the classic on phase transitions and the renormalization group; and Zinn-Justin ([Zinn-Justin 2002](#bib-zinn2002quantum)) is the encyclopedic reference for the functional formalism.
## 3. Large width at initialization: the Gaussian process limit
The next few sections are all about different limits in neural networks, in which some parameter(s) become large. Such limits are where physical/statistical analyses really shine. Studying them teaches some basic lessons about how we should (and should not) initialize networks, and how to choose other hyperparameters. It also leads to analytically solvable models that display many properties we see in practice: lazy vs. rich learning, double descent, phase transitions, scaling laws for optimal learning. For further pedagogical references on much of this material, we recommend ([Roberts et al. 2022](#bib-roberts2022principles); [Spiliopoulos et al. 2025](#bib-spiliopoulos2025foundations); [Ringel et al. 2025](#bib-ringel2025applications)) and the online lectures ([Hanin 2026](#bib-hanin2026mdl)).
We'll begin in this section with a classic result ([Neal 1996](#bib-neal1996bayesian); [Lee et al. 2018](#bib-lee2018deep); [Yang 2019](#bib-yang2019tensorprograms1); [Hanin 2023](#bib-hanin2023random)) about increasing the width of a network at fixed depth and training-sample size.
Modern neural networks are more sophisticated than the simple MLP's discussed here, but they share many of the same features. Layer norms solve the problem of trainability at large depth, but the analysis of rich/lazy regimes, feature learning, and hyperparameter scaling is still relevant. See *e.g.* the discussion of $\mu P$ in [Section 5.5](#55-what-is-this-framework-good-for).
### 3.1 Notation
We'll use the following notation for neural networks and training data.
- $\theta\in \R^{P}$ trainable network parameters (weights, biases)
- $(x^{\mu},y^{\mu})_{\mu=1}^{n}$ training data samples, with $x^{\mu}\in \R^{d}$ and $y^{\mu}\in \R$ (unless otherwise stated)
- $f_{\theta}:\R^{d}\to \R$ network function
For MLP's in particular (also called FCN's = fully connected networks)
- $\{z_{i}^{(\ell)}\}_{i=1}^{N_\ell}$ pre-activations in layer $\ell\in \{0,...,L+1\}$
- $L=$ number of hidden layers (depth); $N_{\ell} =$ width of $\ell$'th layer; $N_{0}=d$ and $N_{L+1}=1$. We will often take constant hidden depth, $N_{1}=...=N_{L}=N$.
- $\phi:\R\to \R$ activation function (assumed the same in all layers except the first)
- network recursion $z_{i}^{(\ell)}= \begin{cases}\frac{1}{\sqrt{N_{\ell-1}}}\sum_{j=1}^{N_{\ell-1}}w^{(\ell)}_{ij}\,\phi(z^{(\ell-1)}_{j})+ b^{(\ell)}_{i}&\ell\geq 2 \\[.2cm] \frac{1}{\sqrt{d}}\sum_{j=1}^{d}w^{(1)}_{ij}\,x_{j}+ b^{(1)}_{i}&\ell = 1\end{cases}$
with $z^{(0)}= x$ and $z^{(L+1)}= f$; parameters are $\theta = \{w^{(\ell)},b^{(\ell)}\}_{\ell=1}^{L+1}$ . The explicit factors of $1/\sqrt{N_{\ell-1}}$ are a normalization convention, explained in [Section 3.2](#32-normalizing-the-network): with them in place, all weights and biases can be taken of order one.
### 3.2 Normalizing the network
Consider a fully connected network of depth $L$ with layer widths $N_{\ell}$ as in [Section 3.1](#31-notation). Suppose we keep $L$ fixed but make the $N_{\ell}$ very large. Suppose also that the network is initialized randomly, with all weights and biases drawn independently from normal distributions of *order-one* variance,
with $\sigma_{w},\sigma_{b}$ fixed as the widths grow. (Here and throughout, $\CN(0,\sigma)$ denotes a Gaussian of *variance* $\sigma$. Letting the variances depend on the layer changes nothing below.) The question is how to normalize the sums in the network recursion so that preactivations stay of order one from layer to layer. From
$$
z_{i}^{(\ell)}= \frac{1}{\sqrt{N_{\ell-1}}}\, w_{i}^{(\ell)}\cdot \phi(z^{(\ell-1)}) + b_{i}^{(\ell)}
$$
(where we've written $w_{i}=(w_{i,1},...,w_{i,N_{\ell-1}})$ and $z^{(\ell-1)}= (z^{(\ell-1)}_{1},...,z^{(\ell-1)}_{N_{\ell-1}})$ as vectors) we get
$$
\begin{aligned}\text{Var}(z_{i}^{(\ell)})&= \frac{1}{N_{\ell-1}}\sum_{j=1}^{N_{\ell-1}}\text{Var}\big(w_{ij}^{(\ell)}\big)\, \E\Big[\phi\big(z_{j}^{(\ell-1)}\big)^{2}\Big] \;+\; \text{Var}\big(b_{i}^{(\ell)}\big) \notag \\&= \sigma_{w} \cdot \frac{1}{N_{\ell-1}}\sum_{j=1}^{N_{\ell-1}}\Big( \text{Var}\big(\phi(z_{j}^{(\ell-1)})\big)+\E\big[\phi(z_{j}^{(\ell-1)})\big]^{2}\Big) +\sigma_{b}\,.\end{aligned}
$$
Here we've conditioned on the previous layer, and used that for independent mean-zero $X$ and arbitrary $Y$, $\text{Var}(XY) = \text{Var}(X)\,\E[Y^{2}]$, along with additivity of variance for independent summands. The average in the second line involves only the activation function and the order-one preactivations of the previous layer, so it stays of order one as $N_{\ell-1}\to\infty$. (We'll sharpen this estimate in [Section 3.7](#37-edge-of-chaos).) Thus
The prefactor $1/\sqrt{N_{\ell-1}}$, the square root of the number of inputs (the "fan-in") feeding each neuron, is exactly what keeps the preactivations of order one at large width. We call this the *fan-in normalization*; it goes back to LeCun and collaborators in the 90's ([LeCun et al. 1998](#bib-lecun1998efficient)). In [Section 4](#4-the-ntk-and-large-width-during-training) it will also be called the *NTK normalization*, and in [Section 5](#5-mean-field-scaling) we will meet the alternative *mean-field normalization*, in which the output layer carries $1/N$ instead of $1/\sqrt{N}$.
**Remark (other conventions).** Common practice, and the convention of much of the literature, is to put no explicit prefactors in the network and instead to initialize with width-dependent variances,
At initialization the two conventions describe *exactly* the same distribution of network functions: [(3.5)](#init-fanin-var) is a change of variables on parameter space, and the network function is the same function of $x$. Everything in this section, which concerns initialization only, is insensitive to the choice. During training, however, the two conventions are *not* equivalent: gradient descent is not invariant under rescaling coordinates on parameter space, as we will explain in [Section 4.3](#43-constant-ntk-at-large-width). Keeping every parameter of order one, with the factors of $N$ displayed explicitly in the network function, makes the comparison between different scaling limits transparent, and it is the convention we use throughout these notes. When reading the literature, expect the factors of $N$ to sit in the variances instead, and translate with [(3.5)](#init-fanin-var).
### 3.3 Gaussian processes
A *Gaussian process* $f:\R^{m}\to \R$ with mean function $\mu$ and covariance kernel $K$ is a random function such that on any $p$ inputs $(x^{1},...,x^{p})\in (\R^{m})^{p}$ the vector of function values is Gaussian
$$
\begin{pmatrix}f(x^{1}) \\ \vdots \\ f(x^{p})\end{pmatrix} \sim \CN \left( \begin{pmatrix}\mu(x^{1}) \\ \vdots \\ \mu(x^{p})\end{pmatrix}, \begin{pmatrix}K(x^{1},x^{1})&\cdots&K(x^{1},x^{p}) \\ \vdots&\ddots&\vdots \\ K(x^{p},x^{1})&\cdots&K(x^{p},x^{p})\end{pmatrix} \right)\,.
$$
Here $\mu:\R^{m}\to \R$ and $K:\R^{m}\times \R^{m}\to \R$ are fixed, deterministic functions that control the mean and covariance of any tuple of values of $f$. The function $K$, often called the kernel of the process, is symmetric: $K(x,x')=K(x',x)$.
Let `$dx$' denote a measure on $\R^{m}$ (it need not be the usual Lebesgue measure). If $K$, as integral kernel, has an inverse $K^{-1}$ satisfying
$$
\int dx' K^{-1}(x,x') K(x',x'') = \delta(x-x'')
$$
then there's also a "physics definition" of a Gaussian process. It's not entirely rigorous, as it invokes multi-dimensional functional integrals. Mathematically, it's best thought of as a heuristic tool, which comes with a package of axiomatic manipulations whose actual results can (and should) all be fully justified by other means. The manipulations are all quite intuitive, though, and familiar to quantum field theorists, so we mention this formalism here. The idea is that, just like a single Gaussian variable `$z$' has a probability measure
$$
P(z) = \frac{1}{Z}e^{-\frac{1}{2}K^{-1}(z-\mu)^2}\,,\qquad Z = \int dz\,e^{-\frac{1}{2}K^{-1}(z-\mu)^2}= \sqrt{2\pi K}\,,
$$
a Gaussian process $f$ has a (non-rigorous) measure
$$
\begin{aligned}P[f]&= \frac{1}{Z}e^{\textstyle-\frac{1}{2} \int dx\,dx'\, (f(x)-\mu(x)) K^{-1}(x,x')(f(x')-\mu(x'))}\,, \notag \\ Z&\; \text{`='}\int Df\, e^{\textstyle-\frac{1}{2} \int dx\,dx'\, (f(x)-\mu(x)) K^{-1}(x,x')(f(x')-\mu(x'))}\,.\end{aligned}
$$
If the kernel and the measure theory on the input space $\R^{m}$ are nice enough to allow $K$ to be diagonalized, with a discrete basis of orthonormal eigenfunctions $\{\phi_{k}(x)\}_{k=1}^{\infty}$ and eigenvalues $\lambda_{1}\geq \lambda_{2}\geq\ldots > 0$ (so that $K^{-1}$ has eigenvalues $1/\lambda_{k}$), then $P[f]$ can be understood as an infinite-dimensional Gaussian. Writing $f(x)-\mu(x) = \sum_{k=1}^{\infty} a_{k} \phi_{k}(x)$,
$$
P[f] = \frac{1}{Z}e^{\textstyle -\frac{1}{2}\sum_{k=1}^\infty a_k^2/\lambda_k}\,,\qquad Z = \prod_{k=1}^{\infty} (2\pi \lambda_{k})^{1/2}\,.
$$
One virtue of the physics formalism is that it organizes expectation values of products of $f(x)$'s. Let $\mu=0$ for simplicity (to generalize to any $\mu$, just shift $f\mapsto f-\mu$ below). Let $J:\R^{m}\to \R$ be an arbitrary (nice...) function and define
The sum is over partitions of $\{1,...,p\}$ into pairs, and vanishes if $p$ is odd. Eqn. [(3.14)](#gp-prods) is called Wick's Theorem in QFT, Isserlis's Theorem in probability (and it's correct, despite the non-rigorous derivation here).
**Exercise 3.1 (Wick's theorem).** If you haven't seen this before, derive [(3.12)](#zj-eval) by completing the square in [(3.11)](#zj), and then derive [(3.14)](#gp-prods) by differentiating. Check the first two cases explicitly: $\E[f(x^{1})f(x^{2})] = K(x^{1},x^{2})$, and $\E[f(x^{1})f(x^{2})f(x^{3})f(x^{4})] = K^{12}K^{34}+K^{13}K^{24}+K^{14}K^{23}$ (with $K^{ij}=K(x^{i},x^{j})$). How many terms appear at $p=6$?
Another virtue of the physics formalism is that it neatly connects *higher cumulants* to *deviations from Gaussianity*. In general, the $p$-th cumulant $\kappa_{p}(X_{1},...,X_{p})$ of random variables $X_{1},...,X_{p}$ is defined recursively via
$$
\E(X_{1}\cdots X_{p}) = \sum_{\text{partition $\Pi$ of $\{1,...,p\}$}}\;\prod_{\text{blocks $\{i_{1},...,i_{k}\}\subset \Pi$}}\kappa_{k}(X_{i_1},...,X_{i_k})
$$
For example, $\kappa_{1}(X)=\E(X)$, $\kappa_{2}(X_{1},X_{2})=\text{Cov}(X_{1},X_{2}) = \E(X_{1}X_{2})-\E(X_{1})E(X_{2})$. In physics, the $\kappa_{p}$'s are called connected correlation functions; they measure the part of a joint expectation that's not captured by expectations of subsets of the variables involved. For a Gaussian process,
$$
\kappa_{2}(f(x^{1}),f(x^{2}))=K(x^{1},x^{2})\,;\qquad \kappa_{p}(f(x^{1}),...,f(x^{p})) = 0\quad\text{for all $p>2$}\,,
$$
as can easily be seen directly from [(3.14)](#gp-prods), where joint expectations for $p>2$ all factor into products of pairs. Thus, nonvanishing higher $\kappa_{p}$'s are a measure of non-Gaussianity.
This appears in the functional-integral formalism as follows. Suppose that $f:\R^{m}\to \R$ is a random function that's governed by a functional probability distribution
$$
P[f] = \frac{1}{Z}\,e^{\textstyle-S[f]}\,,
$$
for some functional $S:\{\text{functions }f\}\to \R$. (In physics $S$ is called the action.) Then just as in [(3.11)](#zj), we can introduce the generating function
$$
\begin{aligned}Z[J]&:= \int Df\, e^{\textstyle-S[f]+\int dx\, J(x)f(x)}\notag \\ \Rightarrow\quad \E\big[f(x^{1})...f(x^{p})\big]&= \frac{1}{Z[0]}\frac{\delta}{\delta J(x^{1})}...\frac{\delta}{\delta J(x^{p})}Z[J]\bigg|_{J=0}\end{aligned}
$$
The cumulants appear naturally as derivatives of $\log Z[J]$, otherwise known as the (sourced) free energy:
$$
\kappa_{p}(f(x^{1}),...,f(x^{p})) = \frac{\delta}{\delta J(x^{1})}...\frac{\delta}{\delta J(x^{p})}\log Z[J] \bigg|_{J=0}\,,
$$
which in turn means (formally) that
$$
Z[J] = \exp\sum_{p=1}^{\infty} \frac{1}{p!}\int dx^{1}...dx^{p}\, \kappa_{p}(f(x^{1}),...,f(x^{p}))\,J(x^{1})...J(x^{p})\,,
$$
a generating function of cumulants. The original $S[f]$ is the inverse Legendre transform of $\log Z[J]$ (by definition of $Z[J]$). Expanded as a formal series in the cumulants (assuming they are small), and assuming $\E[f(x)]=0$ for simplicity, it takes the form
$$
\begin{aligned}S[f] =&\; \frac{1}{2}\int dx^{1}\, dx^{2}\, \kappa_{2}^{-1}(x^{1},x^{2}) f(x^{1})f(x^{2}) \notag \\&- \int dx^{1} dx^{2} dx^{3}\, dy^{1} dy^{2} dy^{3}\; \kappa_{2}^{-1}(x^{1},y^{1})\,\kappa_{2}^{-1}(x^{2},y^{2}) \notag \\&\hspace{1in}\times\, \kappa_{2}^{-1}(x^{3},y^{3})\,\kappa_{3}(y^{1},y^{2},y^{3})\,f(x^{1})f(x^{2})f(x^{3}) \notag \\&- \int dx^{1} dx^{2} dx^{3} dx^{4}\, (\kappa_{2}^{-1})^{4}\big[ \kappa_{4} - \kappa_{3}\kappa_{2}^{-1}\kappa_{3}\;\;\text{terms}\big] \notag \\&\hspace{1in}\times\, f(x^{1})f(x^{2})f(x^{3})f(x^{4})\; - \ldots\end{aligned}
$$
Physicists might recognize this as a sum of Feynman diagrams, with cumulants giving the vertices. A key structure is that the $p$-th order term in the effective action contains only $\kappa_{p}$ and lower cumulants. In particular, $f$ is a Gaussian process if and only if $S[f]$ is quadratic; and small deviations from Gaussianity, measured by cumulants, provide the higher-order terms in $S[f]$.
### 3.4 The NNGP correspondence
A classic result in the theory of neural networks is that, at fixed depth and large width, with the fan-in normalization of [Section 3.2](#32-normalizing-the-network), the network function at initialization converges to a Gaussian process. The theorem appeared in the 90's ([Neal 1996](#bib-neal1996bayesian)) for a single hidden layer, and was generalized more recently to deeper networks and different architectures, *cf.* ([Lee et al. 2018](#bib-lee2018deep); [Yang 2019](#bib-yang2019tensorprograms1); [Hanin 2023](#bib-hanin2023random)).
In this theorem, the source of stochasticity in the network function comes from the random initializations of network parameters. One can think of randomly initializing a swarm of networks and asking what sorts of values the network function outputs. The result is significant because of its universality, and because it feeds into (and greatly simplifies) further analyses of training.
The derivation proceeds via induction. We'll just sketch it here, for MLP's.
[^2]
For simplicity we'll assume that biases have been set to zero, and that all weights have the same order-one variance $\sigma_{w}$, as in [(3.1)](#init-o1), with the fan-in normalization of [Section 3.2](#32-normalizing-the-network).
In the first network layer, the preactivations are
$$
z_{i}^{(1)}(x) = \frac{1}{\sqrt{d}}\sum_{j=1}^{d} w_{ij}^{(1)}x_{j}\,,\qquad w_{ij}^{(1)}\sim \CN(0,\sigma_{w})\,.
$$
At fixed $x$, each $z_{i}^{(1)}(x)$ is a sum of $d$ i.i.d. Gaussian variables (with coefficients $x_{j}$), and thus is Gaussian itself, with $\E[z_{i}^{(1)}] = 0$ and
$$
\text{Cov}\big(z_{i}^{(1)}(x),z_{i}^{(1)}(x')\big) = \frac{1}{d}\sum_{j,j'=1}^{d} \underbrace{\text{Cov}(w_{ij},w_{ij'})}_{\delta_{jj'} \sigma_w}x_{j} x'_{j'}= \sigma_{w}\frac{x\cdot x'}{d}\,,
$$
thus $z_{i}^{(1)}$ is a Gaussian process (GP) with
$$
z_{i}^{(1)}\sim \CN(0,K^{(1)})\,,\qquad K^{(1)}(x,x') = \sigma_{w}\frac{x\cdot x'}{d}\,.
$$
In the second layer, we've got
Thus we again get a Gaussian process. The problem is that if we don't condition on the preactivations, the kernel of this GP is not deterministic, but becomes a random variable itself: $\phi(z^{(1)})\cdot \phi(z^{(1)}{}')$. This is where the infinite-width limit comes in. Heuristically, the central limit theorem implies that the sum of products of (non-Gaussian!) random variables on the LHS of [(3.25)](#z2-rec) will become Gaussian as $N_{1}\to\infty$. More concretely, we can calculate and bound the variance
$$
\begin{aligned}\text{Var}\big(K^{(2)}(z^{(1)},z^{(1)}{}')\big)&= \frac{\sigma_{w}^{2}}{N_{1}^{2}}\text{Var}\big(\phi(z^{(1)})\cdot \phi(z^{(1)}{}')\big) \\&= \frac{\sigma_{w}^{2}}{N_{1}^{2}}N_{1} \cdot \text{Var}\big(\phi(z_{1}^{(1)}) \phi(z_{1}^{(1)}{}')\big) \to 0 \quad\text{as $N_{1}\to\infty$}\end{aligned}
$$
as long as the first-layer preactivations are order one (independent of $N_{1}$), and the activation function is reasonable (any standard choice will be). Thus, as $N_{1}\to\infty$, we can replace the stochastic kernel $K^{(2)}$ by its expectation,
$$
\bar K^{(2)}(x,x') := \sigma_{w} \, \E_{z_1^{(1)}(x),z_1^{(1)}(x')\sim \CN(0,K^{(1)})}[ \phi(z_{1}^{(1)}(x)) \phi(z_{1}^{(1)}(x')) ]
$$
Each subsequent layer is the same. At infinite width, the kernels become deterministic, and obey a recursion
with $\bar K^{L+1}(x,x')$ ultimately governing a GP for the network function $f(x)$.
**Remark (exact kernels).** For a few special activation functions, the Gaussian expectation in [(3.30)](#kbar-rec) can be evaluated in closed form, giving explicit NNGP kernels. Write $K_{11}=\bar K^{(\ell-1)}(x,x)$, $K_{22}=\bar K^{(\ell-1)}(x',x')$ and $K_{12}=\bar K^{(\ell-1)}(x,x')$ for the previous layer's kernel values. For ReLU, $\phi(u)=\max(u,0)$, one finds the *arc-cosine kernel* of Cho and Saul ([Cho & Saul 2009](#bib-cho2009kernel)),
Iterating either formula through the layers gives the exact infinite-width kernel of a deep ReLU or erf network in a few lines of code. (Deriving [(3.31)](#relu-kernel)–[(3.32)](#erf-kernel) is a pleasant but somewhat lengthy Gaussian-integral workout.)
**Exercise 3.2 (The NNGP correspondence, empirically (lab)).** Fix two unit-norm inputs $x, x'\in\R^{d}$ with, say, $d=8$ and $x\cdot x' = 0.5$. For widths $N \in \{3, 30, 300\}$, sample $M = 10^{4}$ independent initializations of the one-hidden-layer erf network $f(x) = \frac{1}{\sqrt{N}}\sum_{i=1}^{N} c_{i}\, \mathrm{erf}\big(a_{i}\cdot x/\sqrt{d}\big)$ with $c_{i}, a_{ij}\sim\CN(0,1)$ (that is, $\sigma_{w} = 1$).
**(a)** Make a histogram of the values $f(x)$ across the ensemble, and overlay the Gaussian density with the variance $\bar K^{(2)}(x,x)$ predicted by [(3.32)](#erf-kernel); watch the visibly non-Gaussian $N=3$ histogram converge.
**(b)** Make a scatter plot of the pairs $\big(f(x), f(x')\big)$ and compare the sample $2\times2$ covariance with the predicted kernel matrix.
**(c)** For each initialization, record the stochastic kernel $K^{(2)}$ of [(3.26)](#k2-rec), and verify that its variance across initializations decays like $1/N$ (that's the mechanism by which the kernel becomes deterministic).
### 3.5 NNGP via functional integrals
One can also analyze large width at initialization via functional integrals. This just reorganizes the above steps. We'll illustrate how it works for a network with one hidden layer and zero biases, $N_{2}=1,\,N_{1}=N,\,N_{0}=d$,
$$
f(x) = \frac{1}{\sqrt{N}}\sum_{i=1}^{N} c_{i}\,\phi\Big(\frac{1}{\sqrt{d}}\,a_{i}\cdot x\Big)\,,\qquad c_{i}\in \R,\;a_{i},x\in \R^{d}\,.
$$
The initial distribution of weights, all of variance $\sigma_{w}$, is
$$
P(c,a) = \frac{1}{Z}\,e^{\textstyle -\frac{1}{2\sigma_{w}}|c|^2-\frac{1}{2\sigma_{w}}|a|^2}\,.
$$
To push this forward to a distribution on the first-layer preactivations $z_{i}=\frac{1}{\sqrt{d}}a_{i}\cdot x$, we integrate against a delta function:
Here the $z_{i}(x)$ are *arbitrary* functions of $x$, and are only tied to the network values $\frac{1}{\sqrt{d}}a_{i}\cdot x$ via the delta functions. Eqn. [(3.35)](#pcz) represents a functional probability density on the functions $z_{i}$. To evaluate it, introduce a second set of "Lagrange multiplier" functions $\tilde z_{i}$, and write
$$
\delta\Big(z_{i}(x)-\tfrac{1}{\sqrt{d}}a_{i}\cdot x\Big) = \int D\tilde z_{i}\, e^{\textstyle i \int dx\, \tilde z_i(x) (z_i(x)-\frac{1}{\sqrt{d}}a_i\cdot x)}\,,
$$
and then do the exact Gaussian integral over $a$'s, and then over $\tilde z$'s:
$$
\begin{aligned}P(c,z)&= \frac{1}{Z}\int da\,D\tilde z\, \exp\bigg[\textstyle -\frac{1}{2\sigma_{w}}|c|^{2}-\sum_{i} \frac{1}{2\sigma_{w}}|a_{i}|^{2} \notag \\&\hspace{1.5in}\textstyle + i \sum_{i} \int dx\, \tilde z_{i}(x)(z_{i}(x)-\frac{1}{\sqrt{d}}a_{i}\cdot x)\bigg] \\&= \frac{1}{Z}\int da\,D\tilde z\, \exp\bigg[\textstyle -\frac{1}{2\sigma_{w}}|c|^{2}+\sum_{i}\Big[- \frac{1}{2\sigma_{w}}\big(a_{i}+\frac{i\sigma_{w}}{\sqrt{d}}\int dx\,\tilde z_{i}(x) x\big)^{2} \notag \\&\hspace{.4in}\textstyle-\frac{\sigma_{w}}{2d}\int dx\,dx'\, \tilde z_{i}(x) (x\cdot x') \tilde z_{i}(x') + i \int dx\, \tilde z_{i}(x)z_{i}(x)\Big]\bigg] \\&= \frac{1}{Z'}\int D\tilde z\,\exp\bigg[\textstyle -\frac{1}{2\sigma_{w}}|c|^{2}-\frac{\sigma_{w}}{2d}\sum_{i}\int dx\,dx'\, \tilde z_{i}(x) (x\cdot x') \tilde z_{i}(x') \notag \\&\hspace{1.5in}\textstyle + i \sum_{i}\int dx\, \tilde z_{i}(x)z_{i}(x)\bigg] \\&= \frac{1}{Z''}\exp\bigg[\textstyle -\frac{1}{2\sigma_{w}}|c|^{2} - \frac{1}{2\sigma_{w}}\sum_{i} \int dx\,dx'\, z_{i}(x) \big(\frac{x\cdot x'}{d}\big)^{-1}z_{i}(x')\bigg]\,.\end{aligned}
$$
We clearly see the first-layer GP kernel $K^{(1)}(x,x') = \frac{\sigma_{w}}{d}x\cdot x'$ appearing.
To push-forward to the network outputs $f = \frac{1}{\sqrt{N}}c\cdot \phi(z)$, we just do this again:
$$
\begin{aligned}P(f)&= \int dc\,Dz\, P(c,z)\, \delta\Big(f(x)-\tfrac{1}{\sqrt{N}}c\cdot \phi(z(x))\Big) \\&= \frac{1}{Z''}\int dc\,Dz\,D\tilde f\, \exp\bigg[\textstyle -\frac{1}{2\sigma_{w}}|c|^{2} \notag \\&\hspace{1in}\textstyle - \frac{1}{2\sigma_{w}}\sum_{i} \int dx\,dx'\, z_{i}(x) \big(\frac{x\cdot x'}{d}\big)^{-1}z_{i}(x') \notag \\&\hspace{1in}\textstyle + i\int dx \tilde f(x)(f(x)-\frac{1}{\sqrt{N}}c\cdot \phi(z(x)))\bigg] \\&\hspace{-.2in}= \frac{1}{Z'''}\int Dz\,D\tilde f\, \exp\bigg[\textstyle - \frac{1}{2\sigma_{w}}\int dx\,dx'\, \big(\frac{x\cdot x'}{d}\big)^{-1}z(x)\cdot z(x') \notag \\&\hspace{.3in}\textstyle - \frac{\sigma_{w}}{2N}\int dx\,dx' \tilde f(x)[\phi(z(x))\cdot \phi(z(x'))]\tilde f(x') +i\int dx\, \tilde f(x) f(x)\bigg]\end{aligned}
$$
The problem is that it's now impossible to do the integral over the $z$'s exactly, because (for general activation function $\phi$) the action is no longer quadratic in them. However, as $N\to \infty$, we can treat $\exp\big[- \frac{\sigma_{w}}{2N}\int dx\,dx' \tilde f(x)[\phi(z(x))\cdot \phi(z(x'))]\tilde f(x')\big]$ as a perturbation, and approximate the *exponent* $[\phi(z(x))\cdot \phi(z(x'))]$ by its expectation value under the Gaussian measure $\exp\big[- \frac{1}{2\sigma_{w}}\int dx\,dx'\, \big(\frac{x\cdot x'}{d}\big)^{-1}z(x)\cdot z(x')\big]$ (with kernel $K^{(1)}$) on the $z$'s. To be clear, the approximation that's happening is of the form
$$
\E \,e^{\frac{1}{N} X}= e^{\frac{1}{N} \E X}+ O(\text{Var}(X)/N^{2})
$$
Up to $O(1/N)$ corrections, that produces the deterministic kernel of the previous section,
$$
\bar K^{(2)}(x,x') = \sigma_{w}\, \E_{z\sim \CN(0,K^{(1)})}\big[ \phi(z(x))\, \phi(z(x'))\big]\,,
$$
and we get
$$
\begin{aligned}P(f)&= \frac{1}{Z''''}\int D\tilde f\, \exp\bigg[\textstyle -\frac{1}{2}\int dx\, dx'\, \tilde f(x) \bar K^{(2)}(x,x') \tilde f(x') \notag \\&\hspace{1.5in}\textstyle + i\int dx\, \tilde f(x) f(x)\bigg] + O\Big(\frac{1}{N}\Big) \\&= \frac{1}{Z'''''}\exp\bigg[\textstyle -\frac{1}{2}\int dx\, dx'\, f(x) \bar K^{(2)}(x,x')^{-1}f(x')\bigg] +O\Big(\frac{1}{N}\Big) \,.\end{aligned}
$$
### 3.6 L/N
The inductive argument of [Section 3.4](#34-the-nngp-correspondence) does not hold at infinite $L$. There are two potential problems.
First, a factor of $\sigma_{w}$ appears with every layer. These factors will accumulate, and can lead either to vanishing kernels (which also turns out to imply vanishing gradients during training), in which case the network will collapse; or to unbounded kernels (and increasingly large gradients), in which case the network will be chaotic and again untrainable. We'll revisit this in [Section 3.7](#37-edge-of-chaos). With careful tuning, the issue can be avoided. For example, in a deep linear network (DLN) with $\phi(z)=z$, one can just take $\sigma_{w}=1$, known as LeCun initialization ([LeCun et al. 1998](#bib-lecun1998efficient)). In a ReLU network, the correct choice instead is $\sigma_{w}=2$, known as He initialization ([He et al. 2015](#bib-he2015delving)).
Even with $\sigma_{w}$ tuned, according to the activation function, so that on average the GP kernels remain order-one from layer to layer, non-Gaussianities can still accumulate. The upshot is that it is not width alone, say of typical size $N$, but rather the ratio $L/N$ of depth-to-width that controls the effective behavior of a network — including its departure from pure Gaussianity at large $N$.
This behavior is analyzed in Chapters 3-4 of ([Roberts et al. 2022](#bib-roberts2022principles)) using functional-integral methods. In particular, Chapter 3 considers the simple example of DLN's, and produces exact recursion relations for higher cumulants of preactivations in each layer. It is shown that in a network of constant hidden width $N$, the leading non-Gaussianity is
$$
\begin{aligned}\kappa_{4}^{(L+1)}&= \Big[ \Big(1+\frac{2}{N}\Big)^{L}-1\Big](\kappa_{2}^{(L+1)})^{2} \\&= \Big[(e^{2L/N}-1)-\frac{1}{N}\frac{2L}{N}e^{2L/N}+O(N^{-2})\Big](\kappa_{2}^{(L+1)})^{2}\,,\end{aligned}
$$
where the expansion on the RHS is at finite $L/N$ and $L,N\to\infty$. For small $L/N$, the non-Gaussianity remains small but finite, of size $L/N$. For large $L/N$, $\kappa_{4}^{(L+1)}$ blows up, as do all the other cumulants. The network outputs will fluctuate chaotically from one instantiation to another, and the network again becomes untrainable.
Here's a more direct analysis, in a special case, following ([Hanin 2026](#bib-hanin2026mdl), Lecture 2). Consider a DLN of length $L$ with input and hidden layers all of width $N$:
$$
\begin{aligned}z^{(\ell)}&= \frac{1}{\sqrt{N}}\,w^{(\ell)}z^{(\ell-1)}\,,\qquad z^{(0)}=x\,, \\ f(x)&= z^{(L+1)}=N^{-\frac{L+1}{2}}\,w^{(L+1)}\cdots w^{(1)}x\,,\end{aligned}
$$
with $w_{ij}^{(\ell)}\sim \CN(0,1)$, i.e. $\sigma_{w}=1$ (and input dimension $d=N$). Conditional on $z^{(\ell-1)}$, the vector $z^{(\ell)}$ is Gaussian $z^{(\ell)}\sim \CN\big(0,\frac{1}{N}I_{N} |z^{(\ell-1)}|^{2}\big)$, so its norm satisfies the recursion
$$
|z^{(\ell)}|^{2} \sim |z^{(\ell-1)}|^{2} \frac{\chi^{2}_{N}}{N}\,.
$$
(Recall that $\chi^{2}_{N}$ is the distribution of the sum of squares of $N$ i.i.d. standard Gaussian variables.) Let's normalize the input so that $|x|=1$. Then the recursion implies that
$$
|f(x)|^{2} \overset{d}{=}\prod_{k=1}^{L+1}\frac{\chi^{2}_{N}}{N}= \exp\sum_{k=1}^{L+1}\big(\log \chi^{2}_{N}-\log N\big)\,,
$$
with each factor of $\chi^{2}_{N}$ independent. At large $N$, we have $\log \chi_{N}^{2}-\log N = \CN(-\frac{1}{N},\frac{2}{N})+O(\frac{1}{N^{2}})$, so at large $N,L$ with $L/N$ fixed
$$
|f(x)|^{2} \overset{d}= e^{\textstyle \CN(-\frac{L}{N},2\frac{L}{N})}\,,
$$
up to higher $O(\frac{1}{N^{2}}$ corrections. Thus, as $L/N$ grows, the expectation of $\E|f(x)|\to 0$ but the variance and higher cumulants blow up.
The lesson once again is
- $N\to \infty$, $L/N=0$: Gaussian process, no correlations across layers (and, as we'll see in [Section 4](#4-the-ntk-and-large-width-during-training), no "feature learning").
- $N,L\to \infty$, $L/N$ small: controlled correlations, trainable network (and, as we'll see, feature learning during training). In open-weight LLM's today, $L/N\approx 0.005$–$0.01$; see [Exercise 3.3](#ex-ln-llms).
- $N,L\to\infty$, $L/N$ large: chaotic and untrainable.
**Exercise 3.3 (★) (Depth-to-width ratios in the wild).** Look up the depth $L$ (number of transformer blocks) and the width $N$ (model/embedding dimension) for a few families of open-weight language models — *e.g.* GPT-2, the Llama family, the Qwen family — and compute $L/N$ for each. How much does the ratio vary across families, across sizes within a family, and across the years? Is everyone training in the "controlled correlations" window?
### 3.7 Edge of chaos
We briefly remark that the analysis of how to tune initialization variances $\sigma_{w},\sigma_{b}$ (etc.) leading to stable network behavior over a large number of layers (the first point in [Section 3.6](#36-ln)) is sometimes called finding the "edge of chaos."
The term dates back to dynamical systems theory from the 80's and 90's, referring to the boundary between ordered and chaotic phases. It was adapted to neural networks around 2016 ([Poole et al. 2016](#bib-poole2016exponential); [Schoenholz et al. 2017](#bib-schoenholz2017deep)). The idea was to treat the progression through layers of a network as a dynamical system, and to ask: given two nearby inputs $x,x'$, for what range of hyperparameters ($\sigma_{w},\sigma_{b}$, etc.) and to what depth do they remain nearby, without either a) collapsing; or (b) blowing up chaotically (with no correlation). The analysis roughly amounts to making sure the kernel recursion of [(3.26)](#k2-rec), [(3.30)](#kbar-rec) leads to kernels that neither collapse nor blow up. Doing it carefully, theoretically, can produce elementary trainable networks (with no layer norms!) over 10,000 layers deep ([Xiao et al. 2018](#bib-xiao2018dynamical)).
Order/chaos phase diagram for a tanh network, with shading indicating the critical depth to which a network is trainable.
For MLP's with ReLu activations, the edge of chaos is just the one point $\sigma_{w}=2,\,\sigma_{b}=0$. Moving slightly away from this point, one can estimate the "critical depth" $L_{c}$ at which the network either collapses or blows up – see ([Schoenholz et al. 2017](#bib-schoenholz2017deep)). For MLP's with other activation functions, the edge of chaos is more interesting. For $\phi=\tanh$, the phase diagram is shown in [Figure 5](#fig-eoc).
The quantity `$\chi_{1}$' that controls the critical depth ($L_{c}= -1/\log \chi_{1}$) is the derivative of the layer-to-layer kernel recursion at its fixed point, which measures how fast nearby inputs separate. Importantly, the same quantity $\chi_{1}$ also controls the size of gradients as they back-prop through a network *during training*. There's a direct link between behavior at initialization and behavior during training. We come to training next.
## 4. The NTK and large width during training
In this section, we consider training networks that are randomly initialized as in [Section 3.2](#32-normalizing-the-network), and their behavior at large width. The lesson will be that as width$\to\infty$ these networks effectively behave like linear models: their parameters $\theta$ never stray far from their values $\theta_{0}$ at initialization, hence
Under suitable conditions, the network itself also converges exponentially fast to a "memorizing network" that fits the training data perfectly. Heuristically, this happens because small changes in many parameters, as in [(4.1)](#f-lin), can add up to an order-one change in the network values themselves. Nevertheless, the network never experiences phase transitions or feature learning — we'll make this more precise below. The large-width regime (at finite depth and training data) is thus known as the *lazy* regime, or the NTK regime. (For a self-contained review of models that are linear in their parameters, and of why training them amounts to kernel ridge regression, see [Appendix B](#b-from-linear-regression-to-kernel-regression).)
Although this regime does not describe the feature learning that makes practical networks interesting, it is important to understand for at least two reasons: it is an honest, analytically solvable model of training; and it is a precise lesson in what must be *avoided* if one wants genuine feature learning — a lesson we act on in [Section 5](#5-mean-field-scaling).
We also emphasize the connection between [Section 3](#3-large-width-at-initialization-the-gaussian-process-limit) and the current analysis. Both involve the *same* large-width limit. At initialization, it leads to a Gaussian process. During training, we'll see that it leads to a linear model with lazy learning, converging once again to a GP at late training times.
### 4.1 Gradient flow and the NTK
Throughout this section, we approximate SGD with continuous-time *gradient flow* (*cf.* [Section 2.1](#21-the-bayesian-model-of-learning)). This involves two approximations. First we update on the entire training dataset $\{(x^{\mu},y^{\mu})\}_{\mu=1}^{n}$ at each step, rather than on batches, so that SGD becomes GD:
$$
\theta_{k+1}- \theta_{k} = - \eta_{GD}\, \nabla_{\theta} \CL_{\theta}(\{x^{\mu},y^{\mu}\}) \big|_{\theta=\theta_k}
$$
where
is the full empirical loss.
Second, we scale down the learning rate $\eta_{GD}=\lambda^{-1}\eta$ and scale up the number of training steps, setting $k = \lfloor \lambda t]$ where $t\geq 0$ is a continuous "time" variable. As $\lambda\to\infty$, the gradient descent converges to a gradient flow
$$
\dot\theta := \frac{d}{dt}\theta = - \eta\, \nabla_{\theta} \CL_{\theta}\,.
$$
All theorems in this section can be extended to SGD, adding additional error terms to bound the difference between SGD and GF. (See for example ([Spiliopoulos et al. 2025](#bib-spiliopoulos2025foundations), Chapter 19).)
#### 4.1.1 Geometry of gradient flows
Suppose that $X$ is a smooth Riemannian manifold, with nondegenerate metric $g = g_{ab}dx^{a} dx^{b}$, inverse metric $g^{-1}= g^{ab}\frac{\partial}{\partial x^{a}}\frac{\partial}{\partial x^{b}}$, and a function $h:X\to \R$. Then one can define gradient flow with respect to $(h,g)$ as
Note that this equation does not make sense (it's not coordinate-free) without an inverse-metric (a bilinear form on each cotangent space $T^{*}X$). The inverse-metric is used to relate the differential $dh = (\partial_{b} h)_{b=1}^{\text{dim} X}$ on the RHS to a *vector field* $g^{-1}dh$, which is the right sort of geometric object to induce a flow on $X$.
In neural networks, $X=\R^{P}$ is the space of parameters $\theta$ and $h = \CL$ is the loss function. One always assumes a constant Euclidean metric in parameter space $g = \sum_{i=1}^{P} (d\theta^{i})^{2} = \delta_{ij}d\theta^{i} d\theta^{j}$, in order to define SGD/GD/GF/etc. This is the source of some of the simplest inductive biases of trained models, and also why one expects "features" to have some linear structure. Eqn. [(4.5)](#gf) with the metric looks like
$$
\dot\theta^{i} = -\eta\, \delta^{ij}\, \partial_{j} \CL_{\theta}\,.
$$
Now return to the general setting of Riemannian manifolds. Suppose we have a smooth map $f : X \to Y$ between two manifolds, and moreover that $h:X\to \R$ is the pull-back of a function $h^{Y}:Y\to \R$, meaning $h(x) = h^{Y}(f(x))$. Let $\{f^{\alpha}\}_{\alpha=1}^{\text{dim}\, Y}$ be the components of the map $f$, in local coordinates $y^{\alpha}$ on $Y$. Then the inverse-metric on $X$ induces an inverse-metric on $Y$ (called its push-forward):
$$
f_{*}(g^{-1}) = f_{*}\bigg[g^{ab}\frac{\partial}{\partial x^{a}}\frac{\partial}{\partial x^{b}}\bigg] = \underbrace{g^{ab} \frac{\partial f^{\alpha}}{\partial x^{a}} \frac{\partial f^{\beta}}{\partial x^{b}}}_{g^{\alpha \beta}}\frac{\partial}{\partial y^{\alpha}}\frac{\partial}{\partial y^{\beta}}\,.
$$
Schematically, $g^{\alpha\beta}= (J g^{-1}J^{T})$, where $J^{\alpha}{}_{a} = \frac{\partial f^{\alpha}}{\partial x^{a}}$ is the Jacobian matrix of $f$. Then the gradient flow on $X$ induces a gradient flow on $Y$, given by
$$
\dot y^{\alpha} = -\eta\, g^{\alpha\beta}\frac{\partial h^{Y}}{\partial y^{\beta}}\,.
$$
Let's apply this to neural networks, where $f:\R^{P}\to \R^{n}$ is the function from parameters to network outputs on the training data $(x^{\mu},y^{\mu})$, *i.e.* $f^{\mu}(\theta) = f_{\theta}(x^{\mu})$ in our standard notation. The loss function is indeed a pullback of a function on $\R^{n}$, since it depends on parameters $\theta$ only through the network values on training data,
$$
\CL_{\theta} = \CL^{(f)}(f(\theta))\,,\qquad \CL^{(f)}(f) := \sum_{\mu=1}^{n} \ell( f^{\mu},y^{\mu})\,.
$$
The push-forward of the gradient flow in parameter space becomes
The inverse-metric $\Theta^{\mu\nu}$ on the space of function values is called the "neural tangent kernel" or NTK. The name was coined in ([Jacot et al. 2018](#bib-jacot2018neural)); it also appeared in ([Du et al. 2019](#bib-du2019gradient)) under a different name. Here we are emphasizing that it's a completely natural geometric object, for any neural network (in fact, for any ML model).
We mention one final bit of geometry. Let $\delta_{\mu\nu}$ be a standard Euclidean metric in the space of function values, with inverse $\delta^{\mu\nu}$. The gradient of the loss $\delta^{\mu\nu}\pd \CL^{(f)}/\pd f^{\nu}$ that appears in [(4.10)](#ntk) (contracted with an inverse metric) is sometimes called the *residual*, and denoted $\Delta^{\mu}$ or $r^{\mu}$. For MSE loss, we've got $\CL^{(f)}= \frac{1}{2}\sum_{\mu=1}^{n} (f^{\mu}-y^{\mu})^{2} = \frac{1}{2}\delta_{\mu\nu}(f^{\mu}-y^{\mu})(f^{\nu}-y^{\nu})$. Then the residual is just the difference between function values and truth for the training data:
$$
\Delta^{\mu} := \delta^{\mu\nu}\pd \CL^{(f)}/\pd f^{\nu} = f^{\mu}-y^{\mu}\,.
$$
Noting that $y^{\mu}$ are constants, the gradient-flow equation can be recast as
Those who have studied gauge theory may recognize this as an equation for *parallel transport*, with respect to an $n\times n$ matrix-valued connection $\eta\,\Theta^{\mu}{}_{\nu}$. If one knows the dependence of $\Theta$ on time (which in general is far from trivial — see the exercise below) the solution is a time-ordered exponential
$$
\Delta^{\mu}(t) = \bigg(P\exp\bigg[-\eta\int \Theta(t)\,dt\bigg]\bigg)^{\mu}_{\;\nu}\Delta^{\nu}(0)\,.
$$
**Exercise 4.1 (The NTK from the chain rule).** Derive the flow equation [(4.10)](#ntk) for the function values $f^{\mu}$ directly, using only the chain rule — no geometry required. (The geometric derivation explains *why* an object like $\Theta^{\mu\nu}$ had to appear; the chain rule shows it appears in three lines.)
**Exercise 4.2 (How fast does the kernel move?).** Show that the time derivative of the NTK satisfies
$$
\frac{d}{dt}\Theta^{\mu\nu}= -\eta\,\big[ \nabla_{\theta} f^{\mu} \cdot \nabla_{\theta}^{2} f^{\nu} \cdot \nabla_{\theta} f^{\rho}\big] \frac{\partial \CL^{(f)}}{\partial f^{\rho}}+ (\mu\leftrightarrow\nu)\,,
$$
where $\nabla_{\theta}^{2} f^{\nu} = \nabla_{\theta}^{2} f_{\theta}(x^{\nu})$ is the Hessian. Bounding the Hessian ultimately leads to the constant-NTK theorem in [Section 4.3](#43-constant-ntk-at-large-width) below.
### 4.2 Features
The Jacobian matrix $J^{\mu}{}_{i}=\pd f^{\mu}/\pd \theta^{i}$ measures how function values on the training data depend on small variations of the network weights. Assuming there are more parameters than training data, its singular value decomposition looks like
$$
J = USV^{T}\,,\qquad S = \bigg[\text{diag}(s_{1},...,s_{n})\;\; \mathbf{0}_{n\times (P-n)}\bigg]\,.
$$
The right singular vectors (columns of $V$) define natural modes in parameter space, while the left singular vectors (columns of $U$) are the natural modes in function-value space — where "natural" means with respect to how network outputs depend on weights. The $n$ left singular vectors (together with their singular values $s_{i}$) can be thought of, abstractly, as the *features* that a network has learned. A common definition of feature *learning* is that these singular vectors+values change during training.
The analysis of features can equivalently be performed via the NTK. As a matrix, we've got
$$
\Theta = JJ^{T} = US^{2}U^{T}\,,
$$
with eigenvectors = the feature vectors, and eigenvalues = $\{s_{i}^{2}\}$ . Thus feature learning happens when the NTK changes during training.
**Remark.** A closely related "dual" object that appears in other parts of learning theory is the Fisher information metric on parameter space $\R^{P}$. Its empirical version is
$$
\mathcal{F} := J^{T}J = f^{*}(\delta)\,,
$$
the pull-back (via the network function) of constant Euclidean metric on the space of function values on the training data. In coordinates, $\mathcal{F}_{ij}=\frac{\pd f^{\mu}}{\pd \theta^{i}}\delta_{\mu\nu}\frac{\pd f^{\mu}}{\pd \theta^{j}}$.
The Fisher metric has eigenvalues $(s_{1}^{2},...,s_{n}^{2},\overbrace{0,...,0}^{P-n})$, and is automatically degenerate when $P>n$ (more parameters than training data points). This is one of the starting points of singular learning theory.
### 4.3 Constant NTK at large width
Now we get to the main theoretical result and simplification of this section. Suppose we've got an MLP of fixed depth $L$, and fixed training data size $n$, and take the width of all hidden layers to infinity. We'll assume the hidden widths all scale in a similar way $N_{L},N_{L-1},...,N_{1}\approx N\to \infty$, and we initialize the network the standard way from [Section 3.2](#32-normalizing-the-network). We want to show that for almost all samples of training data,
- The NTK $\Theta^{\mu\nu}$ remains constant throughout training, up to corrections that vanish as $N\to\infty$.
- The network itself linearizes: at any training time $t$, $f_{\theta(t)}(x) =f_{\theta_0}(x) + (\theta(t)-\theta_{0})\cdot\nabla_{\theta} f_{\theta_0}(x)$
A rigorous derivation is somewhat technical, and we refer the reader to ([Du et al. 2019](#bib-du2019gradient); [Jacot et al. 2018](#bib-jacot2018neural)), ([Roberts et al. 2022](#bib-roberts2022principles), Ch. 11), ([Spiliopoulos et al. 2025](#bib-spiliopoulos2025foundations), Ch. 20) for complete arguments, explained from a variety of different perspectives. We'll just sketch the key ideas behind the result. See also ([Hanin 2026](#bib-hanin2026mdl), Lecture 3).
We'll simplify life by setting biases to zero, so the network becomes
$$
f_{\theta}(x) = z^{(L+1)}(x)\,,\quad z^{(\ell)}_{i} = \frac{1}{\sqrt{N_{\ell-1}}}\sum_{j} w_{ij}^{(\ell)}\phi(z_{j}^{(\ell-1)})\,,\quad w^{(\ell)}_{ij}\sim\CN(0,\sigma_{w})\,,
$$
in the fan-in normalization of [Section 3.2](#32-normalizing-the-network), with order-one weights.
*Size of the NTK.* The first task is to estimate how large the NTK is at large width, because that decides how the learning rate must scale. By the chain rule, the gradient with respect to the layer-$\ell$ weights is an outer product,
$$
\frac{\pd f}{\pd w^{(\ell)}_{ij}}= \frac{1}{\sqrt{N_{\ell-1}}}\;\beta^{(\ell)}_{i}\;\phi(z^{(\ell-1)}_{j})\,,\qquad \beta^{(\ell)}_{i} := \frac{\pd f}{\pd z^{(\ell)}_{i}}\,,
$$
where the backpropagated sensitivities $\beta^{(\ell)}$ obey $\beta^{(L+1)}= 1$ and
$$
\beta^{(\ell)}_{i} = \phi'(z^{(\ell)}_{i})\,\frac{1}{\sqrt{N_{\ell}}}\sum_{k} w^{(\ell+1)}_{ki}\,\beta^{(\ell+1)}_{k}\,.
$$
The NTK is a sum over layers,
The second bracket is an average of $N_{\ell-1}$ terms. At initialization and large width, the law of large numbers replaces it by $\bar K^{(\ell)}_{\mu\nu}/\sigma_{w}$, the Gaussian-process kernel of [(3.30)](#kbar-rec) evaluated on the two inputs, which is of order one. (For $\ell=1$ the bracket is exactly $x^{\mu}\cdot x^{\nu}/d$, which is $\bar K^{(1)}_{\mu\nu}/\sigma_{w}$.) The first bracket is handled by the same recursion run backwards. Define
$$
\dot{\bar K}^{(\ell)}_{\mu\nu}= \dot{\bar K}^{(\ell)}(x^{\mu},x^{\nu}):= \E\big[ \phi'(z^{(\ell-1)}_{i}(x^{\mu}))\,\phi'(z^{(\ell-1)}_{i}(x^{\nu}))\big] \qquad \text{(for any fixed $i$)}\,,
$$
which is of order one for reasonable activation functions. Using the independence at initialization of the weights $w^{(\ell+1)}$ from everything below them, one finds $B^{(\ell)}_{\mu\nu}\approx \sigma_{w}\,\dot{\bar K}^{(\ell+1)}_{\mu\nu}\,B^{(\ell+1)}_{\mu\nu}$ with $B^{(L+1)}_{\mu\nu}=1$, and hence
This is a sum of $L+1$ terms, each of order one, with no factors of $N$ anywhere: *in the fan-in normalization the NTK is of order one at large width*, and every layer contributes to it at the same order. (Equation [(4.23)](#ntk-init) is the standard recursion for the infinite-width NTK ([Jacot et al. 2018](#bib-jacot2018neural)).) By [(4.10)](#ntk), an order-one learning rate $\eta$ then moves the function values $f^{\mu}$ by an order-one amount per unit time, which is what we want, and no rescaling of $\eta$ with width is needed in this convention. For the one-hidden-layer network $f = \frac{1}{\sqrt{N}}\sum_{i} c_{i}\phi(a_{i}\cdot x)$ one can check this by hand:
an average of $N$ order-one terms. (Here and below we absorb the first-layer factor $1/\sqrt{d}$ into the inputs $x$, since $d$ is fixed.)
**How to scale the learning rate.** The rule just used is general, and we will use it again in [Section 5](#5-mean-field-scaling). What must stay of order one as $N\to\infty$ is the *motion of the network values*, $\frac{d}{dt}f^{\mu} = -\eta\,\Theta^{\mu\nu}\,\pd\CL^{(f)}/\pd f^{\nu}$, since those are what fit the data. So: work out how the NTK scales with $N$, and choose $\eta$ so that $\eta\,\Theta = O(1)$. In the fan-in normalization $\Theta = O(1)$ and $\eta = O(1)$. In the mean-field normalization of [Section 5](#5-mean-field-scaling), $\Theta = O(1/N)$ and $\eta$ must grow like $N$.
Other conventions for the parameters are changes of variables, and they can be compared through the same rule. Suppose one puts no prefactors in the network and initializes $\tilde w^{(\ell)}= w^{(\ell)}/\sqrt{N_{\ell-1}}\sim\CN(0,\sigma_{w}/N_{\ell-1})$ as in [(3.5)](#init-fanin-var). Gradient flow with respect to the Euclidean metric $\sum(d\tilde w)^{2}$ is *not* the same flow as gradient flow with respect to $\sum(dw)^{2}$. Under the change of variables, the inverse metric $\delta^{ij}$ of [Section 4.1.1](#411-geometry-of-gradient-flows) picks up a factor $N_{\ell-1}$ in every layer-$\ell$ block, so that in the original coordinates the flow reads $\dot w^{(\ell)}= -\eta\,N_{\ell-1}\,\pd\CL/\pd w^{(\ell)}$. The NTK computed in the $\tilde w$ coordinates is correspondingly larger, $\tilde\Theta = \sum_{\ell} N_{\ell-1}\,\Theta_{(\ell)}$ in terms of the layer contributions of [(4.21)](#ntk-layers), and an order-one motion of $f$ then requires $\tilde\eta\sim 1/N$. A rescaling of the parameters is thus equivalent to a layer-dependent rescaling of the learning rate, and two conventions define the same training dynamics only if their learning rates are matched layer by layer. With a *single* learning rate $\tilde\eta\sim1/N$ in the $\tilde w$ convention, for instance, the hidden layers move as in the fan-in normalization, but the first layer, whose fan-in $d$ does not grow with $N$, is effectively frozen. When comparing learning rates across papers, first translate the parametrizations into one another; for a systematic treatment, see ([Yang & Hu 2021](#bib-yang2021tensor)).
**Caution**: the estimate [(4.23)](#ntk-init) is at fixed depth $L$ and fixed data size $n$. The NTK grows like $L$ with depth, and the operator norm of the $n\times n$ matrix $\Theta$ grows with $n$, so for deep networks or large datasets the learning rate must shrink accordingly, and the constant-NTK statements below need $n$ and $L$ to stay small compared with $N$.
*The Hessian is small.* To see that the NTK is approximately constant during training, we consider its time derivative, which by [Exercise 4.2](#ex-ntk-dot) is bounded by the operator norm of the Hessian of the network function,
The Hessian is a very large but very sparse matrix, and in the fan-in normalization its operator norm is small. We'll sketch how this works for one hidden layer ($L=1$) with one-dimensional output, following ([Hanin 2026](#bib-hanin2026mdl), Lecture 3). Writing
The operator norm of a block-diagonal matrix is the largest operator norm of its blocks. Each block has fixed size $(d+1)\times(d+1)$, contains only order-one quantities (at least at initialization), and carries an explicit factor $1/\sqrt{N}$. Thus
It's reasonable to expect that this bound continues to hold during training, and with a bit more work it can be proven (though we won't do it here). It also holds for networks of arbitrary fixed depth. A more careful analysis, keeping track of cancellations among the fluctuating terms, improves the resulting bound on the kernel drift to $O(1/N)$; the crude $O(1/\sqrt{N})$ is enough for the argument.
*Putting the pieces together.* First, the parameters do not move far in total. From $\dot\theta = -\eta\sum_{\mu}\Delta^{\mu}\nabla_{\theta} f^{\mu}$,
$$
|\dot\theta|^{2} = \eta^{2}\Big|\sum_{\mu=1}^{n}\Delta^{\mu}\nabla_{\theta} f^{\mu}\Big|^{2} \leq \eta^{2}\,|\Delta|^{2}\,|\nabla_{\theta} f|^{2} \leq \eta^{2}\,n\,|\Delta|^{2}\,\|\Theta\| = O(1)\,,
$$
where $|\nabla_{\theta} f|^{2} = \sum_{\mu}\Theta^{\mu\mu}= \mathrm{tr}\,\Theta$. So by the mean value theorem, $\theta(t)-\theta_{0} = t\,\dot\theta(\bar t)$ for some $\bar t\in[0,t]$, and
$$
|\theta(t)-\theta_{0}| = O(1)\qquad\text{for any fixed }t\,.
$$
An order-one total displacement, shared among $P\approx Nd$ parameters, means that each individual weight moves only by $O(1/\sqrt{N})$: relative to its order-one initial value, every weight barely changes. Second, the kernel freezes. By [(4.25)](#ntkdot-bound) and [(4.28)](#hessian-bound),
so the order-one kernel $\Theta$ becomes constant as $N\to\infty$. Third, the network linearizes. By Taylor's theorem with remainder,
$$
f_{\theta}(x)-f_{\theta_0}(x) = \nabla_{\theta} f(x)\big|_{\theta_0}\cdot (\theta-\theta_{0}) +\frac{1}{2}(\theta-\theta_{0})\cdot \nabla_{\theta}^{2} f(x)\big|_{\bar\theta}\cdot (\theta-\theta_{0})
$$
for some $\bar\theta$ on the segment between $\theta_{0}$ and $\theta$. Setting $\theta=\theta(t)$ and using [(4.30)](#theta-bound) and [(4.28)](#hessian-bound), the quadratic term is bounded by
$$
\Big|\frac{1}{2}(\theta-\theta_{0})\cdot \nabla_{\theta}^{2} f(x)\big|_{\bar\theta}\cdot (\theta-\theta_{0})\Big| \leq \frac{1}{2}\|\nabla_{\theta}^{2} f\|\; |\theta-\theta_{0}|^{2} = O\Big( \frac{1}{\sqrt{N}}\Big)\cdot O(1)\,,
$$
so the linear approximation becomes exact as $N\to\infty$, even though the linear term, and with it the network function, moves by an order-one amount. In words: each of the $N$ neurons has $O(1/\sqrt{N})$ leverage on the output, so an order-one change in $f$ can be assembled from $O(1/\sqrt{N})$ changes in every weight, and changes that small leave the features $\nabla_{\theta} f$ untouched.
### 4.4 Lazy learning
When the NTK is constant during training — confusingly often called the *NTK regime* — a network undergoes what is known as "lazy learning." It is not learning features, since the eigenvalues $s_{\mu}^{2}$ and eigenvectors $u_{\mu}$ of the NTK are (by construction) constant. Nevertheless, the network function converges to its training data exponentially fast: from [(4.12)](#gf-delta) with constant $\Theta$, we have
where we recall that $\Delta^{\mu}(t) = f^{\mu}-y^{\mu} = f_{\theta(t)}(x^{\mu})-y^{\mu}$. Writing the residuals in an eigenbasis for the NTK, $\Delta=\hat \Delta^{\mu}u_{\mu}$ we see that
Each feature is learned at a characteristic timescale $t_{\mu} \approx 1/(\eta\, s_{\mu}^{2})$ .
**Exercise 4.3 (Lazy predictions on test points).** Derive that on an arbitrary test point $x$, not necessarily part of the training data,
where $\Theta(x,x^{\mu}) = \nabla_{\theta} f(x)\cdot \nabla_{\theta} f(x^{\mu})$ is the NTK for a test-train pair (which remains constant during training, just like the train-train NTK $\Theta^{\mu\nu}$).
From [Exercise 4.3](#ex-ntk-test), we see that at the end of training $t\to\infty$ the network obeys
$$
f_{\theta_\infty}(x) = f_{\theta_0}(x)- \Theta(x,x^{\mu}) \Theta^{-1}_{\mu\nu}(f_{\theta_0}(x^{\nu})-y^{\nu})\,.
$$
This is an affine-linear transformation of the Gaussian process $f_{\theta_0}(x)$ at initialization, so it is again a Gaussian process. It's straightforward to check that the new mean and covariance are
where $K^{xx'}=K(x,x')$, $\Theta^{x\mu}=\Theta(x,x^{\mu})$, etc.
### 4.5 Summary
So far we've seen that in the fan-in normalization ([Section 3.2](#32-normalizing-the-network)), with order-one weights and an order-one learning rate $\eta$, at large width $N$ we've got
$$
\begin{aligned}\|\Theta\|\sim |\nabla_{\theta} f|^{2} \sim 1\,,\qquad \|\nabla_{\theta}^{2} f\|&\sim \frac{1}{\sqrt{N}}\,,\qquad |f_{\theta}-f_{\theta_0}|\sim 1\,, \\ |\theta-\theta_{0}|&\sim 1 \;\;\text{(i.e. }N^{-1/2}\text{ per weight)}\,,\end{aligned}
$$
which forces $\Theta$ to be constant and the network to linearize as a function of its parameters.
All this breaks down if
- Depth $L$ is large as well, with $L/N$ nonzero. (See ([Roberts et al. 2022](#bib-roberts2022principles), Ch. 11).)
- Data samples $n$ become large, with $n/N$ nonzero.
- Learning rate becomes large.
Both $L$ and $n$ large are interesting cases for modern neural networks!
**Exercise 4.4 (★) (Numerics: the lazy network (lab)).** This lab, and its continuation in [Exercise 5.3](#ex-mf-lab), trains one small network in both the fan-in and the mean-field normalization and tests the claims of [Sections 4](#4-the-ntk-and-large-width-during-training) and [5](#5-mean-field-scaling) numerically. The setup runs in seconds, in about sixty lines of NumPy; reproduce it, then modify it.
*The starter model.* Inputs are one-dimensional: $n=30$ equispaced training points $x^{\mu}\in[-2,2]$, with targets from a two-neuron teacher, $y(x) = 2\tanh(3x) - \frac{3}{2}\tanh(x)$. The student is a one-hidden-layer $\tanh$ network
trained by full-batch gradient descent on $\CL = \frac{1}{2}\sum_{\mu}(f(x^{\mu})-y^{\mu})^{2}$, with residuals $\Delta^{\mu} = f(x^{\mu})-y^{\mu}$. Here take the fan-in exponent $\gamma = \frac{1}{2}$, learning rate $\eta = 0.02$, and $4\times10^{4}$ steps. The empirical NTK to monitor is the $n\times n$ matrix
**(a)** *It fits.* Train and check that the loss becomes small. (The network is overparametrized; fitting is not the interesting part.)
**(b)** *Initialization.* Over about $200$ seeds, histogram $f_{0}(x)$ at a fixed input, say $x=1$. You should see an order-one Gaussian, the Gaussian process of [Section 3.4](#34-the-nngp-correspondence), and the sample covariance across two inputs should match the $\tanh$ GP kernel.
**(c)** *The kernel freezes.* For $N\in\{30,100,300,1000\}$ (a few seeds each), record $\|\Theta_{\rm end}-\Theta_{0}\|/\|\Theta_{0}\|$ and plot it against $N$ on log–log axes. Expected: it decays, with slope between $-\frac{1}{2}$ (the guaranteed bound [(4.31)](#ntk-freeze)) and $-1$ (the typical measurement). Also record the relative parameter displacement $|\theta_{\rm end}-\theta_{0}|/|\theta_{0}|$ and compare with [(4.30)](#theta-bound).
**(d)** *Mode by mode.* Diagonalize $\Theta_{0}$ and plot $|(U^{T}\Delta)_{k}(t)|$ on a semilog axis for a few $k$: straight lines with slopes $-\eta s_{k}^{2}$, the features being learned at $t_{k} = 1/(\eta s_{k}^{2})$, largest eigenvalue first, as in [(4.35)](#d-converge). Compare with the Bayesian version, [Exercise 4.5(b)](#ex-kernel-ridge-modes).
**(e)** *Prediction.* Compare the trained function on fresh test points with the frozen-kernel prediction [(4.36)](#ntk-learn).
### 4.6 Bayesian learning at large width: kernel learning
So far, we have considered MLP's at large width, trained via gradient flow (as an approximation to SGD/GD). We found that their NTK remains constant and that their late-time distribution is a Gaussian process given by [(4.38)](#ntk-gp).
Suppose that we consider the same large-width network, but "train" Bayesian inference instead, on the empirical MSE loss $\CL(\theta) = n L_{n}(\theta) = \frac{1}{2}\sum_{\mu=1}^{n} |f(x^{\mu})-y^{\mu}|^{2}$. (Recall from [Section 2.1](#21-the-bayesian-model-of-learning) that, practically speaking, this means training via Langevin dynamics.) Then we expect a Bayesian posterior
$$
P(\theta) = \frac{1}{Z}e^{-\beta \CL(\theta)}\, \pi(\theta)\,,\qquad \pi(\theta) \propto e^{\textstyle-\frac{1}{2\sigma_{w}}|\theta|^2}\,,
$$
where the prior $\pi(\theta)$ is the Gaussian of order-one variance from [(3.1)](#init-o1), and the network carries the fan-in normalization of [Section 3.2](#32-normalizing-the-network).
We'd like to push the distribution forward to the network function $f_{\theta}(x)$. At large-width, the push-forward of the prior is our familiar Gaussian process from [Section 3.4](#34-the-nngp-correspondence); expressed via a functional integral as in [Section 3.5](#35-nngp-via-functional-integrals) it becomes
$$
\pi[f] = \frac{1}{Z'}e^{\textstyle -\frac{1}{2}\int dx\,dx'\, f(x) K^{-1}(x,x')f(x')}
$$
with $K(x,x')=\bar K^{(L+1)}(x,x')$ the final-layer GP kernel. Since the loss is explicitly a function of network values, we immediately get the late-time trained distribution
$$
P[f] = \frac{1}{Z''}\,e^{\textstyle -\frac{1}{2}\Big[\beta \sum_{\mu=1}^n(f(x^\mu)-y^\mu)^2 + \int dx\,dx'\, f(x) K^{-1}(x,x')f(x') \Big] }
$$
Note that the action here is quadratic in function values, so this is still a Gaussian process!
**Exercise 4.5 (★) (Bayesian lazy learning: kernel ridge regression from the Gaussian posterior).** Lazy gradient flow at large width solved itself in [Section 4.3](#43-constant-ntk-at-large-width): the NTK froze, and the $k$-th kernel eigenmode of the residual decayed as $e^{-\eta s_k^2t}$, so that it is learned at the characteristic time $t_{k} = 1/(\eta s_{k}^{2})$. Here you derive the Bayesian-learning version of the same statement, with $\beta$ in the role of training time. Restricted to the training values $\mathbf{f} = (f(x^{1}),\ldots,f(x^{n}))$, the prior is Gaussian with covariance the kernel matrix $\hat K^{\mu\nu}:= K(x^{\mu},x^{\nu})$, so the posterior on training values is
$$
p(\mathbf{f})\;\propto\;\exp\Big[-\tfrac{1}{2}\,\mathbf{f}^{T}\hat K^{-1}\mathbf{f} - \tfrac{\beta}{2}\,|\mathbf{f}-\mathbf{y}|^{2}\Big]\,.
$$
**(a)** *Posterior mean.* The exponent is quadratic, so the posterior is Gaussian and everything follows from completing the square. Show that the posterior mean of the training values is $\boldsymbol\mu = \hat K(\hat K+\beta^{-1}I_{n})^{-1}\mathbf{y}$.
**(b)** *Eigenmodes.* Diagonalize the kernel, $\hat Ku_{k} = \lambda_{k}u_{k}$, and show that mode by mode
Conclude that the eigenmodes of the GP kernel at initialization are learned one by one as $\beta$ grows, mode $k$ switching on at the characteristic scale $\beta_{k} = 1/\lambda_{k}$, large-eigenvalue modes first (*spectral bias*).
**(c)** *Compare with gradient flow.* Place this alongside [(4.35)](#d-converge), $\hat\Delta^{k}(t) = e^{-\eta s_k^2t}\hat\Delta^{k}(0)$, that is $t_{k} = 1/(\eta s_{k}^{2})$. What is the same, and what differs? Two things to notice: the correspondence $\beta\leftrightarrow\eta t$ between the amount of data (or inverse noise) and training time; and the different shapes of the two filters, the power law $1/(1+\beta\lambda_{k})$ in $\beta$ versus the exponential $e^{-\eta s_k^2t}$ in $t$ (compare [Exercise B.5](#ex-app-filters)). Also: the Bayesian calculation uses the GP kernel $K$, the dynamical one the NTK $\Theta$, two different frozen kernels for the same network.
**(d)** *Prediction.* Show that the full posterior over functions is again a Gaussian process,
$$
P[f] \propto e^{\textstyle -\frac{1}{2} \int dx\,dx'\, [f(x)-\mu(x)] \Sigma(x,x')^{-1} [f(x')-\mu(x')] }\,,
$$
with mean and covariance
$$
\begin{aligned}\mu(x)&= \sum_{\mu,\nu}K(x,x^{\mu})\big(\hat K+\beta^{-1}I_{n}\big)^{-1}_{\mu\nu}y^{\nu}\,, \\ \Sigma(x,x')&= K(x,x')-\sum_{\mu,\nu}K(x,x^{\mu})\big(\hat K+\beta^{-1}I_{n}\big)^{-1}_{\mu\nu}K(x^{\nu},x')\,.\end{aligned}
$$
The posterior mean is *kernel ridge regression* with ridge $\beta^{-1}$ ([Appendix B.4](#b4-kernel-ridge-regression)), and $\beta\to\infty$ interpolates the data exactly. [*Hint:* the exponent is quadratic in $f$; complete the square, and use the Woodbury identity to relate the two natural forms of the answer.]
These equations embody what's sometimes known as *kernel regression*, with *ridge parameter* $\beta^{-1}=T$. ([Appendix B.4](#b4-kernel-ridge-regression) derives the same formula for a finite linear model, with no Bayesian input.) If the training samples are sufficiently dense, the mean can be an approximated by an integral in terms of the untrained kernel
$$
\mu(x) \approx \int dx'\,dx''\, K(x,x')\big(K+(n\beta)^{-1}\text{id}\big)^{-1}(x',x'')\,y(x'')\,.
$$
Writing the data "truths" $y(x) = \sum y_{k} \varphi_{k}(x)$ in an eigenbasis for $K$, with ordered eigenvalues $\lambda_{1}\geq\lambda_{2}\geq...$, we find
$$
\mu(x) = \sum_{k=1}^{\infty} \frac{\lambda_{k}}{\lambda_{k}+(n\beta)^{-1}}y_{k}\varphi_{k}(x) \approx \sum_{k\;\text{s.t.}\, \lambda_k\gtrsim(n\beta)^{-1}}y_{k}\varphi_{k}(x)\,.
$$
This is the Bayesian version of lazy learning. For given $n,\beta$, the "modes" of the data distribution with $\lambda_{k} \gtrsim (n\beta)^{-1}$ are learned by the network, and the rest are projected out.
Recall that in the Bayesian context, the rescaled inverse-temperature $n\beta$ should play a role somewhat analogous to *training time* in SGD/GD. We can see this explicitly here. For example, in SGD/GD lazy learning, the characteristic time to learn an *NTK mode* $s_{\mu}^{2}$ is $t\approx 1/(\eta\,s_{\mu}^{2})$ (see [(4.35)](#d-converge)). Correspondingly, in Bayesian lazy learning, the minimum inverse-temperature required to learn a *GP mode* is $\beta n\approx 1/\lambda_{k}$.
## 5. Mean-field scaling
[Sections 3](#3-large-width-at-initialization-the-gaussian-process-limit) and [4](#4-the-ntk-and-large-width-during-training) analyzed one and the same large-width limit: the fan-in normalization ([Section 3.2](#32-normalizing-the-network)) with order-one weights and an order-one learning rate, and width $N\to\infty$ at fixed depth and data. At initialization the network converges to a Gaussian process; during training the NTK freezes, the network linearizes, and (by the definition of [Section 4.2](#42-features)) no features are learned. It is tempting to conclude that infinite width *implies* laziness. The lesson of this section is that it does not, and that laziness was a consequence of a choice of *parametrization*, rather than of width. Changing the normalization of the output layer by a single extra factor of $\sqrt{N}$ produces a different, equally well-defined infinite-width limit in which the network starts at zero, individual weights move an order-one distance, the NTK evolves in time, and features are genuinely learned.
This alternative limit is called the *mean-field limit*. The name comes by analogy to a classic setup in statistical physics: an interacting particle system with weak all-to-all couplings, whose $N\to\infty$ behavior is governed by a self-consistent "mean field." We'll first explain what physicists mean by a mean-field theory ([Section 5.1](#51-mean-field-theories-in-physics)), then run the analogy in detail for one-hidden-layer networks ([Section 5.2](#52-mean-field-scaling-for-one-hidden-layer-networks)), see how the NTK story of [Section 4](#4-the-ntk-and-large-width-during-training) changes ([Section 5.3](#53-the-ntk-is-no-longer-constant-feature-learning)). Finally we'll discuss some practical applications, including interpolating scalings, $\mu$P, and what this framework is actually used for ([Sections 5.4–5.5](#54-the-scaling-dial-and-maximal-update-parameterization)).
The mean-field limit of one-hidden-layer networks was constructed independently, and essentially simultaneously, by four groups ([Mei et al. 2018](#bib-mei2018mean); [Chizat & Bach 2018](#bib-chizat2018global); [Rotskoff & Vanden-Eijnden 2022](#bib-rotskoff2022trainability); [Sirignano & Spiliopoulos 2020](#bib-sirignano2020mean)); our treatment draws on all of them. Pedagogical accounts may be found in ([Spiliopoulos et al. 2025](#bib-spiliopoulos2025foundations), Ch. 20) and ([Hanin 2026](#bib-hanin2026mdl), Lecture 4).
### 5.1 Mean-field theories in physics
"Mean-field theory" is one of the oldest and simplest approximation schemes in statistical mechanics: *replace the fluctuating environment felt by each degree of freedom with its average (a macroscopic order parameter), then demand self-consistency upon recalculating this average*. For most systems this is an uncontrolled approximation, though often a remarkably successful one *cf.* ([Tong 2012](#bib-tongsft), Chapter 1). There is, however, a family of models — those with weak, all-to-all couplings — for which it becomes exact as $N\to\infty$; and networks in mean-field scaling will turn out to be precisely such a model.
#### 5.1.1 Warm-up: the Curie–Weiss magnet
The Curie–Weiss magnet is the Ising model of [Section 2.4](#24-a-physics-example-the-ising-model) on the *fully connected* graph: every spin couples to every other, with the coupling scaled down to $J/N$ so that the total force on any one spin, and the total energy, stay of order $N^{0}$ and $N^{1}$ respectively. Setting $J=1$, the energy [(2.19)](#ising-e) becomes
$$
E_{N}(s) = -\frac{1}{2N}\sum_{i,j=1}^{N} s_{i} s_{j} - B \sum_{i=1}^{N} s_{i} = -N\Big(\frac{m^{2}}{2}+Bm\Big)\,,\qquad m := \frac{1}{N}\sum_{i} s_{i}\,.
$$
The second equality is exact, not an approximation: on the complete graph the energy depends on the spins only through the magnetization. Compare the mean-field energy [(2.22)](#ising-e-mf) of the lattice model, with $dJ\to 1$. Everything that was approximate in [Section 2.4](#24-a-physics-example-the-ising-model) is therefore exact here. The partition function reduces to the one-dimensional integral [(2.25)](#ising-action) with the same entropy $h(m)$, the saddle-point equation is
$$
m = \tanh\big(\beta(m+B)\big)\,,
$$
and the model has the second-order transition of [Figure 4](#fig-ising-secondorder) at $\beta_{c} = 1$ and the first-order transition of [Figure 2](#fig-ising-wells) in $B$, exactly as analyzed in [Section 2.5](#25-phases-and-phase-transitions).
The reason mean-field theory is exact here is worth stating in a form that carries over to networks. The field felt by a single spin is $-\pd E_{N}/\pd s_{i} = m + B$: spin $i$ does not care about any individual other spin, only about the empirical average $m$, the "mean field." The classic argument then has three features:
1. *Each degree of freedom couples to the others only through an average.* (True by construction here.)
2. *The average barely fluctuates.* If the spins are roughly independent, the central limit theorem makes the fluctuations of $m$ of size $O(1/\sqrt{N})$; as $N\to\infty$ each spin sees a *deterministic* effective field $m+B$.
3. *The average is fixed self-consistently.* A single spin in the field $m+B$ has $\E(s_{i}) = \tanh(\beta(m+B))$; since $m$ was *defined* as the average spin, consistency requires [(5.2)](#cw-sc).
For a finite-dimensional lattice, where each spin has a few strongly coupled neighbors, steps 1–2 are approximations. For the $1/N$ all-to-all model they become exact as $N\to\infty$: the mean field is an average of $N$ weakly correlated terms, and the law of large numbers kills its fluctuations. Networks in mean-field scaling will turn out to be models of exactly this kind.
**Exercise 5.1 (Curie–Weiss by fixed-point iteration).** Solve [(5.2)](#cw-sc) numerically by the fixed-point iteration $m^{(k+1)}= \tanh\big(\beta(m^{(k)}+B)\big)$, for $\beta$ below and above $1$ and for a few values of $B$, and watch the iteration converge. Which fixed points are stable? Keep this iteration in mind: we will meet its infinite-dimensional cousin below.
#### 5.1.2 The general shape of a mean-field model
Nothing in the story above needed binary spins. More abstractly, one can take $N$ "particles" $\theta_{i}\in\R^{p}$ with an extensive energy of the form
with a single-particle potential $V$ and a symmetric pair interaction $U$, again with the characteristic $1/N$ coupling. The force felt by particle $\theta_{i}$ is $-\nabla_{\theta_i}E_{N}$. More so, for particles dominated by drag (highly viscous motion), the dynamics at zero temperature are governed by familiar gradient-descent equations
where $\nabla_{1}$ differentiates the first slot of $U$. Once more each particle feels only an empirical average of forces.
It is worth packaging this observation in a way that will carry over verbatim to networks. Define the *empirical measure*
a probability measure on $\R^{p}$ recording where the collection of $N$ particles sit. Then [(5.4)](#mf-flow) says that every particle descends one and the same effective single-particle potential, shaped by where everyone currently is:
The particles are coupled *only through* $\rho_{N}$ — the exact analogue of "every spin feels only $m$." The same three-step argument as before then suggests that as $N\to\infty$ the empirical measure becomes deterministic and self-consistent. We will carry this out carefully, for networks, in the next subsection.
### 5.2 Mean-field scaling for one-hidden-layer networks
Return to the one-hidden-layer network [(4.26)](#l1-net) of [Sections 3.5](#35-nngp-via-functional-integrals) and [4.3](#43-constant-ntk-at-large-width), but now normalized with $1/N$ out front instead of $1/\sqrt{N}$, with the same *order-one* weights:
with $c_{i}\in\R$, $a_{i},x\in\R^{d}$ (as in [(4.26)](#l1-net), the first-layer factor $1/\sqrt{d}$ has been absorbed into $x$). Compared with the fan-in normalization of [Sections 3–4](#3-large-width-at-initialization-the-gaussian-process-limit), mean-field scaling suppresses the output by *one extra factor of $1/\sqrt{N}$*. The function class is identical — only the coordinates on it have changed. However, gradient flow is *not* invariant under reparametrization: as emphasized in [Section 4.1.1](#411-geometry-of-gradient-flows) and in the callout of [Section 4.3](#43-constant-ntk-at-large-width), gradient descent presupposes the Euclidean metric on parameter space, and rescaling coordinates rescales that metric. Two parametrizations of the same function class can (and here will) flow along completely different trajectories.
For a mean-field analysis, our $N$ particles are going to be the $N$ hidden neurons, with parameters
$$
\theta_{i} = (c_{i},a_{i}) \in \R^{d+1}\,,\qquad i=1,...,N\,.
$$
Their values are encoded in the empirical measure [(5.5)](#mf-rho), and the network [(5.7)](#mf-f) depends on its parameters only through $\rho_{N}$,
We think of the network value $f_{\rho}$ as a linear functional of $\rho$. Neurons are viewed as exchangeable particles, and the measure is their state.
#### 5.2.1 Initialization
At initialization, the summands $c_{i}\phi(a_{i}\cdot x)$ in [(5.7)](#mf-f) are i.i.d. with mean zero (since $c_{i}$ is independent of $a_{i}$ and $\E[c_{i}]=0$). With the $1/N$ prefactor, the law of large numbers gives, for each fixed $x$,
$$
f_{\theta_0}(x) = \underbrace{\E\big[c\,\phi(a\cdot x)\big]}_{=\,0}+ \,O(1/\sqrt{N}) \;\xrightarrow{\;N\to\infty\;}\; 0\,.
$$
Contrast this with [Section 3.4](#34-the-nngp-correspondence): with the $1/\sqrt{N}$ normalization, the same sum sits exactly at the scale of its fluctuations, and the CLT produces an order-one Gaussian process. Mean-field scaling instead sits at the scale of the *mean*.
The network begins training as (essentially) the zero function; there is no random-function prior to speak of, and whatever the trained network becomes is built by training rather than inherited from initialization.
#### 5.2.2 The loss is a mean-field energy
Let's train on MSE loss over the data $(x^{\mu},y^{\mu})_{\mu=1}^{n}$, as in [Section 4](#4-the-ntk-and-large-width-during-training):
$$
\CL(\theta) = \frac{1}{2}\sum_{\mu=1}^{n} \big(f_{\theta}(x^{\mu})-y^{\mu}\big)^{2}\,.
$$
Inserting [(5.7)](#mf-f) and expanding the square, we get
Comparing [(5.12)](#mf-vu) with [(5.3)](#mf-h), we see that up to an overall factor of $N$, the loss is precisely a mean-field energy. In particular, $E_{N}:=N\CL$ has the canonical normalization from [(5.3)](#mf-h) and is extensive, while the loss itself can be interpreted as the *intensive* energy per neuron.
We can also identify physical roles for the two potentials! $V$ rewards neurons for aligning their output function $c\,\phi(a\cdot{}\cdot)$ with the data, playing the role of the external field $B$; while $U$ is a data-weighted interaction between *pairs* of neurons that penalizes redundant features, playing the role of the spin-spin coupling.
#### 5.2.3 Equations of motion and the mean field
Because $\pd f/\pd \theta_{i} = O(1/N)$ in this parametrization, raw gradients of $\CL$ are $O(1/N)$, and gradient flow at fixed learning rate would freeze as $N\to\infty$. We must therefore scale the learning rate *up*,
$$
\eta = N\,\tilde\eta\,,
$$
with $\tilde\eta$ fixed. (Equivalently, we could run vanilla gradient flow on the extensive energy $E_{N}=N\CL$.) This is the rule of [Section 4.3](#43-constant-ntk-at-large-width) at work: the NTK, computed below in [(5.24)](#mf-ntk), is $O(1/N)$, so keeping $\frac{d}{dt}f^{\mu} = -\eta\,\Theta^{\mu\nu}\Delta^{\nu}$ of order one requires $\eta\propto N$. Note the contrast with [Section 4.3](#43-constant-ntk-at-large-width): there the NTK was of order one and the learning rate needed no rescaling; here the NTK *shrinks* like $1/N$ and we scale up. Recall also the caveat from that section about changes of variables. The mean-field network could equally be written without a prefactor, $f = \sum_{i}\tilde c_{i}\phi(a_{i}\cdot x)$ with $\tilde c_{i} = c_{i}/N\sim\CN(0,1/N^{2})$; but then the *same* dynamics [(5.15)](#mf-eom) requires a learning rate $\tilde\eta/N$ for the readout weights $\tilde c_{i}$ and $N\tilde\eta$ for the hidden weights $a_{i}$ — one learning rate per layer, differing by a factor $N^{2}$. Keeping all weights of order one is what lets a single $\eta$ do the job, and it is the reason we adopted that convention in [Section 3.2](#32-normalizing-the-network).
With [(5.14)](#mf-lr), the flow $\dot\theta = -\eta\nabla_{\theta}\CL$ becomes, per neuron,
where $\Delta^{\mu} = f_{\rho_N}(x^{\mu})-y^{\mu}$ is the residual of [Section 4.1](#41-gradient-flow-and-the-ntk). Equivalently, in the potential form [(5.6)](#mf-psi),
with $\Delta_{\rho}^{\mu} = f_{\rho}(x^{\mu})-y^{\mu}$.
As anticipated, the velocity of neuron $i$ depends on its own state $(c_{i},a_{i})$, and on all the other neurons only through the average $f_{\rho_N}= \frac{1}{N}\sum_{j} c_{j}\phi(a_{j}\cdot{}\cdot)$, via the residual $\Delta^{\mu}$. The network function $f$ is the mean field here! We arrive at a dictionary:
| Curie–Weiss | one-hidden-layer network, mean-field scaling |
| --- | --- |
| spin $s_{i}\in\{\pm1\}$ | neuron $\theta_{i}=(c_{i},a_{i})\in\R^{d+1}$ |
| magnetization $m=\frac{1}{N}\sum_{j}s_{j}$ | network function $f_{\rho}=\frac{1}{N}\sum_{j}c_{j}\phi(a_{j}\cdot{}\cdot)$ |
| external field $B$ | data term $V$ (the targets $y^{\mu}$) |
| pair coupling $-\tfrac{1}{N}s_{i}s_{j}$ | neuron-neuron interaction $\tfrac{1}{N}U(\theta_{i},\theta_{j})$ |
| effective field $m+B$ | residual $\Delta^{\mu}_{\rho}= f_{\rho}(x^{\mu})-y^{\mu}$ |
| extensive energy $E_{N}$ | $N\CL$ (the loss is intensive) |
| zero-temperature dynamics | gradient flow |
| Langevin at temperature $T$ | noisy training, [Section 2.1](#21-the-bayesian-model-of-learning) |
| $m = \tanh(\beta(m+B))$ | self-consistent Gibbs measure, [Equation 5.22](#mf-gibbs) below |
**Exercise 5.2 (Structure of the mean-field equations).** **(a)** Derive [(5.12)](#mf-vu)–[(5.13)](#mf-vu-def) and [(5.15)](#mf-eom), and check that [(5.16)](#mf-eom-psi) reproduces them, using the symmetry of $U$.
**(b)** (*Balancedness.*) For a positively homogeneous activation such as ReLU, $\phi'(z)z = \phi(z)$. Show that in this case [(5.15)](#mf-eom) conserves $c_{i}^{2} - |a_{i}|^{2}$ for every neuron, exactly, along the flow. The two layers of a given neuron therefore remain coupled for all time. (On the other hand, distinct neurons will decouple from each other, as we'll see momentarily.)
**(c)** Observe that for MSE loss $\dot c_{i}$ in [(5.15)](#mf-eom) does not depend on $c_{i}$ (only on $a_{i}$ and the mean field). Where does the linearity of $f_{\rho}$ in $c$ enter?
#### 5.2.4 From particles to measures: the continuity equation
We now pass to the $N\to\infty$ limit, in distinct steps.
*Step 1: rewrite the $N$ coupled ODEs as one equation for the empirical measure.* (This remains exact at finite $N$.) Suppose particles move deterministically with some velocity field, $\dot\theta_{t} = v_{t}(\theta_{t})$, and let $\rho_{t}$ denote the law of one such particle.
The expectation of a smooth test function $g$ obeys
$$
\begin{aligned}\frac{d}{dt}\int g\,d\rho_{t}&= \frac{d}{dt}\,\E\big[g(\theta_{t})\big] = \E\big[\nabla g(\theta_{t})\cdot v_{t}(\theta_{t})\big] \\&= \int \nabla g\cdot v_{t}\, d\rho_{t} = -\int g\;\nabla\!\cdot\!(\rho_{t} v_{t})\,,\end{aligned}
$$
integrating by parts in the last step. Since $g$ was arbitrary,
This is known as the *continuity equation*, and it simply expresses conservation of probability. Applied to our system, with $v_{t} = -\tilde\eta\,\nabla\Psi(\cdot\,;\rho_{N,t})$ from [(5.16)](#mf-eom-psi), the empirical measure satisfies (weakly, since it is a sum of delta functions)
An important consequence of the mean-field ($1/N$, all-to-all) coupling here is that [(5.19)](#mf-cont-n) is a closed equation in $\rho_{N}$.
*Step 2: the limit.* At finite $N$, $\rho_{N,t}$ is a random, atomic measure. (Atomic means it's a sum of delta functions; it's random due to random initialization.) As $N\to\infty$ with i.i.d. initialization $\theta_{i}\sim\rho_{0}$, the law of large numbers suggests — and the theorems of ([Mei et al. 2018](#bib-mei2018mean); [Chizat & Bach 2018](#bib-chizat2018global); [Rotskoff & Vanden-Eijnden 2022](#bib-rotskoff2022trainability); [Sirignano & Spiliopoulos 2020](#bib-sirignano2020mean)) confirm — that $\rho_{N,t}$ converges (weakly) to $\rho_{t}$, a deterministic curve of smooth probability measures solving the limiting PDE
There are quantitative $O(1/\sqrt{N})$ error bounds at any fixed time horizon.
#### 5.2.5 The McKean–Vlasov equation and self-consistency
In this context, the analogue of the Curie–Weiss self-consistency [(5.2)](#cw-sc), namely $m=\tanh(\beta(m+B))$, is called the McKean–Vlasov equation. To get it, focus on a single representative neuron $\theta_{i}$. As $N\to \infty$, we can replace the force on it by the ensemble average:
The constraint $\rho_{t} = \mathrm{Law}(\bar\theta_{t})$ is the entire content of the equation: the particle is driven by a field built from a distribution, and that distribution is *its own law*.
(Probabilists call [(5.21)](#mf-mv) a *nonlinear* Markov process, nonlinear not in $\theta$ but in the law.) The McKean–Vlasov equation is equivalent to the continuity equation; the particle trajectories of [(5.21)](#mf-mv) are the characteristic curves of the PDE [(5.20)](#mf-pde).
We remark on a few structural features:
- *Propagation of chaos.* Since neurons interact only through a quantity ($\rho_{N}$, equivalently $f_{\rho_N}$) that becomes deterministic as $N\to\infty$, they become asymptotically *independent*: any fixed collection of $k$ neurons converges to $k$ i.i.d. copies of [(5.21)](#mf-mv), with cross-neuron correlations of size $O(1/N)$. Thus, independence at initialization *propagates* instead of degrading ([Sznitman 1991](#bib-sznitman1991topics)).
- *Gradient-flow structure.* Define the $N\to\infty$ loss as a functional of the measure, $\CL[\rho] = \frac{1}{2}\sum_{\mu} (f_{\rho}(x^{\mu})-y^{\mu})^{2}$. Since $f_{\rho}$ is linear in $\rho$, this functional is convex in $\rho$.
This ultimately leads to global-convergence theorems. For example, it was shown in ([Chizat & Bach 2018](#bib-chizat2018global)) that for homogeneous activations (e.g. ReLU), if training converges at all, it converges to a *global* minimizer of the convex functional $\CL[\rho]$.
- *Temperature.* If one trains instead with the Langevin dynamics of [Section 2.1](#21-the-bayesian-model-of-learning), Eqn. [(2.6)](#training-options), at temperature $T=\beta^{-1}$, one finds convergence to a stationary distribution that satisfies ([Mei et al. 2018](#bib-mei2018mean))
This is another version of the "$m=\tanh(\beta(m+B))$" consistency condition, now for noisy training.
Mean-field training for a $d=1$ example: $f_{\rho}(x)=\frac{1}{N}\sum_{i}c_{i}\tanh(a_{i}x)$, trained by the rescaled gradient flow [(5.15)](#mf-eom) on $n=30$ samples of the two-neuron teacher $y(x)=2\tanh(3x)-\tfrac{3}{2}\tanh(x)$; $\tau=\tilde\eta\, t$ is rescaled training time. *Left:* the empirical measure $\rho_{N,t}$ in the $(a,c)$ plane for $N=300$ (dots, colored by time; thin gray lines trace $40$ of the particle trajectories), drawn over the density of an $N=30{,}000$ run at $\tau=20$ (shaded contours), standing in for the limiting measure $\rho_{t}$. Individual neurons move order-one distances, and the measure condenses onto one-dimensional filaments — the configurations realizing the teacher. (The point symmetry $(a,c)\mapsto(-a,-c)$ reflects the oddness of $\tanh$.) *Right:* the mean field itself, $f_{\rho_t}(x)$, converging to the target (dashed; circles mark the training data).
How would one actually *solve* the self-consistent system [(5.21)](#mf-mv)? The same way one solves [(5.2)](#cw-sc) ([Exercise 5.1](#ex-cw)): freeze a guessed environment $(\rho_{t})$, solve the now-ordinary ODE, compute the law of the solution, feed it back, iterate. On short time horizons this map is a contraction, which is how existence, uniqueness, and propagation of chaos are all proven in one stroke ([Sznitman 1991](#bib-sznitman1991topics)). Alternatively, solve the PDE [(5.20)](#mf-pde) on a grid, which is feasible (only) for small $d$.
In [Figure 6](#fig-mf-transport) we instead solve the system by simply training a network with $d=1$. In the left panel, each dot is a neuron. Training is a flow of the cloud of neurons, and the cloud condenses onto the low-dimensional structures in weight space that realize the target — learning its features. [Exercise 5.3](#ex-mf-lab) reproduces this figure, and sets it beside the frozen cloud of the lazy network.
**Exercise 5.3 (★) (Numerics: the mean-field network (lab)).** Continue the lab of [Exercise 4.4](#ex-ntk-lab): same data, same teacher $y(x) = 2\tanh(3x)-\frac{3}{2}\tanh(x)$, same student [(4.40)](#lab-student), but now with the mean-field exponent $\gamma = 1$ and the rescaled updates [(5.15)](#mf-eom),
$$
\begin{aligned}c_{i}&\mathrel{-}= \Delta\tau\sum_{\mu}\Delta^{\mu}\tanh(a_{i}x^{\mu})\,, \\ a_{i}&\mathrel{-}= \Delta\tau\,c_{i}\sum_{\mu}\Delta^{\mu}\big(1-\tanh^{2}(a_{i}x^{\mu})\big)x^{\mu}\,,\end{aligned}
$$
with $\Delta\tau = 5\times10^{-4}$ for $4\times10^{4}$ steps, so that the rescaled time $\tau = \tilde\eta\, t$ runs to $20$. Note that there is no $1/N$ in these updates: that is the required $O(N)$ learning rate [(5.14)](#mf-lr), already built in ($\eta = N\,\Delta\tau$). Monitor the NTK [(4.41)](#lab-ntk) with $\gamma=1$.
**(a)** *It fits.* Train and check that the loss becomes small, if more slowly in rescaled time than the lazy run.
**(b)** *Initialization.* Over about $200$ seeds, histogram $f_{0}(x)$ at $x=1$ for several $N$: the distribution collapses to $0$ with width $\sim1/\sqrt{N}$. Check the scaling.
**(c)** *The kernel moves.* For $N\in\{30,100,300,1000\}$, record $\|\Theta_{\rm end}-\Theta_{0}\|/\|\Theta_{0}\|$ and plot it together with the lazy curve of [Exercise 4.4(c)](#ex-ntk-lab-freeze): one plot, two curves, the whole dichotomy.
**(d)** *Watch the particles.* Scatter-plot the neuron cloud $(a_{i},c_{i})$ at $\tau\approx0,\,0.4,\,2,\,20$ for both runs. Lazy: the cloud is static (each weight moves $O(1/\sqrt{N})$). Mean-field: the cloud transports and condenses onto thin filaments, the configurations realizing the teacher, as in [Figure 6](#fig-mf-transport). Repeat at $N=3000$: the pattern is the same, because you are looking at a deterministic limiting measure.
**(e)** *Now play.* The starter teacher is itself a two-neuron $\tanh$ network, so the student's features can align *exactly* with the teacher's. Is that feature learning, or feature copying? Try teachers outside the model class and watch what the cloud does instead: $y(x) = \sin(2x)$; $y(x) = |x|-1$; a ReLU student on the $\tanh$ teacher (and vice versa); or move to $d=2$ with $y(x) = \tanh(x_{1}+x_{2})$ and plot the cloud of directions $a_{i}/|a_{i}|$. Which aspects of (c)–(d) are universal, and which were special to the aligned case?
### 5.3 The NTK is no longer constant: feature learning
[Section 4.3](#43-constant-ntk-at-large-width) showed that in the fan-in normalization the NTK freezes as $N\to\infty$, and [Section 4.2](#42-features) defined feature learning as change of the NTK during training. Let's redo that computation in mean-field scaling. From [(5.7)](#mf-f), $\pd f/\pd c_{i} = \frac{1}{N}\phi(a_{i}\cdot x)$ and $\pd f/\pd a_{i} = \frac{1}{N}c_{i}\phi'(a_{i}\cdot x)\,x$, so the NTK of [Section 4.1](#41-gradient-flow-and-the-ntk) is
with the per-neuron kernel
$$
k_{\theta}(x,x') := \phi(a\cdot x)\phi(a\cdot x') + c^{2}\,\phi'(a\cdot x)\phi'(a\cdot x')\; x\cdot x'\,.
$$
Thus $\Theta = O(1/N)$, recovering the previous observation that gradients are small, which was the reason for the rescaling [(5.14)](#mf-lr). The estimates of [Section 4.3](#43-constant-ntk-at-large-width) also explain why the kernel is now free to move. The Hessian blocks [(4.27)](#l1-hessian) carry $1/N$ instead of $1/\sqrt{N}$, so $\|\nabla_{\theta}^{2} f\| = O(1/N)$; but the learning rate is $\eta = N\tilde\eta$, so the bound [(4.25)](#ntkdot-bound) gives for the relative drift $\|\dot\Theta\|/\|\Theta\| \lesssim \eta\,\|\nabla^{2}_{\theta} f\|\,|\Delta| = O(1)$. The bound no longer forces the kernel to freeze, and the equations of motion [(5.15)](#mf-eom) will show that it does move. As $N\to \infty$, the meaningful object is $N\Theta = \int k_{\theta}\, d\rho_{N}$, and the function-space dynamics [(4.10)](#ntk) now becomes a kernel gradient flow *with a time-dependent kernel*:
*At initialization* the LLN makes $N\Theta \to \E_{\rho_0}[k_{\theta}]$, a deterministic kernel — the *same* limiting kernel [(4.24)](#ntk-l1) as in the fan-in normalization. The two limits start from the same place.
*During training*, however, they part ways. In [Section 4.3](#43-constant-ntk-at-large-width), each parameter moved only $O(1/\sqrt{N})$, so $\rho_{N}$ stayed within $O(1/\sqrt{N})$ of $\rho_{0}$ and the kernel was pinned, [(4.31)](#ntk-freeze). Here, the equations of motion [(5.15)](#mf-eom) move each particle an *order-one* distance. (Each neuron has only $O(1/N)$ leverage, so fitting the data *forces* the average neuron to move a distance $O(1)$.) Thus $\rho_{t}$ transports macroscopically, and by [(5.26)](#mf-kernelflow) the kernel is dragged along at order one as well.
### 5.4 The scaling dial and Maximal Update Parameterization
It's quite instructive to look at an entire one-parameter family of scalings, which has fan-in/NTK scaling and mean-field scaling as its endpoints. The family can be encoded via a prefactor
where $c_{i}, a_{ij}\sim \CN(0,1)$. For each $\gamma$, as $N\to\infty$, the learning rate must be simultaneously be scaled as $\eta = N^{2\gamma-1}\tilde\eta$ so that the function-space dynamics stays $O(1)$ (one can derive that from $\Theta\sim N^{1-2\gamma}$). Repeating the estimates above gives (see [Exercise 5.4](#ex-gamma)):
| | $f_{\theta_0}$ at init | per-neuron motion | limit |
| --- | --- | --- | --- |
| $\gamma=\tfrac{1}{2}$ | $O(1)$ (GP, CLT) | $O(N^{-1/2})$ | lazy: [Sections 3](#3-large-width-at-initialization-the-gaussian-process-limit) and [4](#4-the-ntk-and-large-width-during-training) |
| $\tfrac{1}{2}<\gamma<1$ | $O(N^{\frac{1}{2}-\gamma})\to 0$ | $O(N^{\gamma-1})\to 0$ | lazy (frozen kernel, zero init) |
| $\gamma=1$ | $O(N^{-1/2})\to 0$ | $O(1)$ | mean-field: this section |
Notably, every intermediate $\gamma$ collapses onto the kernel limit: the initialization GP is scaled away, but the neurons still move vanishingly little, so training is again linear in an (appropriately rescaled) frozen kernel ([Chizat et al. 2019](#bib-chizat2019lazy); [Mei et al. 2018](#bib-mei2018mean)). Mean-field, $\gamma=1$, is the *unique* scaling in this family with both a well-defined limit and order-one feature movement; it's the "maximal update" that one can take without blowing up. This uniqueness can be generalized from one hidden layer to arbitrary depth (where consistency forces *layer-dependent* learning-rate and initialization scalings), and is known as the Maximal Update Parametrization, $\mu$P, in Yang and Hu's line of work on tensor programs ([Yang & Hu 2021](#bib-yang2021tensor)). In particular, mean-field scaling is $\mu$P for $L=1$.
**Exercise 5.4 (Lazy vs. rich along the dial).** For the family [(5.27)](#gamma-f):
**(a)** show $\Theta = O(N^{1-2\gamma})$ and hence that $\eta = N^{2\gamma-1}\tilde\eta$ is the right co-scaling;
**(b)** show the initialization output is $O(N^{1/2-\gamma})$ and per-neuron velocities are $O(N^{\gamma-1})$;
**(c)** conclude that the (normalized) kernel drift over any fixed training horizon is $O(N^{\gamma-1})$, vanishing for every $\gamma<1$;
**(d)** reconcile $\gamma=\frac{1}{2}$ with the bounds of [Section 4.3](#43-constant-ntk-at-large-width). Optional: relate $\gamma$ to the output-scale parameter $\alpha$ of ([Chizat et al. 2019](#bib-chizat2019lazy)), where multiplying any model by $\alpha\to\infty$ produces lazy training.
### 5.5 What is this framework good for?
Nobody solves the McKean–Vlasov PDE [(5.20)](#mf-pde) for a real network! There are practical payoffs, though, of three different kinds: *scaling prescriptions*, *mechanism identification*, and *reduced descriptions* that downstream theories build on.
Hyperparameter transfer ($\mu$P) is one key application. The mean-field analysis is, fundamentally, a theory of which joint scalings of initialization and learning rate keep the training dynamics nondegenerate and learning features as $N\to\infty$. Its deep-network generalization ([Yang & Hu 2021](#bib-yang2021tensor)) turns this into an engineering recipe, $\mu$Transfer ([Yang et al. 2022](#bib-yang2022tensor)): parametrize the network so that its infinite-width limit is the feature-learning one, and then optimal hyperparameters become approximately width-independent, so one can tune a small proxy model and *transfer* the settings to the large one. In the original demonstration this retuned GPT-3 6.7B at a few percent of pretraining cost and outperformed the published model; variants of this parametrization are now standard in large-scale training.
The dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization) also lets one interpolate between regimes on real architectures, and measure where feature learning actually pays off. Kernel-regime behavior tends to be competitive on small datasets, while on CIFAR-scale image tasks trained networks beat their own frozen NTKs by large margins ([Geiger et al. 2020](#bib-geiger2020disentangling)). Diagnostics from this circle of ideas (kernel distance/alignment during training, e.g. the silent-alignment effect ([Atanasov et al. 2022](#bib-atanasov2022silent))) are now routine tools for asking "is my network actually learning features, or just doing kernel regression?"
The mean-field limit here holds data fixed while sending $N\to\infty$, which unfortunately does not yet make it a quantitative theory of generalization at realistic data/width ratios. The dynamical mean-field theory of [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft) ([Bordelon & Pehlevan 2022](#bib-bordelon2022self)) will rectify this somewhat; it is a field-theoretic refinement of the current analysis, whose extensions do reproduce compute-optimal neural scaling laws, in particular at large dataset size ([Bordelon et al. 2024](#bib-bordelon2024dynamical)).
### 5.6 Grokking as a lazy-to-rich transition: a reading exercise
We finish this section with two guided readings that apply various concepts introduced above to one of the most interesting and fundamental phenomenons in learning: grokking ([Power et al. 2022](#bib-power2022grokking)). Grokking refers to the following behavior. One trains a small network on an algorithmic dataset — the canonical example is addition mod $P$, presented as pairs $(n,m)$ with label $(n+m)\bmod P$, with only a fraction of all $P^{2}$ facts shown. The network reaches perfect *training* accuracy quickly, while *test* accuracy sits at chance: it has simply memorized. Then, orders of magnitude later in training and long after the train loss has flattened, test accuracy jumps abruptly to near $100\%$. Such delayed, abrupt generalization, with emergent capabilities, demands a mechanism.
Theory has produced (at least) two quantitative explanations, and this subsection and the next walk through one each: the *dynamical* lens of Kumar–Bordelon–Gershman–Pehlevan ([Kumar et al. 2024](#bib-kumar2024grokking)) (this subsection) — grokking as a delayed transition from lazy to rich dynamics — and the *equilibrium* lens of Rubin–Seroussi–Ringel ([Rubin et al. 2024](#bib-rubin2024grokking)) ([Section 5.7](#57-bayesian-grokking-a-first-order-transition-alternative-reading)) — grokking as a first-order phase transition with metastable lag. The two walkthroughs are parallel in structure and use largely disjoint toolkits from these notes. The dynamical one presumes comfort with [Sections 4](#4-the-ntk-and-large-width-during-training) and [5.4](#54-the-scaling-dial-and-maximal-update-parameterization); the Bayesian one is a denser read, and combines scaling analysis with a mean-field approximation.
The work of Kumar, Bordelon, Gershman, and Pehlevan ([Kumar et al. 2024](#bib-kumar2024grokking)) invokes no statistical-mechanics ensemble at all. There is only gradient descent, and the lazy/rich dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization), *unfolding in time*. The paper claims that a network groks when 1) it first fits the training data through its cheap, fast channel (kernel regression with the frozen NTK of its initialization), and that lazy fit memorizes without generalizing because the initial kernel is misaligned with the task; then 2) on a slower timescale controlled by the network's laziness, the features finally move (the kernel rotates, as in [Section 5.3](#53-the-ntk-is-no-longer-constant-feature-learning)), a generalizing solution is found, and test loss falls long after train loss did. Two beautiful features of this work are that its clearest examples grok under vanilla gradient descent with zero weight decay, and that the parameter norm increases during grokking. Both are direct counterexamples to earlier accounts of grokking built on weight decay and shrinking norms.
The paper is written almost entirely in the language these notes have already built ([Sections 4](#4-the-ntk-and-large-width-during-training) and [5](#5-mean-field-scaling)), so the dialect box below is short. Read the introduction, the two-layer polynomial-regression model, and the experiments (modular arithmetic, MNIST, a one-layer transformer); skip the dynamical mean-field asymptotics on a first pass — that machinery belongs to [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft). Then work through [Exercise 5.5](#ex-grok-lazy).
A few translations of terminology from ([Kumar et al. 2024](#bib-kumar2024grokking)):
- *Lazy / linearized dynamics*: the regime of [Section 4.3](#43-constant-ntk-at-large-width) — the network moves only in the affine subspace [(4.1)](#f-lin) spanned by its initial gradients, so training is kernel regression with the frozen initial NTK ([Sections 4.1](#41-gradient-flow-and-the-ntk) and [4.6](#46-bayesian-learning-at-large-width-kernel-learning)).
- *Laziness parameter $\alpha$*: an overall multiplier on the network output, $f \mapsto \alpha f$ (with loss and learning rate co-scaled). This is the continuous version of the $\gamma$-dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization), in the form introduced by ([Chizat et al. 2019](#bib-chizat2019lazy)): $\alpha\to\infty$ is fully lazy, small $\alpha$ is rich.
- *Kernel–task (mis)alignment*: the language of [Section 4.2](#42-features) — how much of the target lies in the top eigenmodes of the initial NTK. Misaligned $=$ the target's weight sits in small-eigenvalue modes, which the mode-filtering formula of [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning) says are not learnable from the available data.
- *Sufficient statistics*: a handful of scalar functions of the weights whose closed dynamics determine the test loss — the dynamical analogue of the order parameters of [Section 2](#2-statistical-mechanics-of-learning).
**Exercise 5.5 (Grokking as a lazy-to-rich transition).** Read ([Kumar et al. 2024](#bib-kumar2024grokking)) as assigned above, and answer the following, citing where in the paper each answer lives. (Parts (a)–(e) have exact counterparts in [Exercise 5.6](#ex-grok), the other walkthrough; part (f) is the shared comparison question.)
**(a)** *What regime is the analysis set in, and what dial controls it?* There is no thermodynamic limit and no equilibrium ensemble here, so what plays the role that a "scaling limit" plays elsewhere in this section? How is the laziness parameter $\alpha$ implemented, and how does it relate to the $\gamma$-dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization)? [Section 5.2](#52-mean-field-scaling-for-one-hidden-layer-networks)'s scaling prescription had two halves, an output scale and a learning-rate scale. Which half is doing the work here? (The other walkthrough meets the opposite half: see [Exercise 5.6(a)](#ex-grok-limit).)
**(b)** *What plays the role of the order parameter?* What quantity distinguishes the memorized state from the generalizing one, and how would you measure it in a real training run? Compare with the Fourier overlap $\Phi(w)$ that plays this role in the other walkthrough ([Exercise 5.6(b)](#ex-grok-order)): what information do the two quantities share?
**(c)** *Necessary conditions.* The paper states three conditions for grokking to occur. Find them, and translate each into the language of these notes ([Sections 4.6](#46-bayesian-learning-at-large-width-kernel-learning), [4.3](#43-constant-ntk-at-large-width) and [5.4](#54-the-scaling-dial-and-maximal-update-parameterization)).
**(d)** *What parameter is tuned?* What makes grokking more dramatic, and what makes it vanish entirely? What does the theory predict in the two extremes $\alpha\to\infty$ and $\alpha$ small? Contrast with the control parameters of [Exercise 5.6(d)](#ex-grok-dial).
**(e)** *Where does the delay come from?* Both theories must explain the gap between the train-loss drop and the test-loss jump. Here there is no metastability and no noise — so what produces two well-separated timescales in a single deterministic flow? What does each theory *require* to produce a delay (noise? weight decay? neither?), and how do the paper's zero-weight-decay, growing-norm experiments bear on earlier explanations of grokking?
**(f)** *Synthesis.* Are the equilibrium mechanism of [Section 5.7](#57-bayesian-grokking-a-first-order-transition-alternative-reading) and the dynamical mechanism of this section incompatible, or two corners of one phenomenon? Design at least one experiment whose outcome would favor one over the other — consider, e.g., how the grokking delay should depend on training noise at fixed hyperparameters, or how the delay should fluctuate from seed to seed, under each mechanism.
### 5.7 Bayesian grokking: a first-order transition (alternative reading)
The alternative walkthrough takes the opposite starting point. Delayed, abrupt generalization should by now also trigger a reflex from [Section 2.5](#25-phases-and-phase-transitions): competing branches of a free energy, a first-order transition, and a metastable lag, exactly as dissected for the Ising perceptron in [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning). The paper ([Rubin et al. 2024](#bib-rubin2024grokking)) makes this reflex precise, and it is a fitting capstone for these notes: it assembles, in one real (if small) deep-learning phenomenon, nearly every tool we have developed: the Gibbs posterior of noisy training ([Section 2.1](#21-the-bayesian-model-of-learning)), GP kernels and kernel regression ([Sections 3.3](#33-gaussian-processes) and [4.6](#46-bayesian-learning-at-large-width-kernel-learning)), mean-field scaling and self-consistency (this section), and the anatomy of first-order transitions ([Section 2.5](#25-phases-and-phase-transitions)). It is, however, a genuinely dense paper, written in an idiosyncratic dialect of the "adaptive kernel" school ([Seroussi et al. 2023](#bib-seroussi2023separation)). This subsection is therefore structured as a reading exercise in the style of [Section 1](#1-intro-reading-learning-mechanics): the goal is to *recognize the structures*, not to verify the computations. Read the paper's introduction, setup, and modular-arithmetic results (*do not* read the appendices on a first pass), then work through [Exercise 5.6](#ex-grok).
*Warm-up.* Before any statistical mechanics, it helps to know concretely what the network is trying to become. Gromov ([Gromov 2023](#bib-gromov2023grokking)) writes down *exact* weights for a two-layer network with quadratic activation that solve modular addition: the first-layer weights are plane waves on $\mathbb{Z}_{P}$ — entries $\cos(2\pi k n/P)$, $\sin(2\pi k n/P)$ at various frequencies $k$ — and the quadratic activation's product-to-sum identities convert them into the required functions of $(n+m) \bmod P$. Skim ([Gromov 2023](#bib-gromov2023grokking)) just long enough to believe this. The lesson is that "the Fourier coefficients of the first-layer weights" becomes the obvious coordinate system for measuring how much structure a network has actually learned.
Some translations from the paper's language into that of these notes:
- *Langevin training at equilibrium*: as in [Section 2.1](#21-the-bayesian-model-of-learning) verbatim: noisy gradient descent with weight decay; at late times, the Gibbs posterior; the weight decay $\gamma$ sets the Gaussian prior and the noise level sets the temperature.
- *(Adaptive) kernel*: the GP kernels of [Section 3.3](#33-gaussian-processes), but evaluated in the *posterior* rather than the prior, self-consistently: the kernel determines the posterior over hidden weights, which determines the kernel. This is a finite-temperature sibling of the self-consistent Gibbs measure [(5.22)](#mf-gibbs).
- *Equivalent kernel (EK) approximation*: the dense-data approximation of [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning), where sums over training points are replaced by integrals over the input distribution. This is a cousin of an annealed approximation.
- *Mean-field scaling*: the readout starts with variance $\sigma_{a}^{2}/N$, and $\sigma_{a}^{2}\to0$ as $N\to\infty$; how this relates to the $1/N$ of [(5.7)](#mf-f) is [Exercise 5.6(a)](#ex-grok-limit).
- *GFL ("Gaussian feature learning") phase*: refers to a disordered phase, in which the single-neuron posterior is a Gaussian, centered at zero overlap with the teacher's structure; the kernel adapts only perturbatively.
- *Mixed phase*: a finite fraction of neurons shift to nonzero overlap (here: onto Fourier modes of the task) while the rest remain Gaussian; this *coexistence* is a signal of a first-order transition, cf. [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning).
**Exercise 5.6 (Grokking as a first-order phase transition).** Read ([Rubin et al. 2024](#bib-rubin2024grokking)) as assigned above, and answer the following, citing where in the paper each answer lives. (Parts (a)–(e) have exact counterparts in [Exercise 5.5](#ex-grok-lazy), the other walkthrough; its part (f) is the shared comparison question.)
**(a)** *What scaling limit is taken?* Look for where the output-layer weight decay is fixed, and for the phrase "mean-field scaling." How does the paper's scaling relate to the parametrization [(5.7)](#mf-f) and to the $\gamma$-dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization)? Careful — this is subtler than it looks: in an *equilibrium* (Bayesian) framework there is no learning rate, so which part of [Section 5.2](#52-mean-field-scaling-for-one-hidden-layer-networks)'s scaling prescription survives, and in what combination?
**(b)** *What is the order parameter* for the modular-addition task? What does the posterior distribution of a *single neuron's* input weights look like on each side of the transition, and which figure of [Section 5.2](#52-mean-field-scaling-for-one-hidden-layer-networks) does the ordered side resemble?
**(c)** *What are the main computational techniques?* List them, and match each to the section of these notes where its ancestor appears.
**(d)** *What parameter is tuned* to cross the phase transition in the paper's phase diagrams?
**(e)** *Where did the time go?* Grokking is a statement about training *time*, yet the paper computes a static equilibrium ensemble. How does a first-order transition in an equilibrium free energy produce *delayed, abrupt* generalization in a training run? Compare, line by line, with the metastability discussion for the Ising perceptron in [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning).
## 6. Non-equilibrium physics, MSRJD, and DMFT
So far our statistical-mechanics toolkit has been an equilibrium one: Gibbs measures, partition functions, and free energies ([Section 2](#2-statistical-mechanics-of-learning)), together with limits in which training simplifies to kernel regression ([Section 4](#4-the-ntk-and-large-width-during-training)) or to a deterministic flow of measures ([Section 5](#5-mean-field-scaling)). This section introduces the formalism that confronts the *dynamics* of learning head-on. The strategy is borrowed from non-equilibrium statistical physics: write the stochastic dynamics itself as a path integral — the Martin–Siggia–Rose–Janssen–De Dominicis (MSRJD) construction — then average over the randomness (noise, initialization, and, crucially, the training data, treated as quenched disorder exactly as in [Section 2.3](#23-partition-functions-free-energy-and-cumulants)), and finally extract closed, self-consistent equations for a small set of *two-time* order parameters: correlation functions $C(t,t')$ and response functions $R(t,t')$. The result is known as dynamical mean-field theory (DMFT), and it is the third and final incarnation of "mean field" in these notes: the statics of [Section 2](#2-statistical-mechanics-of-learning) determined scalars self-consistently, the McKean–Vlasov limit of [Section 5](#5-mean-field-scaling) determined a measure, and DMFT determines whole functions of two times.
The payoffs in this section are a complete worked application — double descent as a *dynamical phase transition* in linear regression ([Section 6.2](#62-double-descent-as-a-dynamical-phase-transition)) — and a guided tour of the modern DMFT program for deep learning ([Section 6.3](#63-dynamical-mean-field-theory-dmft)), whose scaling-law application closes the notes in [Section 7](#7-scaling-laws). Throughout, we follow the spirit of Altland and Simons ([Altland & Simons 2010](#bib-altlandcondensedmatterfield2010)), which motivates the formalism in condensed-matter language with a good balance of detail and intuition, and we borrow freely from two pedagogical tutorials written for neuroscientists, ([Chow & Buice 2012](#bib-chowpathintegralmethods2012)) and ([Stapmanns et al. 2020](#bib-stapmannsselfconsistentformulationsstochastic2020)).
### 6.1 From stochastic dynamics to path integrals: the MSRJD formalism
When we model a system with fluctuations — or deliberately add noise to a deterministic model, as in the Langevin dynamics of [Section 2.1](#21-the-bayesian-model-of-learning) — differential equations become stochastic, and the difficulty of solving them jumps. A **stochastic differential equation**, at first order for example, takes the form
$$
\dv{x}{t}= f(x) + \xi\,,
$$
where $\xi(t)$ is a (Gaussian) white noise process — additive noise, $\xi = \xi(t)$, as opposed to multiplicative noise $\xi=\xi(x,t)$ — and $x=x(t)$ is now a random trajectory. The noise is entirely defined by the first and second moments,
$$
\expval{\xi(t)}_{\xi} = 0, \hspace{20pt}\expval{\xi(t)\xi(t')}_{\xi} = D \delta(t-t')\,,
$$
where the expectations are with respect to the noise ensemble[^3]. Obtaining closed-form solutions for the time evolution of the probability distribution for $x(t)$ is usually challenging. These evolutions can be attacked by the traditional routes — the Langevin and Fokker–Planck descriptions — for which the standard references are Risken ([Risken 1996](#bib-risken1996fokker)), Gardiner ([Gardiner 2004](#bib-gardiner2004handbook)), and van Kampen ([van Kampen 2007](#bib-vankampen2007stochastic)). The objects of interest here are instead the moments of $x(t)$ and the probability density function $p(x,t).$ Methods developed in statistical physics (in particular those used in perturbative quantum field theory and classical non-equilibrium systems) using path integrals provide a way to calculate these statistical objects. The path integral, or functional integral, approach is convenient for perturbative calculations (expanding around an exactly solvable theory with a small parameter), but also permits importing other tools like renormalization group methods or other explicitly non-perturbative approaches when there is no small parameter, or the expansion breaks down (as they often do).
The conceptual recipe of the MSRJD formalism is short, and worth internalizing before the derivation:
- write down a *discrete* SDE, fixing a convention (Ito) that enforces causality;
- express the transition probability between time steps as a delta function imposing the update rule;
- assemble the product over infinitesimal steps into a path integral;
- write each delta function in Fourier representation, introducing an auxiliary *response field*;
- average over the noise, then add sources and differentiate to extract statistics.
It turns out that the historic approach of Onsager and Machlup ([Onsager & Machlup 1953](#bib-onsager-machlup-1953)) is the "Lagrangian" version of this "Hamiltonian"-like approach; the two are related by a Legendre transform, or equivalently by integrating out the auxiliary field. That auxiliary field is aptly named the **response field**: as we will see below, its correlation functions compute how the system *responds* to external perturbations, and it is the genuinely new ingredient of the non-equilibrium formalism.
We now give the derivation, quickly but completely. The standard first step is to discretize time for a continuous equation, or simply start with a discrete update rule. Under this discretization $x(t) \to x_{i}$, where $i = 1, ..., N,$ so that,
$$
x_{i} = x_{i-1}+ \Delta t \Big(f(x_{i-1}) + \xi_{i-1}\Big)\,,
$$
where $\Delta t$ is the temporal spacing. We can rearrange this so that the RHS is zero and write the transition probability between steps as delta function, and the product of these representing a trajectory solution. Defining a solution of the stochastic differential equation as $x[\xi]$, an expectation value of an observable $\mathcal{A}[x]$ is formally,
$$
\expval{\mathcal{A} [x]}_{\xi}= \int \mathcal{D}x \mathcal{A}[x] \expval{\delta(x- x[\xi])}_{\xi} = \int \mathcal{D}x \mathcal{A}[x] \abs{\frac{\delta X}{\delta x}}\delta(X)\,.
$$
We used the shorthand $X_{i} \equiv x_{i} - x_{i-1}- \Delta t \big[ f(x_{i-1}) + \xi_{i-1}\big]$ for the constraint imposed at each step. The measure here, is a typical functional path integral measure $\mathcal{D}x = \prod_{i} d x_{i}$ with $\delta(x- x[\xi])$ similarly factorized across steps. In particular we also picked a first order finite difference and the convention that the noise is evaluated at the previous time-step, which is called Ito discretization or convention. Many common processes are manifestly Ito, like SGD for example, where the noise drawn at a batch doesn't affect the weights at that iteration. There are other conventions, the most common alternative is Stratonovich (which is more like a midpoint rule), however, Ito has a unit (volume-preserving) functional determinant and other schemes require careful treatment to avoid anomalous terms from the determinant. In our case the matrix is triangular with a unit diagonal. Substituting in the definition of $X$,
$$
\expval{\mathcal{A}[x]}_{\xi} = \int \mathcal{D}x\, \mathcal{A}[x]\, \expval{\delta \big(\partial_t x - f(x) - \xi\big)}_{\xi}\, .
$$
Finally, we can use the Fourier representation of the delta, introducing an auxiliary variable $\tilde{x}$, which in the continuum is
$$
\expval{\mathcal{A}[x]}_{\xi} = \int \mathcal{D}x\, \mathcal{D}\tilde{x}\; \mathcal{A}[x]\; e^{\textstyle\, i \int dt\, \tilde{x}\, (\partial_t x - f(x) - \xi)}\,.
$$
The final step is the payoff: the noise now appears *linearly* in the exponent, so the Gaussian average over $\xi$ can be done exactly, $\big\langle e^{-i\int dt\, \tilde x\, \xi}\big\rangle_{\xi} = e^{-\frac{D}{2}\int dt\, \tilde x^2}$, leaving the **MSRJD functional integral**
Four structural comments, each of which does real work later:
- *The partition function is one.* Setting $\mathcal{A}=1$, the integral is just $\int \mathcal{D}x$ of a normalized probability: $Z_{\rm dyn}=1$, with no sources. Unlike in equilibrium, all the physics lives in correlation functions, obtained by coupling sources to both $x$ and $\tilde x$ and differentiating — compare the $Z_{\rm dyn}=1$ starting point of [Section 6.2](#62-double-descent-as-a-dynamical-phase-transition).
- *The response field earns its name.* Perturbing the drift, $f \to f + h(t)$, changes the action by $i\int \tilde x\, h$, so
Correlators of $\tilde x$ compute responses to kicks. Causality is built in: $R(t,t')=0$ for $tt'$ (with $T = D/2$ in our normalization) ([Kubo 1957](#bib-kubostatisticalmechanicaltheoryirreversible1957); [Kubo 1966](#bib-kubofluctuationdissipationtheorem1966)). A measured *violation* of the FDT is therefore the standard diagnostic that a system is genuinely out of equilibrium — the warning light used to great effect in [Section 6.2](#62-double-descent-as-a-dynamical-phase-transition).
- *Onsager–Machlup.* The action [(6.7)](#msrjd-action) is quadratic in $\tilde x$, so the response field can be integrated out exactly, giving the probability of a trajectory as $P[x] \propto e^{-\frac{1}{2D}\int dt\, (\partial_t x - f(x))^2}$: the Lagrangian form promised above, in which the most probable path extremizes a classical action.
**Exercise 6.1 (The Ornstein–Uhlenbeck process and the FDT).** **(a)** Fill in the noise average leading to [(6.7)](#msrjd-action).
**(b)** For the linear drift $f(x) = -k x$ (the Ornstein–Uhlenbeck process), compute the stationary correlation and response functions, either by Fourier transforming the quadratic action or by solving the SDE directly: $C(\tau) = \frac{D}{2k}e^{-k|\tau|}$ and $R(\tau) = \theta(\tau)\, e^{-k\tau}$.
**(c)** Verify the fluctuation–dissipation theorem $T R(\tau) = -\pd_{\tau} C(\tau)$ for $\tau>0$, with $T = D/2$.
**(d)** Explain, from the discrete-time construction, why $\expval{\tilde x(t) \tilde x(t')}= 0$ and why $R$ is causal.
### 6.2 Double descent as a dynamical phase transition
As a complete worked application of the formalism, we describe recent work of Li and Goldenfeld ([Li & Goldenfeld 2026](#bib-librokenergodicityviolation2026)), which uses dynamical mean-field theory to analyze the double-descent phenomenon and classifies it as a *dynamical phase transition*, all within the simplest possible setting: a teacher–student model of linear regression. Even this very simple model requires a long calculation, but the paper is quite readable despite the overhead. The headline: in the under-parametrized regime the optimization process is ergodic and relaxes to a finite-error solution, while the over-parametrized regime is a *non-ergodic* phase with genuinely non-equilibrium dynamics — an initial memorization stage leaves a persistent memory of the initial conditions and confines the dynamics to a subspace around the initialization. The diagnostic is the *response function* $R(t,t')$ of [Section 6.1](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism) — how a perturbation at time $t'$ affects the model's state at time $t$; the same object appears as a Green function in quantum mechanics, signal analysis, and information theory. In the under-parametrized regime the response function decays: no memory. In the over-parametrized regime it saturates at a finite value: persistent memory of initialization. Consistently, the first regime obeys the fluctuation–dissipation theorem ([Kubo 1957](#bib-kubostatisticalmechanicaltheoryirreversible1957); [Kubo 1966](#bib-kubofluctuationdissipationtheorem1966); [Goldenfeld 2018](#bib-goldenfeldlecturesphasetransitions2018a)), while in the second the FDT breaks — the standard warning light, flagged in [Section 6.1](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism), that a system has gone out of equilibrium. The model here is quite simple, and the authors propose that the frameworks here may generalize across architectures, with different universality classes based on the variations of the task, proposing random feature models with multiplicative (state-dependent) noise as an extension, with the hopes of eventually attaining an RG framework for interpreting phase transitions in neural networks.
#### 6.2.1 Model definition and formalism
**Notation.** We transcribe the paper into the conventions of these notes: $P$ trainable parameters, $n$ training samples, learning rate $\eta$, teacher weights $w^{*}$, student weights $\hat w$ (with response fields $\tilde w$), and the *overparametrization ratio*
$$
\nu := P/n
$$
as control parameter — the reciprocal of [Section 2](#2-statistical-mechanics-of-learning)'s load $\alpha = n/N$, so that the over-parametrized phase is $\nu>1$. For readers consulting the paper, the dictionary is:
| these notes | Li–Goldenfeld ([Li & Goldenfeld 2026](#bib-librokenergodicityviolation2026)) |
| --- | --- |
| $P$ (parameters) | $d$ |
| $n$ (samples) | $N$ |
| $\nu = P/n$ | $\alpha = d/N$ |
| $w^{*}$, $\hat w$, $\tilde w$ | $\beta$, $\hat\beta$, $\tilde\beta$ |
| $\eta$ (learning rate) | $\zeta$ |
The batch-inclusion probability $\gamma$, the ridge $\lambda$, and the label-noise variance $D$ are as in the paper.
Explicitly, the student and teacher share the same architecture
$$
y^{\mu}=\frac{1}{\sqrt{P}}\sum_{i=1}^{P} w^{*}_{i} x_{i}^{\mu}+\epsilon^{\mu}
$$
where $P$ is the number of model parameters, $\mu=1, \ldots, n$ indexes the dataset with inputs $x^{\mu}$, the $w^{*}_{i}$ are the teacher's weights, and $\epsilon^{\mu}$ is label noise. The thermodynamic limit takes $n \to \infty$ and $P \to \infty$ with the overparametrization ratio $\nu = P/n$ — the paper's **complexity parameter** — held *fixed*. The three main regimes are now definable as under-parametrized: $\nu<1$, near-criticality: $\nu \approx 1$, and over-parametrized: $\nu>1$. We note their optimization procedure, which is SGD with both regularization controlled by $\lambda$ (set to zero for most of the work) and batching controlled by the binary indicator variable $s_{\mu}$ which is drawn for each point in the dataset with probability $\gamma$,
From here they construct the MSRJD integral exactly as in [Section 6.1](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism); taking functional derivatives and averaging over the disorder (the dataset, the initialization, the batch selection, and the label noise) yields the correlation and response functions that are the main objects of study. We skip most of the details, but it is worth displaying the starting point — the dynamical partition function for gradient flow, equal to one just as promised in [Section 6.1](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism):
We want to highlight the general structure of this equation as a path integral over an action with sources coupled to both the "physical" parameter fields and the response fields,
$$
\begin{aligned}Z_{\operatorname{d y n}}[\mathbf{j}, \tilde{\mathbf{j}}]=\Big\langle\int \mathcal{D}\hat{\boldsymbol{w}}(t) \mathcal{D}{\tilde{\boldsymbol{w}}}(t) \exp\Big(S_{\operatorname{d y n}}&+\sum_{i=1}^{P} \int \tilde{j}_{i}(t) \hat{w}_{i}(t) d t \\&+\sum_{i=1}^{n} \int j_{i}(t) \tilde{w}_{i}(t) d t\Big)\Big\rangle,\end{aligned}
$$
which is a very common object to see in statistical field theory! In general $S_{\text{dyn}}$ will hold all of the microscopic details of the theory, and this is a unique case with a fairly tractable model. In general we are not so lucky. The expectation value denoted by angled brackets is averaging over the input dataset $\{\mathbf{x}^{\mu}\}_{\mu=1}^{n}$, the initial condition $\hat{\boldsymbol{w}}(0)$, the stochastic sampling $s_{\mu}(t)$, and the added label noise $\{\epsilon^{\mu}\}_{\mu=1}^{n}$. For comedic relief we also write their full dynamical function integrated over quenched disorder:
$$
\begin{aligned}Z_{\operatorname{d y n}}&\propto \int \mathcal{D}\boldsymbol{\mathcal{B}}(t) \int \mathcal{D}\mathcal{Q}(t,t^{\prime}) \exp \bigg(- i \tilde{\boldsymbol{w}}\partial_{t}{\hat{\boldsymbol{w}}}- \lambda P\, Q_{2}(t,t) \\&-P\, Q_{1}\cdot \hat{Q}_{1} -P\, Q_{2} \cdot \hat{Q}_{2}-P\, Q_{3}\cdot \hat{Q}_{3} -P\, Q_{4}\cdot \hat{Q}_{4} \\&-\frac{n}{2}\log\operatorname{det}{\left(\begin{array}{cc}\gamma Q_2^\top - \gamma Q_4+\delta &\gamma^2 Q_3 \\ Q_1+D& \gamma Q_2-\gamma Q_4^\top+\delta\end{array}\right)}+\tilde{\boldsymbol{j}}\hat{\boldsymbol{w}}+ \boldsymbol{j}\tilde{\boldsymbol{w}}\\&+\int\!\!\int dt\, dt^{\prime}\sum_{k}\Big[\big( \hat{w}_{k}(t) - w^{*}_{k} \big) \hat{Q}_{1}(t,t^{\prime}) \big( \hat{w}_{k}(t^{\prime}) - w^{*}_{k} \big) \\&\hspace{1.2in}+ i\hat{w}_{k}(t) \hat{Q}_{2}(t,t^{\prime}) \tilde{w}_{k}(t^{\prime}) + i\tilde{w}_{k}(t) \hat{Q}_{3}(t,t^{\prime})\, i\tilde{w}_{k}(t^{\prime})\Big] \\&+\int dt \sum_{k} i\tilde{w}_{k}(t)\hat{Q}_{4}(t)w^{*}_{k}\bigg),\end{aligned}
$$
where $\mathcal{D}\boldsymbol{\mathcal{B}}(t)$ abbreviates the measure $\mathcal{D}\hat{\boldsymbol{w}}(t)\, \mathcal{D}\tilde{\boldsymbol{w}}(t)$ over the dynamical fields, and the new measure $\mathcal{D}\mathcal{Q}(t,t^{\prime})$ is
$$
\begin{aligned}\mathcal{D}\mathcal{Q}(t,t^{\prime}) = \left(\frac{P}{2\pi}\right)^{4} \mathcal{D}Q_{1}(t, t^{\prime})\, \mathcal{D}\hat{Q}_{1}(t, t^{\prime})\,&\mathcal{D}Q_{2}(t, t^{\prime})\, \mathcal{D}\hat{Q}_{2}(t, t^{\prime})\, \\ \times\, \mathcal{D}Q_{3}(t, t^{\prime})\, \mathcal{D}\hat{Q}_{3}(t, t^{\prime})\,&\mathcal{D}Q_{4}(t)\, \mathcal{D}\hat{Q}_{4}(t)\,.\end{aligned}
$$
Each of the $\mathcal{Q}$ variables represents a particular correlation function between two variables.
#### 6.2.2 Overview of results via training and generalization error
They define the generalization error as a mean squared error (MSE)
$$
\varepsilon_{g}(t) = \langle \left(\frac{1}{\sqrt{P}}\sum_{i=1}^{P} x_{i}^{0}[\hat{w}_{i}(t)-w^{*}_{i}]-\epsilon^{0}\right)^{2} \rangle
$$
After taking the thermodynamic limit, the generalization error is expressed in terms of the equal time correlator $C(t,t)$ determined from the DMFT equations and $D$, the noise variance:
$$
\varepsilon_{g}(t) = C(t, t)+D
$$
They find two forms of this equation and propose a third, a scaling relation near the critical point. At $\lambda=0$ for $\nu<1$, the asymptotic error is:
Note the pole as criticality is approached. The quantity $\chi=\int_{0}^{\infty}R(\tau)\,d\tau$ is the **susceptibility**; it is well-defined only when the response function decays, i.e. only in equilibrium. Its appearance here signals that in this phase the response to perturbations of the initial conditions really does decay: asymptotically, the system carries no memory of its initialization.
The *non-ergodic phase* with persistent memory is captured by the asymptotic generalization error depending explicitly on the retention of the initial condition $C(0,0)$ and still showing sensitivity to $\nu \to 1$ analogous to the $\nu<1$ regime near criticality:
where fully, $C(0,0) = \langle \hat{w}^{2}(0)\rangle + r$. The term $\langle \hat{w}^{2}(0)\rangle$ is the variance of the initialization and $r$ is the mean-squared norm of the ground-truth parameters per dimension. Note that here, the response function does not decay, because we carry memory of the initial conditions. It saturates at $R(\tau \to \infty)=\frac{\nu-1}{\nu}$, which they call the *ergodicity-breaking degree*, quantifying the persistent memory of a perturbation applied to the initialization.
They also show similar trends with the training error, namely for $\nu<1$ the error is asymptotically finite, as the model is still attempting to approximate the solution, while for $\nu>1$ the training error vanishes as initialization-dependent memorization ensures the model always performs accurately on the training data. They write training error as
$$
\lim_{t \to \infty}\varepsilon_{t}(t)=(\gamma \chi+1)^{-2}\,\varepsilon_{g}(t\to\infty)\,.
$$
For $\nu<1$ it simplifies to $D(1-\nu)$, while again, for $\nu>1$, the system is no longer at equilibrium and the susceptibility $\chi=\int_{0}^{\infty}R(\tau)d\tau$ diverges, which drives the training error to zero.
To characterize the critical behavior of the model near $\nu \approx 1$, they calculate the critical exponents $(\sigma,\phi)$ and propose a scaling function. The first part is relatively straightforward, where from Eqs. [(6.18)](#alpha-1ge) and [(6.19)](#alpha-1ge), one finds $\varepsilon_{g}(t\to\infty)$ diverges from both sides $\nu\to 1^{\pm}$. This suggests the scaling $\varepsilon_{g}(t\to\infty)\sim |1-\nu|^{-\sigma}$ with an exponent $\sigma = 1$ in DMFT. We remind the reader that these are asymptotic long-time steady-state results and both training and generalization error are time-dependent. Since the saturation time is $\nu$-dependent, a general feature of phase transitions is a dynamic scaling data collapse should be present near the transition ([Goldenfeld 2018](#bib-goldenfeldlecturesphasetransitions2018a)). They make this argument and propose a scaling ansatz of form:
where they introduce a new critical exponent $\phi$ and write the universal scaling function $F(z)$, and $\gamma$ is still the probability of samples being included. They make two observations about the form of $F$ near and at the critical point. First, near the point, $\nu \to 1^{\pm}$ corresponds to $z\to 0$, they argue the singularity in front must be absorbed into $F$ and it takes the asymptotic form $F(z) \sim z^{\sigma/\phi}$. Second, at $\nu = 1$, since the generalization error grows in time instead of saturating, $\epsilon_{g} (t) \sim (\gamma t )^{\sigma / \phi}$. The exponent $\phi$ therefore characterizes this growth in training time. Notably, the form of Eq. [(6.21)](#udc) implies a characteristic time scale, which they denote $\tau(\nu) = (\gamma |1-\nu|^{\phi})^{-1}$ to describe convergence in the scaling of the error. Since the error grows, we would expect this to diverge as $\nu \to 1$, which it does, and is an iconic signature of a phase transition, called **critical slowing down**. In MC calculations for lattice systems for example, this shows up in the autocorrelation time increasing, and requiring more time to get uncorrelated samples as we approach the critical point. A great deal of energy goes into designing methods to circumvent this slowing down when attempting to measure precise values of critical exponents using MC techniques, and motivates the choice of entirely different paradigms like RG flow or conformal bootstrap where applicable ([Goldenfeld 2018](#bib-goldenfeldlecturesphasetransitions2018a)).
**Exercise 6.2 (Double descent in the lab).** Simulate the model of this subsection: teacher $w^{*} \sim \CN(0, r\, I_{P})$, inputs $x^{\mu} \sim \CN(0, I_{P})$, labels $y^{\mu} = w^{*}\cdot x^{\mu}/\sqrt{P}+ \epsilon^{\mu}$ with noise variance $D$; student $\hat w$ trained by the update [(6.11)](#eq2) with ridge $\lambda = 0$ and full batches ($\gamma = 1$), from $\hat w(0)\sim\CN(0,\sigma_{0}^{2} I_{P})$. Take $P = 200$ and $\nu = P/n \in \{0.5,\, 0.9,\, 1.1,\, 2\}$.
**(a)** Plot training and generalization error against time, and verify the asymptotics [(6.18)](#alpha-1ge) and [(6.19)](#alpha-1ge) — including the divergence of *both* branches as $\nu\to1$: double descent, resolved in time.
**(b)** Probe ergodicity breaking: run two clones from initializations differing by a small random $\delta\hat w(0)$, and track the normalized separation $|\Delta \hat w(t)|^{2}/|\Delta\hat w(0)|^{2}$. Verify that it decays to zero for $\nu<1$ but saturates for $\nu>1$ — at the predicted plateau $(\nu-1)/\nu$, the paper's ergodicity-breaking degree, measured in twenty lines of code.
**(c)** Near $\nu=1$, extract the relaxation time of $\varepsilon_{g}(t)$ and check for the critical slowing down $\tau \sim |1-\nu|^{-\phi}$ of Eq. [(6.21)](#udc).
### 6.3 Dynamical mean-field theory (DMFT)
Dynamical mean-field theory, of which [Section 6.2](#62-double-descent-as-a-dynamical-phase-transition) was one instance, is a general strategy for reducing the number of effective degrees of freedom in an interacting dynamical system: one follows a *single* representative degree of freedom, moving in a fluctuating background that stands in for all the others — and one demands that the statistics of the background be determined, self-consistently, by the statistics of the representative trajectory itself. This is the same trichotomy as every other mean field in these notes (couple through averages; argue fluctuations are controlled; close self-consistently), upgraded from scalars and measures to two-time functions. An excellent intuition-building tutorial in a biophysics context is ([Blumenthal 2025](#bib-blumenthalbuildingintuitiondynamical2025)), and ([Cugliandolo 2023](#bib-cugliandolorecentapplicationsdynamical2023)) is a broad modern review. In quantum many-body theory the same idea has been enormously effective for otherwise intractable systems ([Georges et al. 1996](#bib-georgesdynamicalmeanfieldtheory1996); [Kotliar et al. 2006](#bib-kotliardmft)). Below we first clarify what *self-consistency* means in a dynamical setting, then survey the recent deep-learning literature; a natural technical entry point for the latter is the out-of-equilibrium DMFT of the perceptron ([Agoritsas et al. 2018](#bib-agoritsas)).
#### 6.3.1 Self-consistency
The term *self-consistency* describes the demand that the predictions of a theory agree with the approximations used to obtain them. A simple example, taken from ([Stapmanns et al. 2020](#bib-stapmannsselfconsistentformulationsstochastic2020)), makes the issue vivid (and is worked in detail in [Exercise 6.3](#ex-selfconsistency)). Consider the stochastic process $d x (t) = f (x(t))\, d t + dW(t)$, with $W$ a Wiener process (so that $\dot W$ is the white noise of [Section 6.1](#61-from-stochastic-dynamics-to-path-integrals-the-msrjd-formalism)), and suppose $f(x_{0}) = 0$ with $f'(x_{0})<0$: a stable fixed point of the noiseless dynamics. One can expand around it, writing $x(t) = x_{0} + \delta x(t)$, to obtain the approximation
$$
\dv{\,}{t}\delta x(t) = f'(x_{0}) \delta x(t) + \xi(t) + \mathcal{O}(\delta^{2} x(t))\,,
$$
where $\xi$ is also centered gaussian white noise. The mean and variance of these fluctuations can be computed straightforwardly via Fourier transform, and in the temporal basis the variance is $\expval{\delta x^2}= \frac{D}{- 2 f'(x_{0})}$, however using these results with the Taylor expansion of $f(x)$ to second order gives a correction to the mean
$$
\expval{x}\approx x_{0} + \frac{1}{2}\frac{f''(x_{0})}{f'(x_{0})}\frac{D}{2 f'(x_{0})}
$$
which is *inconsistent* with the value $x_{0}$ assumed for the mean when constructing the approximation: the noise, filtered through the curvature of $f$, shifts the very point we expanded around.
**Exercise 6.3 (Self-consistency in the smallest possible example).** **(a)** Derive the stationary variance $\expval{\delta x^2}= D/(-2f'(x_{0}))$ of the linearized dynamics (compare [Exercise 6.1(b)](#ex-msrjd-ou)).
**(b)** Carrying the Taylor expansion of $f$ to second order, derive the quoted correction to the mean.
**(c)** Propose the self-consistent repair: expand instead around a point $x_{*}$, chosen so that the *effective* drift vanishes, $f(x_{*}) + \tfrac{1}{2}f''(x_{*})\expval{\delta x^2}_{*} = 0$, with $\expval{\delta x^2}_{*} = D/(-2f'(x_{*}))$ determined simultaneously. Note the structure: two coupled equations, each an input to the other — a miniature of the Curie–Weiss condition of [Section 5.1](#51-mean-field-theories-in-physics), and of the DMFT equations to come.
A framework for resolving this, and also attempting to enforce that the small parameter we expand about is actually the expansion parameter involves using *n-particle irreducible effective actions* (nPI). A full discussion of this is well beyond the scope of this document, but is covered very well in a pedagogical introduction ([Berges 2004](#bib-bergesintroductionnonequilibriumquantum2004)), where the main idea is enforcing self-consistency via increasingly multi-local sources that generate n-point correlation functions. The first section emphasizes that performing standard perturbation theory with an anharmonic oscillator yields so-called *secular* terms which are powers of time $t$ that appear next to the actual expansion parameter $\epsilon$, and using the nPI framework effectively resums these terms properly at each order in $\epsilon$. The 1PI approach is directly related to a local approximation with mean field theory, while in the 2PI approach, the object determined self consistently is the 2-point function, which in the local approximation is related to the local propagator and is the basis of dynamical mean field theory.
#### 6.3.2 Recent applications in deep learning
Fortunately, one rarely needs the general nPI machinery in practice. The working route — structurally the same trick by which large-$N$ models in physics (*e.g.* the $G$–$\Sigma$ formulation of the SYK model) are solved — is to insert bilocal *collective fields* for the correlation and response functions via delta functions, exactly as the $Q$'s of [Section 6.2](#62-double-descent-as-a-dynamical-phase-transition), and then to evaluate their path integral by saddle point in a suitable $N\to\infty$ limit. This is how the deep-learning DMFT of ([Bordelon & Pehlevan 2022](#bib-bordelon2022self); [Bordelon et al. 2024](#bib-bordelon2024dynamical)) is constructed.
A brief, opinionated tour of the literature, oldest first.
*Spin-glass origins.* The mean-field dynamics of disordered systems is where all of this machinery was forged. Amit, Gutfreund, and Sompolinsky ([Amit et al. 1985](#bib-amitspinglassmodelsneural1985)) brought spin-glass methods into neural networks (in the guise of Hopfield-style associative memories), and Sompolinsky, Crisanti, and Sommers ([Sompolinsky et al. 1988](#bib-sompolinskychaosrandomneural1988)) derived the DMFT of random recurrent networks, discovering a transition to chaos. The structural point to retain: an ordered system like the Ising magnet has one global coupling $J\sum_{\langle ij\rangle}\sigma_{i}\sigma_{j}$ and is already hard; a *disordered* system has random couplings $\sum_{ij}J_{ij}\sigma_{i}\sigma_{j}$, and it is precisely the all-to-all randomness that makes mean-field methods exact — the same structural feature that data-dependent couplings give learning problems, starting with the perceptron of [Section 2](#2-statistical-mechanics-of-learning).
*The Bordelon–Pehlevan program.* For deep networks trained by gradient descent at fixed dataset, the DMFT of ([Bordelon & Pehlevan 2022](#bib-bordelon2022self)) takes the infinite-width limit and finds that the order parameters are precisely the *feature kernels* $\Phi(t,s)$ and *gradient kernels* $G(t,s)$ — two-time generalizations of the NNGP kernel and NTK of [Sections 3–4](#3-large-width-at-initialization-the-gaussian-process-limit). The lazy regime of [Section 4](#4-the-ntk-and-large-width-during-training) is recovered when these kernels freeze; rich, feature-learning dynamics is their genuine evolution, computed self-consistently. (Fair warning: the notation is heavy — response fields and conjugate order parameters abound — but the structure is exactly MSRJD plus a saddle point.) The program now includes finite-width kernel fluctuations ([Bordelon & Pehlevan 2023](#bib-bordelondynamicsfinitewidth2023)), hyperparameter transfer across *depth* ([Bordelon et al. 2023](#bib-bordelondepthwisehyperparametertransfer2023)), and the dynamical theory of neural scaling laws ([Bordelon et al. 2024](#bib-bordelon2024dynamical)) that we meet again in [Section 7.9](#79-all-the-knobs-at-once-a-dynamical-theory).
*The Mignacco–Urbani line.* A complementary series of works treats the *data* as quenched disorder at proportional scalings of dimension and sample size, in the spirit of [Section 2](#2-statistical-mechanics-of-learning): DMFT for Gaussian-mixture classification ([Mignacco et al. 2021](#bib-mignaccodynamicalmeanfieldtheory2021)), and a DMFT characterization of the effective noise of SGD at late times ([Mignacco & Urbani 2022](#bib-mignaccoeffectivenoisestochastic2022)) — the dynamical refinement of the batch-noise discussion of [Section 2.1](#21-the-bayesian-model-of-learning).
*Rigor.* What began as a physics heuristic is increasingly a theorem: rigorous derivations of DMFT for first-order methods now exist ([Celentano et al. 2021](#bib-celentano2021high); [Gerbelot et al. 2024](#bib-gerbelot2024rigorous)), with a growing line of extensions — e.g. to adaptive Langevin diffusions, with propagation-of-chaos guarantees and convergence of the linear response ([Fan et al. 2025](#bib-fan2025adaptive)). A recent tour de force in this direction ([Montanari & Urbani 2025](#bib-montanaridynamicaldecouplinggeneralization2025)) solves the dynamics of a two-layer network at scale — a long calculation, but a well-organized paper — finding: the emergence of a slow timescale tied to the growth of the network's complexity (in the Gaussian/Rademacher sense); an inductive bias toward low complexity when the initialization has low complexity; a dynamical *decoupling* of the feature-learning and overfitting regimes; and a non-monotone test error, with a "feature unlearning" regime at very late times.
The connective tissue of this whole section, worth stating once more: the MSRJD path integral is to training dynamics what the Gibbs measure was to Bayesian learning — the object on which every physics technique (saddle points, disorder averages, collective fields, perturbation theory) can be brought to bear.
## 7. Scaling laws
In this final section we return to the theme of [Section 1](#1-intro-reading-learning-mechanics) and discuss the *neural scaling laws*. Over many orders of magnitude, the test loss of large models falls as a clean power law in each resource spent on it: parameters, data, and compute. These laws were found by fitting curves to measurements, and they have proved to be a remarkably simple, robust, and consequential empirical observation, one that any successful model of deep learning should be able to capture. Many models do capture their form, and several go further and make quantitative predictions for the exponents, based on features of the data distribution.
The organizing question of the section is: *where do the exponents come from?* We first assemble the empirical picture ([Sections 7.1–7.2](#71-the-empirical-laws)): the three laws, what each one holds fixed, and the joint laws that govern compute-optimal training, where the Kaplan–Chinchilla discrepancy lived and was resolved. We then turn to theory. [Section 7.3](#73-where-do-the-laws-come-from) states the mechanism common to every proposal — power laws in the data produce power laws in the loss — and [Section 7.4](#74-a-survey-of-theoretical-proposals) surveys the main proposals in one table. The rest of the section works out four rows of that table, in order of increasing machinery: the frequency statistics of skills ([Section 7.5](#75-the-quantization-model-exponents-from-the-statistics-of-skills), run as a derive-it-yourself exercise, with its payoff for "emergent capabilities" in [Section 7.6](#76-emergence-smooth-laws-from-staggered-cliffs)); the geometry of data and the spectra of kernels ([Section 7.7](#77-exponents-from-geometry-and-from-spectra)); and solvable large-$N$ field theories, static and then dynamical ([Sections 7.8–7.9](#78-a-solvable-large-n-field-theory)). We close with the newest measurements, for which there is no theory yet: post-training and test-time scaling ([Section 7.10](#710-beyond-pretraining-post-training-and-test-time-scaling)).
**Notation.** We keep the conventions of these notes: $P$ is the parameter count and $n$ the dataset size. Be warned that the scaling-law literature writes $N$ for our $P$ and $D$ for our $n$ (their $N$ is *not* our width!), so the famous exponents $\alpha_{N},\alpha_{D}$ are our $\alpha_{P},\alpha_{n}$. For reading the sources:
| these notes | scaling-law literature ([Kaplan et al. 2020](#bib-kaplan2020scaling); [Hoffmann et al. 2022](#bib-hoffmann2022training); [Michaud et al. 2023](#bib-michaud2023quantization); [Maloney et al. 2022](#bib-maloney2022solvable)) |
| --- | --- |
| $P$ (parameters) | $N$ |
| $n$ (dataset size) | $D$ |
| $\alpha_{P}$, $\alpha_{n}$ (exponents) | $\alpha_{N}$, $\alpha_{D}$ |
| $m$ (quanta learned, [Section 7.5](#75-the-quantization-model-exponents-from-the-statistics-of-skills)) | $n$ |
Note also that in this section $\alpha$ (no subscript) denotes the Zipf exponent of the quantization model, not the load of [Section 2](#2-statistical-mechanics-of-learning).
### 7.1 The empirical laws
Kaplan et al. ([Kaplan et al. 2020](#bib-kaplan2020scaling)) established, for autoregressive language models, that when a model is not bottlenecked by the other resources its test loss obeys
where $P$ is the number of (non-embedding) parameters, $n$ the dataset size in tokens, and $C$ the training compute in FLOPs. The measured exponents were $\alpha_{P} \approx 0.076$, $\alpha_{n}\approx 0.095$, and $\alpha_{C}\approx 0.050$, with the trends holding over up to seven orders of magnitude. Three points of bookkeeping are usually left implicit. They are worth making explicit before asking where the laws come from.
*Why FLOPs, and why $C \approx 6Pn$.* Compute is measured in floating-point operations because FLOPs are the only measure of training effort that does not depend on how the job is run. Wall-clock time depends on the hardware and the parallelism. The number of training *steps* depends on the batch size: the same $n$ tokens can be one giant batch or a million small ones, and below the "critical batch size" (which ([Kaplan et al. 2020](#bib-kaplan2020scaling)) measure separately) the loss depends only on the number of tokens processed. FLOPs count the arithmetic itself, which is also what energy and money are proportional to. The standard estimate is
$$
C \;\approx\; 6\,P\,n\,.
$$
Per token, the forward pass costs about $2P$ FLOPs: essentially every parameter sits in some matrix multiplication and is used once, for one multiply and one add. The backward pass costs about twice the forward pass, since each layer needs the gradient with respect to its activations and the gradient with respect to its weights, two matrix multiplications of the same size as the forward one. That gives $6P$ FLOPs per token, times $n$ tokens. The $\approx$ absorbs attention (whose per-token cost grows with context length rather than parameter count), embeddings, and optimizer overhead.
*The single-epoch assumption.* The estimate [(7.2)](#six-nd) assumes that every token processed is a fresh token, so that "tokens processed" and "dataset size" are the same number. This was true of the runs on which the laws were measured, which used well under one epoch of the corpus. Conceptually the two are different resources. The *statistical* resource is the number of distinct tokens $n$; the *computational* resource is the number of tokens processed, $T = E\,n$ over $E$ epochs, so that in general $C\approx 6PT$. Fixing $n$ and growing $C$ by repetition buys optimization, not statistics: the loss converges to the $n$-limited floor and then overfits. In the language of [Section 2.1](#21-the-bayesian-model-of-learning), more epochs means larger $\beta$ at a fixed quenched sample, while fresh data changes the quenched free energy itself. Empirically the one-epoch idealization is forgiving: for up to about four epochs, repeated tokens are nearly as valuable as fresh ones, after which their value decays quickly ([Muennighoff et al. 2023](#bib-muennighoff2023scaling)). This matters increasingly, now that frontier training is data-constrained.
*The compute law is not independent.* Since $C\approx 6Pn$ depends only on the other two resources, how can compute ever be the unbottlenecked one? The answer is that $L(C)$ is measured by a different protocol. It is the lower envelope
traced out by training many model sizes and keeping the best loss at each value of compute. On this curve compute is the binding constraint by construction, and "not bottlenecked" means "optimally allocated." This makes $\alpha_{C}$ a derived exponent. Minimizing an additive joint law $L = A\,P^{-\alpha_P}+ B\,n^{-\alpha_n}$ at fixed $Pn$ balances the two terms, and gives
$$
\frac{1}{\alpha_{C}}\;=\; \frac{1}{\alpha_{P}}+ \frac{1}{\alpha_{n}}\,, \qquad P \propto C^{\,\alpha_n/(\alpha_P + \alpha_n)}\,,\quad n \propto C^{\,\alpha_P/(\alpha_P + \alpha_n)}\,.
$$
The harmonic sum makes $\alpha_{C}$ smaller than either single-resource exponent, since a factor of $C$ must be split between growing $P$ and growing $n$. Kaplan's values give $\alpha_{C}\approx 0.042$, against the measured $0.050$: roughly consistent, with the gap traceable to the form of the fit and to early-stopping details.
*Universal form; data-dependent exponents.* The same power-law form holds far beyond language. Henighan et al. ([Henighan et al. 2020](#bib-henighan2020scaling)) trained autoregressive transformers on many kinds of data and found power laws everywhere, but with exponents that depend on the data:
| modality | $\alpha_{P}$ | $\alpha_{C}$ |
| --- | --- | --- |
| language | $0.070$ | $0.048$ |
| images $32{\times}32$ | $0.13$ | $0.10$ |
| images $8{\times}8$ | $0.24$ | $0.19$ |
| video | $0.24$ | $0.14$ |
| math (procedural) | $0.16$ | $0.17$ |
("Math" here means procedurally generated problems, in the style of the DeepMind mathematics dataset. They are emitted by a known program, so the supply is effectively unlimited. Remember this case: it is the one modality whose data distribution is known in closed form.) The exponents vary by a factor of $3$–$5$ across modalities. Architecture, by contrast, barely matters: within reasonable architectures, changing the shape (depth vs. width at fixed size) or even the family moves the constants and the plateaus, not the slopes ([Kaplan et al. 2020](#bib-kaplan2020scaling)). This asymmetry is the central clue of the section: *whatever sets the exponent lives in the data*. [Section 7.7](#77-exponents-from-geometry-and-from-spectra) makes it quantitative.
### 7.2 Joint laws, bottlenecks, and the Chinchilla revision
Each single-resource law assumes the other resources are unlimited. Joint fits describe the crossover to being bottlenecked, with the loss governed by whichever resource is scarcest. Two parametric forms dominate the literature. Kaplan et al. used
while the Chinchilla analysis of Hoffmann et al. ([Hoffmann et al. 2022](#bib-hoffmann2022training)) refit with the additive ansatz
$$
L(P,n) \;\approx\; E + A\, P^{-0.34}+ B\, n^{-0.28}\,.
$$
The additive form is easy to read. $E\approx 1.69$ is the irreducible entropy of text, the noise floor that no amount of scale removes; the two deficits decay independently; and the bottleneck is simply the larger deficit. Both forms reduce to the single-resource laws [(7.1)](#scaling-laws) in the appropriate corners.
The difference between the fits mattered a great deal in practice, because of what they say about the compute-optimal allocation in [(7.3)](#compute-envelope). Kaplan et al.'s analysis gave $P\propto C^{0.73}$ and $n\propto C^{0.27}$: with more compute, grow the model. For two years the field trained ever-larger models on comparatively little data. Chinchilla's fit gives $P\propto C^{0.5}$ and $n\propto C^{0.5}$: grow both equally, or as a rule of thumb about $20$ tokens per parameter. The prediction was tested directly. At the compute budget of DeepMind's 280B-parameter Gopher, the new rule said to train a model four times smaller on four times more data, and the resulting 70B-parameter Chinchilla beat Gopher across the board. (Both models, and the law, are named in DeepMind's animal series.) Chinchilla's allocation is now the accepted baseline. The discrepancy was later dissected by Porian et al. ([Porian et al. 2024](#bib-porian2024resolving)): Kaplan's fit was skewed by three small-scale artifacts — how last-layer and embedding compute were counted, over-long warmup, and optimizer settings that were not retuned with scale — and with these corrected, the two experiments agree. One coda: frontier models now deliberately *over*-train relative to Chinchilla, on many more than $20$ tokens per parameter, because a smaller model at equal loss is cheaper at inference time, a cost the training-compute-optimal analysis does not see.
With the empirical picture assembled, here is what a theory of [(7.1)](#scaling-laws) owes us: (i) why power laws at all, over so many decades; (ii) the values of the exponents (why so small? why do they depend on the modality but not on the architecture?); (iii) relations between exponents (is $\alpha_{n}$ determined by $\alpha_{P}$?); and (iv) when the laws break. The models surveyed below each answer a subset.
### 7.3 Where do the laws come from?
The basic answer, common to every proposal surveyed below, is that the loss curve reads off the coarse distribution of structure in the data. Power laws go in with the data and come out in the loss.
The prototype of "power laws in" is *Zipf's law*: rank the words of a large corpus by frequency, and the frequency of the word of rank $k$ falls off as $f_{k}\propto 1/k$, over about four decades of rank ([Figure 7](#fig-zipf)). It is a scale-free tail of ever rarer structure, and it is not special to words. Sequences of $n$ consecutive words ($n$-grams, the unit of pre-neural language modeling) are Zipf-like at every $n$, with slopes that flatten as $n$ grows. Facts, meaning discrete pieces of world knowledge attested in text, inherit the heavy-tailed popularity of the entities they mention, and a model's ability to answer a factual question tracks the number of training documents that support it ([Kandpal et al. 2023](#bib-kandpal2023large)). The covariance and kernel spectra of text and image data are measured power laws ([Sections 7.7–7.8](#77-exponents-from-geometry-and-from-spectra)). A learner working through such data never runs out of structure to absorb. It moves ever further down the tail, at ever slower returns. That is, qualitatively, why the loss is a power law rather than an exponential (there is structure at all scales) and why it does not simply plateau (the tail does not end).
Zipf's law, illustrated on a synthetic corpus: $2\times10^{6}$ tokens drawn from a $1/k$ distribution over a vocabulary of $5\times10^{4}$ words, with the empirical rank–frequency curve against the $f\propto 1/\mathrm{rank}$ guide. The curve is a consistency check rather than data, since the words were sampled from a Zipf distribution to begin with. Real corpora look the same over roughly four decades, with deviations at both ends: the Zipf–Mandelbrot correction at low rank, and a steeper tail beyond about $10^{4}$ words.
The exponents of these tails differ, and the differences matter, because in every model below the loss exponent is inherited from the data exponent. Pure Zipf is in fact a marginal case. Write the tail of skill frequencies as $p_{k}\propto k^{-(\alpha+1)}$. The quantization model of [Section 7.5](#75-the-quantization-model-exponents-from-the-statistics-of-skills) gives a loss exponent equal to $\alpha$, and exact Zipf ($f_{k}\propto 1/k$, that is, $\alpha = 0$) is the boundary at which the tail sum decays only logarithmically and power-law scaling degenerates. Read this way, the measured $\alpha_{P}\approx 0.076$ for language says that its effective skill distribution is $p_{k}\sim k^{-1.08}$, barely steeper than pure Zipf. Language scaling is slow because its structure sits just above the marginal $1/k$ tail. Image data has steeper tails and scales faster; compare the modality table of [Section 7.1](#71-the-empirical-laws).
Where does the theory stand? The single-resource laws have been derived from this mechanism by many models of quite different flavor, and a few models derive multi-resource laws; the survey below keeps score. Three things are missing. First, no theory computes the data's exponents from first principles. Every model converts measured data statistics into loss curves, and the statistics themselves are an empirical input. This is true even for procedurally generated math, where the generating program is known exactly but the computation has not been carried out. Second, there are no joint laws with genuine feature learning. The best dynamical derivations live in lazy, random-feature settings ([Section 7.9](#79-all-the-knobs-at-once-a-dynamical-theory)), while in solvable cases feature learning is known to improve exponents for structured targets, so the corrections need not be small. Third, the bridge from the smooth loss to capabilities is missing. Benchmark performance can jump within a decade of scale while the loss glides smoothly ([Wei et al. 2022](#bib-wei2022emergent)); this has been contested as an artifact of the metric ([Schaeffer et al. 2023](#bib-schaeffer2023emergent)), and the quantization model of [Section 7.6](#76-emergence-smooth-laws-from-staggered-cliffs) gives it a sharper reading. A fourth gap is newer: the post-training and test-time scaling laws of [Section 7.10](#710-beyond-pretraining-post-training-and-test-time-scaling) have no theory at all.
### 7.4 A survey of theoretical proposals
The table below organizes the most prominent derivations of neural scaling laws by four questions: what models the *data* (equivalently, what a "feature" is and how features are distributed); what models the *network*, if anything; whether actual *exponents* are predicted, and in terms of what; and whether *multi-resource* laws are obtained.
| **proposal** | **data model / "features"** | **NN framework** | **exponents?** | **multi-resource?** |
| --- | --- | --- | --- | --- |
| quantization; Michaud et al. ([Michaud et al. 2023](#bib-michaud2023quantization)) | discrete skill "quanta," Zipf-used, learned in order | none (a counting argument) | in terms of the Zipf exponent $\alpha$; plus the relation $\tfrac{1}{\alpha_n}= \tfrac{1}{\alpha_P}+ 1$ | separate $P$, $n$, $S$ laws |
| data manifold; Sharma–Kaplan ([Sharma & Kaplan 2020](#bib-sharma2020neural)) | smooth target on a $d$-dimensional manifold | piecewise-linear interpolation (function class only; no training) | yes: $\alpha_{P}\approx 4/d$ | no ($P$ only) |
| percolation; Brill ([Brill 2024](#bib-brill2024scaling)) | latent subtask clusters at criticality, with power-law cluster sizes | none: a generic nonparametric learner $+$ regression estimates | from percolation critical exponents | two regimes; not joint |
| kernel eigenlearning; Bahri et al. ([Bahri et al. 2024](#bib-bahri2021explaining)), Bordelon–Canatar–Pehlevan ([Bordelon et al. 2020](#bib-bordelon2020spectrum); [Canatar et al. 2021](#bib-canatar2021spectral); [Simon et al. 2023](#bib-simon2023eigenlearning)) | power-law kernel spectra $+$ target coefficients | **frozen kernel**: infinite-width NTK, or GP/Bayesian | in terms of the spectral exponents | partially: resolution- vs. variance-limited regimes |
| solvable field theory; Maloney–Roberts–Sully ([Maloney et al. 2022](#bib-maloney2022solvable)) | latent power-law covariance, extended by random features | linear model on **fixed random features** (frozen kernel); large-$N$ RMT, statics | in terms of the data's spectral exponent | yes: joint $L(P,n)$, with duality $P \leftrightarrow n$ |
| DMFT; Bordelon–Atanasov–Pehlevan ([Bordelon et al. 2024](#bib-bordelon2024dynamical)) | power-law (source/capacity) spectra | DMFT of gradient-flow training on **fixed random features** (frozen kernel; dynamics) | yes, including the time exponent | fully joint $L(t,P,n)$ $+$ compute frontier |
Five remarks for reading the table. First, everyone agrees on the origin. In every row the loss exponent is inherited from a power-law tail in the data, whether of skill frequencies, cluster sizes, manifold volume, or covariance spectra. The rows differ in what the tail is a tail *of*, and in how much network machinery is retained. Second, and this deserves emphasis: *every row that models the network at all does so with a frozen kernel*. The eigenlearning row works with the infinite-width NTK or GP kernel of [Sections 4.3](#43-constant-ntk-at-large-width) and [4.6](#46-bayesian-learning-at-large-width-kernel-learning); the two field-theory rows train only the readout of a fixed set of random features, which is a frozen kernel by construction ([Exercise B.6(b)](#ex-app-tangent-readout)). The manifold row models only the network's function class, as a count of linear pieces, and never trains it; the quantization and percolation rows contain no network. So none of the current derivations involves feature learning, the very mechanism that [Sections 5–6](#5-mean-field-scaling) argued is what makes finite networks differ from kernels. Third, equilibrium routes (Bayesian GP learning curves and their finite-width refinements) and dynamical routes (gradient flow, DMFT) give matching exponents wherever they overlap. This is the "two limits of one system" correspondence of [Sections 2.1](#21-the-bayesian-model-of-learning) and [6](#6-non-equilibrium-physics-msrjd-and-dmft) at work again. Fourth, only the two field-theory rows currently predict joint laws, and only the dynamical one prices training time, deriving a Chinchilla-like compute-optimal frontier rather than assuming it. Fifth, no row computes the data's own exponent. The closest attempt is the percolation proposal, which derives power-law cluster-size statistics from a critical model of the data distribution; its two criticality regimes (power-law-distributed discrete subtasks, and a dominant data manifold) map onto the quantization and manifold rows respectively. We will not develop the percolation model further. The remaining subsections work through the other rows, in order of machinery.
*Do the predicted exponents agree with experiment?* Where the comparison has been made, the answer is encouraging but not yet sharp. The quantization model's authors extracted "quanta" from a small language model by clustering gradients, measured how often each is used in the training distribution, and found a power law whose exponent roughly matches the model's scaling exponent, as the model requires ([Michaud et al. 2023](#bib-michaud2023quantization)); on the other hand, the model's parameter-free relation between $\alpha_{P}$ and $\alpha_{n}$ is violated by the measured ordering of Kaplan's exponents ([Exercise 7.1(e)](#ex-quanta-falsify)). The manifold prediction $\alpha_{P}\approx 4/d$ was confirmed in teacher–student experiments, where $d$ and $\alpha_{P}$ can be measured independently, and tested more loosely on image classifiers and GPT-type language models ([Sharma & Kaplan 2020](#bib-sharma2020neural)); for language, taking $\alpha_{P}\approx0.076$ at face value gives $d\approx 50$ ([Exercise 7.2](#ex-manifold)). The kernel predictions have been tested on standard architectures and datasets, with exponents computed from measured kernel spectra, and the predicted duality between the $P$- and $n$-exponents holds up ([Bahri et al. 2024](#bib-bahri2021explaining)). The two field-theory rows take the data's spectral exponent as input. Their qualitative predictions — exponents independent of architecture, $\alpha_{P}\approx\alpha_{n}$, and equiparameterization as the compute-optimal rule — match the Chinchilla-era measurements, but a quantitative test needs the covariance spectrum of the training data, which is rarely reported for the large runs.
The honest summary is that these theories convert an exponent measured in the *data* into an exponent measured in the *loss*, and the conversion has held up wherever both have been measured.
### 7.5 The quantization model: exponents from the statistics of skills
The quantization model of Michaud, Liu, Girit, and Tegmark ([Michaud et al. 2023](#bib-michaud2023quantization)) is a minimal model of scaling. It has no architecture, no optimizer, and essentially no mathematics beyond a tail sum, and it produces all three laws of [(7.1)](#scaling-laws), plus a sharp point of view on "emergent capabilities." The name is deliberately provocative: the proposal is that what a network learns comes in *discrete units*, dubbed *quanta*.
[^4]
The assumptions, in the form we'll use:
- **QH1 (discreteness).** The computations and knowledge relevant to the prediction task decompose into enumerable *quanta*: basic skills or facts (*e.g.* "carry the one," "Paris follows `the capital of France is'," ...). Each quantum is either fully learned or not learned at all, with no partial credit.
- **QH2 (ordered acquisition).** Networks learn quanta roughly in order of decreasing usefulness, so a given model has learned the first $m$ quanta of a fixed sequence.
- **QH3 (Zipfian use).** Quantum $k$ is *used* on a fraction
of samples ($\alpha>0$ a constant), and on those samples, knowing it reduces the loss by a constant $\Delta$ (independent of $k$).
Consequently, a model that has learned the first $m$ quanta has expected loss
That is the entire model. Everything else is asking how $m$ is limited by each resource.
**Exercise 7.1 (★) (Three scaling laws from three assumptions).** **(a)** Estimate the tail sum in [(7.8)](#quanta-loss) to show $L_{m} - L_{\infty} \propto m^{-\alpha}$.
**(b)** *Parameter-limited (data and training unlimited).* Suppose each quantum costs a constant capacity of $c$ parameters to learn. Derive $\alpha_{P}$ in terms of $\alpha$.
**(c)** *Data-limited (capacity and training unlimited).* Suppose a quantum is learnable only if it is used by at least $\tau$ training samples, for constant $\tau$. Given $n$ i.i.d. samples, which quanta are learnable? Derive $\alpha_{n}$.
**(d)** *Training-limited, single epoch.* Suppose instead that learning quantum $k$ requires $\tau$ gradient updates on batches containing it. After $S$ steps of online training (fresh batches), derive $\alpha_{S}$.
**(e)** *Falsifiability.* Parts (b)–(c) predict a parameter-free relation between $\alpha_{P}$ and $\alpha_{n}$: write it down, and note which of the two exponents must be larger. Compare with the measured values quoted after [(7.1)](#scaling-laws). Something is wrong — what? List the assumptions you used and discuss which are most suspect. (This is a discussion question with no settled answer; see the solution for candidates.)
### 7.6 Emergence: smooth laws from staggered cliffs
The quantization model's deepest payoff is conceptual. Under QH1, the learning of any *single* quantum is a step function: a capability is absent, then abruptly present. Yet the aggregate loss [(7.8)](#quanta-loss) is a smooth power law. There is no contradiction. *A smooth law is a sum of many staggered cliffs*, each weighted by $p_{k}$ and individually invisible in the aggregate, since quantum $k$ contributes only the small fraction $p_{k}$ of the loss. Michaud et al. sharpen this with a genetics metaphor. A test sample is *monogenic* if its loss is governed by a single quantum, so that its loss curve during training is a plateau followed by a cliff; it is *polygenic* if many quanta contribute, in which case its loss falls smoothly. Both kinds of per-token curves are observed in real LLMs. From this perspective, "emergent capabilities" ([Wei et al. 2022](#bib-wei2022emergent)) are neither magic nor mirage ([Schaeffer et al. 2023](#bib-schaeffer2023emergent)): any benchmark dominated by a few quanta *must* show abrupt arrival, while any sufficiently aggregated metric *must* look smooth. Emergence is a property of what you measure, not (necessarily) of a phase transition inside the model.
Multitask sparse parity ([Michaud et al. 2023](#bib-michaud2023quantization)): smooth scaling from staggered cliffs. Setup: $T=48$ subtasks with Zipfian frequencies $p_{i}\propto i^{-(\alpha+1)}$, $\alpha=0.4$; inputs are a one-hot task indicator ($T$ bits) concatenated with $24$ random $\pm1$ bits; the label is the parity of a fixed random $3$-subset of the bits, the subset depending on the task. Model: one-hidden-layer ReLU MLP, width $512$, MSE loss on $\pm1$ labels, trained online (fresh batches of $256$; Adam, learning rate $2\times10^{-3}$, $1.2\times10^{5}$ steps; fixed seed). *Left:* test loss per subtask vs. training step. Each subtask is learned in an abrupt cliff, frequent tasks first. *Right:* the frequency-weighted aggregate $L(S)=\sum_{i}p_{i}L_{i}(S)$ is smooth; the dashed guide is the *idealized* quanta-model slope $-\alpha/(\alpha+1)$ from [Exercise 7.1(d)](#ex-quanta-training). The measured decay is visibly steeper, as discussed in the text.
**Lab: multitask sparse parity.** The entire mechanism fits in a laptop-scale experiment, shown in [Figure 8](#fig-sparse-parity); students are encouraged to reproduce and modify it (all details are in the caption). A network is trained online on $T$ separate parity subtasks with Zipfian frequencies. The input is a one-hot task indicator concatenated with random bits, and the label is the parity of a small task-specific subset of the bits. Parity is a needle-in-a-haystack skill with no partial credit, *i.e.* an honest hand-built quantum. The result: each subtask's test loss is a staggered cliff (frequent tasks fall first), while the frequency-weighted aggregate loss decays smoothly. Note also the honest failure visible in the right panel. The measured aggregate decays *more steeply* than the idealized quanta slope $-\alpha/(\alpha+1)$ of [Exercise 7.1(d)](#ex-quanta-training), because the cliffs bunch up rather than spreading Zipf-like; the rare tasks are learned sooner than the threshold rule $S\,p_{k}\gtrsim\tau$ predicts. The subtasks share a hidden layer and an adaptive optimizer, so they are not learned independently, and assumption (iv) in the solution to [Exercise 7.1(e)](#ex-quanta-falsify) fails.
### 7.7 Exponents from geometry and from spectra
The quantization model takes the Zipf exponent $\alpha$ as *input*. The next two arguments produce exponents as *output*, from properties of the data.
#### 7.7.1 Geometry: the data-manifold argument
Sharma and Kaplan ([Sharma & Kaplan 2020](#bib-sharma2020neural)) propose the simplest continuous mechanism. A network with ReLU activations computes a piecewise-linear function, so regression on a data manifold of intrinsic dimension $d$ is a problem of *piecewise-linear interpolation*. A network with $P$ parameters can afford $\sim P$ linear pieces; tiling a $d$-dimensional manifold with $P$ pieces gives spacing $s \sim P^{-1/d}$; a linear patch approximates a smooth target to accuracy $\sim s^{2}$ (Taylor); and MSE is the *square* of the error. Hence
Exponents are small because data manifolds are high-dimensional. The formula is testable, since $d$ can be measured with intrinsic-dimension estimators, which is exactly what ([Sharma & Kaplan 2020](#bib-sharma2020neural)) do.
**Exercise 7.2 (Exponents from the data manifold).** **(a)** Fill in the derivation of [(7.9)](#manifold-law), stating precisely where smoothness of the target is used.
**(b)** How does the answer change for mean-absolute-error loss?
**(c)** Taking the measured $\alpha_{P}\approx0.076$ for language models at face value, estimate the effective dimension of the "manifold of text." Discuss whether the input dimension, the embedding dimension, or something else is being measured.
**(d)** What does this mechanism predict for the *data* exponent $\alpha_{n}$, if $n$ samples are what limits the resolution of the tiling?
#### 7.7.2 Spectra: scaling laws from kernel eigenvalues
The second mechanism uses machinery we already own: the Bayesian kernel regression of [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning), whose mode-filtering formula said that with $n$ samples (and inverse temperature $\beta$), the modes of the GP kernel that get learned are those with $\lambda_{k} \gtrsim (n\beta)^{-1}$. Suppose the kernel spectrum and the target's decomposition are power laws,
$$
\lambda_{k} \sim k^{-a}\quad (a>1)\,,\qquad y_{k}^{2} \sim k^{-b}\quad (b>1)\,,
$$
in a shared eigenbasis $\{\varphi_{k}\}$. The unlearned modes dominate the test error:
giving $L(n) \sim n^{-(b-1)/a}$: *the scaling exponent is a ratio of spectral decay rates*. Everything about the data and the architecture enters through two measurable power laws. This is the starting point of the "eigenlearning" analyses ([Bordelon et al. 2020](#bib-bordelon2020spectrum); [Canatar et al. 2021](#bib-canatar2021spectral); [Simon et al. 2023](#bib-simon2023eigenlearning)), which sharpen [(7.11)](#spectral-law) into precise learning curves (including the overfitting factor whose divergence produces the double-descent peak of [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft)). Bahri et al. ([Bahri et al. 2024](#bib-bahri2021explaining)) organize the possibilities into a clean $2\times2$ taxonomy: *resolution-limited* regimes (exponents from spectra or geometry, as here and in [(7.9)](#manifold-law)) versus *variance-limited* regimes (universal exponent $1$), in each of $P$ and $n$.
**Exercise 7.3 (Spectral scaling laws, and quanta as eigenmodes).** **(a)** Derive [(7.11)](#spectral-law) from the mode-filtering formula of [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning), and check that the crossover modes $k\sim k_{*}$ do not change the exponent.
**(b)** Show that the quantization model's data-scaling law is the special case $b=a$ of [(7.11)](#spectral-law) under the dictionary: quanta $=$ eigenmodes, Zipf law $=$ spectrum, $a = \alpha+1$. (Hint: in the quanta model, $p_{k}$ plays *both* roles, setting the learning threshold *and* the contribution $\Delta\,p_{k}$ to the loss, so the two power laws coincide.)
**(c)** In the quanta model, what would it mean for $b \neq a$? Invent a modification of QH3 that realizes it.
### 7.8 A solvable large-N field theory
Maloney, Roberts, and Sully ([Maloney et al. 2022](#bib-maloney2022solvable)) supply what the arguments above lack: a *microscopic, exactly solvable* model in which the joint law $L(P,n)$ is computed, not modeled. The ingredients are: a generative data model in which observed inputs are (nonlinear functions of) random projections of a latent space, engineered so that the data covariance has a power-law spectrum $\lambda_{k} \sim k^{-(1+\zeta)}$, as measured in real image and text data (the spectral exponent is the paper's $\gamma$; we reserve $\gamma$ for the scaling dial of [Section 5.4](#54-the-scaling-dial-and-maximal-update-parameterization)); a random-feature *linear* student, so we are in the lazy limit of [Section 4](#4-the-ntk-and-large-width-during-training), whose $P$ random features are its trainable parameters, trained on $n$ samples; and the joint thermodynamic limit $P,n\to\infty$ at fixed ratio. The computation proceeds at large $P$ and $n$, where expectations over data and features become Gaussian integrals, the test loss becomes a product of resolvents (Green's functions), and the $1/P$ expansion organizes into planar ('t Hooft) diagrams that can be resummed exactly.
[^5]
*Frozen, but random.* A natural first reaction is that this model should be trivial. The student is a linear model on fixed random features. Its kernel never moves, so we are squarely in the lazy regime of [Section 4](#4-the-ntk-and-large-width-during-training), and in the strict infinite-width lazy limit the loss would not depend on $P$ at all. The $P$-scaling law is a *finite-size effect on the lazy limit*. At finite $P$ the feature Gram matrix is a *random* matrix, correlated with the random sample, and its rank caps how many modes of the data spectrum the student can resolve. All of the physics lives in the disorder average over this randomness, which is what the resolvents and the planar diagrams compute. Frozen does not mean deterministic.
Three notable results:
- *Where exponents come from:* the loss inherits the *data's* spectral exponent, $L \sim \min(P,n)^{-\zeta}$ up to logarithms. Remarkably, generic nonlinear feature maps *extend* the data's power-law spectrum rather than reshaping it. This explains why exponents are roughly architecture-independent and can only be improved by changing the data (or the task).
- *Duality and Chinchilla:* the loss is (to leading order) *symmetric* under $P\leftrightarrow n$. This predicts $\alpha_{P} \approx \alpha_{n}$, consistent with observation, and opposite in spirit to the quanta model's asymmetric prediction ([Exercise 7.1(e)](#ex-quanta-falsify)). It makes *equiparameterization* $P\propto n$ compute-optimal, a first-principles Chinchilla rule.
- *Where laws break:* scaling stops when $\min(P,n)$ exceeds the latent dimensionality of the data model, and the loss plateaus. Power laws are a statement about the resolvable structure in the data, and they end when the structure is exhausted.
### 7.9 All the knobs at once: a dynamical theory
Everything so far is statics, or one knob at a time. The current frontier is the *dynamical* mean-field theory of [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft) applied to scaling. Bordelon, Atanasov, and Pehlevan ([Bordelon et al. 2024](#bib-bordelon2024dynamical)) solve a random-feature model trained by gradient flow, with power-law data spectra, in the joint limit of large width, large data, and long time. They obtain a single self-consistent prediction for $L(t, P, n)$, with the correlation and response functions of [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft) as the order parameters. Each resource, when it is the binding constraint, contributes its own power law, recovering the resolution-limited exponents of [Section 7.7](#77-exponents-from-geometry-and-from-spectra) in the appropriate corners. Moreover, the theory predicts the full *compute-optimal frontier*: minimizing $L$ at fixed compute $C \propto P t$ yields Chinchilla-like allocation rules and the frontier exponent, rather than assuming them.
A fair question: the features here are again frozen and random, exactly the setting of [Section 7.8](#78-a-solvable-large-n-field-theory), so is dynamical mean-field theory overkill? No, for two reasons. First, time is a genuine third resource. Predicting $L(t,P,n)$ and the compute frontier requires the whole training trajectory, and the statics only sees the $t\to\infty$ corner. Second, the naive mode-by-mode solution of the lazy regime, in which each kernel eigenmode of the residual decays as $e^{-\eta\lambda_k t}$ as in [(4.34)](#exp-const), applies to a *deterministic* kernel. Here the kernel's eigenvalues and eigenvectors are random and correlated with the sampled data, so the disorder average couples the modes. The DMFT correlation and response functions are the dynamical generalization of the static resolvents of [Section 7.8](#78-a-solvable-large-n-field-theory), and at $t\to\infty$ the equations collapse back onto the static answer, as they must. What the theory is *not* yet doing is feature learning. The adaptive-kernel version of the formalism exists ([Bordelon & Pehlevan 2022](#bib-bordelon2022self)), but it has not been pushed through power-law data to closed-form scaling laws; these are the "joint laws with genuine feature learning" missing from [Section 7.3](#73-where-do-the-laws-come-from). Since feature learning is known to improve exponents for structured targets in solvable cases, the correction need not be small.
### 7.10 Beyond pretraining: post-training and test-time scaling
Frontier models are no longer just pretrained. They are fine-tuned on curated demonstrations (SFT), optimized against reward models or verifiable rewards (RLHF, RLVR), and allowed to spend variable compute at inference time. Do scaling laws extend to these stages? Increasingly, yes. But the functional forms change, and the changes are instructive.
*Transfer and fine-tuning* still look Kaplan-like. The value of pretraining for a downstream task can be quantified as the *effective data transferred*, meaning the amount of fine-tuning data that the pretrained initialization is worth, and it scales as a power law in both the fine-tuning dataset size and the model size ([Hernandez et al. 2021](#bib-hernandez2021transfer)).
*Reward-model optimization* introduces a new scaling axis: not data or parameters, but the distance moved from the base policy, $d = \sqrt{\mathrm{KL}(\pi \,\|\, \pi_{0})}$. Gao, Schulman, and Hilton ([Gao et al. 2023](#bib-gao2023reward)) found that the true ("gold") reward follows simple laws in $d$, namely $d\,(a - b\log d)$ for best-of-$n$ sampling and $d\,(a - b\,d)$ for RL, with coefficients that scale smoothly with the size of the reward model. Note the shape: not a power law but a rise and a *fall*. Optimizing a learned proxy eventually degrades the true objective. This is a Goodhart peak, quantified.
*RL on verifiable rewards*, the training recipe of the reasoning models, scales by *saturating sigmoids* rather than power laws. In the largest systematic study to date, performance as a function of RL compute is predictable, but by a fit with an asymptote, an efficiency, and a midpoint, and most recipe choices move the efficiency without changing the asymptote ([Khatri et al. 2025](#bib-khatri2025art)). A physics gloss suggests why the form differs, though it is a hypothesis rather than a result. Pretraining mines an effectively unbounded Zipf tail of new structure, hence power laws. RL post-training mostly elicits and reweights capabilities the base model already contains, on a bounded task distribution. In quanta language, pretraining acquires quanta, while RL sharpens the policy over quanta it already has, and one runs out. *Test-time compute* is the other new axis: the coverage achieved by repeated sampling grows roughly as a power law in the number of samples ([Brown et al. 2024](#bib-brown2024monkeys)), and train-time and test-time compute can be traded against each other along a compute-optimal frontier ([Snell et al. 2024](#bib-snell2024test)).
What does not yet exist is a joint law over all four budgets (pretraining, SFT, RL, inference), and any *theory* at all, even at the level of the quantization model, that derives the sigmoid or the Goodhart peak from a model of the data and skill distribution. The empirical forms are waiting for a statistical-mechanics treatment.
That every ingredient of the pretraining synthesis — Gibbs measures, kernels, mean-field limits, MSRJD actions — appeared earlier in these notes as a physics import is, we hope, a fair summary of the state of the field. A theory of the scaling laws is not finished, but it is recognizably being built, and it is being built out of statistical physics. The post-training laws are younger still: their regularities are only now being settled, and there is no theory of them yet. For a student entering the field, that is a map of where the open territory begins.
## A. The replica method
[Section 2.3](#23-partition-functions-free-energy-and-cumulants) introduced the quenched free energy $-\E_{\rm data}\log Z$ and the annealed approximation [(2.18)](#annealed-swap) that replaces it by $-\log\E_{\rm data}Z$, and [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning) used that approximation to analyze the Ising perceptron. This appendix explains how the quenched average is computed exactly, by the *replica method*, and what the annealed approximation misses. We keep the Ising perceptron of [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning) as the running example, in the notation and normalization of that section: $N$ binary weights $\theta\in\{\pm1\}^{N}$, a binary teacher $\theta^{*}$, Gaussian inputs, $n=\alpha N$ examples, and the partition function [(2.39)](#perc-z). The classic reference for everything here is Engel and Van den Broeck ([Engel & Van den Broeck 2001](#bib-engel2001statistical)).
### A.1 What the annealed approximation gets wrong
Jensen's inequality gives $\log\E Z\geq\E\log Z$, so the annealed free energy $-\log\E Z$ is always a lower bound on the quenched one. To see what drives the gap, write $Z = e^{-N\varphi}$, where the free-energy density $\varphi$ fluctuates from dataset to dataset. Exponentiation turns modest fluctuations of $\varphi$ into enormous fluctuations of $Z$: a dataset whose $\varphi$ is below the typical value by a constant is exponentially rare, but it contributes an exponentially large $Z$, and the two exponentials compete at the same order in $N$. The average $\E Z$ can therefore be dominated by vanishingly rare datasets, while the typical dataset, which is what $\E\log Z$ describes, behaves differently. (This is the familiar behavior of lognormal random variables, whose mean sits far in the upper tail. The mean is a poor summary of a heavy-tailed quantity.)
For the Ising perceptron the direction of the error follows. The datasets that dominate $\E Z$ are those that constrain the student unusually weakly, and weakly constrained students sit at low overlap. The annealed average therefore overweights poorly aligned students, and overestimates the amount of data needed to learn: it puts the transition to perfect generalization at $\alpha_{c}^{\rm ann}\approx1.45$, where the exact calculation sketched below gives $\alpha_{c}\approx1.245$ ([Gy"orgyi 1990](#bib-gyorgyi1990first)). The structure of the transition, two basins exchanging dominance, is the same in both.
**Remark.** Annealing does not always err so gently. Here the model is realizable (a teacher exists) and the fluctuations of $\varphi$ are tame, so annealing shifts a number without changing the picture. In general the annealed computation can fail qualitatively, most famously for models with discrete weights at large load, where it can predict negative entropies (fewer than one consistent student, nonsense for a count) and spurious transition points.
### A.2 Replicas
The full calculation computes $\E_{\rm data}\log Z$ through the identity
$$
\E\log Z = \lim_{Q\to0}\frac{\E Z^{Q} - 1}{Q}\qquad\big(\text{from }Z^{Q} = e^{Q\log Z}= 1 + Q\log Z + O(Q^{2})\big)\,.
$$
The strategy is to compute $\E Z^{Q}$ for *integer* $Q$, where it has a transparent meaning, obtain a formula analytic in $Q$, and continue it to $Q\to0$. For integer $Q$, $Z^{Q}$ is the partition function of $Q$ independent copies, or *replicas*, of the student, all judged against the same dataset. For the Ising perceptron, from [(2.39)](#perc-z),
The data average is easy, because the inputs are Gaussian. For fixed replicas and teacher, each example enters only through the $Q+1$ numbers $u^{a} = \theta^{a}\cdot x^{\mu}/\sqrt{N}$ and $v = \theta^{*}\cdot x^{\mu}/\sqrt{N}$, which are exactly jointly Gaussian (they are linear in $x^{\mu}$), with mean zero and a covariance determined entirely by the *overlap matrix*
The product of step functions asks that every $u^{a}$ have the sign of $v$, and the probability of that event is some function $e^{-e(q,R)}$ of the overlaps alone: a Gaussian orthant probability in $Q+1$ dimensions. Since the examples are independent, the data average of the product over $\mu$ is the $n$-th power, $e^{-n\,e(q,R)}= e^{-\alpha N\,e(q,R)}$. For $Q=1$ the only overlap is $R$, the orthant is a wedge, and $e(R) = -\log(1-\varepsilon(R))$ is the energy term of [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning).
The disorder is gone. In its place, the replicas are coupled through the overlaps $q^{ab}$, and the sum over the $Q$ binary vectors can be organized by the value of the overlaps, exactly as the sum over one binary vector was organized by $R$ in [Exercise 2.3(c)](#ex-ising-count). Writing $e^{N\,h(q,R)}$ for the number of $Q$-tuples of students with prescribed overlaps (for $Q=1$ this is the binomial count of [Exercise 2.3(c)](#ex-ising-count)), and replacing the sum over overlaps by an integral,
This is the replicated version of the one-variable action [(2.44)](#perc-action) of [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning), to which it reduces at $Q=1$. At large $N$ the integral is dominated by the minimum of $s_{Q}$ over the overlaps, and the replica identity [(A.1)](#replica-identity) then gives the quenched free-energy density
$$
\Phi(\alpha) := -\frac{1}{N}\,\E_{\rm data}\log Z \;=\; \lim_{Q\to0}\;\frac{1}{Q}\,\min_{q,R}\,s_{Q}(q,R)\,.
$$
Compare the annealed approximation, which is $\min_{R} s_{1}(R)$: the annealed calculation keeps a *single* replica, and so never sees the coupling $q^{ab}$ between replicas that the shared data induce. That coupling is precisely the information about dataset-to-dataset fluctuations that [Appendix A.1](#a1-what-the-annealed-approximation-gets-wrong) found missing.
*Replica symmetry.* To make [(A.5)](#replica-saddle) concrete one needs $h$ and $e$ as functions of the $Q(Q+1)/2$ overlaps, and a way to continue the result in $Q$. The standard ansatz, *replica symmetry*, is to take $q^{ab}= q$ for all $a\neq b$ and $R^{a} = R$. It is natural, since the saddle-point equations are symmetric under permuting replicas, and it reduces the problem to two scalars. Under it, both $h$ and $e$ become one-dimensional Gaussian integrals: the orthant probability collapses to an integral over a single shared Gaussian "field" of the Gaussian tail function $H(x) = \int_{x}^{\infty}\frac{dz}{\sqrt{2\pi}}e^{-z^2/2}$, and the count $h$ is handled by introducing conjugate variables for the constraints, which reduces it to a one-site problem in the manner of the mean-field Ising model of [Section 2.4](#24-a-physics-example-the-ising-model). The resulting expression is analytic in $Q$ and can be differentiated at $Q=0$. As a caution, symmetric equations *can* have asymmetric solutions, and this happens in physical systems, famously in spin glasses; the assumption must be checked, and for the perceptron problems considered here it holds.
*The result.* For the Ising perceptron, the replica-symmetric calculation ([Gy"orgyi 1990](#bib-gyorgyi1990first)) reproduces the two-basin structure found in [Exercise 2.3](#ex-ising): an interior solution with $R^{*}(\alpha)<1$ that drifts toward the teacher as data accumulate, and the teacher solution $R=1$, which exchange dominance in a first-order transition. Only the number moves, from $\alpha_{c}^{\rm ann}\approx1.45$ to $\alpha_{c}\approx1.245$. A guided walk through the calculation follows.
**Exercise A.1 (The quenched calculation, via replicas).** Work through the following program for the Ising perceptron, in as much detail as your appetite allows; full details are in ([Engel & Van den Broeck 2001](#bib-engel2001statistical), chs. 2 and 7).
**(a)** *Replicate.* For integer $Q$, write $\E Z^{Q}$ as in [(A.2)](#ezn): a $Q$-fold copy of the system, coupled only through the shared data.
**(b)** *Average over data.* Show that for fixed replicas and teacher, $u^{a} = \theta^{a}\cdot x/\sqrt{N}$ and $v = \theta^{*}\cdot x/\sqrt{N}$ are exactly jointly Gaussian with the covariance [(A.3)](#overlap-matrix), and that the data average factorizes across examples, each contributing the same function $e^{-e(q,R)}$ of the overlaps. Check that at $Q=1$, $e(R) = -\log(1-\varepsilon(R))$ with $\varepsilon$ the wedge formula [(2.41)](#eps-of-r).
**(c)** *Push forward to the order parameters.* Insert delta functions fixing $q^{ab}$ and $R^{a}$ (in integral representation, with conjugate variables $\hat q^{ab}$, $\hat R^{a}$), so that the sum over the binary vectors factorizes over sites $i=1,\ldots,N$ and becomes the $N$-th power of a $Q$-spin sum. Conclude that $\E Z^{Q}$ takes the form [(A.4)](#replica-action), and that the annealed action [(2.44)](#perc-action) is its $Q=1$ case.
**(d)** *Replica symmetry.* Take $q^{ab}= q$ for $a\neq b$ and $R^{a} = R$ (and likewise for the conjugates). Show that the orthant probability becomes $e^{-e_{\rm RS}}= 2\int Dt\,H(\cdot)^{Q}$ for an appropriate argument of the tail function $H$, built from $q$, $R$, and a single Gaussian variable $t$, and that the site sum becomes a one-dimensional Gaussian integral of $(2\cosh(\cdot))^{Q}$. Continue both in $Q$ and expand to first order at $Q=0$.
**(e)** *Use the symmetry of Gibbs learning.* A student drawn uniformly from the consistent set is statistically indistinguishable from the teacher, since the teacher is itself a uniformly random binary vector consistent with the data. Argue that this forces $q = R$ at the saddle: the overlap between two Gibbs students equals their overlap with the teacher. One scalar remains, with one conjugate.
**(f)** *Solve.* The saddle-point equations for $R(\alpha)$ are transcendental. Solve them numerically and show that there are two competing solutions, an interior one and $R=1$, whose free energies cross at $\alpha_{c}\approx1.245$. Compare with the annealed $\alpha_{c}^{\rm ann}\approx1.45$ of [Exercise 2.3(g)](#ex-ising-g), and explain the direction of the shift using [Appendix A.1](#a1-what-the-annealed-approximation-gets-wrong).
### A.3 A variant: the spherical perceptron
The replica method is easiest to push all the way to closed form in a variant of the problem with *continuous* weights, which is also the classic case in which the annealed and quenched learning curves can be compared analytically. Keep everything in [Section 2.6](#26-the-ising-perceptron-a-first-order-transition-in-learning) the same, but let the weights be real, $\theta,\theta^{*}\in\R^{N}$. Since $f_{\theta} = \mathrm{sign}(\theta\cdot x)$ is invariant under $\theta\to\lambda\theta$ with $\lambda>0$, remove the redundant scale by putting the weights on a sphere,
$$
|\theta|^{2} = |\theta^{*}|^{2} = N\,,
$$
normalized so that individual weights stay of order one, and let the prior be the uniform measure $\pi(\theta)\,d\theta$ on the sphere. This is the *spherical perceptron*. The wedge formula [(2.41)](#eps-of-r) for the generalization error is unchanged, since its derivation used only that $(\theta\cdot x,\theta^{*}\cdot x)/\sqrt{N}$ is a Gaussian pair with correlation $R$, which holds for any $\theta$ with $|\theta|^{2} = N$. What changes is the entropy: the number of students at overlap $R$ becomes the *volume* of the slice of the sphere at overlap $R$. Decompose $\theta = R\,\theta^{*} + \theta_{\perp}$ with $\theta_{\perp}\perp\theta^{*}$. The constraint [(A.6)](#sph-constraint) gives $|\theta_{\perp}|^{2} = N(1-R^{2})$, so the slice is a sphere of radius $\sqrt{N(1-R^{2})}$ in the $(N-1)$-dimensional orthogonal complement, of volume $\propto\big(N(1-R^{2})\big)^{(N-2)/2}$, and
[Figure 9](#fig-gann) shows $S_{\rm ann}$ for several loads. The entropy–energy competition is the same as for the Ising perceptron, but the outcome is different: the entropy $-\frac{1}{2}\log(1-R^{2})$ diverges as $R\to1$, so the perfectly aligned state is infinitely penalized, the minimum always sits in the interior, and it slides continuously from $R=0$ toward $R=1$ as $\alpha$ grows. There is no phase transition. This is the answer to the closing question of [Exercise 2.3(h)](#ex-ising-h): the binary weights of the Ising perceptron give the aligned state finite entropy, $h(1)=0$, and that is what allows it to compete as a separate basin.
The annealed action $S_{\rm ann}(R)$ of the spherical perceptron, [(A.8)](#annealed-integral), for $\alpha\in\{0,0.25,0.5,0.75,1\}$ (bottom to top at $R=0$), with the minimizer marked on each curve. The entropy pins the minimum at $R=0$ when $\alpha=0$; increasing the load slides it smoothly toward larger overlap, and the entropic wall at $R=1$ keeps it in the interior. One basin at every $\alpha$, moving continuously: no phase transition.
*The learning curve.* The large-$\alpha$ behavior follows by expanding near $R=1$. Writing $R = \cos(\pi\varepsilon)$, so that $\varepsilon$ is the generalization error itself, $1-R^{2}\approx(\pi\varepsilon)^{2}$ and
The generalization error decays as a *power law* in the amount of data, with exponent one: the simplest instance of the scaling laws of [Section 7](#7-scaling-laws). The exact, quenched treatment ([Seung et al. 1992](#bib-seung1992statistical); [Engel & Van den Broeck 2001](#bib-engel2001statistical)) runs through the replica steps of [Exercise A.1](#ex-replica) with the volume entropy [(A.7)](#coarea) in place of the count (here replica symmetry can be justified by the convexity of the region of the sphere consistent with the data), and yields
the same exponent, with a different constant. As in [Appendix A.1](#a1-what-the-annealed-approximation-gets-wrong), the annealed average overweights poorly aligned students and overestimates the error, $1/\alpha>0.625/\alpha$; and since the model is realizable, the exponent survives the approximation even though the constant does not.
**Exercise A.2 (The annealed learning curve of the spherical perceptron).** **(a)** Derive the slice entropy [(A.7)](#coarea), and assemble the annealed action [(A.8)](#annealed-integral) by repeating the steps of [Exercise 2.3(d)](#ex-ising-pushforward).
**(b)** Derive the saddle-point condition for [(A.9)](#annealed-saddle), and show that $R=1$ is never a minimum, in contrast with [Exercise 2.3(f)](#ex-ising-f).
**(c)** Verify the large-$\alpha$ learning curve [(A.11)](#annealed-eps), and show that at small $\alpha$ the annealed theory predicts $R\propto\alpha$. Sketch the full learning curve $\varepsilon(\alpha)$.
**(d)** Which steps of [Exercise A.1](#ex-replica) change when the count of binary students is replaced by the volume [(A.7)](#coarea), and which do not?
## B. From linear regression to kernel regression
This appendix is a self-contained review of models that are *linear in their parameters*: linear regression, regression on fixed nonlinear features, kernel regression, and kernel *ridge* regression. It ends with the observation that a neural network linearized around its initialization — the ansatz [(4.1)](#f-lin) of [Section 4](#4-the-ntk-and-large-width-during-training) — is exactly such a model, so that training it amounts to kernel ridge regression. Nothing here requires large width, gradient-flow asymptotics, or moving features; those are the subjects of [Sections 3–5](#3-large-width-at-initialization-the-gaussian-process-limit), and this appendix supplies the linear algebra they rest on. Readers who found the kernel formulas of [Sections 4.1](#41-gradient-flow-and-the-ntk) and [4.6](#46-bayesian-learning-at-large-width-kernel-learning) appearing out of thin air may want to read it first. [^6]
We use the notation of [Section 3.1](#31-notation): training data $(x^{\mu},y^{\mu})_{\mu=1}^{n}$ with $x^{\mu}\in \R^{d}$ and $y^{\mu}\in\R$; a model $f_{\theta}:\R^{d}\to \R$ with parameters $\theta\in \R^{P}$; and $f^{\mu} := f_{\theta}(x^{\mu})$ for its values on the training inputs. It is convenient to stack the data into matrices,
$$
X := \begin{pmatrix}(x^{1})^{T} \\ \vdots \\ (x^{n})^{T}\end{pmatrix} \in \R^{n\times d}\,,\qquad Y := \begin{pmatrix}y^{1} \\ \vdots \\ y^{n}\end{pmatrix} \in \R^{n}\,,\qquad f := \begin{pmatrix}f^{1} \\ \vdots \\ f^{n}\end{pmatrix} \in \R^{n}\,,
$$
so that $X^{\mu}{}_{j} = x^{\mu}_{j}$ is the $j$-th coordinate of the $\mu$-th input, and the inputs and their labels are stacked in the same way. Throughout, the loss is the MSE of [Section 4.1.1](#411-geometry-of-gradient-flows),
with $\Delta\in\R^{n}$ the residual. Many references (including the notebook this appendix is drawn from) normalize the loss as $\frac{1}{2n}|f-Y|^{2}$ instead; this rescales the ridge parameter $\lambda$ below by a factor of $n$ and the learning rate $\eta$ by a factor of $1/n$, and changes nothing else.
### B.1 Linear regression
The simplest model is linear in the inputs,
so that $P=d$ and the vector of training predictions is $f = X\theta$. (An intercept can be included by appending a constant coordinate $x_{d+1}=1$ to every input, as in [Exercise B.1](#ex-app-lsq) below.) The loss $\CL_{\theta} = \frac{1}{2}|X\theta - Y|^{2}$ is a convex quadratic function of $\theta$, with gradient
This has the shape that every gradient in these notes has: the transpose of the Jacobian $\pd f^{\mu}/\pd\theta^{j}$ — here the constant matrix $X$ — acting on the residual, *cf.* [(4.10)](#ntk). Setting the gradient to zero gives the *normal equations*
$$
X^{T}X\,\theta = X^{T}Y\,.
$$
Since the loss is convex, every solution of [(B.5)](#app-normal) is a global minimizer, and there are no other critical points. What the solutions look like depends on how $P=d$ compares with $n$.
**Fewer parameters than data, $P0$ for $v\neq 0$. There is a unique minimizer,
$$
\theta_{*} = (X^{T}X)^{-1}X^{T}Y\,.
$$
The normal equations say that the residual $\Delta_{*} = X\theta_{*}-Y$ is orthogonal to every column of $X$. Equivalently, the prediction $X\theta_{*}$ is the *orthogonal projection* of $Y$ onto the column space $\mathrm{col}(X)\subset\R^{n}$, the set of all prediction vectors the model can produce ([Figure 10](#fig-lsq)). Fitting is projecting. The projection $X\theta_{*}$ is unique even in cases where $\theta_{*}$ is not.
Least squares, geometrically. The normal equations $X^{T}\Delta_{*}=0$ say that the residual is orthogonal to the column space of $X$, so the fitted prediction vector $X\theta_{*}$ is the orthogonal projection of $Y$ onto $\mathrm{col}(X)$. "Orthogonal" refers to the Euclidean metric $\delta_{\mu\nu}$ on the space $\R^{n}$ of function values on the training set.
It is worth pausing on the word "orthogonal." Least squares minimizes the Euclidean ($L^{2}$) distance between $X\theta$ and $Y$, and the projection is orthogonal with respect to the Euclidean metric $\delta_{\mu\nu}$ on $\R^{n}$ — the same metric on the space of function values that was made explicit in [Section 4.1.1](#411-geometry-of-gradient-flows). This metric is a *choice*, hidden inside the phrase "mean squared error." It declares that an error of a given size counts the same at every data point, and that errors at different data points are independent. Any other positive-definite metric $g_{\mu\nu}$ on $\R^{n}$ defines a different loss $\CL^{g} = \frac{1}{2}g_{\mu\nu}\Delta^{\mu}\Delta^{\nu}$, with a different minimizer
$$
\theta_{*}^{g} = (X^{T}gX)^{-1}X^{T}gY\,,
$$
which is the $g$-orthogonal projection of $Y$ onto $\mathrm{col}(X)$ ("weighted least squares"). Unless $Y$ already lies in $\mathrm{col}(X)$, different metrics give genuinely different fitted functions; see [Exercise B.1(d)](#ex-app-lsq-weighted). Note that $g_{\mu\nu}$ is a metric on the space of function values *indexed by* the data points, not a metric on the input space $\R^{d}$.
**More parameters than data, $P>n$.** If $d>n$, or more generally if the columns of $X$ are linearly dependent, then $X^{T}X$ is singular. The normal equations still have solutions ([Exercise B.1(b)](#ex-app-lsq-posdef)), but now infinitely many: an affine subspace of $\R^{P}$ of dimension $P-\mathrm{rank}(X)$. If moreover $X$ has full row rank $n$, then $Y\in\mathrm{col}(X)$ and *every* solution interpolates the data exactly, $X\theta = Y$, with zero loss. The loss has no way to prefer one of them. Something else must make the choice; that is the subject of [Appendix B.3](#b3-kernel-regression), and it is where kernels come from. Far from being an edge case, $P>n$ is the normal situation for neural networks.
**Exercise B.1 (Least squares by hand).** **(a)** Differentiate $\CL_{\theta} = \frac{1}{2}|X\theta-Y|^{2}$ with respect to a single coordinate $\theta^{j}$, assemble the components into [(B.4)](#app-lsq-grad), and derive the normal equations [(B.5)](#app-normal). Explain geometrically why they say that the residual is orthogonal to every column of $X$.
**(b)** Show that if $X$ has full column rank then $X^{T}X$ is positive definite, so that [(B.6)](#app-lsq) is the unique minimizer. Show that the normal equations are always consistent (*i.e.* $X^{T}Y\in\mathrm{col}(X^{T}X)$), and that if $P>n$ they have infinitely many solutions.
**(c)** Fit the affine function $f(x) = \theta^{1} + \theta^{2} x$ to the three points $x\in\{-1,0,1\}$ with labels $y=x^{2}$; that is, take
$$
X = \begin{pmatrix}1 & -1 \\ 1 & 0 \\ 1 & 1\end{pmatrix}\,,\qquad Y = \begin{pmatrix}1 \\ 0 \\ 1\end{pmatrix}\,.
$$
Find $\theta_{*}$, the predictions, the residual, and the loss. Is the remaining error a failure of optimization or of representation?
**(d)** Repeat (c) with the weighted loss [(B.7)](#app-wls) for $g = \mathrm{diag}(1,4,1)$, which trusts the middle data point four times as much as the others. Check that the fitted function changes.
### B.2 Linear regression with nonlinear features
"Linear model" is a statement about how the parameters enter, not about the shape of the fitted function. Choose, *before looking at any labels*, a *feature map*
$$
\psi:\R^{d}\to \R^{P}\,,\qquad \psi(x) = \big(\psi_{1}(x),\ldots,\psi_{P}(x)\big)\,,
$$
and consider the model
This is a weighted sum of $P$ fixed shapes $\psi_{i}$. It can be as nonlinear in $x$ as the shapes are, but it is linear in $\theta$. Standard choices include monomials in the coordinates of $x$ up to some degree $m$ (giving $P=\binom{d+m}{m}$), Fourier modes, Gaussian bumps $e^{-|x-c_k|^2/2\sigma^2}$ centered at chosen points $c_{k}$, and *random features* $\psi_{i}(x) = \phi(a_{i}\cdot x)$ with the $a_{i}$ drawn at random and then frozen — the last being a one-hidden-layer network whose first layer is never trained.
Stacking the feature vectors of the training inputs into the $n\times P$ *feature matrix*
$$
\Psi^{\mu}{}_{i} := \psi_{i}(x^{\mu})\,,
$$
the training predictions are $f = \Psi\theta$. Note that $\Psi^{\mu}{}_{i} = \pd f^{\mu}/\pd\theta^{i}$: the feature matrix is precisely the Jacobian $J$ of [Section 4.2](#42-features), which for a linear model is a constant. Linear regression is the special case $\psi(x) = x$, $\Psi = X$, and *everything* in [Appendix B.1](#b1-linear-regression) goes through verbatim with $X$ replaced by $\Psi$: the normal equations are $\Psi^{T}\Psi\,\theta = \Psi^{T}Y$; if $P\leq n$ and $\Psi$ has full column rank, the unique minimizer is $\theta_{*} = (\Psi^{T}\Psi)^{-1}\Psi^{T}Y$; and if $P>n$ there are infinitely many interpolants.
Choosing $\psi$ is choosing the model. The same seven data points (black), interpolated exactly by three different fixed feature maps with $P=n=7$: polynomials of degree six, trigonometric functions of period one with three harmonics, and Gaussian bumps centered at the data. All three drive the training loss to zero; they disagree everywhere outside the training range (shaded), and that disagreement was decided before the labels were read.
Replacing $X$ by $\Psi$ changes the algebra not at all, and changes the model completely. [Figure 11](#fig-feature-maps) shows three feature maps, each with $P=n$, fitted to the same seven points. Each interpolates the data exactly, and they disagree everywhere else. Which is right? The data cannot say: the behavior away from the training set was fixed by the choice of $\psi$, which was made before the labels were seen. The feature map is the model's *inductive bias*, made explicit.
One fixed extra feature can make a hard problem linear. Left: two classes of points in the plane, on concentric rings. No line separates them, since the convex hull of each class contains the origin. Right: after appending the single feature $\psi_{3}= c\,(u^{2}+v^{2})$, the rings sit at two different heights, and a horizontal plane $\psi_{3}= c\,r_{0}^{2}$ separates them. Nothing has been learned: we supplied the radial coordinate.
Why would one ever want a very large $P$? [Figure 12](#fig-radial-lift) gives the classic answer. Two classes of points on concentric rings cannot be separated by a line in the plane. Append one fixed feature, $\psi(u,v) = (u,\,v,\,c\,(u^{2}+v^{2}))$, and the classes become separable by a plane in $\R^{3}$: a model that is linear in the lifted coordinates is a circle in the original ones. The price is that $P$ grows fast. All monomials of degree $3$ in $d=1000$ input coordinates already give $P = \binom{1002}{3}\approx 1.7\times 10^{8}$, and Gaussian bumps centered at every point of a continuum give $P=\infty$. Yet there are only $n$ data points, and the fit is determined by $n$ numbers. Nothing about the problem should require a $P\times P$ inverse — so where does $P$ actually enter? Answering that question honestly is the content of the next two sections.
**Exercise B.2 (A three-point feature calculation).** Take the same three points $x\in\{-1,0,1\}$ with labels $y=x^{2}$ as in [Exercise B.1(c)](#ex-app-lsq-affine), but now the feature map $\psi(x) = (1,\,\sqrt{2}\,x,\,x^{2})$.
**(a)** Form the feature matrix $\Psi$ and find coefficients $\theta$ that interpolate all three labels. The representational failure of [Exercise B.1(c)](#ex-app-lsq-affine) is gone: what changed?
**(b)** Explain precisely how the predictor can be nonlinear in $x$ while the regression problem remains linear.
**(c)** If the last feature is replaced by $c\,x^{2}$ for some $c\neq 0$, what changes: the fitted function, its coefficient, or both?
**(d)** Write the inner product $K(x,z) := \psi(x)\cdot\psi(z)$ in closed form, without vector notation. This is a preview of [Appendix B.3](#b3-kernel-regression).
### B.3 Kernel regression
Return to the overparametrized case $P>n$, and assume that $\Psi$ has full row rank $n$ (the generic situation), so that there is a $(P-n)$-dimensional affine family of exact interpolants $\Psi\theta = Y$. To choose one, we need a criterion that the loss does not provide. The standard choice is the interpolant of *minimum norm*,
This introduces a *second* implicit Euclidean metric — now $\delta_{ij}$ on parameter space $\R^{P}$, the metric that also defines gradient descent in [Section 4.1.1](#411-geometry-of-gradient-flows). ([Appendix B.1](#b1-linear-regression) used $\delta_{\mu\nu}$ on the space $\R^{n}$ of function values.) The two metrics play dual roles: the first decides how to fit when the data cannot be matched, the second decides which fit to choose when it can.
To solve [(B.11)](#app-minnorm), introduce Lagrange multipliers $\alpha\in\R^{n}$ for the $n$ constraints and extremize $\frac{1}{2}|\theta|^{2} - \alpha\cdot(\Psi\theta - Y)$. Stationarity in $\theta$ gives
This is the key structural fact: however large $P$ is, the minimum-norm solution lies in the $n$-dimensional subspace of $\R^{P}$ spanned by the feature vectors of the training inputs, and is described by $n$ numbers. (In the language of [Section 4.2](#42-features), $\theta_{*}$ lies in the row space of the Jacobian.) The constraint $\Psi\theta_{*} = Y$ then becomes $\Psi\Psi^{T}\alpha = Y$, an $n\times n$ linear system. Define the *kernel* and its *Gram matrix* on the training set,
The kernel is an inner product of feature vectors: a notion of *similarity* between inputs, fixed once $\psi$ is fixed, and involving no labels. The Gram matrix is symmetric and positive semidefinite, since $c^{T}\hat K c = |\Psi^{T}c|^{2}\geq 0$; it is positive definite exactly when $\Psi$ has full row rank, which we assumed. Then $\alpha = \hat K^{-1}Y$, and the prediction at an arbitrary input $x$ is
This is *kernel regression* (more precisely, kernel interpolation). The feature map has disappeared from both the fit and the prediction. To fit, one needs $K$ on pairs of training points; to predict, one needs $K$ between the test point and the training points; no feature vector is ever constructed. The prediction is a weighted sum of similarities to the training points, with weights $\alpha$ that depend on the labels — but the labels only reweight the training points, they never change which inputs $K$ considers similar.
Compare [(B.14)](#app-kernel-reg) with the late-time mean [(4.38)](#ntk-gp) of a lazily trained network: they are the same formula, with $\hat K$ in place of the NTK. This is no coincidence. For any model linear in its parameters the Jacobian is $J=\Psi$, so the NTK [(4.10)](#ntk) is
constant in time. Moreover, gradient flow finds the minimum-norm solution: $\dot\theta = -\eta\,\Psi^{T}\Delta$ always points along the row space of $\Psi$, so $\theta(t)-\theta_{0}$ stays in that $n$-dimensional subspace and converges to the interpolant nearest to $\theta_{0}$ — for $\theta_{0} = 0$, to $\theta_{*}$ ([Exercise B.3(c)](#ex-app-kernel-flow)). This is the simplest instance of the *implicit bias* of gradient descent, and it is the Euclidean metric on parameter space entering through the back door: the flow is defined by $\delta^{ij}$, and so is the norm that it ends up minimizing.
**The kernel trick.** There are two readings of [(B.14)](#app-kernel-reg). Read from the feature side, it rewrites a $P\times P$ problem as an $n\times n$ one. Read from the kernel side, it says the feature map is dispensable: any function $K(x,x')$ that arises as an inner product of features will do, whether or not one can write the features down. For $x,z\in\R^{2}$ and $\psi(x) = (x_{1}^{2},\,\sqrt{2}\,x_{1}x_{2},\,x_{2}^{2})$, one checks that $\psi(x)\cdot\psi(z) = (x\cdot z)^{2}$: a single dot product, squared, stands in for three features. In $d$ dimensions, $K(x,z) = (x\cdot z)^{m}$ costs $O(d)$ to evaluate and stands in for all $\binom{d+m-1}{m}$ monomials of degree $m$. The Gaussian kernel $K(x,z) = e^{-|x-z|^2/2\sigma^2}$ corresponds to an infinite-dimensional feature map that no one ever writes down. Conversely, any symmetric $K$ all of whose Gram matrices are positive semidefinite arises from some feature map (Mercer's theorem): the eigenfunctions of $K$, weighted by the square roots of its eigenvalues, furnish one (*cf.* [Sections 3.3](#33-gaussian-processes) and [4.6](#46-bayesian-learning-at-large-width-kernel-learning)). This is a trade, not a free lunch. The feature route scales with $P$ and the kernel route with $n$; with many examples and few features, the explicit route is the cheap one.
**Exercise B.3 (Minimum-norm interpolation and the kernel).** **(a)** Carry out the Lagrange-multiplier computation leading to [(B.12)](#app-representer)–[(B.14)](#app-kernel-reg). Alternatively, use the singular value decomposition $\Psi = USV^{T}$ of [Section 4.2](#42-features) to show that the minimum-norm interpolant is $\theta_{*} = \Psi^{T}(\Psi\Psi^{T})^{-1}Y$, and write it in terms of $U,S,V$.
**(b)** Show that $c^{T}\hat Kc = |\Psi^{T}c|^{2}\geq 0$ for all $c\in\R^{n}$. Interpret a zero eigenvalue of $\hat K$ as a linear combination of training examples that vanishes in feature space. What does it take for $\hat K$ to be invertible?
**(c)** Show that under gradient flow $\dot\theta = -\eta\,\Psi^{T}(\Psi\theta - Y)$ the displacement $\theta(t)-\theta_{0}$ stays in the row space of $\Psi$ for all $t$; derive $\Delta(t) = e^{-\eta t\hat K}\Delta(0)$; and conclude that, for invertible $\hat K$, the flow converges to the interpolant closest to $\theta_{0}$. Derive [(4.36)](#ntk-learn) for this model, with $\Theta = \hat K$.
**(d)** Verify that $\psi(x)\cdot\psi(z) = (x\cdot z)^{2}$ for $\psi(x) = (x_{1}^{2},\sqrt{2}\,x_{1}x_{2},x_{2}^{2})$, and explain why the $\sqrt{2}$ is not decoration. Which feature map on $\R$ produces the kernel $(1+xz)^{2}$ of [Exercise B.2(d)](#ex-app-features-kernel)? How many features does $(x\cdot z)^{m}$ on $\R^{d}$ stand in for?
### B.4 Kernel ridge regression
Rather than imposing the minimum-norm criterion as a hard constraint, one can add it to the loss as a penalty. For $\lambda>0$, *ridge regression* minimizes
The two terms measure distance to the data with $\delta_{\mu\nu}$ and distance to the origin with $\delta_{ij}$: ridge regression is a weighted compromise between the two metrics of [Appendices B.1](#b1-linear-regression) and [B.3](#b3-kernel-regression), with $\lambda$ setting the exchange rate. (In the training of neural networks the penalty is called *weight decay*.) The penalty does two jobs at once. When $P>n$ it selects one solution out of infinitely many; and for any $P$ it makes the problem well posed even when the features are redundant or the labels are noisy. Note that the penalty is not invariant under rescaling features: it decides what counts as a "large" coefficient in the coordinates $\psi$ that were supplied.
The gradient is $\Psi^{T}(\Psi\theta - Y) + \lambda\theta$, so the stationarity condition is
The matrix on the left is positive definite for every $\lambda>0$, whatever the rank of $\Psi$, since $v^{T}(\Psi^{T}\Psi+\lambda I_{P})v = |\Psi v|^{2}+\lambda|v|^{2}>0$. So there is always a unique minimizer,
at the price of a $P\times P$ inverse. To move the inverse to the $n\times n$ side, expand and regroup the product
$$
\Psi^{T}(\Psi\Psi^{T} + \lambda I_{n}) = \Psi^{T}\Psi\Psi^{T} + \lambda\Psi^{T} = (\Psi^{T}\Psi + \lambda I_{P})\Psi^{T}\,.
$$
Both bracketed matrices are positive definite, so multiplying by their inverses on the appropriate sides gives the *push-through identity*
Nothing was cancelled — $\Psi^{T}$ is rectangular and has no inverse — and that is precisely why the identity holds. Substituting into [(B.18)](#app-ridge-p), the ridge solution again lies in the span of the training feature vectors, $\theta_{\lambda} = \Psi^{T}\alpha_{\lambda}$, now with $\alpha_{\lambda} = (\hat K+\lambda I_{n})^{-1}Y$, and the prediction at a new input is
$$
f_{\lambda}(x) = \sum_{\mu,\nu=1}^{n} K(x,x^{\mu})\,\big(\hat K + \lambda I_{n}\big)^{-1}_{\mu\nu}\,y^{\nu}\,.
$$
This is *kernel ridge regression*, the central formula of this appendix. As $\lambda\to 0$ it reduces to the minimum-norm interpolant [(B.14)](#app-kernel-reg) when $\hat K$ is invertible, and in general to the least-squares solution of minimum norm ([Exercise B.4(c)](#ex-app-ridge-limit)). As $\lambda\to\infty$, $\theta_{\lambda}\to 0$.
Formula [(B.21)](#app-krr) is exactly the posterior mean of [Exercise 4.5](#ex-kernel-ridge), with ridge parameter $\lambda = \beta^{-1}$ — there derived for the infinite-width Gaussian process, here for a finite linear model. The Bayesian reading is worth making explicit. With a Gaussian likelihood $\propto e^{-\frac{\beta}{2}|\Psi\theta-Y|^2}$ and a Gaussian prior $\theta\sim\CN(0,\sigma^{2} I_{P})$, the loss [(B.16)](#app-ridge-loss) is $-\beta^{-1}\log$ of the posterior density, up to a constant, with $\lambda = 1/(\beta\sigma^{2})$. Ridge regression is thus maximum a posteriori estimation. The prior pushed forward to $f_{\theta}(x) = \theta\cdot\psi(x)$ is a Gaussian process with kernel $\sigma^{2}K(x,x')$, a finite-$P$ instance of [Section 3.3](#33-gaussian-processes), and the posterior mean of [Exercise 4.5](#ex-kernel-ridge) for this process reproduces [(B.21)](#app-krr).
**Spectral view.** With the singular value decomposition $\Psi = USV^{T}$ of [Section 4.2](#42-features), the Gram matrix is $\hat K = US^{2}U^{T}$, with eigenvalues $s_{a}^{2}$ and orthonormal eigenvectors $u_{a}\in\R^{n}$: the "features" of [Section 4.2](#42-features), *i.e.* the natural modes in the space of function values. The ridge predictions on the training set are $f_{\lambda} = \Psi\theta_{\lambda} = \hat K(\hat K+\lambda I_{n})^{-1}Y$, or, mode by mode,
Ridge regression is a *spectral filter*: modes of the data with $s_{a}^{2}\gg\lambda$ are fit, modes with $s_{a}^{2}\ll\lambda$ are suppressed, and the crossover is at $s_{a}^{2}\approx\lambda$. Compare gradient flow from $\theta_{0} = 0$ on the unpenalized loss, whose training predictions are, by [(4.35)](#d-converge),
$$
f(t) = \sum_{a} \big(1-e^{-\eta t s_a^2}\big)\,(u_{a}\cdot Y)\,u_{a}\,.
$$
This is also a spectral filter, with $1/\eta t$ playing the role of $\lambda$: *early stopping* regularizes much as ridge does. But the two filters have different shapes ([Figure 13](#fig-filters)), and no choice of $\lambda$ reproduces a given stopping time on every mode at once ([Exercise B.5](#ex-app-filters)). This is the precise content of the remark closing [Section 4.6](#46-bayesian-learning-at-large-width-kernel-learning), that the inverse temperature of Bayesian learning and the training time of gradient descent play analogous roles.
Ridge regression and early-stopped gradient flow as spectral filters: the fraction of a data mode $u_{a}\cdot Y$ that is fit, as a function of the mode's kernel eigenvalue $s_{a}^{2}$. The two are tuned to agree on the mode with $s^{2}=1$, and disagree elsewhere; the dots mark the two-mode example of [Exercise B.5](#ex-app-filters).
**Exercise B.4 (Ridge regression, two ways).** **(a)** Derive [(B.17)](#app-ridge-normal) and [(B.18)](#app-ridge-p), and prove that $\Psi^{T}\Psi+\lambda I_{P}$ is invertible even when $\Psi$ is rank deficient.
**(b)** Prove the push-through identity [(B.20)](#app-pushthrough) and deduce [(B.21)](#app-krr). Explain why "cancelling $\Psi^{T}$" is not a valid step, and why the identity is nevertheless true. (It is a special case of $A(BA+I)^{-1}= (AB+I)^{-1}A$.)
**(c)** What happens as $\lambda\to\infty$? Using the SVD, show that $\theta_{\lambda}$ applies the scalar filter $s/(s^{2}+\lambda)$ to each singular value $s$ of $\Psi$, and hence that as $\lambda\to 0$ it converges to the least-squares solution of minimum norm. When does this limit interpolate the labels exactly?
**(d)** Return to $x\in\{-1,0,1\}$, $y=x^{2}$ and the kernel $K(x,z) = (1+xz)^{2}$ of [Exercise B.2](#ex-app-features). With $\lambda = 1$, verify that
$$
\hat K = \begin{pmatrix}4 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 4\end{pmatrix}\,,\qquad \alpha_{\lambda} = \frac{1}{4}\begin{pmatrix}1 \\ -1 \\ 1\end{pmatrix}\,,
$$
and predict at $z=2$ using only kernel evaluations. Compare with the interpolant of [Exercise B.2(a)](#ex-app-features-interp), which predicts $z^{2}=4$ there, and explain the direction of the discrepancy.
**(e)** The feature route solves a $P\times P$ system and the kernel route an $n\times n$ one. When is each attractive?
**Exercise B.5 (Early stopping vs. ridge, on two data points).** Take $n=2$ training points with Gram matrix and labels
$$
\hat K = \begin{pmatrix}2 & 1 \\ 1 & 2\end{pmatrix}\,,\qquad Y = \begin{pmatrix}2 \\ 0\end{pmatrix}\,,
$$
and initial predictions $f(0) = 0$.
**(a)** Verify the eigenpairs $u_{\pm} = \frac{1}{\sqrt{2}}(1,\pm 1)^{T}$ with $s_{\pm}^{2} = 3, 1$, and decompose $Y$ into these modes.
**(b)** Using [(B.23)](#app-flow-filter), find the training predictions $f(t)$ under gradient flow, and evaluate them at $\eta t = \log 2$. Why do both the eigenvalue $s_{a}^{2}$ and the initial projection $u_{a}\cdot Y$ matter for what is visible during training?
**(c)** Using [(B.22)](#app-ridge-filter), find the ridge predictions for $\lambda = 1$. The two filters agree on the mode with $s^{2} = 1$; show that the ridge parameter matching gradient flow at time $t$ on a mode with eigenvalue $s^{2}$ is $\lambda_{\rm eff}(s^{2},t) = s^{2}/(e^{\eta t s^2}-1)$, which depends on the mode — so early stopping is not ridge regression.
**(d)** State precisely what "spectral bias" means for this finite, fixed-kernel calculation. What additional structure would be needed to turn it into a claim about Fourier frequencies, about robustness to label noise, or about test error?
### B.5 Linearized networks are kernel ridge regressors
Now let $f_{\theta}$ be a neural network, with $\theta\in\R^{P}$ the vector of *all* its weights and biases. It is not linear in $\theta$. But its linearization around the initialization $\theta_{0}$, the ansatz [(4.1)](#f-lin),
is a model of the form [(B.9)](#app-featmodel). Compare term by term. The parameters are the displacement $\delta\theta := \theta-\theta_{0}\in\R^{P}$. The feature map is
the *tangent features*: one number per parameter, the sensitivity of the output at $x$ to that parameter, evaluated at initialization and then frozen. And the constant offset $f_{\theta_0}(x)$ is absorbed by shifting the labels to
The offset matters: the linearized model is fit to the *residuals at initialization*, not to the raw labels. The feature matrix [(B.10)](#app-psi) is the Jacobian at initialization, $\Psi = J_{0}$, and the Gram matrix [(B.13)](#app-kernel) is
the neural tangent kernel [(4.10)](#ntk) at initialization, extended to arbitrary pairs of inputs as in [Exercise 4.3](#ex-ntk-test). Everything in [Appendices B.3](#b3-kernel-regression) and [B.4](#b4-kernel-ridge-regression) now applies verbatim. In particular, ridge regression on the linearized network with penalty $\frac{\lambda}{2}|\theta-\theta_{0}|^{2}$ — weight decay toward the initialization — gives, by [(B.21)](#app-krr),
At $\lambda\to 0$ this is the $t\to\infty$ limit of [(4.36)](#ntk-learn), the minimum-norm interpolant of the residuals; and the finite-time solution [(4.36)](#ntk-learn) is its early-stopped cousin, in the sense of [(B.23)](#app-flow-filter). The kernel trick is not a convenience here but a necessity: $P$ is the parameter count of the network, often $10^{7}$ or more, while $\Theta_{0}$ is $n\times n$.
Two remarks on what is and is not being claimed. First, for the linearized model $f^{\rm lin}$ everything above is exact, and its kernel is constant by construction. For the actual network, [(B.24)](#app-flin) is an ansatz. [Section 4.3](#43-constant-ntk-at-large-width) shows that it becomes exact in the large-width limit at fixed depth and data; whether and when it is a good approximation otherwise — and what happens when the tangent features move — is the subject of [Sections 5–6](#5-mean-field-scaling), not of this appendix. Second, the initialization $\theta_{0}$ is random, so $f_{\theta_0}$, the tangent features, and the kernel $\Theta_{0}$ are all random too, and [(B.28)](#app-lin-krr) is the trained function for one draw. Averaging over draws, at large width where $\Theta_{0}$ becomes deterministic, gives the Gaussian process [(4.38)](#ntk-gp).
**Exercise B.6 (Tangent features of a one-hidden-layer network).** Consider $f_{\theta}(x) = \sum_{i=1}^{N} c_{i}\,\phi(a_{i}\cdot x)$ with $\theta = (c_{i}, a_{i})_{i=1}^{N}$, so that $P = N(d+1)$.
**(a)** Compute the tangent features $\nabla_{\theta} f_{\theta_0}(x)$, and show that the tangent kernel is
$$
\Theta(x,x') = \sum_{i=1}^{N} \phi(a_{i}\cdot x)\,\phi(a_{i}\cdot x') + \sum_{i=1}^{N} c_{i}^{2}\,\phi'(a_{i}\cdot x)\,\phi'(a_{i}\cdot x')\,(x\cdot x')\,.
$$
**(b)** Suppose only the output weights $c$ are trained, with the $a_{i}$ frozen at their random initial values. Show that this is exactly a fixed-feature model [(B.9)](#app-featmodel) with random features, that its linearization is exact, and that its kernel is the first term above. Which term do you get by training only the $a_{i}$?
**(c)** Show that gradient flow with weight decay toward the initialization, $\dot\theta = -\eta\big[\nabla_{\theta} \CL + \lambda(\theta-\theta_{0})\big]$, applied to the linearized model $f^{\rm lin}$, converges to the ridge solution [(B.28)](#app-lin-krr) for every $\lambda>0$ and from any starting point.
**(d)** Check that [(B.28)](#app-lin-krr) at $\lambda\to 0$ coincides with [(4.36)](#ntk-learn) at $t\to\infty$.
## C. Solutions to exercises
Worked solutions to the exercises, in the order in which the exercises appear. (Instructors: uncommenting `\solutionsfalse` in the preamble removes this appendix's contents, and all other solution boxes, from the compiled PDF.)
## References
Elisabeth Agoritsas, Giulio Biroli, Pierfrancesco Urbani, and Francesco Zamponi (2018). *Out-of-equilibrium dynamical mean-field equations for the perceptron model*. Journal of Physics A: Mathematical and Theoretical.
Alexander Altland and Ben D. Simons (2010). *Condensed Matter Field Theory*. Cambridge University Press.
Daniel J. Amit, Hanoch Gutfreund, and H. Sompolinsky (1985). *Spin-Glass Models of Neural Networks*. Physical Review A.
Alexander Atanasov, Blake Bordelon, and Cengiz Pehlevan (2022). *Neural networks as kernel learners: the silent alignment effect*. International Conference on Learning Representations (ICLR).
Yasaman Bahri, Jonathan Kadmon, Jeffrey Pennington, Samuel S. Schoenholz, Jascha Sohl-Dickstein, and Surya Ganguli (2020). *Statistical mechanics of deep learning*. Annual Review of Condensed Matter Physics.
Yasaman Bahri, Ethan Dyer, Jared Kaplan, Jaehoon Lee, and Utkarsh Sharma (2024). *Explaining neural scaling laws*. Proceedings of the National Academy of Sciences.
J"urgen Berges (2004). *Introduction to Nonequilibrium Quantum Field Theory*. AIP Conference Proceedings.
Emmy Blumenthal (2025). [*Building Intuition for Dynamical Mean-Field Theory: A Simple Model and the Cavity Method*](https://arxiv.org/abs/2507.16654). arXiv.
Blake Bordelon and Cengiz Pehlevan (2022). *Self-consistent dynamical field theory of kernel evolution in wide neural networks*. Advances in Neural Information Processing Systems (NeurIPS).
Blake Bordelon and Cengiz Pehlevan (2023). [*Dynamics of Finite Width Kernel and Prediction Fluctuations in Mean Field Neural Networks*](https://arxiv.org/abs/2304.03408). arXiv.
Blake Bordelon, Abdulkadir Canatar, and Cengiz Pehlevan (2020). *Spectrum dependent learning curves in kernel regression and wide neural networks*. International Conference on Machine Learning (ICML).
Blake Bordelon, Lorenzo Noci, Mufan Bill Li, Boris Hanin, and Cengiz Pehlevan (2023). [*Depthwise Hyperparameter Transfer in Residual Networks: Dynamics and Scaling Limit*](https://arxiv.org/abs/2309.16620). arXiv.
Blake Bordelon, Alexander Atanasov, and Cengiz Pehlevan (2024). *A dynamical model of neural scaling laws*. International Conference on Machine Learning (ICML).
Ari Brill (2024). *Neural scaling laws rooted in the data distribution*.
Bradley Brown, Jordan Juravsky, Ryan Ehrlich, Ronald Clark, Quoc V. Le, Christopher R'e, and Azalia Mirhoseini (2024). *Large language monkeys: scaling inference compute with repeated sampling*.
Abdulkadir Canatar, Blake Bordelon, and Cengiz Pehlevan (2021). *Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks*. Nature Communications.
Michael Celentano, Chen Cheng, and Andrea Montanari (2021). *The high-dimensional asymptotics of first order methods with random data*.
Zhongtian Chen, Edmund Lau, Jake Mendel, Susan Wei, and Daniel Murfet (2023). [*Dynamical versus Bayesian Phase Transitions in a Toy Model of Superposition*](https://arxiv.org/abs/2310.06301). arXiv preprint arXiv:2310.06301.
L'ena"ic Chizat and Francis Bach (2018). *On the global convergence of gradient descent for over-parameterized models using optimal transport*. Advances in Neural Information Processing Systems (NeurIPS).
L'ena"ic Chizat, Edouard Oyallon, and Francis Bach (2019). *On lazy training in differentiable programming*. Advances in Neural Information Processing Systems (NeurIPS).
Youngmin Cho and Lawrence K. Saul (2009). *Kernel methods for deep learning*. Advances in Neural Information Processing Systems (NeurIPS).
Carson C. Chow and Michael A. Buice (2012). [*Path Integral Methods for Stochastic Differential Equations*](https://arxiv.org/abs/1009.5966). arXiv.
Leticia F. Cugliandolo (2023). [*Recent Applications of Dynamical Mean-Field Methods*](https://arxiv.org/abs/2305.01229). arXiv.
Simon S. Du, Xiyu Zhai, Barnab'as P'oczos, and Aarti Singh (2019). [*Gradient Descent Provably Optimizes Over-parameterized Neural Networks*](https://arxiv.org/abs/1810.02054). International Conference on Learning Representations (ICLR).
Nelson Elhage, Tristan Hume, Catherine Olsson, Nicholas Schiefer, Tom Henighan, Shauna Kravec, Zac Hatfield-Dodds, Robert Lasenby, Dawn Drain, Carol Chen, Roger Grosse, Sam McCandlish, Jared Kaplan, Dario Amodei, Martin Wattenberg, and Christopher Olah (2022). [*Toy Models of Superposition*](https://arxiv.org/abs/2209.10652). arXiv preprint arXiv:2209.10652.
Andreas Engel and Christian Van den Broeck (2001). *Statistical Mechanics of Learning*. Cambridge University Press.
Zhou Fan, Justin Ko, Bruno Loureiro, Yue M. Lu, and Yandi Shen (2025). *Dynamical mean-field analysis of adaptive Langevin diffusions: propagation-of-chaos and convergence of the linear response*.
Leo Gao, John Schulman, and Jacob Hilton (2023). *Scaling laws for reward model overoptimization*. International Conference on Machine Learning (ICML).
Crispin W. Gardiner (2004). *Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences*. Springer.
Mario Geiger, Stefano Spigler, Arthur Jacot, and Matthieu Wyart (2020). *Disentangling feature and lazy training in deep neural networks*. Journal of Statistical Mechanics: Theory and Experiment.
Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg (1996). *Dynamical Mean-Field Theory of Strongly Correlated Fermion Systems and the Limit of Infinite Dimensions*. Reviews of Modern Physics.
C'edric Gerbelot, Emanuele Troiani, Francesca Mignacco, Florent Krzakala, and Lenka Zdeborov'a (2024). *Rigorous dynamical mean-field theory for stochastic gradient descent methods*. SIAM Journal on Mathematics of Data Science.
Nigel Goldenfeld (2018). *Lectures On Phase Transitions And The Renormalization Group*. CRC Press.
G'eza Gy"orgyi (1990). *First-order transition to perfect generalization in a neural network with binary synapses*. Physical Review A.
Boris Hanin (2023). [*Random Neural Networks in the Infinite Width Limit as Gaussian Processes*](https://arxiv.org/abs/2107.01562). The Annals of Applied Probability.
Boris Hanin (2026). *Special Lectures on the Mathematics of Deep Learning*.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun (2015). *Delving deep into rectifiers: surpassing human-level performance on ImageNet classification*. International Conference on Computer Vision (ICCV).
Tom Henighan, Jared Kaplan, Mor Katz, Mark Chen, Christopher Hesse, Jacob Jackson, Heewoo Jun, Tom B. Brown, Prafulla Dhariwal, Scott Gray, and others (2020). *Scaling laws for autoregressive generative modeling*.
Danny Hernandez, Jared Kaplan, Tom Henighan, and Sam McCandlish (2021). *Scaling laws for transfer*.
Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, and others (2022). *Training compute-optimal large language models*. Advances in Neural Information Processing Systems (NeurIPS).
Arthur Jacot, Franck Gabriel, and Cl'ement Hongler (2018). *Neural tangent kernel: convergence and generalization in neural networks*. Advances in Neural Information Processing Systems (NeurIPS).
Nikhil Kandpal, Haikang Deng, Adam Roberts, Eric Wallace, and Colin Raffel (2023). *Large language models struggle to learn long-tail knowledge*. International Conference on Machine Learning (ICML).
Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei (2020). *Scaling laws for neural language models*. arXiv preprint arXiv:2001.08361.
Devvrit Khatri, Lovish Madaan, Rishabh Tiwari, Rachit Bansal, Sai Surya Duvvuri, Manzil Zaheer, Inderjit S. Dhillon, David Brandfonbrener, and Rishabh Agarwal (2025). *The art of scaling reinforcement learning compute for LLMs*.
G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti (2006). [*Electronic structure calculations with dynamical mean-field theory*](https://link.aps.org/doi/10.1103/RevModPhys.78.865). Rev. Mod. Phys.
Ryogo Kubo (1957). *Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems*. Journal of the Physical Society of Japan.
R. Kubo (1966). *The Fluctuation-Dissipation Theorem*. Reports on Progress in Physics.
Tanishq Kumar, Blake Bordelon, Samuel J. Gershman, and Cengiz Pehlevan (2024). *Grokking as the transition from lazy to rich training dynamics*.
T. Lancaster and S. Blundell (2014). *Quantum Field Theory for the Gifted Amateur*. OUP Oxford.
Yann LeCun, L'eon Bottou, Genevieve B. Orr, and Klaus-Robert M"uller (1998). *Efficient BackProp*. Neural Networks: Tricks of the Trade.
Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S. Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein (2018). *Deep neural networks as Gaussian processes*. International Conference on Learning Representations (ICLR).
Chan Li and Nigel Goldenfeld (2026). [*Broken Ergodicity and the Violation of the Fluctuation-Dissipation Theorem Lead to Generalization Beyond Overfitting in Machine Learning*](https://arxiv.org/abs/2607.04135). arXiv.
J. Lin (2025). *Feature identification via the empirical NTK*.
Alexander Maloney, Daniel A. Roberts, and James Sully (2022). *A solvable model of neural scaling laws*. arXiv preprint arXiv:2210.16859.
Song Mei, Andrea Montanari, and Phan-Minh Nguyen (2018). *A mean field view of the landscape of two-layer neural networks*. Proceedings of the National Academy of Sciences.
Eric J. Michaud, Ziming Liu, Uzay Girit, and Max Tegmark (2023). *The quantization model of neural scaling*. Advances in Neural Information Processing Systems (NeurIPS).
Francesca Mignacco and Pierfrancesco Urbani (2022). [*The Effective Noise of Stochastic Gradient Descent*](https://arxiv.org/abs/2112.10852). Journal of Statistical Mechanics: Theory and Experiment.
Francesca Mignacco, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborov'a (2021). *Dynamical Mean-Field Theory for Stochastic Gradient Descent in Gaussian Mixture Classification**. Journal of Statistical Mechanics: Theory and Experiment.
Andrea Montanari and Pierfrancesco Urbani (2025). [*Dynamical Decoupling of Generalization and Overfitting in Large Two-Layer Networks*](https://arxiv.org/abs/2502.21269). arXiv.
Niklas Muennighoff, Alexander M. Rush, Boaz Barak, Teven Le Scao, Aleksandra Piktus, Nouamane Tazi, Sampo Pyysalo, Thomas Wolf, and Colin Raffel (2023). *Scaling data-constrained language models*. Advances in Neural Information Processing Systems (NeurIPS).
Radford M. Neal (1996). *Bayesian Learning for Neural Networks*. Springer.
Lars Onsager and Stefan Machlup (1953). *Fluctuations and irreversible processes*. Physical Review.
R. K. Pathria and Paul D. Beale (2011). *Statistical Mechanics*. Elsevier/Academic Press.
Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli (2016). *Exponential expressivity in deep neural networks through transient chaos*. Advances in Neural Information Processing Systems (NeurIPS).
Tomer Porian, Mitchell Wortsman, Jenia Jitsev, Ludwig Schmidt, and Yair Carmon (2024). *Resolving discrepancies in compute-optimal scaling of language models*. Advances in Neural Information Processing Systems (NeurIPS).
Alethea Power, Yuri Burda, Harri Edwards, Igor Babuschkin, and Vedant Misra (2022). *Grokking: generalization beyond overfitting on small algorithmic datasets*.
Charles Renshaw-Whitman (2026). *Iliad Physics Learning: lectures and exercise notebooks for the August 2026 Iliad Intensive*.
Zohar Ringel, Noa Rubin, Edo Mor, Moritz Helias, and Inbar Seroussi (2025). [*Applications of Statistical Field Theory in Deep Learning*](https://arxiv.org/abs/2502.18553). Now Publishers.
Hannes Risken (1996). *The Fokker--Planck Equation: Methods of Solution and Applications*. Springer.
Daniel A. Roberts, Sho Yaida, and Boris Hanin (2022). *The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks*. Cambridge University Press.
Grant M. Rotskoff and Eric Vanden-Eijnden (2022). *Trainability and accuracy of artificial neural networks: an interacting particle system approach*. Communications on Pure and Applied Mathematics.
Noa Rubin, Inbar Seroussi, and Zohar Ringel (2024). [*Grokking as a First Order Phase Transition in Two Layer Networks*](https://arxiv.org/abs/2310.03789). The Twelfth International Conference on Learning Representations (ICLR).
Rylan Schaeffer, Brando Miranda, and Sanmi Koyejo (2023). *Are emergent abilities of large language models a mirage?*. Advances in Neural Information Processing Systems (NeurIPS).
Samuel S. Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein (2017). *Deep information propagation*. International Conference on Learning Representations (ICLR).
Inbar Seroussi, Gadi Naveh, and Zohar Ringel (2023). *Separation of scales and a thermodynamic description of feature learning in some CNNs*. Nature Communications.
H. Sebastian Seung, Haim Sompolinsky, and Naftali Tishby (1992). *Statistical mechanics of learning from examples*. Physical Review A.
Utkarsh Sharma and Jared Kaplan (2020). *A neural scaling law from the dimension of the data manifold*.
James B. Simon, Madeline Dickens, Dhruva Karkada, and Michael R. DeWeese (2023). *The eigenlearning framework: a conservation law perspective on kernel ridge regression and wide neural networks*. Transactions on Machine Learning Research (TMLR).
Jamie Simon, Daniel Kunin, Alexander Atanasov, Jeremy Cohen, Arthur Jacot, Eric J. Michaud, Enric Boix-Adser\`a, Nikhil Ghosh, Mason Kamb, Berkan Ottlik, Blake Bordelon, Florentin Guth, Dhruva Karkada, and Joseph Turnbull (2026). *There Will Be a Scientific Theory of Deep Learning*. arXiv preprint arXiv:2604.21691.
Justin Sirignano and Konstantinos Spiliopoulos (2020). *Mean field analysis of neural networks: a law of large numbers*. SIAM Journal on Applied Mathematics.
Charlie Snell, Jaehoon Lee, Kelvin Xu, and Aviral Kumar (2024). *Scaling LLM test-time compute optimally can be more effective than scaling model parameters*.
H. Sompolinsky, A. Crisanti, and H. J. Sommers (1988). *Chaos in Random Neural Networks*. Physical Review Letters.
Konstantinos Spiliopoulos, Richard B. Sowers, and Justin Sirignano (2025). *Mathematical Foundations of Deep Learning Models and Algorithms*. American Mathematical Society.
Jonas Stapmanns, Tobias K"uhn, David Dahmen, Thomas Luu, Carsten Honerkamp, and Moritz Helias (2020). [*Self-Consistent Formulations for Stochastic Nonlinear Neuronal Dynamics*](https://arxiv.org/abs/1812.09345). Physical Review E.
Alain-Sol Sznitman (1991). *Topics in propagation of chaos*. 'Ecole d''Et'e de Probabilit'es de Saint-Flour XIX --- 1989.
David Tong (2006). [*Lectures on Quantum Field Theory*](http://cam.ac.uk).
David Tong (2012). *Statistical Field Theory*. University of Cambridge.
David Tong (2012). [*Lectures on Kinetic Theory*](https://www.damtp.cam.ac.uk/user/tong/kintheory/). University of Cambridge.
David Tong (2012). *Statistical Physics*. University of Cambridge.
Nicolaas G. van Kampen (2007). *Stochastic Processes in Physics and Chemistry*. North-Holland.
Jason Wei, Yi Tay, Rishi Bommasani, Colin Raffel, Barret Zoph, Sebastian Borgeaud, Dani Yogatama, Maarten Bosma, Denny Zhou, Donald Metzler, and others (2022). *Emergent abilities of large language models*. Transactions on Machine Learning Research.
Christopher K. I. Williams (1997). *Computing with infinite networks*. Advances in Neural Information Processing Systems (NeurIPS).
Lechao Xiao, Yasaman Bahri, Jascha Sohl-Dickstein, Samuel S. Schoenholz, and Jeffrey Pennington (2018). *Dynamical isometry and a mean field theory of CNNs: how to train 10,000-layer vanilla convolutional neural networks*. International Conference on Machine Learning (ICML).
Greg Yang and Edward J. Hu (2021). *Tensor programs IV: feature learning in infinite-width neural networks*. International Conference on Machine Learning (ICML).
Greg Yang (2019). [*Tensor Programs I: Wide Feedforward or Recurrent Neural Networks of Any Architecture are Gaussian Processes*](https://arxiv.org/abs/1910.12478). Advances in Neural Information Processing Systems 32 (NeurIPS 2019).
Greg Yang, Edward J. Hu, Igor Babuschkin, Szymon Sidor, Xiaodong Liu, David Farhi, Nick Ryder, Jakub Pachocki, Weizhu Chen, and Jianfeng Gao (2022). *Tensor programs V: tuning large neural networks via zero-shot hyperparameter transfer*. Advances in Neural Information Processing Systems (NeurIPS).
Lenka Zdeborov'a and Florent Krzakala (2016). *Statistical physics of inference: thresholds and algorithms*. Advances in Physics.
J. Zinn-Justin (2002). *Quantum Field Theory and Critical Phenomena*. Clarendon Press.
[^1]: A physicist's free energy is $-\beta^{-1}\log Z$; we absorb the factor of $\beta$ throughout, so that $F$ is dimensionless.
[^2]: For transformers, the limit analogous to "large width" that leads to a Gaussian process is a large number of attention heads.
[^3]: As a point of notation, $\eta$ and $\zeta$ are also commonly used for noise, and we will distinguish between discrete and continuum delta functions as $\delta_{ij}$ and $\delta(t-t')$ respectively.
[^4]: Michaud's blog post [ericjmichaud.com/quanta](https://ericjmichaud.com/quanta/) is an excellent companion to the paper and can substitute for it in this exercise. But attempt [Exercise 7.1](#ex-quanta) *before* reading either; the derivations are short, and they are more instructive to find than to follow.
[^5]: This is also a working example for [Section 6](#6-non-equilibrium-physics-msrjd-and-dmft)'s toolkit: the same model can be solved dynamically, and its double-descent peak, at the interpolation threshold $P\sim n$, is regularized by ridge or early stopping, connecting back to the divergence of the overfitting factor mentioned in [Section 7.7](#77-exponents-from-geometry-and-from-spectra).
[^6]: This appendix adapts a lecture and an exercise notebook prepared by Charles Renshaw-Whitman for the Iliad Intensive, translated into the notation of [Section 3.1](#31-notation).