Physics of Deep Learning
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.
By Aleksander Czejdo, Tudor Dimofte, Brianna Grado-White, Charles Renshaw-Whitman
Contents
- -1. Prerequisites
- 0. A one-day teaching plan
- 1. Intro reading: learning mechanics
- 2. Statistical mechanics of learning
- 3. Large width at initialization: the Gaussian process limit
- 4. The NTK and large width during training
- 5. Mean-field scaling
- 5.1 Mean-field theories in physics
- 5.2 Mean-field scaling for one-hidden-layer networks
- 5.3 The NTK is no longer constant: feature learning
- 5.4 The scaling dial and Maximal Update Parameterization
- 5.5 What is this framework good for?
- 5.6 Grokking as a lazy-to-rich transition: a reading exercise
- 5.7 Bayesian grokking: a first-order transition (alternative reading)
- 6. Non-equilibrium physics, MSRJD, and DMFT
- 7. Scaling laws
- 7.1 The empirical laws
- 7.2 Joint laws, bottlenecks, and the Chinchilla revision
- 7.3 Where do the laws come from?
- 7.4 A survey of theoretical proposals
- 7.5 The quantization model: exponents from the statistics of skills
- 7.6 Emergence: smooth laws from staggered cliffs
- 7.7 Exponents from geometry and from spectra
- 7.8 A solvable large-N field theory
- 7.9 All the knobs at once: a dynamical theory
- 7.10 Beyond pretraining: post-training and test-time scaling
- A. The replica method
- B. From linear regression to kernel regression
- 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 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. 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 |
| 11:00–11:15 | Break | |
| 11:15–12:00 | Exercises, Part I (choose one) | Exercises 2.1, 2.2 and 2.3 |
| 12:00–12:30 | Lecture, Part IIa: large width at initialization; the NTK | Sections 3.1–3.4, Section 3.6, Sections 4.1–4.2 |
| 12:30–1:30 | Lunch | |
| 1:30–2:00 | Lecture, Part IIb: large width during training; mean-field scaling | Sections 4.3–4.4, Sections 5.2–5.4 |
| 2:00–3:00 | Exercises or lab, Part II (choose one) | Exercises 3.3, 4.5, 4.4 and 5.3 |
| 3:00–3:45 | Lecture, Part III: neural scaling laws, in theory and practice | Sections 7.1–7.5 |
| 3:45–4:15 | Break | |
| 4:15–5:00 | Exercise, Part III | Exercise 7.1 |
| 5:00–5:30 | Reading on "learning mechanics" | Section 1 |
| 5:30–5:50 | Discussion in small groups | the prompts of Section 1 |
| 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; Tong 2012), the physics behind the lecture; Engel and Van den Broeck (Engel & Van den Broeck 2001) for the perceptron tradition; and Gyorgyi (Gy"orgyi 1990) for the Ising perceptron transition. For Part II: Hanin's lectures (Hanin 2026) and Roberts, Yaida, and Hanin (Roberts et al. 2022) on large-width limits and the expansion; Ringel et al. (Ringel et al. 2025) on statistical field theory in deep learning; Spiliopoulos, Sowers, and Sirignano (Spiliopoulos et al. 2025) for rigorous NTK and mean-field limits; Lin (Lin 2025) for experiments identifying features through the empirical NTK; and the two accounts of grokking (Kumar et al. 2024; Rubin et al. 2024), which were not lectured but are the natural follow-up reading (Sections 5.6 and 5.7). For Part III: Kaplan et al. (Kaplan et al. 2020) and Henighan et al. (Henighan et al. 2020) for the empirical laws; Hoffmann et al. (Hoffmann et al. 2022) and Porian et al. (Porian et al. 2024) for Chinchilla and the reconciliation with Kaplan; Michaud et al. (Michaud et al. 2023) for the quantization model; and Maloney, Roberts, and Sully (Maloney et al. 2022) and Bordelon, Atanasov, and Pehlevan (Bordelon et al. 2024) 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); Appendix B is adapted from it.
1. Intro reading: learning mechanics
The review (Simon et al. 2026) 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), and think about the following:
- Rank the physics analogies in Table 1 of the reading (solvable models hydrogen atom; limits thermodynamic limit; empirical laws Kepler and Planck; hyperparameters Reynolds number; universality critical phenomena). Which are the most compelling? Which do you think are most useful in understanding deep learning?
- 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?
- 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); modern reviews include (Zdeborov'a & Krzakala 2016; Bahri et al. 2020); and an excellent reference for the physics itself is David Tong's lecture course on statistical field theory (Tong 2012), which we follow closely in Sections 2.4–2.5.
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). We then switch to physics. We introduce partition functions, free energies, and cumulants (Section 2.3), 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). 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). Section 2.7 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 with parameters , trained on samples to minimize a loss
a sum over samples of a per-sample loss . In practice, the network is initialized at random, , and trained with some version of gradient descent,
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
Here the initialization distribution plays the role of the Bayesian prior, and is interpreted as the likelihood of the training data given the parameters (Section 2.2 explains why this is reasonable). The positive constant , 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 . It is convenient to write the loss as times an average over samples,
The empirical loss stays finite as , converging to the population loss , the expectation over the true data distribution. Thus for large the posterior is
and the sample size always accompanies the inverse temperature, in the combination . Increasing the number of samples is the same as lowering the temperature. If, moreover, samples are drawn fresh at each training step, then is also a rough measure of training time. Altogether: more data more training larger colder. We will often write for short, but it is worth remembering that the effective grows with .
Bayes vs. gradient descent. We do not train networks by sampling from (2.3); 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, (more samples, longer training) (lower temperature). As we will see, sweeping a parameter like 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), 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.
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 is Gaussian noise with mean and covariance , is the learning rate, and 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 ), the distribution converges at late times to the Gibbs distribution
(See (Tong 2012, Chapter 3) for a physics-focused discussion of this convergence.) A prior can be included by training on the modified loss , which reproduces (2.3). Thus one can engineer the Bayesian posterior by Langevin training, and this is how theoretical predictions derived from (2.3) are usually tested. Note that the Langevin route gives a different Bayes-to-gradient connection from the first: here training time at fixed gives the posterior at , while at fixed gives gradient flow.
Remark. SGD itself is more complicated than Langevin dynamics. The covariance 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 (gradient flow) versus Section 4.6 (Bayesian); for toy models of grokking in Section 5 (Rubin et al. 2024); and, for toy models of superposition (Elhage et al. 2022) through the lens of singular learning theory, in (Chen et al. 2023).
2.2 Why the exponential of a loss is a likelihood
The Bayesian model treats 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, with . Then the probability density of a label given the input and the parameters is
and for 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 over classes, typically the softmax of its output values. The cross-entropy loss is , so that
exactly, with : 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, is the likelihood and (2.3) at is the textbook Bayesian posterior. A value gives a tempered posterior. It arises naturally when the noise level is not known (as in (2.9)), and it is also a useful theoretical knob: as the interpretation of above shows, temperature and sample size are interchangeable in the Bayesian model. There is also a useful limiting case. Some losses are hard constraints, if the prediction is right and if it is wrong. Then is a step function, the posterior is uniform on the parameters that fit every training sample, and no appears at all. This is the zero-temperature limit, and it is the setting of the Ising perceptron in Section 2.6.
2.3 Partition functions, free energy, and cumulants
The Bayesian posterior (2.3) 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 | microscopic state of a system |
| (e.g. the spins of electrons) | |
| loss | energy , or Euclidean action |
| posterior | Boltzmann distribution at temperature |
| number of samples , or training time | inverse temperature |
| #parameters | #degrees of freedom |
To include the prior , simply redefine the energy as , so that . In what follows we write for whatever sits in the exponent, and for the number of degrees of freedom, as physicists do; for a network, .
The thermodynamic limit. Statistical mechanics is about systems with very many degrees of freedom, and one studies their behavior as : the thermodynamic limit. In this limit one distinguishes extensive quantities, which grow in proportion to (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 parameters and samples, the loss (2.1) is extensive in and the empirical loss 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),
is called the partition function. In Bayesian terms it is the marginal likelihood of the data. It depends on and on any other hyperparameters of the problem. Its logarithm defines the free energy1 and the free energy density,
There are three reasons to care about rather than . First, is extensive and intensive: the free energy density has a limit as , while itself is exponentially large or small. Second, " is minimized": at large the integral (2.11) is dominated by the regions of parameter space with the smallest local free energy, so the system's macroscopic state is found by minimizing . Section 2.4 makes this precise. Third, derivatives of compute expectation values, as follows.
Sources and cumulants. Suppose we care about the expectation values of some observables , functions of the parameters. Add source terms to the exponent,
Differentiating under the integral,
and in general the -th derivative is the -th cumulant,
the mean, the covariance, the third central moment, and from 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, 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: already couples to , 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 , an observable , and the sourced partition function , with . All expectations are in , that is, at .
(a) Show that and .
(b) The third cumulant is defined by continuing the pattern, . Compute it in terms of the moments , , , and show that it equals the third centered moment .
(c) Show that for a Gaussian , is exactly quadratic in , so that all cumulants with vanish. Compute 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).
Solution
(a) Differentiating under the integral, , which at is ; since , this is . Once more, .
(b) Using on gives . Expanding gives the same expression.
(c) For (or more generally when is Gaussian under ), , so is quadratic and every derivative beyond the second vanishes. One more derivative of the pattern in (b), with , gives ; for a Gaussian , so while the fourth centered moment is not.
(d) , so ; , so ; and by the same computation as (b) with , so .
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 samples were drawn, and at large 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 is usually hard, and the annealed approximation swaps the order,
which is an upper bound by Jensen's inequality and often a good approximation. We will use it in Section 2.6. The exact treatment goes through the replica method, sketched in Appendix A; the classic reference is (Engel & Van den Broeck 2001).
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 , on the most famous model in statistical mechanics. For the physics of this example, well beyond what we cover here, see (Tong 2012, Chapter 1).
The model. Let be a graph with vertices and undirected edges ; typically is a regular lattice in dimensions, and every vertex has the same number of neighbors ( for a hypercubic lattice). At each vertex sits a spin . A state of the system is a choice of all spins, so there are states, and the energy of a state is
The coupling rewards neighboring spins for aligning with each other, and the external magnetic field rewards spins for aligning with it. At inverse temperature the Boltzmann distribution is with . Dimensionally, only the combinations and 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 -dimensional microscopic state for the purpose of the questions we want to ask. The question here is whether the spins align ( close to , an ordered phase) or not ( near , a disordered phase), and how this depends on , , and .
Mean-field approximation. We would like to push the distribution over states forward to a distribution over the single variable . 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 by ,
(The factor compensates for counting each edge twice.) The energy becomes a function of 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 and quantitatively right for (Tong 2012, 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 as the prototype of mean-field limits of neural networks.
Entropy, and the reduction to one variable. With the energy depending only on , the sum over states can be organized by the value of :
since a state with magnetization has exactly spins up. Stirling's formula gives
the entropy per site at magnetization : the log of the number of microscopic states compatible with the macroscopic value . It is largest at , where and almost all of the states live, and it vanishes at , where exactly one state lives. Replacing the sum over by an integral (the spacing is ), the partition function becomes a one-dimensional integral,
The function is called the effective action, or the free-energy density at fixed : energy minus temperature times entropy, per site. We have reduced the problem from 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 the integrand is exponentially peaked at the minima of , and the integral is dominated by a window of width around the lowest one. This is the saddle-point approximation (also called Laplace's method), and it gives
The free energy is minimized, in exactly the sense promised in Section 2.3. The minimizer is the magnetization the system actually displays. Setting and using ,
This transcendental equation is best solved graphically (Figure 1). Depending on , , and it has one solution, or three; when there are three, the outer two are minima of and the middle one is a maximum.
Graphical solution of the mean-field equation (2.27): intersections of with the diagonal. Left: high temperature, one solution at . Middle: low temperature at , three solutions, of which is a maximum of and are degenerate minima. Right: low temperature with , three solutions, with the positive one the global minimum.
Self-consistency (an aside). There is another route to (2.27), which is the form in which mean-field theory usually appears and which we will use in Section 5. Look at the effective field felt by a single spin :
A single spin in this field has partition function and mean . If the approximation is to be consistent, the average that comes out must equal the average we put in, , and that is (2.27). 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 , which control where the probability concentrates, as , , and are varied. Taylor-expanding the entropy, , 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 : the coupling beats the temperature, the quadratic coefficient in (2.29) is negative, and is a double well (Figure 2). The system wants to be ordered, with neighboring spins aligned; but aligned which way? The field decides. For the global minimum is near , for near , and at the two minima are exactly degenerate. As is swept through , the magnetization the system displays jumps from to . This is a first-order phase transition.
The action of (2.25) in the ordered phase (), for a negative, zero, and positive external field. The two minima tilt with ; the global minimum (orange) switches sides at . 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) is a derivative of the free energy with respect to the field:
At finite , the two minima of are two competing phases, and is a probability-weighted average of their magnetizations. As the saddle-point approximation puts all the weight on the lower minimum, so develops a discontinuity at , and itself develops a kink (Figure 3, 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
is a peak whose height grows like and whose width shrinks like : a nascent delta function (Figure 3, right). Genuine non-analyticity can occur only at : at finite the partition function is a finite sum of analytic functions of , so everything is smooth, and what one sees is a rapid crossover that sharpens as grows. Watching the peak in grow with , rather than saturate, is the standard numerical diagnostic that separates a true transition from a crossover. Exercise 2.2 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).
Fingerprints of a first-order transition, in the mean-field Ising model at , computed from the one-dimensional integral (2.25) at finite . Left: the expected magnetization as a function of the field is a sigmoid of width , sharpening into a step. Right: the susceptibility is a peak of height and width .
Metastability and hysteresis. The suppressed minimum in Figure 2 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 . In practice the system stays in the metastable phase until the field is pushed well past , and the transition observed dynamically lags the thermodynamic one. This is hysteresis, familiar from real magnets. It is the reason the " training time" analogy of Section 2.1 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.
A second-order transition. An -th order phase transition is a discontinuity in the -th derivative of the free-energy density as . The Ising model has a second-order transition too, at fixed , as the temperature crosses (Figure 4). Here nothing crosses anything: the single basin at flattens and splits into two symmetric ones. The quadratic coefficient in (2.29) changes sign at , and the minima sit at
so the order parameter moves continuously but not smoothly; and are continuous, and the discontinuity sits in (the heat capacity). The point is a critical point. Below it the two ordered states are exact mirror images, and the system must pick one: this is spontaneous symmetry breaking, the symmetry of the energy at being broken by the state. The exponent in (2.32) 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 in (Onsager), and there is no transition at all in . Mean-field exponents become exact for . Understanding why is the subject of the renormalization group, for which see (Tong 2012).
The second-order transition of the mean-field Ising model at . Left: the action for temperatures above, at, and below ; a single basin flattens and splits into two symmetric ones. Right: the resulting magnetization , continuous at but with the square-root singularity (2.32) (dashed).
Back to learning. In the Bayesian model of Section 2.1, the natural control parameter is , the amount of data or training. The Ising analogue is sweeping at fixed , 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 is the cleanest example. Before turning to it, the following exercise establishes what every first-order transition looks like at finite .
Exercise 2.2 (★) (Anatomy of a first-order transition). Consider any system whose partition function has the form , depending on some hyperparameter . Here could be the inverse temperature, an external field, or the load of Section 2.6; we have absorbed into to keep the setup general. Suppose that has exactly two isolated minima, at and . Then at large the integral is dominated by two competing phases. Letting be small neighborhoods of ,
with the free-energy densities of the two basins. Suppose further that as crosses a critical value the two minimum values cross: with , where a prime denotes .
(a) Occupation. Show that the probability of finding the system in basin is
and that, expanding to first order in ,
a logistic function of . Over what window of does the system switch allegiance?
(b) First derivative. Show from (2.33) that
the probability-weighted average of the branch slopes. Sketch it at finite and at : a sigmoid of height sharpening into a jump. The kink in itself is the defining discontinuity of a first-order transition. Show also, directly from the integral, that : the observable that jumps is the one conjugate to the control parameter. When it is the energy, and the jump is the latent heat; for the Ising model with the field, it is the magnetization.
(c) Second derivative. Using (verify this), show that
Sketch this curve. Near , where , the first term is a peak of height and width , which converges to a delta function as . Show also that when does not depend on , . (For this is the heat capacity; for the Ising field it is the susceptibility.)
(d) Plot it. Take the minimal model , and plot , , and for a few values of . (Use a numerically stable .) Watching the peak grow with rather than saturate is the standard diagnostic separating a true transition from a smooth crossover.
Solution
(a) is the fraction of contributed by basin . With near the crossing, . The switch happens over the window , which shrinks with system size.
(b) Differentiating gives exactly; near the slopes are . Directly from the integral, : the control parameter is itself a source, coupling to . Sketch: at finite , is a sigmoid interpolating from to over the window of (a); as , , continuous with a kink at , and becomes a step of height .
(c) and . Then . Sketch: a smooth background from the terms plus a peak at of height and width , the nascent delta function of the kink. For the variance identity: with independent of , differentiating once more gives , so , which is why the curve is a peak and not a dip. The two-branch formula exhibits the leading piece of this variance: in one basin and in the other, and the variance of that two-point mixture is exactly .
(d) See Figure 14. The peak in the right panel grows in proportion to ; a crossover's bump would saturate.
Anatomy of a first-order transition at finite , in the minimal two-branch model , , with the control parameter written as . Left: the branch free energies cross at ; the limiting free-energy density is their minimum, which has a kink. Middle: at finite , is the probability-weighted average of the branch slopes (2.36), a sigmoid of width . Right: the susceptibility (2.37) peaks at with height and width .
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); see (Engel & Van den Broeck 2001, Chapter 7). It is presented as a walk-through exercise; solutions are in Appendix C, but try each part before reading them.
Setup. A teacher is a fixed vector of binary weights . Inputs are Gaussian, , and the teacher labels them by
A student is another binary vector , guessing labels the same way, . This model is called the Ising perceptron, because its weights take the values of Ising spins. We observe examples with , and do Bayesian learning at zero temperature with a uniform prior, in the sense of Section 2.2: the per-sample loss is if the student classifies the example correctly and otherwise, so the posterior is uniform on the students that classify every example correctly. Since teacher and student agree on exactly when ,
with for and otherwise. With the factor , is the fraction of all students consistent with the data. Both and will be taken large, at fixed ratio: the load
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 and the alignment variables . Show that : the overlap is the magnetization of the spins . Learning the teacher means magnetizing them to .
(b) Generalization error. Show that a student at overlap misclassifies a fresh Gaussian input with probability
[Hint: project onto the plane spanned by and ; the projection is an isotropic two-dimensional Gaussian, and the two decision boundaries are diameters meeting at angle . When do the signs disagree?]
(c) Entropy: the same as the magnet. Show that the number of students at overlap is , so that by Stirling's formula
identical to the magnet's entropy (2.24) with .
(d) Push forward to one variable. Average over datasets. The examples are independent, and each is classified correctly by a student at overlap with probability , so
Trade the sum over for an integral over using (c), and show that at fixed load and large ,
This is the exact analogue of the magnet's (2.25): an entropy term counting states, and a "field" term rewarding alignment. The aligning field is now the data, and its strength is the load . (Averaging rather than is the annealed approximation (2.18). It preserves all the qualitative structure here and shifts the numbers only slightly; Appendix A explains how to do the exact calculation.)
(e) Plot for several values of and locate its critical points. Do you expect a phase transition? At which ?
(f) The teacher basin. is the perfectly aligned state, where the student is the teacher. Show that for every (zero entropy, but zero data cost), and that is always a local minimum: writing , deviating from the teacher gains entropy but costs data-fit (expand ), and . 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 , is minimized at with : the disordered cloud of typical students dominates. As grows, this interior minimum drifts toward and its value rises. Numerically, find the load at which the two basins exchange dominance, . You should find , with the interior minimum still at , that is, generalization error . Conclude that at the typical student jumps from error to zero error: a first-order transition to perfect generalization. Then connect to Exercise 2.2 with as the control parameter: what do and look like at finite , 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.
(h) The dictionary. Complete the table:
| Ising magnet | Ising perceptron |
|---|---|
| spins | alignments |
| magnetization | overlap |
| entropy | entropy (the same function) |
| aligning field | ? |
| knob swept across the transition | ? |
| ordered phase | ? |
| jump in as a function of | ? |
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 , , where the slice at overlap has volume ? Appendix A.3 works this case out.)
Solution
(a) since . The map is a bijection of the hypercube (a "gauge transformation"), so the uniform measure goes to the uniform measure, and the perceptron's order parameter is literally a magnetization.
(b) Both labels depend only on the projection of onto , which is an isotropic two-dimensional Gaussian. The decision boundaries are two diameters at angle ; the signs disagree exactly on the double wedge between them, of angular fraction . Note that the answer does not depend on , and holds for any with , binary or not.
(c) A student at overlap agrees with the teacher on exactly coordinates (), giving the binomial count; Stirling gives , the binary entropy at magnetization . Note and : the latter is finite, because exactly one state sits at .
(d) Independence over examples gives the -th power. Trading the sum for an integral, , and ; collecting exponents gives as stated. (The normalizes at .)
(e) The curves are shown in Figure 15. At , has a single minimum at , the bulk of the hypercube. As grows, the interior minimum drifts toward larger and its value rises (note ), while stays pinned at . For large enough there are two competing local minima separated by a barrier: exactly the two-basin structure of Exercise 2.2. One should therefore expect a first-order transition at the load where the interior minimum's value crosses , near from the plot; at the interior minimum survives only as a metastable dip above . No transition is expected at small , where the interior minimum is the only relevant one and moves smoothly.
The annealed action of the Ising perceptron, (2.44), for several loads. Interior local minima are marked with circles; the teacher state , where exactly, with squares. The two basins exchange dominance near .
(f) At : and , so for all . Near with : , while gives data cost . Since as , : a basin around the teacher at every . It is protected by the singular cost of the first errors, which is what a smooth well cannot provide.
(g) Minimizing over the interior numerically (Figure 16, left), the interior branch value crosses at , with minimizer and . The order parameter jumps from to (middle panel). With as the control parameter of Exercise 2.2: has a kink at ; at finite , is a sigmoid interpolating between the branch slopes and , over a window of width (right panel), and is a peak of height . The role of the latent heat is played by the jump in , the per-example log-likelihood released when the student snaps onto the teacher. Past the interior minimum survives as a local minimum until the spinodal , where it merges with the barrier and disappears. In the window a local learning dynamics (flipping a few weights at a time to improve alignment) can stay stuck at even though perfect generalization is thermodynamically dominant: the information needed to identify the teacher is present in the data before the dynamics can find it, and generalization arrives abruptly and later than the data first permits. This is the metastability of Section 2.5, and a prototype of the delayed transitions of Section 5.7. (The exact, quenched treatment (Gy"orgyi 1990; Engel & Van den Broeck 2001) moves the transition to but changes nothing structurally.)
The Ising perceptron's phase transition. Left: the free-energy densities of the interior basin, , and of the teacher basin, , cross at ; the limiting is their minimum, with a kink. Middle: the order parameter jumps from to ; the dashed continuation is the metastable interior minimum, which disappears at the spinodal . Right: at finite , , computed from the integral (2.44), is a sigmoid sharpening into the step of height .
(h) Field the data term, the per-example log-likelihood , with strength ; the knob swept the load , the amount of data; the ordered phase perfect generalization ; the jump in as a function of the discontinuous drop of from to as a function of . The structural difference: in the magnet the two competing minima are symmetry partners, while here they are an entropy-dominated basin and a data-dominated basin. The ingredient that makes the sharp transition possible is being finite: discreteness gives the perfectly aligned state a finite entropy, so it can compete as a separate basin. For spherical weights the slice entropy diverges to as , the perfectly aligned state is infinitely penalized, the minimum always sits in the interior and slides continuously with , and there is no transition; the generalization error then decays smoothly, as in the exact treatment (Appendix A.3).
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 , partition functions, free energies. For a field theorist, a partition function with sources,
is a moment generating functional; its logarithm generates connected correlation functions (cumulants); the action 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 is played variously by the weights of a network, by an order parameter such as the magnetization or the overlap of Sections 2.4 and 2.6, or — most usefully — by the network function itself, with the input playing the role of the spacetime coordinate. Functional integrals over the network function are developed carefully in Sections 3.3 and 3.5, 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 : the weights are real and positive, as after the standard continuation
The role of — the parameter controlling the size of fluctuations around saddle points — is played here by temperature, or by at large width, or by at large sample size.
When we eventually study the dynamics of training as a field theory in Section 6, 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; Tong 2012; Tong 2012; Tong 2006) are a friendly introduction to statistical physics, statistical field theory, kinetic theory, and QFT respectively; Lancaster and Blundell (Lancaster & Blundell 2014) is a gentle conceptual route into QFT; Altland and Simons (Altland & Simons 2010) bridges toward condensed-matter field theory (and contains the Langevin/MSRJD material used in Section 6); Pathria (Pathria & Beale 2011) is a standard graduate statistical mechanics text; Goldenfeld (Goldenfeld 2018) is the classic on phase transitions and the renormalization group; and Zinn-Justin (Zinn-Justin 2002) 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; Spiliopoulos et al. 2025; Ringel et al. 2025) and the online lectures (Hanin 2026).
We'll begin in this section with a classic result (Neal 1996; Lee et al. 2018; Yang 2019; Hanin 2023) 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 in Section 5.5.
3.1 Notation
We'll use the following notation for neural networks and training data.
- trainable network parameters (weights, biases)
- training data samples, with and (unless otherwise stated)
- network function
For MLP's in particular (also called FCN's = fully connected networks)
- pre-activations in layer
- number of hidden layers (depth); width of 'th layer; and . We will often take constant hidden depth, .
- activation function (assumed the same in all layers except the first)
- network recursion with and ; parameters are . The explicit factors of are a normalization convention, explained in Section 3.2: 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 with layer widths as in Section 3.1. Suppose we keep fixed but make the 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 fixed as the widths grow. (Here and throughout, denotes a Gaussian of variance . 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
(where we've written and as vectors) we get
Here we've conditioned on the previous layer, and used that for independent mean-zero and arbitrary , , 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 . (We'll sharpen this estimate in Section 3.7.) Thus
The prefactor , 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). In Section 4 it will also be called the NTK normalization, and in Section 5 we will meet the alternative mean-field normalization, in which the output layer carries instead of .
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) is a change of variables on parameter space, and the network function is the same function of . 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. Keeping every parameter of order one, with the factors of 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 to sit in the variances instead, and translate with (3.5).
3.3 Gaussian processes
A Gaussian process with mean function and covariance kernel is a random function such that on any inputs the vector of function values is Gaussian
Here and are fixed, deterministic functions that control the mean and covariance of any tuple of values of . The function , often called the kernel of the process, is symmetric: .
Let `' denote a measure on (it need not be the usual Lebesgue measure). If , as integral kernel, has an inverse satisfying
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 `' has a probability measure
a Gaussian process has a (non-rigorous) measure
If the kernel and the measure theory on the input space are nice enough to allow to be diagonalized, with a discrete basis of orthonormal eigenfunctions and eigenvalues (so that has eigenvalues ), then can be understood as an infinite-dimensional Gaussian. Writing ,
One virtue of the physics formalism is that it organizes expectation values of products of 's. Let for simplicity (to generalize to any , just shift below). Let be an arbitrary (nice...) function and define
The integral evaluates formally to
Then expectation values can be expressed via derivatives, which by (3.12) become
The sum is over partitions of into pairs, and vanishes if is odd. Eqn. (3.14) 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) by completing the square in (3.11), and then derive (3.14) by differentiating. Check the first two cases explicitly: , and (with ). How many terms appear at ?
Solution
Complete the square: shifting (schematically, ) decouples the source, leaving . Each derivative acting on the exponential either brings down a factor (which must later be hit by another derivative, else it dies at ) or closes a pair by hitting a previously produced factor, contributing . Every surviving term therefore corresponds to a perfect pairing of . For this gives ; for , the three pairings ; for there are terms; in general .
Another virtue of the physics formalism is that it neatly connects higher cumulants to deviations from Gaussianity. In general, the -th cumulant of random variables is defined recursively via
For example, , . In physics, the '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,
as can easily be seen directly from (3.14), where joint expectations for all factor into products of pairs. Thus, nonvanishing higher 's are a measure of non-Gaussianity.
This appears in the functional-integral formalism as follows. Suppose that is a random function that's governed by a functional probability distribution
for some functional . (In physics is called the action.) Then just as in (3.11), we can introduce the generating function
The cumulants appear naturally as derivatives of , otherwise known as the (sourced) free energy:
which in turn means (formally) that
a generating function of cumulants. The original is the inverse Legendre transform of (by definition of ). Expanded as a formal series in the cumulants (assuming they are small), and assuming for simplicity, it takes the form
Physicists might recognize this as a sum of Feynman diagrams, with cumulants giving the vertices. A key structure is that the -th order term in the effective action contains only and lower cumulants. In particular, is a Gaussian process if and only if is quadratic; and small deviations from Gaussianity, measured by cumulants, provide the higher-order terms in .
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, the network function at initialization converges to a Gaussian process. The theorem appeared in the 90's (Neal 1996) for a single hidden layer, and was generalized more recently to deeper networks and different architectures, cf. (Lee et al. 2018; Yang 2019; Hanin 2023).
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.
For simplicity we'll assume that biases have been set to zero, and that all weights have the same order-one variance , as in (3.1), with the fan-in normalization of Section 3.2.
In the first network layer, the preactivations are
At fixed , each is a sum of i.i.d. Gaussian variables (with coefficients ), and thus is Gaussian itself, with and
thus is a Gaussian process (GP) with
In the second layer, we've got
If we condition on the preactivations of the first layer, the same argument gives
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: . 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) will become Gaussian as . More concretely, we can calculate and bound the variance
as long as the first-layer preactivations are order one (independent of ), and the activation function is reasonable (any standard choice will be). Thus, as , we can replace the stochastic kernel by its expectation,
Each subsequent layer is the same. At infinite width, the kernels become deterministic, and obey a recursion
with ultimately governing a GP for the network function .
Remark (exact kernels). For a few special activation functions, the Gaussian expectation in (3.30) can be evaluated in closed form, giving explicit NNGP kernels. Write , and for the previous layer's kernel values. For ReLU, , one finds the arc-cosine kernel of Cho and Saul (Cho & Saul 2009),
and for the sigmoid-like activation Williams (Williams 1997) found
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)–(3.32) is a pleasant but somewhat lengthy Gaussian-integral workout.)
Exercise 3.2 (The NNGP correspondence, empirically (lab)). Fix two unit-norm inputs with, say, and . For widths , sample independent initializations of the one-hidden-layer erf network with (that is, ).
(a) Make a histogram of the values across the ensemble, and overlay the Gaussian density with the variance predicted by (3.32); watch the visibly non-Gaussian histogram converge.
(b) Make a scatter plot of the pairs and compare the sample covariance with the predicted kernel matrix.
(c) For each initialization, record the stochastic kernel of (3.26), and verify that its variance across initializations decays like (that's the mechanism by which the kernel becomes deterministic).
Solution
Everything is closed-form thanks to the erf kernel (3.32), which is the point of choosing erf. Expected observations:
(a) at the histogram has visibly heavy tails and excess kurtosis (the output is a sum of only three i.i.d. non-Gaussian terms); at deviations are still detectable; at the histogram is indistinguishable from the predicted Gaussian at this sample size.
(b) The sample covariance matches the kernel matrix built from (3.32) to within the Monte Carlo error .
(c) A log–log plot of against has slope , as in the variance computation of Section 3.4. The whole lab is 20 lines of NumPy; the one subtlety is to histogram over initializations at fixed inputs, not over inputs.
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, ,
The initial distribution of weights, all of variance , is
To push this forward to a distribution on the first-layer preactivations , we integrate against a delta function:
Here the are arbitrary functions of , and are only tied to the network values via the delta functions. Eqn. (3.35) represents a functional probability density on the functions . To evaluate it, introduce a second set of "Lagrange multiplier" functions , and write
and then do the exact Gaussian integral over 's, and then over 's:
We clearly see the first-layer GP kernel appearing.
To push-forward to the network outputs , we just do this again:
The problem is that it's now impossible to do the integral over the 's exactly, because (for general activation function ) the action is no longer quadratic in them. However, as , we can treat as a perturbation, and approximate the exponent by its expectation value under the Gaussian measure (with kernel ) on the 's. To be clear, the approximation that's happening is of the form
Up to corrections, that produces the deterministic kernel of the previous section,
and we get
3.6 L/N
The inductive argument of Section 3.4 does not hold at infinite . There are two potential problems.
First, a factor of 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. With careful tuning, the issue can be avoided. For example, in a deep linear network (DLN) with , one can just take , known as LeCun initialization (LeCun et al. 1998). In a ReLU network, the correct choice instead is , known as He initialization (He et al. 2015).
Even with 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 , but rather the ratio of depth-to-width that controls the effective behavior of a network — including its departure from pure Gaussianity at large .
This behavior is analyzed in Chapters 3-4 of (Roberts et al. 2022) 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 , the leading non-Gaussianity is
where the expansion on the RHS is at finite and . For small , the non-Gaussianity remains small but finite, of size . For large , 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, Lecture 2). Consider a DLN of length with input and hidden layers all of width :
with , i.e. (and input dimension ). Conditional on , the vector is Gaussian , so its norm satisfies the recursion
(Recall that is the distribution of the sum of squares of i.i.d. standard Gaussian variables.) Let's normalize the input so that . Then the recursion implies that
with each factor of independent. At large , we have , so at large with fixed
up to higher corrections. Thus, as grows, the expectation of but the variance and higher cumulants blow up.
The lesson once again is
- , : Gaussian process, no correlations across layers (and, as we'll see in Section 4, no "feature learning").
- , small: controlled correlations, trainable network (and, as we'll see, feature learning during training). In open-weight LLM's today, –; see Exercise 3.3.
- , large: chaotic and untrainable.
Exercise 3.3 (★) (Depth-to-width ratios in the wild). Look up the depth (number of transformer blocks) and the width (model/embedding dimension) for a few families of open-weight language models — e.g. GPT-2, the Llama family, the Qwen family — and compute 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?
Solution
Representative numbers at the time of writing (take the width to be the model/embedding dimension): GPT-2 small has and GPT-2 XL ; Llama-2-7B has ; Llama-2-70B, ; Llama-3-405B, . The striking finding is how stable the ratio is: modern open-weight families cluster around – across two orders of magnitude in parameter count, comfortably inside the "controlled correlations" window, with the earliest (GPT-2-era) models running noticeably deeper for their width. Whether this stability reflects tuned wisdom or copied convention is a fair discussion question.
3.7 Edge of chaos
We briefly remark that the analysis of how to tune initialization variances (etc.) leading to stable network behavior over a large number of layers (the first point in Section 3.6) 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; Schoenholz et al. 2017). The idea was to treat the progression through layers of a network as a dynamical system, and to ask: given two nearby inputs , for what range of hyperparameters (, 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), (3.30) 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).
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 . Moving slightly away from this point, one can estimate the "critical depth" at which the network either collapses or blows up – see (Schoenholz et al. 2017). For MLP's with other activation functions, the edge of chaos is more interesting. For , the phase diagram is shown in Figure 5.
The quantity `' that controls the critical depth () 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 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, and their behavior at large width. The lesson will be that as width these networks effectively behave like linear models: their parameters never stray far from their values 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), 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.)
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.
We also emphasize the connection between Section 3 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). This involves two approximations. First we update on the entire training dataset at each step, rather than on batches, so that SGD becomes GD:
where
is the full empirical loss.
Second, we scale down the learning rate and scale up the number of training steps, setting where is a continuous "time" variable. As , the gradient descent converges to a gradient flow
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, Chapter 19).)
4.1.1 Geometry of gradient flows
Suppose that is a smooth Riemannian manifold, with nondegenerate metric , inverse metric , and a function . Then one can define gradient flow with respect to 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 ). The inverse-metric is used to relate the differential on the RHS to a vector field , which is the right sort of geometric object to induce a flow on .
In neural networks, is the space of parameters and is the loss function. One always assumes a constant Euclidean metric in parameter space , 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) with the metric looks like
Now return to the general setting of Riemannian manifolds. Suppose we have a smooth map between two manifolds, and moreover that is the pull-back of a function , meaning . Let be the components of the map , in local coordinates on . Then the inverse-metric on induces an inverse-metric on (called its push-forward):
Schematically, , where is the Jacobian matrix of . Then the gradient flow on induces a gradient flow on , given by
Let's apply this to neural networks, where is the function from parameters to network outputs on the training data , i.e. in our standard notation. The loss function is indeed a pullback of a function on , since it depends on parameters only through the network values on training data,
The push-forward of the gradient flow in parameter space becomes
The inverse-metric on the space of function values is called the "neural tangent kernel" or NTK. The name was coined in (Jacot et al. 2018); it also appeared in (Du et al. 2019) 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 be a standard Euclidean metric in the space of function values, with inverse . The gradient of the loss that appears in (4.10) (contracted with an inverse metric) is sometimes called the residual, and denoted or . For MSE loss, we've got . Then the residual is just the difference between function values and truth for the training data:
Noting that 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 matrix-valued connection . If one knows the dependence of on time (which in general is far from trivial — see the exercise below) the solution is a time-ordered exponential
Exercise 4.1 (The NTK from the chain rule). Derive the flow equation (4.10) for the function values directly, using only the chain rule — no geometry required. (The geometric derivation explains why an object like had to appear; the chain rule shows it appears in three lines.)
Solution
By the chain rule and the gradient-flow equation,
using in the middle step that depends on only through the function values .
Exercise 4.2 (How fast does the kernel move?). Show that the time derivative of the NTK satisfies
where is the Hessian. Bounding the Hessian ultimately leads to the constant-NTK theorem in Section 4.3 below.
Solution
Differentiate along the flow: . Substituting gives the stated formula. Taking operator norms, : the kernel can only move as fast as the Hessian of the network function (not of the loss) allows — which is why bounding that Hessian at large width freezes the kernel.
4.2 Features
The Jacobian matrix 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
The right singular vectors (columns of ) define natural modes in parameter space, while the left singular vectors (columns of ) are the natural modes in function-value space — where "natural" means with respect to how network outputs depend on weights. The left singular vectors (together with their singular values ) 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
with eigenvectors = the feature vectors, and eigenvalues = . 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 . Its empirical version is
the pull-back (via the network function) of constant Euclidean metric on the space of function values on the training data. In coordinates, .
The Fisher metric has eigenvalues , and is automatically degenerate when (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 , and fixed training data size , and take the width of all hidden layers to infinity. We'll assume the hidden widths all scale in a similar way , and we initialize the network the standard way from Section 3.2. We want to show that for almost all samples of training data,
- The NTK remains constant throughout training, up to corrections that vanish as .
- The network itself linearizes: at any training time ,
A rigorous derivation is somewhat technical, and we refer the reader to (Du et al. 2019; Jacot et al. 2018), (Roberts et al. 2022, Ch. 11), (Spiliopoulos et al. 2025, 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, Lecture 3).
We'll simplify life by setting biases to zero, so the network becomes
in the fan-in normalization of Section 3.2, 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- weights is an outer product,
where the backpropagated sensitivities obey and
The NTK is a sum over layers,
The second bracket is an average of terms. At initialization and large width, the law of large numbers replaces it by , the Gaussian-process kernel of (3.30) evaluated on the two inputs, which is of order one. (For the bracket is exactly , which is .) The first bracket is handled by the same recursion run backwards. Define
which is of order one for reasonable activation functions. Using the independence at initialization of the weights from everything below them, one finds with , and hence
This is a sum of terms, each of order one, with no factors of 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) is the standard recursion for the infinite-width NTK (Jacot et al. 2018).) By (4.10), an order-one learning rate then moves the function values by an order-one amount per unit time, which is what we want, and no rescaling of with width is needed in this convention. For the one-hidden-layer network one can check this by hand:
an average of order-one terms. (Here and below we absorb the first-layer factor into the inputs , since is fixed.)
How to scale the learning rate. The rule just used is general, and we will use it again in Section 5. What must stay of order one as is the motion of the network values, , since those are what fit the data. So: work out how the NTK scales with , and choose so that . In the fan-in normalization and . In the mean-field normalization of Section 5, and must grow like .
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 as in (3.5). Gradient flow with respect to the Euclidean metric is not the same flow as gradient flow with respect to . Under the change of variables, the inverse metric of Section 4.1.1 picks up a factor in every layer- block, so that in the original coordinates the flow reads . The NTK computed in the coordinates is correspondingly larger, in terms of the layer contributions of (4.21), and an order-one motion of then requires . 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 in the convention, for instance, the hidden layers move as in the fan-in normalization, but the first layer, whose fan-in does not grow with , is effectively frozen. When comparing learning rates across papers, first translate the parametrizations into one another; for a systematic treatment, see (Yang & Hu 2021).
Caution: the estimate (4.23) is at fixed depth and fixed data size . The NTK grows like with depth, and the operator norm of the matrix grows with , so for deep networks or large datasets the learning rate must shrink accordingly, and the constant-NTK statements below need and to stay small compared with .
The Hessian is small. To see that the NTK is approximately constant during training, we consider its time derivative, which by Exercise 4.2 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 () with one-dimensional output, following (Hanin 2026, Lecture 3). Writing
the Hessian splits into blocks, one per hidden neuron,
The operator norm of a block-diagonal matrix is the largest operator norm of its blocks. Each block has fixed size , contains only order-one quantities (at least at initialization), and carries an explicit factor . 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 ; the crude is enough for the argument.
Putting the pieces together. First, the parameters do not move far in total. From ,
where . So by the mean value theorem, for some , and
An order-one total displacement, shared among parameters, means that each individual weight moves only by : relative to its order-one initial value, every weight barely changes. Second, the kernel freezes. By (4.25) and (4.28),
so the order-one kernel becomes constant as . Third, the network linearizes. By Taylor's theorem with remainder,
for some on the segment between and . Setting and using (4.30) and (4.28), the quadratic term is bounded by
so the linear approximation becomes exact as , even though the linear term, and with it the network function, moves by an order-one amount. In words: each of the neurons has leverage on the output, so an order-one change in can be assembled from changes in every weight, and changes that small leave the features 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 and eigenvectors of the NTK are (by construction) constant. Nevertheless, the network function converges to its training data exponentially fast: from (4.12) with constant , we have
where we recall that . Writing the residuals in an eigenbasis for the NTK, we see that
Each feature is learned at a characteristic timescale .
Exercise 4.3 (Lazy predictions on test points). Derive that on an arbitrary test point , not necessarily part of the training data,
where is the NTK for a test-train pair (which remains constant during training, just like the train-train NTK ).
Solution
From Exercise 4.3, we see that at the end of training the network obeys
This is an affine-linear transformation of the Gaussian process at initialization, so it is again a Gaussian process. It's straightforward to check that the new mean and covariance are
where , , etc.
4.5 Summary
So far we've seen that in the fan-in normalization (Section 3.2), with order-one weights and an order-one learning rate , at large width we've got
which forces to be constant and the network to linearize as a function of its parameters.
All this breaks down if
- Depth is large as well, with nonzero. (See (Roberts et al. 2022, Ch. 11).)
- Data samples become large, with nonzero.
- Learning rate becomes large.
Both and 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, trains one small network in both the fan-in and the mean-field normalization and tests the claims of Sections 4 and 5 numerically. The setup runs in seconds, in about sixty lines of NumPy; reproduce it, then modify it.
The starter model. Inputs are one-dimensional: equispaced training points , with targets from a two-neuron teacher, . The student is a one-hidden-layer network
trained by full-batch gradient descent on , with residuals . Here take the fan-in exponent , learning rate , and steps. The empirical NTK to monitor is the 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 seeds, histogram at a fixed input, say . You should see an order-one Gaussian, the Gaussian process of Section 3.4, and the sample covariance across two inputs should match the GP kernel.
(c) The kernel freezes. For (a few seeds each), record and plot it against on log–log axes. Expected: it decays, with slope between (the guaranteed bound (4.31)) and (the typical measurement). Also record the relative parameter displacement and compare with (4.30).
(d) Mode by mode. Diagonalize and plot on a semilog axis for a few : straight lines with slopes , the features being learned at , largest eigenvalue first, as in (4.35). Compare with the Bayesian version, Exercise 4.5(b).
(e) Prediction. Compare the trained function on fresh test points with the frozen-kernel prediction (4.36).
Solution
Reference results for the setup as specified (single seeds; small- numbers fluctuate).
(a) Final loss –. Keep for stability; here , so is safe.
(b) At , , and the sample covariance across two inputs matches the GP kernel.
(c) Measured at : decaying, with a log–log slope between and ; average a few seeds to clean up the trend. The relative parameter displacement decays like , since the total displacement is by (4.30) while . Compute from per-example gradients (an Jacobian contracted with itself), in float64.
(d) The residual modes are straight lines on a semilog plot, with slopes , until they hit the small nonlinear floor; the smallest-eigenvalue modes are visibly not yet learned at the end of training. This is spectral bias on display, and the dynamical twin of Exercise 4.5(b).
(e) At the trained function on test points agrees with (4.36) to within a few percent; the agreement degrades as shrinks, which is feature learning beginning to stir. The continuation is Exercise 5.3.
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).
Suppose that we consider the same large-width network, but "train" Bayesian inference instead, on the empirical MSE loss . (Recall from Section 2.1 that, practically speaking, this means training via Langevin dynamics.) Then we expect a Bayesian posterior
where the prior is the Gaussian of order-one variance from (3.1), and the network carries the fan-in normalization of Section 3.2.
We'd like to push the distribution forward to the network function . At large-width, the push-forward of the prior is our familiar Gaussian process from Section 3.4; expressed via a functional integral as in Section 3.5 it becomes
with the final-layer GP kernel. Since the loss is explicitly a function of network values, we immediately get the late-time trained distribution
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: the NTK froze, and the -th kernel eigenmode of the residual decayed as , so that it is learned at the characteristic time . Here you derive the Bayesian-learning version of the same statement, with in the role of training time. Restricted to the training values , the prior is Gaussian with covariance the kernel matrix , so the posterior on training values is
(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 .
(b) Eigenmodes. Diagonalize the kernel, , and show that mode by mode
Conclude that the eigenmodes of the GP kernel at initialization are learned one by one as grows, mode switching on at the characteristic scale , large-eigenvalue modes first (spectral bias).
(c) Compare with gradient flow. Place this alongside (4.35), , that is . What is the same, and what differs? Two things to notice: the correspondence between the amount of data (or inverse noise) and training time; and the different shapes of the two filters, the power law in versus the exponential in (compare Exercise B.5). Also: the Bayesian calculation uses the GP kernel , the dynamical one the NTK , two different frozen kernels for the same network.
(d) Prediction. Show that the full posterior over functions is again a Gaussian process,
with mean and covariance
The posterior mean is kernel ridge regression with ridge (Appendix B.4), and interpolates the data exactly. [Hint: the exponent is quadratic in ; complete the square, and use the Woodbury identity to relate the two natural forms of the answer.]
Solution
(a) Completing the square, the exponent is with , so the posterior is Gaussian with covariance and mean (multiply through by ).
(b) In the eigenbasis all matrices are diagonal, giving (4.46). The filter is for and for : a switch at . On a axis it is a sigmoid of universal shape, pleasant to compare with the logistic occupations of Section 2.5, though nothing here is a phase transition: the switches are smooth crossovers, one per mode.
(c) The same: frozen-kernel linear structure, mode-by-mode learning, spectral bias with scale , and the dictionary " training time" of Section 2.1 made quantitative as . Different: the suppression of the unlearned part is , a power law in , versus , exponential in (the two filters of Figure 13); and the kernels differ. Bayesian learning of the frozen network is governed by the GP kernel , a property of the function distribution at initialization, while gradient flow is governed by the NTK , a property of the parametrization. Both are fixed, data-independent kernels: in neither picture does the lazy network adapt its features.
(d) Varying the exponent of with respect to gives the stationarity equation , that is . Evaluating at the training points gives a closed linear system whose solution is (a), and substituting back yields the stated . The covariance is the inverse of the quadratic form, , which the Woodbury identity converts into the stated form. As this pins : noiseless GP interpolation; finite is ridge regularization with ridge .
These equations embody what's sometimes known as kernel regression, with ridge parameter . (Appendix B.4 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
Writing the data "truths" in an eigenbasis for , with ordered eigenvalues , we find
This is the Bayesian version of lazy learning. For given , the "modes" of the data distribution with are learned by the network, and the rest are projected out.
Recall that in the Bayesian context, the rescaled inverse-temperature 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 is (see (4.35)). Correspondingly, in Bayesian lazy learning, the minimum inverse-temperature required to learn a GP mode is .
5. Mean-field scaling
Sections 3 and 4 analyzed one and the same large-width limit: the fan-in normalization (Section 3.2) with order-one weights and an order-one learning rate, and width 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) 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 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 behavior is governed by a self-consistent "mean field." We'll first explain what physicists mean by a mean-field theory (Section 5.1), then run the analogy in detail for one-hidden-layer networks (Section 5.2), see how the NTK story of Section 4 changes (Section 5.3). Finally we'll discuss some practical applications, including interpolating scalings, P, and what this framework is actually used for (Sections 5.4–5.5).
The mean-field limit of one-hidden-layer networks was constructed independently, and essentially simultaneously, by four groups (Mei et al. 2018; Chizat & Bach 2018; Rotskoff & Vanden-Eijnden 2022; Sirignano & Spiliopoulos 2020); our treatment draws on all of them. Pedagogical accounts may be found in (Spiliopoulos et al. 2025, Ch. 20) and (Hanin 2026, 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, Chapter 1). There is, however, a family of models — those with weak, all-to-all couplings — for which it becomes exact as ; 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 on the fully connected graph: every spin couples to every other, with the coupling scaled down to so that the total force on any one spin, and the total energy, stay of order and respectively. Setting , the energy (2.19) becomes
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) of the lattice model, with . Everything that was approximate in Section 2.4 is therefore exact here. The partition function reduces to the one-dimensional integral (2.25) with the same entropy , the saddle-point equation is
and the model has the second-order transition of Figure 4 at and the first-order transition of Figure 2 in , exactly as analyzed in Section 2.5.
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 : spin does not care about any individual other spin, only about the empirical average , the "mean field." The classic argument then has three features:
- Each degree of freedom couples to the others only through an average. (True by construction here.)
- The average barely fluctuates. If the spins are roughly independent, the central limit theorem makes the fluctuations of of size ; as each spin sees a deterministic effective field .
- The average is fixed self-consistently. A single spin in the field has ; since was defined as the average spin, consistency requires (5.2).
For a finite-dimensional lattice, where each spin has a few strongly coupled neighbors, steps 1–2 are approximations. For the all-to-all model they become exact as : the mean field is an average of 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) numerically by the fixed-point iteration , for below and above and for a few values of , and watch the iteration converge. Which fixed points are stable? Keep this iteration in mind: we will meet its infinite-dimensional cousin below.
Solution
The iteration converges geometrically to a stable fixed point, at rate set by the slope of the right-hand side there. For the only fixed point is and it is stable. For the fixed point has slope and is unstable: any seed away from it flows to one of the two nonzero fixed points , which are stable, and a small selects the one with the sign of (the other survives as a stable fixed point too, the metastable state of Section 2.5). Its infinite-dimensional cousin is the fixed-point iteration for the McKean–Vlasov equation in Section 5.2.
5.1.2 The general shape of a mean-field model
Nothing in the story above needed binary spins. More abstractly, one can take "particles" with an extensive energy of the form
with a single-particle potential and a symmetric pair interaction , again with the characteristic coupling. The force felt by particle is . More so, for particles dominated by drag (highly viscous motion), the dynamics at zero temperature are governed by familiar gradient-descent equations
where differentiates the first slot of . 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 recording where the collection of particles sit. Then (5.4) 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 — the exact analogue of "every spin feels only ." The same three-step argument as before then suggests that as 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) of Sections 3.5 and 4.3, but now normalized with out front instead of , with the same order-one weights:
with , (as in (4.26), the first-layer factor has been absorbed into ). Compared with the fan-in normalization of Sections 3–4, mean-field scaling suppresses the output by one extra factor of . 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 and in the callout of Section 4.3, 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 particles are going to be the hidden neurons, with parameters
Their values are encoded in the empirical measure (5.5), and the network (5.7) depends on its parameters only through ,
We think of the network value as a linear functional of . Neurons are viewed as exchangeable particles, and the measure is their state.
5.2.1 Initialization
At initialization, the summands in (5.7) are i.i.d. with mean zero (since is independent of and ). With the prefactor, the law of large numbers gives, for each fixed ,
Contrast this with Section 3.4: with the 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 , as in Section 4:
Inserting (5.7) and expanding the square, we get
with
Comparing (5.12) with (5.3), we see that up to an overall factor of , the loss is precisely a mean-field energy. In particular, has the canonical normalization from (5.3) 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! rewards neurons for aligning their output function with the data, playing the role of the external field ; while 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 in this parametrization, raw gradients of are , and gradient flow at fixed learning rate would freeze as . We must therefore scale the learning rate up,
with fixed. (Equivalently, we could run vanilla gradient flow on the extensive energy .) This is the rule of Section 4.3 at work: the NTK, computed below in (5.24), is , so keeping of order one requires . Note the contrast with Section 4.3: there the NTK was of order one and the learning rate needed no rescaling; here the NTK shrinks like 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, with ; but then the same dynamics (5.15) requires a learning rate for the readout weights and for the hidden weights — one learning rate per layer, differing by a factor . Keeping all weights of order one is what lets a single do the job, and it is the reason we adopted that convention in Section 3.2.
With (5.14), the flow becomes, per neuron,
where is the residual of Section 4.1. Equivalently, in the potential form (5.6),
with .
As anticipated, the velocity of neuron depends on its own state , and on all the other neurons only through the average , via the residual . The network function is the mean field here! We arrive at a dictionary:
| Curie–Weiss | one-hidden-layer network, mean-field scaling |
|---|---|
| spin | neuron |
| magnetization | network function |
| external field | data term (the targets ) |
| pair coupling | neuron-neuron interaction |
| effective field | residual |
| extensive energy | (the loss is intensive) |
| zero-temperature dynamics | gradient flow |
| Langevin at temperature | noisy training, Section 2.1 |
| self-consistent Gibbs measure, Equation 5.22 below |
Exercise 5.2 (Structure of the mean-field equations). (a) Derive (5.12)–(5.13) and (5.15), and check that (5.16) reproduces them, using the symmetry of .
(b) (Balancedness.) For a positively homogeneous activation such as ReLU, . Show that in this case (5.15) conserves 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 in (5.15) does not depend on (only on and the mean field). Where does the linearity of in enter?
Solution
(a) Expanding the square in the loss and collecting single-sum and double-sum terms gives (5.12)–(5.13) directly. For the equations of motion, and ; the chain rule plus gives (5.15). Writing the same gradient through , the symmetry ensures the double sum contributes twice the one-sided derivative, matching (5.16).
(b) From (5.15): and . For positively homogeneous , , so the two right-hand sides coincide and .
(c) Because is linear in the output weights, contains no ; the only route by which influences is through the mean field — i.e. through the population average, weighted by . This linearity is also what makes the loss a quadratic (rather than higher-order) polynomial in each particle, i.e. what limits the interaction to be pairwise.
5.2.4 From particles to measures: the continuity equation
We now pass to the limit, in distinct steps.
Step 1: rewrite the coupled ODEs as one equation for the empirical measure. (This remains exact at finite .) Suppose particles move deterministically with some velocity field, , and let denote the law of one such particle.
The expectation of a smooth test function obeys
integrating by parts in the last step. Since was arbitrary,
This is known as the continuity equation, and it simply expresses conservation of probability. Applied to our system, with from (5.16), the empirical measure satisfies (weakly, since it is a sum of delta functions)
An important consequence of the mean-field (, all-to-all) coupling here is that (5.19) is a closed equation in .
Step 2: the limit. At finite , is a random, atomic measure. (Atomic means it's a sum of delta functions; it's random due to random initialization.) As with i.i.d. initialization , the law of large numbers suggests — and the theorems of (Mei et al. 2018; Chizat & Bach 2018; Rotskoff & Vanden-Eijnden 2022; Sirignano & Spiliopoulos 2020) confirm — that converges (weakly) to , a deterministic curve of smooth probability measures solving the limiting PDE
There are quantitative 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), namely , is called the McKean–Vlasov equation. To get it, focus on a single representative neuron . As , we can replace the force on it by the ensemble average:
The constraint 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) a nonlinear Markov process, nonlinear not in but in the law.) The McKean–Vlasov equation is equivalent to the continuity equation; the particle trajectories of (5.21) are the characteristic curves of the PDE (5.20).
We remark on a few structural features:
- Propagation of chaos. Since neurons interact only through a quantity (, equivalently ) that becomes deterministic as , they become asymptotically independent: any fixed collection of neurons converges to i.i.d. copies of (5.21), with cross-neuron correlations of size . Thus, independence at initialization propagates instead of degrading (Sznitman 1991).
- Gradient-flow structure. Define the loss as a functional of the measure, . Since is linear in , this functional is convex in .
This ultimately leads to global-convergence theorems. For example, it was shown in (Chizat & Bach 2018) that for homogeneous activations (e.g. ReLU), if training converges at all, it converges to a global minimizer of the convex functional .
- Temperature. If one trains instead with the Langevin dynamics of Section 2.1, Eqn. (2.6), at temperature , one finds convergence to a stationary distribution that satisfies (Mei et al. 2018)
This is another version of the "" consistency condition, now for noisy training.
Mean-field training for a example: , trained by the rescaled gradient flow (5.15) on samples of the two-neuron teacher ; is rescaled training time. Left: the empirical measure in the plane for (dots, colored by time; thin gray lines trace of the particle trajectories), drawn over the density of an run at (shaded contours), standing in for the limiting measure . Individual neurons move order-one distances, and the measure condenses onto one-dimensional filaments — the configurations realizing the teacher. (The point symmetry reflects the oddness of .) Right: the mean field itself, , converging to the target (dashed; circles mark the training data).
How would one actually solve the self-consistent system (5.21)? The same way one solves (5.2) (Exercise 5.1): freeze a guessed environment , 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). Alternatively, solve the PDE (5.20) on a grid, which is feasible (only) for small .
In Figure 6 we instead solve the system by simply training a network with . 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 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: same data, same teacher , same student (4.40), but now with the mean-field exponent and the rescaled updates (5.15),
with for steps, so that the rescaled time runs to . Note that there is no in these updates: that is the required learning rate (5.14), already built in (). Monitor the NTK (4.41) with .
(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 seeds, histogram at for several : the distribution collapses to with width . Check the scaling.
(c) The kernel moves. For , record and plot it together with the lazy curve of Exercise 4.4(c): one plot, two curves, the whole dichotomy.
(d) Watch the particles. Scatter-plot the neuron cloud at for both runs. Lazy: the cloud is static (each weight moves ). Mean-field: the cloud transports and condenses onto thin filaments, the configurations realizing the teacher, as in Figure 6. Repeat at : 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 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: ; ; a ReLU student on the teacher (and vice versa); or move to with and plot the cloud of directions . Which aspects of (c)–(d) are universal, and which were special to the aligned case?
Solution
Reference results (single seeds).
(a) Final loss at and still slowly falling; the mean-field run is measured in rescaled time, so run longer to go lower.
(b) At , , with standard deviation falling like across : the Gaussian process of Section 3 survives only as the fluctuation around zero.
(c) Measured at : order one, with no decay, against the decaying lazy curve. Use the same formula with the appropriate in both runs.
(d) The mean-field cloud reproduces Figure 6: neurons stream along order-one trajectories and pile onto one-dimensional filaments through (teacher directions, matched output weights), with the point symmetry from the oddness of . The lazy cloud is visually identical at and (displacements ).
(e) For or the loss still falls (more slowly), and the cloud condenses onto different, task-shaped structures, with more filaments spread over a range of slopes , since no finite set of teacher neurons exists to copy; the kernel still moves by . That is the cleaner demonstration that mean-field training builds features rather than merely locating the teacher's. Mismatched activations behave similarly. The lazy runs, by contrast, change in only one way, the quality of the fit (the random features either happen to span the target well or do not); the kernel and the cloud never adapt.
5.3 The NTK is no longer constant: feature learning
Section 4.3 showed that in the fan-in normalization the NTK freezes as , and Section 4.2 defined feature learning as change of the NTK during training. Let's redo that computation in mean-field scaling. From (5.7), and , so the NTK of Section 4.1 is
with the per-neuron kernel
Thus , recovering the previous observation that gradients are small, which was the reason for the rescaling (5.14). The estimates of Section 4.3 also explain why the kernel is now free to move. The Hessian blocks (4.27) carry instead of , so ; but the learning rate is , so the bound (4.25) gives for the relative drift . The bound no longer forces the kernel to freeze, and the equations of motion (5.15) will show that it does move. As , the meaningful object is , and the function-space dynamics (4.10) now becomes a kernel gradient flow with a time-dependent kernel:
At initialization the LLN makes , a deterministic kernel — the same limiting kernel (4.24) as in the fan-in normalization. The two limits start from the same place.
During training, however, they part ways. In Section 4.3, each parameter moved only , so stayed within of and the kernel was pinned, (4.31). Here, the equations of motion (5.15) move each particle an order-one distance. (Each neuron has only leverage, so fitting the data forces the average neuron to move a distance .) Thus transports macroscopically, and by (5.26) 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 . For each , as , the learning rate must be simultaneously be scaled as so that the function-space dynamics stays (one can derive that from ). Repeating the estimates above gives (see Exercise 5.4):
| at init | per-neuron motion | limit | |
|---|---|---|---|
| (GP, CLT) | lazy: Sections 3 and 4 | ||
| lazy (frozen kernel, zero init) | |||
| mean-field: this section |
Notably, every intermediate 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; Mei et al. 2018). Mean-field, , 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, P, in Yang and Hu's line of work on tensor programs (Yang & Hu 2021). In particular, mean-field scaling is P for .
Exercise 5.4 (Lazy vs. rich along the dial). For the family (5.27):
(a) show and hence that is the right co-scaling;
(b) show the initialization output is and per-neuron velocities are ;
(c) conclude that the (normalized) kernel drift over any fixed training horizon is , vanishing for every ;
(d) reconcile with the bounds of Section 4.3. Optional: relate to the output-scale parameter of (Chizat et al. 2019), where multiplying any model by produces lazy training.
Solution
(a) and similarly for , so ; demanding forces .
(b) At initialization is times a sum of i.i.d. mean-zero terms , so . Per-neuron velocities: , and likewise for .
(c) The normalized kernel is an empirical average over particles, so over any fixed horizon its drift is of the order of the per-particle displacement, — vanishing for every , order one exactly at .
(d) At : per-neuron motion , matching (4.30) (an total displacement shared among neurons); initialization output , the GP of Section 3. For the optional part: relative to the mean-field parametrization (5.7), with , so in the language of (Chizat et al. 2019) every sends the laziness with width — lazy — while keeps of order one.
5.5 What is this framework good for?
Nobody solves the McKean–Vlasov PDE (5.20) 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 (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 . Its deep-network generalization (Yang & Hu 2021) turns this into an engineering recipe, Transfer (Yang et al. 2022): 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 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). Diagnostics from this circle of ideas (kernel distance/alignment during training, e.g. the silent-alignment effect (Atanasov et al. 2022)) 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 , 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 (Bordelon & Pehlevan 2022) 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).
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). Grokking refers to the following behavior. One trains a small network on an algorithmic dataset — the canonical example is addition mod , presented as pairs with label , with only a fraction of all 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 . 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) (this subsection) — grokking as a delayed transition from lazy to rich dynamics — and the equilibrium lens of Rubin–Seroussi–Ringel (Rubin et al. 2024) (Section 5.7) — 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 and 5.4; 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) invokes no statistical-mechanics ensemble at all. There is only gradient descent, and the lazy/rich dial of Section 5.4, 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), 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 and 5), 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. Then work through Exercise 5.5.
A few translations of terminology from (Kumar et al. 2024):
- Lazy / linearized dynamics: the regime of Section 4.3 — the network moves only in the affine subspace (4.1) spanned by its initial gradients, so training is kernel regression with the frozen initial NTK (Sections 4.1 and 4.6).
- Laziness parameter : an overall multiplier on the network output, (with loss and learning rate co-scaled). This is the continuous version of the -dial of Section 5.4, in the form introduced by (Chizat et al. 2019): is fully lazy, small is rich.
- Kernel–task (mis)alignment: the language of Section 4.2 — 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 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.
Exercise 5.5 (Grokking as a lazy-to-rich transition). Read (Kumar et al. 2024) 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, 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 implemented, and how does it relate to the -dial of Section 5.4? Section 5.2'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).)
(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 that plays this role in the other walkthrough (Exercise 5.6(b)): 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, 4.3 and 5.4).
(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 and small? Contrast with the control parameters of Exercise 5.6(d).
(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 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.
Solution
(a) Finite network, deterministic gradient descent/flow, no limit required: the analysis is organized around the linearization of Section 4.3 and departures from it. The dial is an output multiplier (with appropriate co-scaling of the loss/learning rate), exactly the scale parameter of lazy training (Chizat et al. 2019): at large an change of requires only an change of weights, so features barely move — a continuous, finite- version of the -dial, where played the role of (Exercise 5.4). And this paper is precisely about the half of the scaling story that the Bayesian treatment of Exercise 5.6(a) could not see: there, only the induced function prior survived and the learning-rate/dynamics half dropped out; here there is no prior at all, and everything is the dynamics half — the two walkthroughs literally partition Section 5.2's prescription between them.
(b) The alignment between the (current) NTK and the task — equivalently the kernel's distance from its initialization, or the paper's low-dimensional sufficient statistics: in the memorized state the kernel is still (close to) its initialization and the fit lives in misaligned modes; in the generalizing state the kernel's top eigenspaces have rotated onto the target. Measurable in any run by tracking the empirical NTK (Section 4.2). The Fourier overlap of Exercise 5.6(b) carries the same information — features aligned with the task's structure — read off the weights and their posterior rather than off the kernel; for modular addition, a rotated kernel and condensed Fourier modes are two descriptions of the same event.
(c) (i) The top eigenvectors of the initial NTK are misaligned with the labels: in the language of Section 4.6, the target's weight sits in modes with small eigenvalues, so the lazy (kernel-regression) solution fits the training points without generalizing — memorization, in the precise mode-filtering sense. (ii) The dataset is in a Goldilocks window: large enough that a generalizing solution exists (cf. the critical sample size of Exercise 5.6(d)), small enough that train and test loss do not simply track each other. (iii) Training starts lazy (large , or large initialization scale): feature learning is initially suppressed, as in Section 4.3, so the kernel fit happens first.
(d) The laziness (and the initialization scale, which acts similarly): increasing it stretches the delay between the lazy fit and feature learning, making grokking more dramatic; decreasing it makes the network rich from the start, so train and test fall together and grokking vanishes. In the extreme the network is a pure kernel machine forever: by condition (i) it never generalizes, and the grokking time diverges. Contrast Exercise 5.6(d): there the knobs were sample size and temperature moving the system across a phase boundary; here the knob does not change where the system ends up so much as how long the misaligned detour lasts.
(e) From a separation of timescales built into the flow: fitting the train set within the frozen-kernel subspace is fast (the exponential convergence of Section 4, eq. (4.34)), while rotating the kernel is slow — suppressed by the laziness, since feature motion is -small per unit of function change. Train loss bottoms out on the fast clock; test loss must wait for the slow one. Nothing stochastic is needed: no noise, no weight decay, no barrier crossing — and that is the force of the paper's counterexamples, since grokking with zero weight decay and growing weight norm rules out explanations in which regularization-driven norm shrinkage causes the transition. The equilibrium mechanism of Section 5.7, by contrast, needs temperature: its delay is the escape from a metastable free-energy basin, which pure gradient flow would never make.
(f) They live in different corners and need not compete: Rubin et al. describe noisy (Langevin/Bayesian) training near equilibrium, Kumar et al. deterministic training far from it; both say the cheap channel is used first — kernel fit or metastable memorizing phase — and features arrive late. Distinguishing experiments:
- noise dependence — an activated (metastable-escape) delay should shorten rapidly with increasing training noise, roughly Arrhenius-like, while a timescale-separation delay survives at zero noise and depends on it only weakly; measure grokking time vs. Langevin temperature at fixed , .
- Seed-to-seed statistics — nucleation out of a metastable state is a rare-event process with broad (roughly exponential) waiting-time distributions across seeds; a deterministic mechanism predicts narrow, initialization-scale-set delays.
- The -dial itself — the lazy-to-rich account makes a sharp prediction that the delay is monotonically stretched by with grokking eliminated at small ; the equilibrium account ties the delay to the free-energy landscape at fixed ensemble, not to . Running these on modular addition is a very feasible project.
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: 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. The paper (Rubin et al. 2024) 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), GP kernels and kernel regression (Sections 3.3 and 4.6), mean-field scaling and self-consistency (this section), and the anatomy of first-order transitions (Section 2.5). It is, however, a genuinely dense paper, written in an idiosyncratic dialect of the "adaptive kernel" school (Seroussi et al. 2023). This subsection is therefore structured as a reading exercise in the style of Section 1: 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.
Warm-up. Before any statistical mechanics, it helps to know concretely what the network is trying to become. Gromov (Gromov 2023) writes down exact weights for a two-layer network with quadratic activation that solve modular addition: the first-layer weights are plane waves on — entries , at various frequencies — and the quadratic activation's product-to-sum identities convert them into the required functions of . Skim (Gromov 2023) 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 verbatim: noisy gradient descent with weight decay; at late times, the Gibbs posterior; the weight decay sets the Gaussian prior and the noise level sets the temperature.
- (Adaptive) kernel: the GP kernels of Section 3.3, 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).
- Equivalent kernel (EK) approximation: the dense-data approximation of Section 4.6, 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 , and as ; how this relates to the of (5.7) is Exercise 5.6(a).
- 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.
Exercise 5.6 (Grokking as a first-order phase transition). Read (Rubin et al. 2024) 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, 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) and to the -dial of Section 5.4? Careful — this is subtler than it looks: in an equilibrium (Bayesian) framework there is no learning rate, so which part of Section 5.2'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 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.
Solution
(a) A Bayesian two-layer network whose readout prior (set by the output-layer weight decay) has variance ; the limit taken is together with — the paper's "mean-field scaling" — plus the EK (dense-data) approximation. The subtlety: at equilibrium, only the induced prior on functions matters. For the linear readout, a prefactor in and a prior variance for enter the function-space prior only through the product (a Gaussian change of variables). So the mean-field parametrization (5.7) — prefactor , order-one weights — corresponds to bare readout variance , while variance at fixed would be precisely the standard GP scaling of Section 3. The mean-field regime is the additional suppression . Meanwhile the other half of Section 5.2's prescription — the learning-rate rescaling (5.14) — has no equilibrium counterpart at all: the learning rate drops out of the stationary distribution (Section 2.1). In the language of the -dial: gradient-flow scalings are classified by (prefactor, initialization, learning rate) jointly, but their Bayesian shadows are classified by the single combination .
(b) The overlap of a hidden neuron's input weights with the teacher-relevant directions — for modular addition, the Fourier coefficients of on (the warm-up (Gromov 2023) says which ones matter). The sharper statement is about the single-neuron posterior : in the GFL phase it is a single Gaussian centered at ; in the mixed phase it becomes a mixture, with a finite fraction of neurons condensed at and the rest still Gaussian. The order parameter jumps discontinuously: first order. The ordered side is a finite-temperature cousin of Figure 6, where the measure over neurons condenses onto the low-dimensional structures realizing the target.
(c) (i) Noisy training Gibbs posterior (Section 2.1); (ii) mean-field decoupling of the posterior into readout hidden-layer factors, with neurons coupled only through averages (this section); (iii) the self-consistent "adaptive" kernel , whose defining equation has exactly the structure of (5.22) — the kernel depends on a posterior that depends on the kernel; (iv) saddle-point evaluation (Section 2.4); (v) a variational Gaussian-mixture ansatz, needed precisely because a first-order transition means several competing saddles; (vi) the EK approximation for the data average (Section 4.6). Note what is absent: in the main text the data averaging is done by EK smoothing, not by the replica trick of Appendix A.2.
(d) Primarily the sample size at fixed temperature: below a critical number of training facts the disordered (memorizing) phase is dominant and the network never generalizes; above it, the mixed phase wins. (The transition line lives in the plane, so it can equally be crossed by tuning the noise level or weight decay at fixed — but the grokking narrative crosses it in .)
(e) In equilibrium there is no time; the free-energy crossing itself is instantaneous. The time-dependence of grokking is the dynamical lag of a first-order transition: Langevin training equilibrates quickly within the metastable disordered basin — which is enough to drive the train loss to zero (memorization is easy there) while test accuracy stays at chance — and only much later escapes to the thermodynamically dominant mixed phase, whereupon test accuracy jumps. This is precisely the Ising-perceptron story of Section 2.6: the information needed to generalize is present in the data well before the dynamics finds the ordered state, so capability arrives "abruptly and later than it should."
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), together with limits in which training simplifies to kernel regression (Section 4) or to a deterministic flow of measures (Section 5). 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), and finally extract closed, self-consistent equations for a small set of two-time order parameters: correlation functions and response functions . 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 determined scalars self-consistently, the McKean–Vlasov limit of Section 5 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) — and a guided tour of the modern DMFT program for deep learning (Section 6.3), whose scaling-law application closes the notes in Section 7. Throughout, we follow the spirit of Altland and Simons (Altland & Simons 2010), 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) and (Stapmanns et al. 2020).
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 — differential equations become stochastic, and the difficulty of solving them jumps. A stochastic differential equation, at first order for example, takes the form
where is a (Gaussian) white noise process — additive noise, , as opposed to multiplicative noise — and is now a random trajectory. The noise is entirely defined by the first and second moments,
where the expectations are with respect to the noise ensemble3. Obtaining closed-form solutions for the time evolution of the probability distribution for 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), Gardiner (Gardiner 2004), and van Kampen (van Kampen 2007). The objects of interest here are instead the moments of and the probability density function 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) 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 , where so that,
where 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 , an expectation value of an observable is formally,
We used the shorthand for the constraint imposed at each step. The measure here, is a typical functional path integral measure with 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 ,
Finally, we can use the Fourier representation of the delta, introducing an auxiliary variable , which in the continuum is
The final step is the payoff: the noise now appears linearly in the exponent, so the Gaussian average over can be done exactly, , leaving the MSRJD functional integral
Four structural comments, each of which does real work later:
- The partition function is one. Setting , the integral is just of a normalized probability: , with no sources. Unlike in equilibrium, all the physics lives in correlation functions, obtained by coupling sources to both and and differentiating — compare the starting point of Section 6.2.
- The response field earns its name. Perturbing the drift, , changes the action by , so
Correlators of compute responses to kicks. Causality is built in: for (the Ito discretization is what guarantees this, along with ).
- Two-time observables. The natural objects of the formalism are the correlation function and the response . In equilibrium they are not independent: the fluctuation–dissipation theorem (FDT) locks them together, for (with in our normalization) (Kubo 1957; Kubo 1966). 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.
- Onsager–Machlup. The action (6.7) is quadratic in , so the response field can be integrated out exactly, giving the probability of a trajectory as : 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).
(b) For the linear drift (the Ornstein–Uhlenbeck process), compute the stationary correlation and response functions, either by Fourier transforming the quadratic action or by solving the SDE directly: and .
(c) Verify the fluctuation–dissipation theorem for , with .
(d) Explain, from the discrete-time construction, why and why is causal.
Solution
(a) The white-noise measure is Gaussian with , so , and collecting terms in the exponent gives (6.7).
(b) The linear SDE is solved by (stationary state). Then with . Adding a perturbation to the drift shifts the solution by , so .
(c) For : , i.e. with .
(d) In the discrete construction each delta function propagates information strictly forward in time (Ito: the noise at step affects , never earlier variables). Diagrammatically, is a retarded propagator; would be a closed loop of retarded propagators, which vanishes — equivalently, "a response to no perturbation is not an observable." Causality of is the same statement.
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), 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 of Section 6.1 — how a perturbation at time affects the model's state at time ; 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; Kubo 1966; Goldenfeld 2018), while in the second the FDT breaks — the standard warning light, flagged in Section 6.1, 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: trainable parameters, training samples, learning rate , teacher weights , student weights (with response fields ), and the overparametrization ratio
as control parameter — the reciprocal of Section 2's load , so that the over-parametrized phase is . For readers consulting the paper, the dictionary is:
| these notes | Li–Goldenfeld (Li & Goldenfeld 2026) |
|---|---|
| (parameters) | |
| (samples) | |
| , , | , , |
| (learning rate) |
The batch-inclusion probability , the ridge , and the label-noise variance are as in the paper.
Explicitly, the student and teacher share the same architecture
where is the number of model parameters, indexes the dataset with inputs , the are the teacher's weights, and is label noise. The thermodynamic limit takes and with the overparametrization ratio — the paper's complexity parameter — held fixed. The three main regimes are now definable as under-parametrized: , near-criticality: , and over-parametrized: . We note their optimization procedure, which is SGD with both regularization controlled by (set to zero for most of the work) and batching controlled by the binary indicator variable which is drawn for each point in the dataset with probability ,
From here they construct the MSRJD integral exactly as in Section 6.1; 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:
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,
which is a very common object to see in statistical field theory! In general 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 , the initial condition , the stochastic sampling , and the added label noise . For comedic relief we also write their full dynamical function integrated over quenched disorder:
where abbreviates the measure over the dynamical fields, and the new measure is
Each of the 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)
After taking the thermodynamic limit, the generalization error is expressed in terms of the equal time correlator determined from the DMFT equations and , the noise variance:
They find two forms of this equation and propose a third, a scaling relation near the critical point. At for , the asymptotic error is:
Note the pole as criticality is approached. The quantity 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 and still showing sensitivity to analogous to the regime near criticality:
where fully, . The term is the variance of the initialization and 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 , 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 the error is asymptotically finite, as the model is still attempting to approximate the solution, while for the training error vanishes as initialization-dependent memorization ensures the model always performs accurately on the training data. They write training error as
For it simplifies to , while again, for , the system is no longer at equilibrium and the susceptibility diverges, which drives the training error to zero.
To characterize the critical behavior of the model near , they calculate the critical exponents and propose a scaling function. The first part is relatively straightforward, where from Eqs. (6.18) and (6.19), one finds diverges from both sides . This suggests the scaling with an exponent 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 -dependent, a general feature of phase transitions is a dynamic scaling data collapse should be present near the transition (Goldenfeld 2018). They make this argument and propose a scaling ansatz of form:
where they introduce a new critical exponent and write the universal scaling function , and is still the probability of samples being included. They make two observations about the form of near and at the critical point. First, near the point, corresponds to , they argue the singularity in front must be absorbed into and it takes the asymptotic form . Second, at , since the generalization error grows in time instead of saturating, . The exponent therefore characterizes this growth in training time. Notably, the form of Eq. (6.21) implies a characteristic time scale, which they denote to describe convergence in the scaling of the error. Since the error grows, we would expect this to diverge as , 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).
Exercise 6.2 (Double descent in the lab). Simulate the model of this subsection: teacher , inputs , labels with noise variance ; student trained by the update (6.11) with ridge and full batches (), from . Take and .
(a) Plot training and generalization error against time, and verify the asymptotics (6.18) and (6.19) — including the divergence of both branches as : double descent, resolved in time.
(b) Probe ergodicity breaking: run two clones from initializations differing by a small random , and track the normalized separation . Verify that it decays to zero for but saturates for — at the predicted plateau , the paper's ergodicity-breaking degree, measured in twenty lines of code.
(c) Near , extract the relaxation time of and check for the critical slowing down of Eq. (6.21).
Solution
(a) With full batches and , gradient descent converges for ; keep a safe factor below. The late-time errors match (6.18) and (6.19), and plotting against traces the double-descent peak at from both sides.
(b) The clone separation obeys the linear dynamics : it decays iff the sampled covariance has full rank, i.e. iff (). For the covariance has a null space of dimension , and the projection of onto it survives forever; for an isotropic perturbation the expected surviving fraction is — precisely the paper's ergodicity-breaking degree .
(c) The relaxation time is set by the smallest nonzero eigenvalue of , which by Marchenko–Pastur closes like near threshold, so one expects ; compare with the paper's measured exponent and with the collapse ansatz Eq. (6.21).
6.3 Dynamical mean-field theory (DMFT)
Dynamical mean-field theory, of which Section 6.2 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), and (Cugliandolo 2023) 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; Kotliar et al. 2006). 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).
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), makes the issue vivid (and is worked in detail in Exercise 6.3). Consider the stochastic process , with a Wiener process (so that is the white noise of Section 6.1), and suppose with : a stable fixed point of the noiseless dynamics. One can expand around it, writing , to obtain the approximation
where 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 , however using these results with the Taylor expansion of to second order gives a correction to the mean
which is inconsistent with the value assumed for the mean when constructing the approximation: the noise, filtered through the curvature of , shifts the very point we expanded around.
Exercise 6.3 (Self-consistency in the smallest possible example). (a) Derive the stationary variance of the linearized dynamics (compare Exercise 6.1(b)).
(b) Carrying the Taylor expansion of to second order, derive the quoted correction to the mean.
(c) Propose the self-consistent repair: expand instead around a point , chosen so that the effective drift vanishes, , with 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, and of the DMFT equations to come.
Solution
(a) The linearized equation is an Ornstein–Uhlenbeck process with ; by Exercise 6.1(b) its stationary variance is .
(b) Taking the stationary average of with expanded to second order about : , and using , , the quoted shift — nonzero, contradicting the assumption .
(c) Demand instead that the expansion point satisfy the effective stationarity condition together with . These are two equations determining jointly: the fluctuations shape the effective potential, and the effective potential shapes the fluctuations. Structurally this is again — and DMFT is the same fixed-point demand imposed on entire two-time functions , instead of scalars.
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), 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 that appear next to the actual expansion parameter , and using the nPI framework effectively resums these terms properly at each order in . 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- models in physics (e.g. the – 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 's of Section 6.2, and then to evaluate their path integral by saddle point in a suitable limit. This is how the deep-learning DMFT of (Bordelon & Pehlevan 2022; Bordelon et al. 2024) 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) brought spin-glass methods into neural networks (in the guise of Hopfield-style associative memories), and Sompolinsky, Crisanti, and Sommers (Sompolinsky et al. 1988) 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 and is already hard; a disordered system has random couplings , 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.
The Bordelon–Pehlevan program. For deep networks trained by gradient descent at fixed dataset, the DMFT of (Bordelon & Pehlevan 2022) takes the infinite-width limit and finds that the order parameters are precisely the feature kernels and gradient kernels — two-time generalizations of the NNGP kernel and NTK of Sections 3–4. The lazy regime of Section 4 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), hyperparameter transfer across depth (Bordelon et al. 2023), and the dynamical theory of neural scaling laws (Bordelon et al. 2024) that we meet again in Section 7.9.
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: DMFT for Gaussian-mixture classification (Mignacco et al. 2021), and a DMFT characterization of the effective noise of SGD at late times (Mignacco & Urbani 2022) — the dynamical refinement of the batch-noise discussion of Section 2.1.
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; Gerbelot et al. 2024), 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). A recent tour de force in this direction (Montanari & Urbani 2025) 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 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): 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 states the mechanism common to every proposal — power laws in the data produce power laws in the loss — and Section 7.4 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, run as a derive-it-yourself exercise, with its payoff for "emergent capabilities" in Section 7.6); the geometry of data and the spectra of kernels (Section 7.7); and solvable large- field theories, static and then dynamical (Sections 7.8–7.9). We close with the newest measurements, for which there is no theory yet: post-training and test-time scaling (Section 7.10).
Notation. We keep the conventions of these notes: is the parameter count and the dataset size. Be warned that the scaling-law literature writes for our and for our (their is not our width!), so the famous exponents are our . For reading the sources:
| these notes | scaling-law literature (Kaplan et al. 2020; Hoffmann et al. 2022; Michaud et al. 2023; Maloney et al. 2022) |
|---|---|
| (parameters) | |
| (dataset size) | |
| , (exponents) | , |
| (quanta learned, Section 7.5) |
Note also that in this section (no subscript) denotes the Zipf exponent of the quantization model, not the load of Section 2.
7.1 The empirical laws
Kaplan et al. (Kaplan et al. 2020) established, for autoregressive language models, that when a model is not bottlenecked by the other resources its test loss obeys
where is the number of (non-embedding) parameters, the dataset size in tokens, and the training compute in FLOPs. The measured exponents were , , and , 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 . 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 tokens can be one giant batch or a million small ones, and below the "critical batch size" (which (Kaplan et al. 2020) 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
Per token, the forward pass costs about 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 FLOPs per token, times tokens. The 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) 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 ; the computational resource is the number of tokens processed, over epochs, so that in general . Fixing and growing by repetition buys optimization, not statistics: the loss converges to the -limited floor and then overfits. In the language of Section 2.1, more epochs means larger 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). This matters increasingly, now that frontier training is data-constrained.
The compute law is not independent. Since depends only on the other two resources, how can compute ever be the unbottlenecked one? The answer is that 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 a derived exponent. Minimizing an additive joint law at fixed balances the two terms, and gives
The harmonic sum makes smaller than either single-resource exponent, since a factor of must be split between growing and growing . Kaplan's values give , against the measured : 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) trained autoregressive transformers on many kinds of data and found power laws everywhere, but with exponents that depend on the data:
| modality | ||
|---|---|---|
| language | ||
| images | ||
| images | ||
| video | ||
| math (procedural) |
("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 – 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). This asymmetry is the central clue of the section: whatever sets the exponent lives in the data. Section 7.7 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) refit with the additive ansatz
The additive form is easy to read. 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) 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). Kaplan et al.'s analysis gave and : with more compute, grow the model. For two years the field trained ever-larger models on comparatively little data. Chinchilla's fit gives and : grow both equally, or as a rule of thumb about 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): 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 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) 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 determined by ?); 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 falls off as , over about four decades of rank (Figure 7). It is a scale-free tail of ever rarer structure, and it is not special to words. Sequences of consecutive words (-grams, the unit of pre-neural language modeling) are Zipf-like at every , with slopes that flatten as 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). The covariance and kernel spectra of text and image data are measured power laws (Sections 7.7–7.8). 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: tokens drawn from a distribution over a vocabulary of words, with the empirical rank–frequency curve against the 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 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 . The quantization model of Section 7.5 gives a loss exponent equal to , and exact Zipf (, that is, ) is the boundary at which the tail sum decays only logarithmically and power-law scaling degenerates. Read this way, the measured for language says that its effective skill distribution is , barely steeper than pure Zipf. Language scaling is slow because its structure sits just above the marginal tail. Image data has steeper tails and scales faster; compare the modality table of Section 7.1.
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), 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); this has been contested as an artifact of the metric (Schaeffer et al. 2023), and the quantization model of Section 7.6 gives it a sharper reading. A fourth gap is newer: the post-training and test-time scaling laws of Section 7.10 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) | discrete skill "quanta," Zipf-used, learned in order | none (a counting argument) | in terms of the Zipf exponent ; plus the relation | separate , , laws |
| data manifold; Sharma–Kaplan (Sharma & Kaplan 2020) | smooth target on a -dimensional manifold | piecewise-linear interpolation (function class only; no training) | yes: | no ( only) |
| percolation; Brill (Brill 2024) | 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), Bordelon–Canatar–Pehlevan (Bordelon et al. 2020; Canatar et al. 2021; Simon et al. 2023) | 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) | latent power-law covariance, extended by random features | linear model on fixed random features (frozen kernel); large- RMT, statics | in terms of the data's spectral exponent | yes: joint , with duality |
| DMFT; Bordelon–Atanasov–Pehlevan (Bordelon et al. 2024) | power-law (source/capacity) spectra | DMFT of gradient-flow training on fixed random features (frozen kernel; dynamics) | yes, including the time exponent | fully joint 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 and 4.6; 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)). 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 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 and 6 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); on the other hand, the model's parameter-free relation between and is violated by the measured ordering of Kaplan's exponents (Exercise 7.1(e)). The manifold prediction was confirmed in teacher–student experiments, where and can be measured independently, and tested more loosely on image classifiers and GPT-type language models (Sharma & Kaplan 2020); for language, taking at face value gives (Exercise 7.2). The kernel predictions have been tested on standard architectures and datasets, with exponents computed from measured kernel spectra, and the predicted duality between the - and -exponents holds up (Bahri et al. 2024). The two field-theory rows take the data's spectral exponent as input. Their qualitative predictions — exponents independent of architecture, , 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) 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), 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.
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 quanta of a fixed sequence.
- QH3 (Zipfian use). Quantum is used on a fraction
of samples ( a constant), and on those samples, knowing it reduces the loss by a constant (independent of ).
Consequently, a model that has learned the first quanta has expected loss
That is the entire model. Everything else is asking how is limited by each resource.
Exercise 7.1 (★) (Three scaling laws from three assumptions). (a) Estimate the tail sum in (7.8) to show .
(b) Parameter-limited (data and training unlimited). Suppose each quantum costs a constant capacity of parameters to learn. Derive in terms of .
(c) Data-limited (capacity and training unlimited). Suppose a quantum is learnable only if it is used by at least training samples, for constant . Given i.i.d. samples, which quanta are learnable? Derive .
(d) Training-limited, single epoch. Suppose instead that learning quantum requires gradient updates on batches containing it. After steps of online training (fresh batches), derive .
(e) Falsifiability. Parts (b)–(c) predict a parameter-free relation between and : write it down, and note which of the two exponents must be larger. Compare with the measured values quoted after (7.1). 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.)
Solution
(a) , so . (The same replace-sum-by-integral move as the entropy estimates of Section 2.4.)
(b) quanta fit, so : .
(c) Quantum is seen times, so it is learnable iff , i.e. by (7.7). Then : .
(d) Identical threshold logic with in place of : . Note the model predicts — data and single-epoch training time are interchangeable bottlenecks, a foreshadowing of the joint treatment in Section 7.9.
(e) Eliminating : , hence always. But the measured LLM values have : the predicted ordering is violated, not just the numbers. Candidate culprits, roughly in order of suspicion: (i) the measured comes from multi-epoch training with early stopping, not the clean "capacity-unlimited, data-limited" protocol assumed in (c); (ii) the threshold and per-quantum cost need not be -independent (rare quanta may be intrinsically harder or easier); (iii) constant across quanta is crude; (iv) real capabilities may not decompose into independent quanta at all — they share structure, so learning one changes the cost of others (see the lab of Section 7.6, where exactly this happens). The instructive point is not which repair is right but that the model is sharp enough to be wrong — compare the falsifiability discussion of Section 1.
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) is a smooth power law. There is no contradiction. A smooth law is a sum of many staggered cliffs, each weighted by and individually invisible in the aggregate, since quantum contributes only the small fraction 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) are neither magic nor mirage (Schaeffer et al. 2023): 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): smooth scaling from staggered cliffs. Setup: subtasks with Zipfian frequencies , ; inputs are a one-hot task indicator ( bits) concatenated with random bits; the label is the parity of a fixed random -subset of the bits, the subset depending on the task. Model: one-hidden-layer ReLU MLP, width , MSE loss on labels, trained online (fresh batches of ; Adam, learning rate , 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 is smooth; the dashed guide is the idealized quanta-model slope from Exercise 7.1(d). 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; students are encouraged to reproduce and modify it (all details are in the caption). A network is trained online on 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 of Exercise 7.1(d), because the cliffs bunch up rather than spreading Zipf-like; the rare tasks are learned sooner than the threshold rule 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) fails.
7.7 Exponents from geometry and from spectra
The quantization model takes the Zipf exponent 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) propose the simplest continuous mechanism. A network with ReLU activations computes a piecewise-linear function, so regression on a data manifold of intrinsic dimension is a problem of piecewise-linear interpolation. A network with parameters can afford linear pieces; tiling a -dimensional manifold with pieces gives spacing ; a linear patch approximates a smooth target to accuracy (Taylor); and MSE is the square of the error. Hence
Exponents are small because data manifolds are high-dimensional. The formula is testable, since can be measured with intrinsic-dimension estimators, which is exactly what (Sharma & Kaplan 2020) do.
Exercise 7.2 (Exponents from the data manifold). (a) Fill in the derivation of (7.9), stating precisely where smoothness of the target is used.
(b) How does the answer change for mean-absolute-error loss?
(c) Taking the measured 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 , if samples are what limits the resolution of the tiling?
Solution
(a) As in the text; smoothness (bounded second derivatives) enters in the per-patch error .
(b) Mean-absolute error is the error itself, not its square: .
(c) : an effective, intrinsic dimension of the distribution of text as seen by the model — vastly smaller than the nominal input dimension (vocabulary context), which is the whole point.
(d) If data limits the tiling, samples resolve spacing and : this mechanism predicts equal exponents — closer to what is observed than the quanta prediction, and the same conclusion the field theory of Section 7.8 reaches from spectra.
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, whose mode-filtering formula said that with samples (and inverse temperature ), the modes of the GP kernel that get learned are those with . Suppose the kernel spectrum and the target's decomposition are power laws,
in a shared eigenbasis . The unlearned modes dominate the test error:
giving : 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; Canatar et al. 2021; Simon et al. 2023), which sharpen (7.11) into precise learning curves (including the overfitting factor whose divergence produces the double-descent peak of Section 6). Bahri et al. (Bahri et al. 2024) organize the possibilities into a clean taxonomy: resolution-limited regimes (exponents from spectra or geometry, as here and in (7.9)) versus variance-limited regimes (universal exponent ), in each of and .
Exercise 7.3 (Spectral scaling laws, and quanta as eigenmodes). (a) Derive (7.11) from the mode-filtering formula of Section 4.6, and check that the crossover modes do not change the exponent.
(b) Show that the quantization model's data-scaling law is the special case of (7.11) under the dictionary: quanta eigenmodes, Zipf law spectrum, . (Hint: in the quanta model, plays both roles, setting the learning threshold and the contribution to the loss, so the two power laws coincide.)
(c) In the quanta model, what would it mean for ? Invent a modification of QH3 that realizes it.
Solution
(a) The learned modes contribute for ; the unlearned ones contribute ; the crossover window contributes as well, so the exponent is unchanged. Hence .
(b) With : , i.e. the exponent of Exercise 7.1(c). The threshold is the threshold ; the tail sum of is the tail sum of . Kernel regression is a quantization model whose quanta are eigenmodes.
(c) decouples how often a skill is tested (its loss weight) from how easy it is to learn (its threshold) — e.g. let quantum require occurrences. Any monotone relation between frequency and difficulty generates a two-exponent family, and (7.11) is its universal form.
7.8 A solvable large-N field theory
Maloney, Roberts, and Sully (Maloney et al. 2022) supply what the arguments above lack: a microscopic, exactly solvable model in which the joint law 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 , as measured in real image and text data (the spectral exponent is the paper's ; we reserve for the scaling dial of Section 5.4); a random-feature linear student, so we are in the lazy limit of Section 4, whose random features are its trainable parameters, trained on samples; and the joint thermodynamic limit at fixed ratio. The computation proceeds at large and , where expectations over data and features become Gaussian integrals, the test loss becomes a product of resolvents (Green's functions), and the expansion organizes into planar ('t Hooft) diagrams that can be resummed exactly.
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, and in the strict infinite-width lazy limit the loss would not depend on at all. The -scaling law is a finite-size effect on the lazy limit. At finite 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, 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 . This predicts , consistent with observation, and opposite in spirit to the quanta model's asymmetric prediction (Exercise 7.1(e)). It makes equiparameterization compute-optimal, a first-principles Chinchilla rule.
- Where laws break: scaling stops when 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 applied to scaling. Bordelon, Atanasov, and Pehlevan (Bordelon et al. 2024) 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 , with the correlation and response functions of Section 6 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 in the appropriate corners. Moreover, the theory predicts the full compute-optimal frontier: minimizing at fixed compute 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, so is dynamical mean-field theory overkill? No, for two reasons. First, time is a genuine third resource. Predicting and the compute frontier requires the whole training trajectory, and the statics only sees the corner. Second, the naive mode-by-mode solution of the lazy regime, in which each kernel eigenmode of the residual decays as as in (4.34), 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, and at 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), 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. 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).
Reward-model optimization introduces a new scaling axis: not data or parameters, but the distance moved from the base policy, . Gao, Schulman, and Hilton (Gao et al. 2023) found that the true ("gold") reward follows simple laws in , namely for best-of- sampling and 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). 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), and train-time and test-time compute can be traded against each other along a compute-optimal frontier (Snell et al. 2024).
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 introduced the quenched free energy and the annealed approximation (2.18) that replaces it by , and Section 2.6 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 as the running example, in the notation and normalization of that section: binary weights , a binary teacher , Gaussian inputs, examples, and the partition function (2.39). The classic reference for everything here is Engel and Van den Broeck (Engel & Van den Broeck 2001).
A.1 What the annealed approximation gets wrong
Jensen's inequality gives , so the annealed free energy is always a lower bound on the quenched one. To see what drives the gap, write , where the free-energy density fluctuates from dataset to dataset. Exponentiation turns modest fluctuations of into enormous fluctuations of : a dataset whose is below the typical value by a constant is exponentially rare, but it contributes an exponentially large , and the two exponentials compete at the same order in . The average can therefore be dominated by vanishingly rare datasets, while the typical dataset, which is what 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 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 , where the exact calculation sketched below gives (Gy"orgyi 1990). 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 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 through the identity
The strategy is to compute for integer , where it has a transparent meaning, obtain a formula analytic in , and continue it to . For integer , is the partition function of independent copies, or replicas, of the student, all judged against the same dataset. For the Ising perceptron, from (2.39),
The data average is easy, because the inputs are Gaussian. For fixed replicas and teacher, each example enters only through the numbers and , which are exactly jointly Gaussian (they are linear in ), with mean zero and a covariance determined entirely by the overlap matrix
The product of step functions asks that every have the sign of , and the probability of that event is some function of the overlaps alone: a Gaussian orthant probability in dimensions. Since the examples are independent, the data average of the product over is the -th power, . For the only overlap is , the orthant is a wedge, and is the energy term of Section 2.6.
The disorder is gone. In its place, the replicas are coupled through the overlaps , and the sum over the binary vectors can be organized by the value of the overlaps, exactly as the sum over one binary vector was organized by in Exercise 2.3(c). Writing for the number of -tuples of students with prescribed overlaps (for this is the binomial count of Exercise 2.3(c)), and replacing the sum over overlaps by an integral,
This is the replicated version of the one-variable action (2.44) of Section 2.6, to which it reduces at . At large the integral is dominated by the minimum of over the overlaps, and the replica identity (A.1) then gives the quenched free-energy density
Compare the annealed approximation, which is : the annealed calculation keeps a single replica, and so never sees the coupling between replicas that the shared data induce. That coupling is precisely the information about dataset-to-dataset fluctuations that Appendix A.1 found missing.
Replica symmetry. To make (A.5) concrete one needs and as functions of the overlaps, and a way to continue the result in . The standard ansatz, replica symmetry, is to take for all and . It is natural, since the saddle-point equations are symmetric under permuting replicas, and it reduces the problem to two scalars. Under it, both and become one-dimensional Gaussian integrals: the orthant probability collapses to an integral over a single shared Gaussian "field" of the Gaussian tail function , and the count 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. The resulting expression is analytic in and can be differentiated at . 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) reproduces the two-basin structure found in Exercise 2.3: an interior solution with that drifts toward the teacher as data accumulate, and the teacher solution , which exchange dominance in a first-order transition. Only the number moves, from to . 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, chs. 2 and 7).
(a) Replicate. For integer , write as in (A.2): a -fold copy of the system, coupled only through the shared data.
(b) Average over data. Show that for fixed replicas and teacher, and are exactly jointly Gaussian with the covariance (A.3), and that the data average factorizes across examples, each contributing the same function of the overlaps. Check that at , with the wedge formula (2.41).
(c) Push forward to the order parameters. Insert delta functions fixing and (in integral representation, with conjugate variables , ), so that the sum over the binary vectors factorizes over sites and becomes the -th power of a -spin sum. Conclude that takes the form (A.4), and that the annealed action (2.44) is its case.
(d) Replica symmetry. Take for and (and likewise for the conjugates). Show that the orthant probability becomes for an appropriate argument of the tail function , built from , , and a single Gaussian variable , and that the site sum becomes a one-dimensional Gaussian integral of . Continue both in and expand to first order at .
(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 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 are transcendental. Solve them numerically and show that there are two competing solutions, an interior one and , whose free energies cross at . Compare with the annealed of Exercise 2.3(g), and explain the direction of the shift using Appendix A.1.
Solution
The six steps are their own outline; we record only the checkpoints. (Step 2) Conditioned on the and , the and are linear in the Gaussian vector , hence exactly Gaussian, with the covariances (A.3); independence over examples gives the -th power, so the per-example factor is with , a -dimensional Gaussian orthant probability. At the orthant is the double wedge of Exercise 2.3(b), of probability . (Step 3) With and likewise for , the exponent is a sum over sites of identical terms, so ; the saddle point over the conjugates then defines . At only survives and the site sum is , reproducing the binomial entropy after the saddle in . (Step 4) Under replica symmetry, the correlated Gaussians can be written as with shared and the independent, and correlated with through ; integrating out the gives the orthant probability as , where the argument of is linear in with coefficients built from and . The site sum becomes after decoupling the term by the same trick. Both are analytic in , and as . (Step 5) The teacher is a uniformly random binary vector consistent with every example, exactly as a Gibbs student is, so the two are exchangeable, and at the saddle: . (Step 6) Solving the remaining pair of equations numerically, the interior solution exists for all up to a spinodal and the teacher solution always exists; their free energies cross at (Gy"orgyi 1990), below the annealed . The annealed average is dominated by rare, weakly constraining datasets, which overweight poorly aligned students and so require more data before the teacher basin wins; removing that bias moves the transition earlier.
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 the same, but let the weights be real, . Since is invariant under with , remove the redundant scale by putting the weights on a sphere,
normalized so that individual weights stay of order one, and let the prior be the uniform measure on the sphere. This is the spherical perceptron. The wedge formula (2.41) for the generalization error is unchanged, since its derivation used only that is a Gaussian pair with correlation , which holds for any with . What changes is the entropy: the number of students at overlap becomes the volume of the slice of the sphere at overlap . Decompose with . The constraint (A.6) gives , so the slice is a sphere of radius in the -dimensional orthogonal complement, of volume , and
Repeating the annealed calculation of Exercise 2.3(d) with this entropy gives
and the annealed free-energy density is
Figure 9 shows for several loads. The entropy–energy competition is the same as for the Ising perceptron, but the outcome is different: the entropy diverges as , so the perfectly aligned state is infinitely penalized, the minimum always sits in the interior, and it slides continuously from toward as grows. There is no phase transition. This is the answer to the closing question of Exercise 2.3(h): the binary weights of the Ising perceptron give the aligned state finite entropy, , and that is what allows it to compete as a separate basin.
The annealed action of the spherical perceptron, (A.8), for (bottom to top at ), with the minimizer marked on each curve. The entropy pins the minimum at when ; increasing the load slides it smoothly toward larger overlap, and the entropic wall at keeps it in the interior. One basin at every , moving continuously: no phase transition.
The learning curve. The large- behavior follows by expanding near . Writing , so that is the generalization error itself, and
which is minimized at
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. The exact, quenched treatment (Seung et al. 1992; Engel & Van den Broeck 2001) runs through the replica steps of Exercise A.1 with the volume entropy (A.7) 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, the annealed average overweights poorly aligned students and overestimates the error, ; 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), and assemble the annealed action (A.8) by repeating the steps of Exercise 2.3(d).
(b) Derive the saddle-point condition for (A.9), and show that is never a minimum, in contrast with Exercise 2.3(f).
(c) Verify the large- learning curve (A.11), and show that at small the annealed theory predicts . Sketch the full learning curve .
(d) Which steps of Exercise A.1 change when the count of binary students is replaced by the volume (A.7), and which do not?
Solution
(a) Decompose . The spherical constraint forces , so the slice at fixed is a sphere of radius in the -dimensional orthogonal complement, of volume , giving up to constants. The data factor is as in Exercise 2.3(d), and collecting exponents gives (A.8).
(b) Setting with gives
As the entropy term while the energy term stays finite, so and is never a minimum; the minimizer is always interior.
(c) With and : , so , minimized at . For small the entropic part is and , so , and minimizing gives . The learning curve rises linearly from , bends over, and approaches with the power-law tail .
(d) Steps 1, 2, 5, and 6 are unchanged: the data average sees only the overlaps, and the Gibbs-symmetry argument holds for any prior invariant under the symmetries of the problem. Step 3 changes: the sum over binary vectors becomes an integral over points of the sphere, and the count becomes the log-volume of the set of -tuples with prescribed overlaps, a Gaussian integral with a of the overlap matrix; under replica symmetry it is up to constants, whose derivative reproduces (A.7) at . Step 4 is then the same -function integral. The outcome differs because this entropy diverges at : no competing basin, no transition, and a learning curve (Seung et al. 1992).
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) of Section 4 — 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, and this appendix supplies the linear algebra they rest on. Readers who found the kernel formulas of Sections 4.1 and 4.6 appearing out of thin air may want to read it first. 6
We use the notation of Section 3.1: training data with and ; a model with parameters ; and for its values on the training inputs. It is convenient to stack the data into matrices,
so that is the -th coordinate of the -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,
with the residual. Many references (including the notebook this appendix is drawn from) normalize the loss as instead; this rescales the ridge parameter below by a factor of and the learning rate by a factor of , and changes nothing else.
B.1 Linear regression
The simplest model is linear in the inputs,
so that and the vector of training predictions is . (An intercept can be included by appending a constant coordinate to every input, as in Exercise B.1 below.) The loss is a convex quadratic function of , with gradient
This has the shape that every gradient in these notes has: the transpose of the Jacobian — here the constant matrix — acting on the residual, cf. (4.10). Setting the gradient to zero gives the normal equations
Since the loss is convex, every solution of (B.5) is a global minimizer, and there are no other critical points. What the solutions look like depends on how compares with .
Fewer parameters than data, . If the columns of are linearly independent — the generic situation when — then is positive definite, since for . There is a unique minimizer,
The normal equations say that the residual is orthogonal to every column of . Equivalently, the prediction is the orthogonal projection of onto the column space , the set of all prediction vectors the model can produce (Figure 10). Fitting is projecting. The projection is unique even in cases where is not.
Least squares, geometrically. The normal equations say that the residual is orthogonal to the column space of , so the fitted prediction vector is the orthogonal projection of onto . "Orthogonal" refers to the Euclidean metric on the space of function values on the training set.
It is worth pausing on the word "orthogonal." Least squares minimizes the Euclidean () distance between and , and the projection is orthogonal with respect to the Euclidean metric on — the same metric on the space of function values that was made explicit in Section 4.1.1. 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 on defines a different loss , with a different minimizer
which is the -orthogonal projection of onto ("weighted least squares"). Unless already lies in , different metrics give genuinely different fitted functions; see Exercise B.1(d). Note that is a metric on the space of function values indexed by the data points, not a metric on the input space .
More parameters than data, . If , or more generally if the columns of are linearly dependent, then is singular. The normal equations still have solutions (Exercise B.1(b)), but now infinitely many: an affine subspace of of dimension . If moreover has full row rank , then and every solution interpolates the data exactly, , 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, and it is where kernels come from. Far from being an edge case, is the normal situation for neural networks.
Exercise B.1 (Least squares by hand). (a) Differentiate with respect to a single coordinate , assemble the components into (B.4), and derive the normal equations (B.5). Explain geometrically why they say that the residual is orthogonal to every column of .
(b) Show that if has full column rank then is positive definite, so that (B.6) is the unique minimizer. Show that the normal equations are always consistent (i.e. ), and that if they have infinitely many solutions.
(c) Fit the affine function to the three points with labels ; that is, take
Find , 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) for , which trusts the middle data point four times as much as the others. Check that the fitted function changes.
Solution
(a) . The -th normal equation sets the inner product of the -th column of with the residual to zero, so ; hence is the orthogonal projection of onto .
(b) vanishes only if , which for is impossible when the columns are independent; so is positive definite and invertible. Consistency: (if then ), and for any matrix one has , so . If then , and adding any element of to a solution gives another.
(c) , , so : the constant function . Predictions , residual , loss . The optimization is exact (the problem is convex and was solved in closed form); the error is representational — no affine function agrees with at three points. Exercise B.2 repairs it by adding a feature.
(d) and , so : the constant function , pulled toward the heavily weighted middle label . Residual . A different metric on gives a different fit.
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
and consider the model
This is a weighted sum of fixed shapes . It can be as nonlinear in as the shapes are, but it is linear in . Standard choices include monomials in the coordinates of up to some degree (giving ), Fourier modes, Gaussian bumps centered at chosen points , and random features with the 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 feature matrix
the training predictions are . Note that : the feature matrix is precisely the Jacobian of Section 4.2, which for a linear model is a constant. Linear regression is the special case , , and everything in Appendix B.1 goes through verbatim with replaced by : the normal equations are ; if and has full column rank, the unique minimizer is ; and if there are infinitely many interpolants.
Choosing is choosing the model. The same seven data points (black), interpolated exactly by three different fixed feature maps with : 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 by changes the algebra not at all, and changes the model completely. Figure 11 shows three feature maps, each with , 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 , 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 , the rings sit at two different heights, and a horizontal plane separates them. Nothing has been learned: we supplied the radial coordinate.
Why would one ever want a very large ? Figure 12 gives the classic answer. Two classes of points on concentric rings cannot be separated by a line in the plane. Append one fixed feature, , and the classes become separable by a plane in : a model that is linear in the lifted coordinates is a circle in the original ones. The price is that grows fast. All monomials of degree in input coordinates already give , and Gaussian bumps centered at every point of a continuum give . Yet there are only data points, and the fit is determined by numbers. Nothing about the problem should require a inverse — so where does 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 with labels as in Exercise B.1(c), but now the feature map .
(a) Form the feature matrix and find coefficients that interpolate all three labels. The representational failure of Exercise B.1(c) is gone: what changed?
(b) Explain precisely how the predictor can be nonlinear in while the regression problem remains linear.
(c) If the last feature is replaced by for some , what changes: the fitted function, its coefficient, or both?
(d) Write the inner product in closed form, without vector notation. This is a preview of Appendix B.3.
Solution
(a) , which is invertible (), so exactly one interpolates: , giving and predictions . What changed is that now lies in .
(b) All nonlinearity in sits inside the fixed functions ; depends on linearly, so the loss is quadratic in and the normal equations are linear.
(c) The coefficient becomes , the product is unchanged, and so is the fitted function. (With the ridge penalty of Appendix B.4 the rescaling would change the fit, since is not invariant under it.)
(d) .
B.3 Kernel regression
Return to the overparametrized case , and assume that has full row rank (the generic situation), so that there is a -dimensional affine family of exact interpolants . 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 on parameter space , the metric that also defines gradient descent in Section 4.1.1. (Appendix B.1 used on the space 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), introduce Lagrange multipliers for the constraints and extremize . Stationarity in gives
This is the key structural fact: however large is, the minimum-norm solution lies in the -dimensional subspace of spanned by the feature vectors of the training inputs, and is described by numbers. (In the language of Section 4.2, lies in the row space of the Jacobian.) The constraint then becomes , an 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 is fixed, and involving no labels. The Gram matrix is symmetric and positive semidefinite, since ; it is positive definite exactly when has full row rank, which we assumed. Then , and the prediction at an arbitrary input 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 on pairs of training points; to predict, one needs 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 that depend on the labels — but the labels only reweight the training points, they never change which inputs considers similar.
Compare (B.14) with the late-time mean (4.38) of a lazily trained network: they are the same formula, with in place of the NTK. This is no coincidence. For any model linear in its parameters the Jacobian is , so the NTK (4.10) is
constant in time. Moreover, gradient flow finds the minimum-norm solution: always points along the row space of , so stays in that -dimensional subspace and converges to the interpolant nearest to — for , to (Exercise B.3(c)). 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 , and so is the norm that it ends up minimizing.
The kernel trick. There are two readings of (B.14). Read from the feature side, it rewrites a problem as an one. Read from the kernel side, it says the feature map is dispensable: any function that arises as an inner product of features will do, whether or not one can write the features down. For and , one checks that : a single dot product, squared, stands in for three features. In dimensions, costs to evaluate and stands in for all monomials of degree . The Gaussian kernel corresponds to an infinite-dimensional feature map that no one ever writes down. Conversely, any symmetric all of whose Gram matrices are positive semidefinite arises from some feature map (Mercer's theorem): the eigenfunctions of , weighted by the square roots of its eigenvalues, furnish one (cf. Sections 3.3 and 4.6). This is a trade, not a free lunch. The feature route scales with and the kernel route with ; 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)–(B.14). Alternatively, use the singular value decomposition of Section 4.2 to show that the minimum-norm interpolant is , and write it in terms of .
(b) Show that for all . Interpret a zero eigenvalue of as a linear combination of training examples that vanishes in feature space. What does it take for to be invertible?
(c) Show that under gradient flow the displacement stays in the row space of for all ; derive ; and conclude that, for invertible , the flow converges to the interpolant closest to . Derive (4.36) for this model, with .
(d) Verify that for , and explain why the is not decoration. Which feature map on produces the kernel of Exercise B.2(d)? How many features does on stand in for?
Solution
(a) Stationarity of in gives ; the constraint then reads , so and ; the prediction is . It is a minimum, not merely a critical point, because the objective is strictly convex and the constraint set is affine. In the SVD of Section 4.2, with , one has and , where are the first columns of : only right singular vectors with nonzero singular value (the row space) appear.
(b) . If then : a linear dependence among the training feature vectors (for instance two identical inputs, or simply ). is invertible iff the training feature vectors are linearly independent, which requires .
(c) , the row space of , so stays there. Since is constant, , whence . Integrating, , and is (4.36) with and . At , interpolates, and . Any other interpolant has , so by Pythagoras : the flow converges to the interpolant closest to , and for to .
(d) ; without the the cross term would be and the sum would not be a perfect square. comes from , as in Exercise B.2. By the multinomial theorem , so over the monomials of degree .
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 , ridge regression minimizes
The two terms measure distance to the data with and distance to the origin with : ridge regression is a weighted compromise between the two metrics of Appendices B.1 and B.3, with setting the exchange rate. (In the training of neural networks the penalty is called weight decay.) The penalty does two jobs at once. When it selects one solution out of infinitely many; and for any 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 that were supplied.
The gradient is , so the stationarity condition is
The matrix on the left is positive definite for every , whatever the rank of , since . So there is always a unique minimizer,
at the price of a inverse. To move the inverse to the side, expand and regroup the product
Both bracketed matrices are positive definite, so multiplying by their inverses on the appropriate sides gives the push-through identity
Nothing was cancelled — is rectangular and has no inverse — and that is precisely why the identity holds. Substituting into (B.18), the ridge solution again lies in the span of the training feature vectors, , now with , and the prediction at a new input is
This is kernel ridge regression, the central formula of this appendix. As it reduces to the minimum-norm interpolant (B.14) when is invertible, and in general to the least-squares solution of minimum norm (Exercise B.4(c)). As , .
Formula (B.21) is exactly the posterior mean of Exercise 4.5, with ridge parameter — 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 and a Gaussian prior , the loss (B.16) is of the posterior density, up to a constant, with . Ridge regression is thus maximum a posteriori estimation. The prior pushed forward to is a Gaussian process with kernel , a finite- instance of Section 3.3, and the posterior mean of Exercise 4.5 for this process reproduces (B.21).
Spectral view. With the singular value decomposition of Section 4.2, the Gram matrix is , with eigenvalues and orthonormal eigenvectors : the "features" of Section 4.2, i.e. the natural modes in the space of function values. The ridge predictions on the training set are , or, mode by mode,
Ridge regression is a spectral filter: modes of the data with are fit, modes with are suppressed, and the crossover is at . Compare gradient flow from on the unpenalized loss, whose training predictions are, by (4.35),
This is also a spectral filter, with playing the role of : early stopping regularizes much as ridge does. But the two filters have different shapes (Figure 13), and no choice of reproduces a given stopping time on every mode at once (Exercise B.5). This is the precise content of the remark closing Section 4.6, 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 that is fit, as a function of the mode's kernel eigenvalue . The two are tuned to agree on the mode with , and disagree elsewhere; the dots mark the two-mode example of Exercise B.5.
Exercise B.4 (Ridge regression, two ways). (a) Derive (B.17) and (B.18), and prove that is invertible even when is rank deficient.
(b) Prove the push-through identity (B.20) and deduce (B.21). Explain why "cancelling " is not a valid step, and why the identity is nevertheless true. (It is a special case of .)
(c) What happens as ? Using the SVD, show that applies the scalar filter to each singular value of , and hence that as it converges to the least-squares solution of minimum norm. When does this limit interpolate the labels exactly?
(d) Return to , and the kernel of Exercise B.2. With , verify that
and predict at using only kernel evaluations. Compare with the interpolant of Exercise B.2(a), which predicts there, and explain the direction of the discrepancy.
(e) The feature route solves a system and the kernel route an one. When is each attractive?
Solution
(a) Setting gives (B.17); for , so the matrix is positive definite, hence invertible, for every .
(b) uses only associativity. Multiply on the left by and on the right by , both of which exist by (a). is and, for , has no two-sided inverse, so it cannot be "cancelled"; the identity holds because was regrouped, not divided out. Then and .
(c) As , . With , , so . As , for and for : this is the pseudo-inverse , the least-squares solution of minimum norm (an arbitrary unregularized minimizer need not have minimum norm). It interpolates iff , i.e. iff for every zero singular value.
(d) on gives the displayed . Check . At the kernel vector is , so , versus for the interpolant. The ridge fit is shrunk toward zero: it no longer even interpolates the training labels (its training predictions are ), and the quadratic growth it does fit is damped.
(e) The feature route costs and is attractive when and the features are explicit; the kernel route costs and is attractive when , or when can be evaluated without constructing at all. Neither is free: is fatal for large .
Exercise B.5 (Early stopping vs. ridge, on two data points). Take training points with Gram matrix and labels
and initial predictions .
(a) Verify the eigenpairs with , and decompose into these modes.
(b) Using (B.23), find the training predictions under gradient flow, and evaluate them at . Why do both the eigenvalue and the initial projection matter for what is visible during training?
(c) Using (B.22), find the ridge predictions for . The two filters agree on the mode with ; show that the ridge parameter matching gradient flow at time on a mode with eigenvalue is , 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?
Solution
(a) , so , ; and .
(b) . At , and , giving . The eigenvalue sets a mode's rate, while sets whether it is present and with what amplitude; a fast mode with zero initial amplitude leaves no trace in the trajectory.
(c) Ridge multipliers and , so . The flow multipliers at are and : the two agree on only. Solving gives ; at this is for but for . No single matches both modes.
(d) Here "spectral bias" means exactly this: components of the residual along eigenvectors of the fixed Gram matrix with larger eigenvalue decay faster, at rate . To speak of Fourier frequencies one needs a relation between the and Fourier modes (true, for instance, for translation-invariant kernels on uniform grids). To speak of noise one needs a model of how label noise projects onto the (typically spread evenly across modes, so that it lives mostly in the many slow modes — which is why early stopping and ridge both filter it). To speak of test error one needs a distribution over inputs and the population version of the kernel operator, as in Sections 4.6 and 7.7.
B.5 Linearized networks are kernel ridge regressors
Now let be a neural network, with the vector of all its weights and biases. It is not linear in . But its linearization around the initialization , the ansatz (4.1),
is a model of the form (B.9). Compare term by term. The parameters are the displacement . The feature map is
the tangent features: one number per parameter, the sensitivity of the output at to that parameter, evaluated at initialization and then frozen. And the constant offset 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) is the Jacobian at initialization, , and the Gram matrix (B.13) is
the neural tangent kernel (4.10) at initialization, extended to arbitrary pairs of inputs as in Exercise 4.3. Everything in Appendices B.3 and B.4 now applies verbatim. In particular, ridge regression on the linearized network with penalty — weight decay toward the initialization — gives, by (B.21),
At this is the limit of (4.36), the minimum-norm interpolant of the residuals; and the finite-time solution (4.36) is its early-stopped cousin, in the sense of (B.23). The kernel trick is not a convenience here but a necessity: is the parameter count of the network, often or more, while is .
Two remarks on what is and is not being claimed. First, for the linearized model everything above is exact, and its kernel is constant by construction. For the actual network, (B.24) is an ansatz. Section 4.3 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, not of this appendix. Second, the initialization is random, so , the tangent features, and the kernel are all random too, and (B.28) is the trained function for one draw. Averaging over draws, at large width where becomes deterministic, gives the Gaussian process (4.38).
Exercise B.6 (Tangent features of a one-hidden-layer network). Consider with , so that .
(a) Compute the tangent features , and show that the tangent kernel is
(b) Suppose only the output weights are trained, with the frozen at their random initial values. Show that this is exactly a fixed-feature model (B.9) 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 ?
(c) Show that gradient flow with weight decay toward the initialization, , applied to the linearized model , converges to the ridge solution (B.28) for every and from any starting point.
Solution
(a) and . Taking the inner product of the tangent features at and and summing over in the second block gives the displayed .
(b) With the frozen, with is linear in the trainable parameters , so exactly, , and is the first term: the random-features model. Training only the gives the second term, but there the linearization is not exact, since is nonlinear in .
(c) Writing , the linearized loss is , so with . This is a linear ODE with positive-definite , so converges from any start to , by push-through. Then is (B.28), since has components .
(d) As , ; as in (4.36), . Both give .
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. 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. 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. 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. 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. arXiv.
Leticia F. Cugliandolo (2023). Recent Applications of Dynamical Mean-Field Methods. arXiv.
Simon S. Du, Xiyu Zhai, Barnab'as P'oczos, and Aarti Singh (2019). Gradient Descent Provably Optimizes Over-parameterized Neural Networks. 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. 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.
Andrey Gromov (2023). Grokking modular arithmetic.
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. 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. 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. 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. 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. 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. 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. 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. 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.
David Tong (2012). Statistical Field Theory. University of Cambridge.
David Tong (2012). Lectures on Kinetic Theory. 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. 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.
Footnotes
-
A physicist's free energy is ; we absorb the factor of throughout, so that is dimensionless. ↩
-
For transformers, the limit analogous to "large width" that leads to a Gaussian process is a large number of attention heads. ↩
-
As a point of notation, and are also commonly used for noise, and we will distinguish between discrete and continuum delta functions as and respectively. ↩
-
Michaud's blog post ericjmichaud.com/quanta is an excellent companion to the paper and can substitute for it in this exercise. But attempt Exercise 7.1 before reading either; the derivations are short, and they are more instructive to find than to follow. ↩
-
This is also a working example for Section 6's toolkit: the same model can be solved dynamically, and its double-descent peak, at the interpolation threshold , is regularized by ridge or early stopping, connecting back to the divergence of the overfitting factor mentioned in Section 7.7. ↩
-
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. ↩