Double descent is the principle of least action Congzhou M Sha
Penn Medicine Doylestown Hospital, 595 W State St, Doylestown, PA 18901, United States
arXiv:2609.19076v1 [cs.LG] 16 Sep 2026
Abstract The test error of a model plotted against its number of parameters d falls, peaks when the model can just fit the training data, and falls again, exhibiting the double descent phenomenon. We explain the phenomenon with statistical mechanics. The training trajectory of a stochastic gradient-based method is a particle wandering over the energy landscape of the training loss at an induced temperature T , and a run that has equilibrated visits every parameter vector of a given training loss equally often, the fundamental postulate of statistical mechanics, with probability given by the Boltzmann distribution. Because training starts at an initial point and has only finite time to diffuse, it carries an effective weight decay, which makes every parameter a quadratic degree of freedom. The equipartition theorem then distributes the energy among the d degrees of freedom in shares of T /2, so at a fixed training loss adding parameters lowers the temperature and drives the Boltzmann distribution toward the stationary path. Finally, adding parameters can only lower the L2 norm of the stationary path, so a solution sampled at fixed loss is less likely to be large with increasing d, effectively increasing weight regularization.
1
Introduction
Consider the simple regression task as described by Nakkiran et al. (2019; 2020): fitting 20 noisy samples of a cubic function on the interval [−1, 1] using polynomials with varying numbers of parameters d (Fig. 1A–D). d counts all the coefficients of the polynomial including the constant term, so a polynomial of degree d − 1 has d parameters. d = 2, i.e. linear regression f (x) = θ1 x + θ0 , produces an underfit whereas d = 4, i.e. the cubic f (x) = θ3 x3 + θ2 x2 + θ1 x + θ0 , is the intended ground truth solution. For d = 20, there are exactly as many coefficients as there are points, so exactly one polynomial passes through all of the training data, and it oscillates wildly between the data. Surprisingly, fitting a d = 1000 curve also passes through every data point, yet it tracks the true cubic nearly as well as the d = 4 solution within the training data range, but without the enormous oscillations of the d = 20 polynomial. Plotting the test error as a function of d (Fig. 1E) shows a classical U-shape initially, with a peak at d = n, however there is a long decreasing tail to the right. This double descent phenomenon (Belkin et al., 2019; Nakkiran et al., 2020) appears in deep networks as a function of layer width, network depth, and number of training epochs. Exact treatments exist for linear and random-feature models (Hastie et al., 2022; Mei & Montanari, 2022; Bartlett et al., 2020), and the threshold has been identified with a jamming transition (Spigler et al., 2019; Geiger et al., 2020) and with broken ergodicity (Li & Goldenfeld, 2026). The aim of this work is to provide further elementary physical intuition for the double descent phenomenon, with a pedagogical focus so that it is understandable to undergraduates.
2
The statistical mechanics of the loss function
2.1
Microstates and the fundamental postulate
Statistical mechanics is the science of observing nature at large scales (Reif, 1965). Just as we may know the energy of a gas but are not concerned with the coordinates and velocity of every molecule, we may know the training loss of a model but are not concerned with the specific numerical values of the model’s parameters. 1
A
E
B
C
F
D
G
Figure 1: A–D Legendre-basis fits with d parameters (degree d − 1, constant term included) to n = 20 points yi = 3x3i − 2xi + ξi , ξi ∼ N (0, 0.42 ): least squares for d ≤ n, minimum-norm interpolation for d > n. E Test error of the minimal solution, the least squares fit for d ≤ n and the minimum-norm interpolant for d > n, which peaks at d = n. F Root mean square distance from the true cubic over the range I = [x1 , xn ] of the 1/2 R 1 training inputs, D2 = |I| (f −q)2 with q(x) = 3x3 −2x, of the minimal solution θ ∗ = (A⊤ A+λI)−1 A⊤ y I of Eq. (6), for four values of the weight decay λ. G The square of the smallest singular value of the fitting matrix, σmin (A)2 ; the horizontal lines mark the three regularization levels of F. The amplification of label noise into the fit is governed by 1/(σmin (A)2 + λ), so each curve in F peaks where σmin (A)2 falls below its λ.
Any set of trained weights which produces a low test loss generally makes the user happy. However, we can reason a posteriori about the properties of the minima and how the double descent phenomenon arises, given some assumptions on the definition of training for a sufficient length of time, i.e. ergodicity. Say that we flip ten coins and only know that seven are heads. A microstate in this context is defined as the specific sequence of coin flip results, e.g. HTHTHHHTHH. Since we are given only the information about the 10 total number of heads, multiple sequences (microstates) are possible, i.e. 10 = 7 3 = 120. We may then define the macrostate M of “sequences of 10 coin flips with 7 heads" as the set of these 120 sequences, and the number of heads is known as an observable O. The number of such microstates in a macrostate is known as the multiplicity Ω(M ). The entropy then is defined as S ∼ log Ω (we omit the Boltzmann factor kB ). Fundamental postulate of statistical mechanics. If all that is known about a physical system is its macroscopic observable, then every microstate consistent with that macrostate is equally probable to be the current physical state of the system. In the context of optimization, a microstate is a specific choice of numerical parameters θ for a model M , the observable is the training loss L(M (θ)) which P we will think of as the total energy E of the system, for 1 21 example the least squares objective L(M (θ)) = 2n i (M (xi ; θ) − yi ) , and the macrostate is the set of all parameters for which the model achieves this loss.
1We include an arbitrary factor of 1 to interpret the loss as potential energy. 2
2
2.2
Ergodicity
During model training, stochastic gradient descent (Robbins & Monro, 1951) or another stochastic gradientbased method (Kingma & Ba, 2015; Loshchilov & Hutter, 2019) is used to move from higher to lower values of the loss function. When the training error levels off, training is assumed to have reached an equilibrium. During training, a trajectory of parameters {θ 1 , θ 2 , · · · } are visited. A particular training trajectory {θ i } is called ergodic if it visits a representative sample of microstates, such that computing any macroscopic observable based on that sample yields the same value as if we performed the analytic sample across all solutions. An ergodic trajectory that has settled at a training loss E therefore samples the macrostate at that loss, the level set {θ : L(θ) = E}. In the polynomial fitting example this set is explicit. For d ≥ n the loss can reach E = 0 exactly, the interpolating regime, and the macrostate is the (d − n)-dimensional set of interpolants; for E > 0 it is a contour of dimension d − 1. What the analysis needs is that the trajectory samples this set without preference for any part of it, which we take as an assumption: Ergodicity of model training. Training a randomly-initialized model M with d parameters to a specified tolerance (e.g. L(M (θ)) = E) produces a microstate θ d which is drawn uniformly from the corresponding macrostate. With this assumption in place, there is an analogy between coin flips and model training (Table 1 in Appendix A.1). 2.3
The Boltzmann distribution and equipartition
The fundamental postulate weights every microstate of a fixed energy equally. Training however, does not fix the energy: the loss fluctuates from step to step. Each mini-batch reports a slightly different gradient, the gradient sampling process itself is often drawn from a noisy distribution (Kingma & Ba, 2015; Loshchilov & Hutter, 2019). The most accurate description of training is a probability distribution over energies, and the fundamental postulate determines that this distribution must be the Boltzmann distribution (Appendix A.2): X e−Es /T P (s) = , Z= e−Es /T , (1) Z s where the sum, or integral, runs over all microstates of the system regardless of energy, and Z normalizes the probabilities (Reif, 1965). In the context of optimization, the energy is the training loss. Training for a finite time regularizes the problem, and is equivalent to weight decay of strength λ ∝ 1/t (Krogh & Hertz, 1992; Ali et al., 2019), so we take the energy to be the loss with weight decay, Lλ (θ) = L(θ) + λ2 ∥θ∥2 ,
λ > 0.
(2)
For the least squares problem in Fig. 1, L(θ) = 12 ∥Aθ − y∥2 and A is the fitting matrix of the basis functions at the training inputs Eq. (20). The probability density of observing a candidate parameter θ during training given the observed labels is then h L (θ) i λ P (θ | y) ∝ exp − , (3) T with T determined by the properties of the optimization algorithm. For least squares, the loss function with weight decay is a paraboloid. Expanding the norms and completing the square, Lλ (θ) = 21 θ ⊤ A⊤ Aθ − y ⊤ Aθ + 12 ∥y∥2 + λ2 θ ⊤ θ
(4)
M = A⊤ A + λI = 21 θ ⊤ M θ − y ⊤ Aθ + 12 ∥y∥2 , = 21 (θ − θ ∗ )⊤ M (θ − θ ∗ ) + Emin , θ ∗ = M −1 A⊤ y, 3
(5) Emin = 21 ∥y∥2 − 12 y ⊤ AM −1 A⊤ y,
(6)
the macrostate ME = { : L( ) = E}, every parameter vector at training loss E
a trained model: a point picked at random
the center * , the minimum of the loss
its projection onto the parameter plane, an ellipse
Figure 2: The fundamental postulate. The training loss L(θ) as a landscape over two parameters θ = (θ1 , θ2 ), with the contour L(θ) = E drawn on the surface (orange): it is the macrostate ME , every parameter vector at training loss E, and its projection onto the parameter plane is an ellipse (dashed). By the fundamental postulate, a converged training run limited to a particular fixed energy visits every point of the contour of that energy equally often, so a trained model at that energy is a point picked at random from it (blue). The center θ ∗ (black) is the local minimum. where the third line follows from 12 (θ − θ ∗ )⊤ M (θ − θ ∗ ) = 12 θ ⊤ M θ − θ ⊤ M θ ∗ + 12 θ ∗⊤ M θ ∗ and M θ ∗ = A⊤ y. The matrix M is symmetric positive definite for λ > 0, so the loss has its minimum Emin = Lλ (θ ∗ ) at θ ∗ , the ridge regression fit, and each contour Lλ = E is the ellipsoid (θ − θ ∗ )⊤ M (θ − θ ∗ ) = 2(E − Emin ). Any C ω loss2 has this form to leading order near a minimum θ ∗ , with the Hessian H in place of M 3 (Fig. 2). A key consequence of the Boltzmann distribution is the equipartition theorem. Near a local energy minimum, the loss is quadratic, and under the Boltzmann weight each of the d principal directions carries the same share T /2 of the loss above the minimum (Appendix A.3), so E[L] = Emin +
Td . 2
(7)
Near the minimum, L = Emin + 12 ∥w∥2 in the coordinates w = H 1/2 (θ − θ ∗ ) (Eq. (6) for least squares, with H = M ), and the contour at E = Emin + ∆E is the sphere ∥w∥2 = 2∆E, on which the fundamental postulate makes w uniform. A uniform point on a sphere has E[wi2 ] = 2∆E/d for every coordinate, so each of the d directions carries the same share ∆E/d of the excess loss. For least squares, the data term 12 ∥Aθ − y∥2 changes only along the n directions of the row space of AM −1/2 , which we call seen by the data, and not along the d − n directions of its null space, which we call unseen. The seen directions, which alone determine the fit at the training inputs, therefore carry h i n EE 21 ∥ws ∥2 = ∆E, (8) d 2We require that the loss function be real analytic (C ω ) and not simply smoothly differentiable (C ∞ ) to avoid pathologic functions such as e−1/x , whose Taylor series does not converge to the function itself at x = 0. 3 By definition, at a critical point of a function, its gradient is the zero vector 0, and the point curvature is positive
4
which falls as 1/d. At a fixed training loss E = Emin + ∆E, additional parameters beyond the interpolant d > n therefore result in cooling of the temperature of the training algorithm along the seen directions. 2.4
The path integral over curves
Read as a statement about curves rather than parameters, Eq. (3) is a path integral (Feynman & Hibbs, 1965). The model maps each parameter vector to a curve f (·; θ), so the sum over microstates is a sum over every curve the model can draw, each weighted by e−Lλ /T , and a prediction is the average over all of them, Z Z 2 1 −L(θ)/T −λ∥θ∥2 /2T dθ f (x; θ) e−L(θ)/T e−λ∥θ∥ /2T . (9) Z = dθ e e , ⟨f (x)⟩ = Z Under Eq. (8) the temperature is Teff ∝ ∆E/d, and the variance of the coefficients becomes 2∆E/(λd), shrinking with every added parameter. The principle of least action states that when T → 0, the primary contribution to Z becomes that of the stationary path, where Lλ (θ) is minimized, and thus the exponential is maximized (Feynman & Hibbs (1965)). In general, the Boltzmann distribution is equivalent to the overdamped Langevin equation √ θ̇ = −∇L(θ) − λθ + 2T ξ(t),
(10)
with ξ white noise of unit strength: gradient descent on the loss with weight decay, driven by isotropic noise of temperature T , the stochastic gradient Langevin dynamics of Welling & Teh (2011). The density P (θ, t) of a cloud of runs obeys the Fokker–Planck equation (Risken, 1996), in which the drift −∇Lλ transports probability and the noise diffuses it with coefficient T , ∂t P = ∇ · P ∇Lλ + T ∇2 P = ∇ · J, J = P ∇Lλ + T ∇P. (11) A stationary density with no probability current, J = 0, satisfies T ∇P = −P ∇Lλ , that is ∇ ln P = −
∇Lλ , T
P (θ) =
so
1 −Lλ (θ)/T e , Z
(12)
which is Eq. (3), and it is normalizable exactly when Z of Eq. (9) converges: for d > n this requires λ > 0, since without the weight decay the loss is flat along the d − n unseen directions. For a quadratic loss, e−Lλ /T = e−Emin /T exp[− 21 (θ − θ ∗ )⊤ (H/T )(θ − θ ∗ )] is the Gaussian of mean θ ∗ and covariance T H −1 , with Z = (2πT )d/2 (det H)−1/2 e−Emin /T . For a model linear in its parameters, f (x) = p(x)⊤ θ with p(x) the vector of basis functions at x, the curve is a Gaussian process with mean p(x)⊤ θ ∗ and covariance T p(x)⊤ H −1 p(x′ ) (Rasmussen & Williams, 2006). For least squares with weight decay, H = M and the mean curve is the ridge regression fit f ∗ (x) = p(x)⊤ θ ∗ . Its two limits are the two ends of the path integral. As T → 0 the noise vanishes and the parameters descend to the stationary point ∇L(θ)+λθ = 0. For least squares, L(θ) = 12 ∥Aθ − y∥2 = 12 (Aθ − y)⊤ (Aθ − y) has gradient ∇L(θ) = A⊤ (Aθ − y), so the stationary point solves A⊤ (Aθ − y) + λθ = 0,
that is
(A⊤ A + λI) θ = A⊤ y,
θ = θ ∗ = M −1 A⊤ y,
(13)
mirroring Eq. (6). It is a minimum because the Hessian M = A⊤ A + λI is positive definite for λ > 0. p As T → ∞, we perform a change of variables θ = T /λ η; for least squares ∇L(θ) = A⊤ Aθ − A⊤ y, and Eq. (10) becomes p √ √ η̇ = −(A⊤ A + λI)η + λ/T A⊤ y + 2λ ξ(t) −−−−→ −M η + 2λ ξ(t). (14) T →∞
The T → ∞ limit Eq. (14) is the Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein, 1930), resulting in √ Rt the motion θ(t) − θ ∗ = e−Ht (θ 0 − θ ∗ ) + 2T 0 e−H(t−s) dWs . This solution is Gaussian at every t, with Rt mean θ ∗ + e−Ht (θ 0 − θ ∗ ) and covariance 2T 0 e−2Hv dv = T H −1 (I − e−2Ht ), which tends to the stationary 5
law as t → ∞ (Pavliotis, 2014, Chapter 4). For the value at x, with p̃i (x) the components of p(x) along the eigenvectors of M and µi its eigenvalues, ⟨f (x)2 ⟩t = ⟨f (x)⟩2t + T
X i
p̃i (x)2
1 − e−2µi t , µi
⟨f (x)⟩t = f ∗ (x) + p(x)⊤ e−M t (θ 0 − θ ∗ ).
(15)
At stationarity ⟨f (x)2 ⟩ = f ∗ (x)2 + T p(x)⊤ M −1 p(x). At a finite time the unseen directions, µi = λ, have reached only T (1 − e−2λt )/λ ≈ 2T t of their stationary variance T /λ: the run behaves as if the stiffness were capped at 1/(2t), the λ ∝ 1/t of Eq. (2). Theorem 1 (Finite training time is weight decay, Appendix C). For least squares with weight decay λ ≥ 0, the run of Eq. (10) from θ 0 = 0 has, at time t, 1 Cov[θ(t)] ⪯ 2T M + 2t I
−1
E[θ(t)] = (I − e−M t )M −1 A⊤ y,
,
and the mean agrees with the ridge fit (M + 1t I)−1 A⊤ y at weight decay λ + 1/t along every eigenvector of M up to a factor between 1 and 1.3. By Theorem 1, a run of length t has the fluctuations of a run with weight decay λ + 1/(2t), up to a factor two, and the mean fit of a run with weight decay λ + 1/t, up to a factor 1.3 (Ali et al., 2019); the finite training time caps the size of the fit whatever the conditioning of the data. Late stopping removes the cap and leaves the ridge fit at the explicit λ, whose quality depends on the conditioning of the data: the collapse of σmin (A) at d = n in Fig. 1G and its recovery beyond. 2.5
Measuring the size of the fitted functions R1 2 The Legendre polynomials are orthogonal, −1 Pk Pl dx = 2k+1 δkl (Szegő, 1975, Eq. (4.3.3) with α = β = 0), so that Z 1 X 2 f 2 dx = θk2 ≤ 2∥θ∥2 , (16) 2k + 1 −1 k<d
therefore the least-action curve through the data is the smallest one in the L2 sense. For d ≥ n and distinct inputs, AA⊤ is invertible and as λ → 0 the least-action curve is the interpolant of least norm, θ ∗ = A⊤ (AA⊤ )−1 y (Theorem B.1 in Appendix B), whose energy 12 ∥θ ∗ ∥2 = 12 y ⊤ (AA⊤ )−1 y changes with every added coefficient. We now show that this energy is nonincreasing. Theorem 2 (The energy of the interpolant, Appendix C). Let d ≥ n, let pd = (Pd (xi ))i≤n be the values at the inputs of the basis function added when d becomes d + 1, and let αd = (AA⊤ )−1 y, so that θ ∗ = A⊤ αd . Then 2 (p⊤ d αd ) ∥θ ∗(d+1) ∥2 = ∥θ ∗(d) ∥2 − : (17) ⊤ ⊤ 1 + pd (AA )−1 pd the energy of the interpolant is nonincreasing in d, and strictly decreasing unless p⊤ d αd = 0, which for given inputs holds only on a hyperplane of labels y. Thus, each new basis function with nonzero overlap with the training data lowers the energy.
3
Conclusion
The training trajectory produced by stochastic gradient-based methods can be thought of as a particle wandering over the energy landscape with some induced temperature T , and so can be treated with statistical mechanics. Because real world machine learning algorithms start at some initial point and have only finite time to diffuse in the energy landscape, there is effectively weight decay regularization. By increasing the degrees of freedom in this landscape, the equipartition theorem distributes energy among the quadratic coordinates in multiples of T2 . Increasing the number of degrees of freedom d (which are made quadratic due to weight decay regularization) thus effectively decreases the temperature for a given energy E, driving the Boltzmann distribution toward the stationary path. Furthermore, increasing d can only lower the energy of 6
the stationary path (Theorem 2) and thus the resulting L2 norm Eq. (16), therefore reducing the likelihood at fixed E that the sampled solution has large L2 norm. Thus, increasing d beyond the interpolation threshold is a form of weight decay regularization as well.
AI disclosure The author used Claude (Anthropic; Claude Opus 5 through Claude Code) as a research and writing assistant. Under the author’s direction it drafted and revised text, derived and typeset the theorems and their proofs, wrote the Lean 4 verification of the theorems, wrote the Python scripts that produce every figure, and ran the numerical experiments. The author’s initial idea was to explore double descent through the lens of statistical mechanics, in particular to bound the norms of solutions near local minima, as a function of d. The author noted that there is symmetry among minima for neural networks by permuting neurons in a hidden layer. From there, the author was reminded of the spontaneous symmetry breaking in physics, of the path integral, and finally of the connection between the path integral and the Boltzmann distribution. Conflict of interest statement The author is a 2026–2027 Doximity AI Fellow, and receives items and services of minimal monetary value from Doximity, Inc in exchange for consulting services.
References Alnur Ali, J Zico Kolter, and Ryan J Tibshirani. A continuous-time view of early stopping for least squares regression. In Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS), volume 89 of Proceedings of Machine Learning Research, pp. 1370–1378, 2019. URL https://proceedings.mlr.press/v89/ali19a.html. Peter L Bartlett, Philip M Long, Gábor Lugosi, and Alexander Tsigler. Benign overfitting in linear regression. Proceedings of the National Academy of Sciences, 117(48):30063–30070, 2020. doi: 10.1073/ pnas.1907378117. Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine-learning practice and the classical bias–variance trade-off. Proceedings of the National Academy of Sciences, 116(32): 15849–15854, 2019. doi: 10.1073/pnas.1903070116. Richard P Feynman and Albert R Hibbs. Quantum Mechanics and Path Integrals. McGraw-Hill, New York, 1965. Mario Geiger, Arthur Jacot, Stefano Spigler, Franck Gabriel, Levent Sagun, Stéphane d’Ascoli, Giulio Biroli, Clément Hongler, and Matthieu Wyart. Scaling description of generalization with number of parameters in deep learning. Journal of Statistical Mechanics: Theory and Experiment, 2020(2):023401, 2020. doi: 10.1088/1742-5468/ab633c. Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. Annals of Statistics, 50(2):949–986, 2022. doi: 10.1214/21-AOS2133. Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations (ICLR), 2015. arXiv:1412.6980. Anders Krogh and John A Hertz. A simple weight decay can improve generalization. In Advances in Neural Information Processing Systems, volume 4, pp. 950–957. Morgan Kaufmann, 1992. URL https: //proceedings.neurips.cc/paper/1991/hash/8eefcfdf5990e441f0fb6f3fad709e21-Abstract.html. Chan Li and Nigel Goldenfeld. Broken ergodicity and the violation of the fluctuation-dissipation theorem lead to generalization beyond overfitting in machine learning, 2026. arXiv:2607.04135. Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Conference on Learning Representations (ICLR), 2019. arXiv:1711.05101. 7
Song Mei and Andrea Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75(4):667–766, 2022. doi: 10.1002/cpa.22008. Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever. Deep double descent. Windows on Theory blog, https://windowsontheory.org/2019/12/05/deep-double-descent/, December 2019. Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever. Deep double descent: Where bigger models and more data hurt. In International Conference on Learning Representations (ICLR), 2020. arXiv:1912.02292. Grigorios A Pavliotis. Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations, volume 60 of Texts in Applied Mathematics. Springer, New York, 2014. doi: 10.1007/978-1-4939-1323-7. Carl Edward Rasmussen and Christopher K I Williams. Gaussian Processes for Machine Learning. MIT Press, 2006. doi: 10.7551/mitpress/3206.001.0001. Frederick Reif. Fundamentals of Statistical and Thermal Physics. McGraw-Hill, New York, 1965. Hannes Risken. The Fokker–Planck Equation: Methods of Solution and Applications. Springer, Berlin, 2nd edition, 1996. doi: 10.1007/978-3-642-61544-3. Herbert Robbins and Sutton Monro. A stochastic approximation method. Annals of Mathematical Statistics, 22(3):400–407, 1951. doi: 10.1214/aoms/1177729586. Jack Sherman and Winifred J Morrison. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics, 21(1):124–127, 1950. doi: 10.1214/aoms/1177729893. Stefano Spigler, Mario Geiger, Stéphane d’Ascoli, Levent Sagun, Giulio Biroli, and Matthieu Wyart. A jamming transition from under- to over-parametrization affects generalization in deep learning. Journal of Physics A: Mathematical and Theoretical, 52(47):474001, 2019. doi: 10.1088/1751-8121/ab4c8b. Gábor Szegő. Orthogonal Polynomials, volume 23 of Colloquium Publications. American Mathematical Society, 4th edition, 1975. doi: 10.1090/coll/023. George E Uhlenbeck and Leonard S Ornstein. On the theory of the Brownian motion. Physical Review, 36 (5):823–841, 1930. doi: 10.1103/PhysRev.36.823. Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th International Conference on Machine Learning (ICML), pp. 681–688, 2011. URL https: //icml.cc/2011/papers/398_icmlpaper.pdf.
A
Worked examples and derivations for Section 2.1
A.1
Coin flips
With the fundamental postulate, one may count physical states of a collection of objects, given macrostates defined by a macroscopic observable, such as the total energy E or the average pressure P , and derive laws regarding how these quantities relate with each other, for example the ideal gas law P V = N kB T (Reif, 1965). With the ergodicity assumption of Section 2.2 in place, there is an analogy between the coin flips and model training (Table 1). 8
Ten coin flips
Machine learning model M A specific vector of parameters θ
Multiplicity Ω
One sequence of outcomes, e.g. HTHTHHHTHH Total number of heads All sequences with the observed number of heads, e.g. the 120 sequences with 7 heads 10 = 120 7
Fundamental postulate
Each of the 120 sequences has probability 1/120
Microstate Observable Macrostate
Training loss L(M (θ)) All parameter vectors that achieve the observed loss of E: Θ = {θ : L(M (θ)) = E} Number/density of the set of parameters with loss E All θ ∈ Θ are equally likely to have been reached once training has converged
Table 1: The statistical mechanics analogy between coin flips and machine learning models. A.2
The Boltzmann distribution
Consider a small system S in contact with a large reservoir R, the two together isolated with total energy Etot and able to exchange energy. A microstate of the whole is a pair (microstate of S, microstate of R), and the postulate weights every such pair with energy Etot equally. Fix one microstate s of the system, of energy Es . The number of pairs in which S is in state s is the number of reservoir microstates of energy Etot − Es , that is, the multiplicity ΩR (Etot − Es ). Hence P (s) ∝ ΩR (Etot − Es ) = exp SR (Etot − Es ) . (18) The system is small, so Es ≪ Etot , and the exponent can be expanded to first order: SR (Etot − Es ) = SR (Etot ) − Es ∂SR /∂E + . . . . The first term is a constant. The derivative ∂SR /∂E is a property of the reservoir alone, and we call its reciprocal the temperature T (again ignoring the historical factor of kB ). The result is Eq. (1). With the training loss as the energy, Eq. (3), a model that fits the data less well is not forbidden, only exponentially less probable, by a factor e−∆L/T per unit of extra loss. If, in addition, the coefficients are penalized by λ2 ∥θ∥2 , the penalty adds to the energy, and the weight becomes exp[−Lλ (θ)/T ] with Lλ = L + λ2 ∥θ∥2 ; the factor exp[−λ∥θ∥2 /2T ] is a Gaussian prior of variance s2 = T /λ on each coefficient. Zero temperature, T → 0, concentrates the weight on the minimizers of Lλ . Zero penalty, λ → 0, recovers the initial uniform weighting of microstates. A.3
Equipartition
Pd In a neighborhood about a local energy minimum Emin , the landscape is quadratic, L(θ) = Emin + 12 i=1 µi u2i , in coordinates u along its principal curvatures µi > 0 (Eq. (6) makes this exact for least squares). Under 2 the Boltzmann weight e−L/T the coordinates are independent, each with density ∝ e−µi ui /2T , a Gaussian of 1 variance T /µi , so each term of the sum has average 2 µi · T /µi = T /2 whatever its curvature, and E[L] = Emin +
d X T i=1
2
= Emin +
Td , 2
which is Eq. (7).
B
The Legendre fit in each regime
The fits of Fig. 1 and the matrix A of Section 2.3 are defined as follows. We work in the Legendre basis P0 , . . . , Pd−1 , the polynomials of Rodrigues’ formula (Szegő, 1975, Eq. (4.3.1) with α = β = 0) Pk (x) =
1
dk
2k k! dxk
k x2 − 1 ,
P0 = 1,
P1 = x, 9
P2 = 12 (3x2 − 1),
P3 = 21 (5x3 − 3x), . . . ,
(19)
R1 2 so that Pk has degree exactly k and they are orthogonal on [−1, 1], −1 Pk Pl dx = 2k+1 δkl .4 Every Legendre polynomial is bounded by one on the whole P interval, |Pk (x)| ≤ 1 with Pk (1) = 1, so the columns of A have entries of comparable size. We write f (x) = k θk Pk (x) = p(x)⊤ θ, so that the value f (x) at any input is a linear observable of the microstate θ. The training constraints are Aθ = y, one row per training point and one column per basis function, 1 x1 12 (3x21 − 1) · · · P0 (x1 ) P1 (x1 ) · · · Pd−1 (x1 ) y1 θ0 P (x ) P (x ) · · · P (x ) 1 x 1 2 y2 θ1 1 2 d−1 2 2 0 2 2 (3x2 − 1) · · · , = θ = . , y = . , A= .. .. .. .. .. .. . . . .. . . . . . . . . . . . yn θd−1 P0 (xn ) P1 (xn ) · · · Pd−1 (xn ) 1 xn 12 (3x2n − 1) · · · (20) an n × d matrix, the fitting matrix, whose ith row is p(xi )⊤ , so that (Aθ)i = f (xi ). Theorem B.1 (The fit in each regime). Let the n training inputs be distinct. Then A has rank min(d, n), and the fit that training aims at is linear in the labels, θ = W y for a d × n matrix W that depends only on the inputs; in each regime: d < n:
W = (A⊤ A)−1 A⊤ , least squares. The macrostate M at E = 0 is empty; the training loss has the unique minimizer θ = W y, with Emin = 12 ∥(I − AW )y∥2 > 0 unless y lies in the column space of A;
d = n:
W = A−1 . A is invertible and M = {W y}, the unique interpolant known as the Lagrange polynomial through the n points;
d > n:
W = A⊤ (AA⊤ )−1 . M = θ ∗ + ker A is an affine subspace of dimension d − n, and its member of least norm is θ ∗ = W y, the minimum-norm interpolant, which is where gradient descent from θ = 0 converges.
Proof. The rank. Restrict A to its first min(d, n) columns and the first min(d, n) inputs; this square block is Vik = Pk (xi ). Since Pk has degree exactly k, V = W T with Wik = xki the Vandermonde matrix and T upper triangular with nonzero diagonal, and Y det W = (xj − xi ) ̸= 0 i<j
for distinct inputs. So V is invertible and A has rank min(d, n). d < n. The columns of A are independent, so A⊤ A is invertible, and the range of A is a proper subspace of Rn : Aθ = y has no solution unless y lies in the range, and the macrostate at E = 0 is empty. The strictly convex quadratic 21 ∥Aθ − y∥2 has the unique stationary point A⊤ Aθ = A⊤ y,
θ = (A⊤ A)−1 A⊤ y,
whose residual y − A(A⊤ A)−1 A⊤ y is the projection of y off the range, so Emin is as stated. d = n. A = V is invertible and θ = A−1 y is the one solution, the Lagrange interpolant. d > n. The rows of A are independent, so AA⊤ is invertible, and θ ∗ = A⊤ (AA⊤ )−1 y satisfies Aθ ∗ = y. The solution set is the coset θ ∗ + ker A with dim ker A = d − n. Since θ ∗ lies in the row space of A, which is orthogonal to ker A, it is the solution of least norm. Gradient descent from θ = 0 moves only within the row space, since every gradient A⊤ (Aθ − y) lies there, and converges to it. The minimal solution of Fig. 1 is the fit of this theorem at λ = 0: least squares for d < n, the interpolant for d = n, and the minimum-norm interpolant for d > n. With weight decay, the three formulas merge into the single one of Eq. (6), θ ∗ = (A⊤ A + λI)−1 A⊤ y = A⊤ (AA⊤ + λI)−1 y, which tends to the corresponding case as λ → 0. 4 The monomial, or Taylor, basis {1, x, x2 , . . . } spans the same polynomials but is numerically unstable: the corresponding A is the Vandermonde matrix, whose condition number grows exponentially with d, because the monomials of high degree are all nearly zero on most of the interval and nearly equal to one another near x = 1 and x = −1.
10
C
Proofs
The theorems below, Theorem B.1 and the identity of Eq. (6) are also verified in Lean 4 with Mathlib; the source files accompany the arXiv submission as ancillary files (anc/lean/), and their README.md lists which statement each file proves. For Theorem 1 the Ornstein–Uhlenbeck solution is taken as given and the three scalar inequalities of the proof, including the constant 1.3, are what is verified. Theorem 1 (Finite training time is weight decay). For least squares with weight decay λ ≥ 0, the run of Eq. (10) from θ 0 = 0 has, at time t, 1 Cov[θ(t)] ⪯ 2T M + 2t I
−1
E[θ(t)] = (I − e−M t )M −1 A⊤ y,
,
and the mean agrees with the ridge fit (M + 1t I)−1 A⊤ y at weight decay λ + 1/t along every eigenvector of M up to a factor between 1 and 1.3. Proof. Both statements are diagonal in the eigenbasis of M ; let µ ≥ 0 be an eigenvalue. By the covariance T M −1 (I − e−2M t ) of the Ornstein–Uhlenbeck solution the variance along its eigenvector is T (1 − e−2µt )/µ (equal to 2T t when µ = 0), and with x = 2µt, 1 − e−x 2 2 1 − e−2µt = 2t ≤ 2t = , µ x 1+x µ + 1/(2t) because (1 − e−x )(1+ x) ≤ 2x: for x ≤ 1, 1 − e−x ≤ x and x(1+ x) ≤ 2x; for x ≥ 1, 1 − e−x ≤ 1 and 1+ x ≤ 2x. This is the covariance bound. For the mean, the solution with θ 0 = 0 gives E[θ(t)] = (I − e−M t )θ ∗ = (I − e−M t )M −1 A⊤ y, whose coefficient along the eigenvector is (1 − e−µt )/µ times that of A⊤ y, against 1/(µ + 1/t) for the ridge fit at λ + 1/t. Their ratio is r(x) = (1 − e−x )(1 + x)/x with x = µt, which satisfies 1 ≤ r(x) ≤ 1.3 for x ≥ 0: r → 1 at both ends, r ≥ 1 since (1 − e−x )(1 + x) ≥ x, and its single maximum, at x ≈ 1.79, is 1.299. Theorem 2 (The energy of the interpolant). Let d ≥ n, let pd = (Pd (xi ))i≤n be the values at the inputs of the basis function added when d becomes d + 1, and let αd = (AA⊤ )−1 y, so that θ ∗ = A⊤ αd . Then ∥θ ∗(d+1) ∥2 = ∥θ ∗(d) ∥2 −
2 (p⊤ d αd ) : ⊤ ⊤ 1 + pd (AA )−1 pd
(17)
the energy of the interpolant is nonincreasing in d, and strictly decreasing unless p⊤ d αd = 0, which for given inputs holds only on a hyperplane of labels y. ∗(d) −1 Proof. Write Ad for the fitting matrix with d coefficients and K = Ad A⊤ = A⊤ y, ∥θ ∗(d) ∥2 = d . Since θ dK ⊤ −1 ⊤ −1 ⊤ −1 ⊤ ⊤ y K Ad Ad K y = y K y. Adding the column pd gives Ad+1 Ad+1 = K + pd pd , and the Sherman– Morrison formula (Sherman & Morrison, 1950) −1 (K + pd p⊤ = K −1 − d)
−1 K −1 pd p⊤ dK −1 p 1 + p⊤ d dK
−1 −1 2 −1 gives y ⊤ (K + pd p⊤ y = y ⊤ K −1 y − (p⊤ y) /(1 + p⊤ pd ), which is Eq. (17) with αd = K −1 y. The d) dK dK ⊤ subtracted term is nonnegative, and it vanishes only when pd αd = (K −1 pd )⊤ y = 0, a linear condition on y whose coefficient vector K −1 pd is nonzero whenever pd ̸= 0.
11