# A Dynamical Model of Neural Scaling Laws

Blake Bordelon  Alexander Atanasov  Cengiz Pehlevan

###### Abstract

On a variety of tasks, the performance of neural networks predictably improves
with training time, dataset size and model size across many orders of
magnitude. This phenomenon is known as a neural scaling law. Of fundamental
importance is the compute-optimal scaling law, which reports the performance
as a function of units of compute when choosing model sizes optimally. We
analyze a random feature model trained with gradient descent as a solvable
model of network training and generalization. This reproduces many
observations about neural scaling laws. First, our model makes a prediction
about why the scaling of performance with training time and with model size
have different power law exponents. Consequently, the theory predicts an
asymmetric compute-optimal scaling rule where the number of training steps are
increased faster than model parameters, consistent with recent empirical
observations. Second, it has been observed that early in training, networks
converge to their infinite-width dynamics at a rate $1/\textit{width}$ but at
late time exhibit a rate $\textit{width}^{-c}$, where $c$ depends on the
structure of the architecture and task. We show that our model exhibits this
behavior. Lastly, our theory shows how the gap between training and test loss
can gradually build up over time due to repeated reuse of data.

Scaling Laws, Mean Field Theory, Infinite Width Limits, Deep Ensembles

##  1 Introduction

(a) One-Pass CIFAR-5M Test Dynamics  (b) CIFAR-5M Compute Scaling

(c) Wikitext-100M Train and Test (d) Wikitext-5M Train and Test

Figure 1: Train and test losses (cross-entropy) as a function of training time
$t$ and width $N$. For models trained online, we do not make a distinction
between training and test error because each new batch is drawn fresh and
would have the same loss in expectation as an independent test set. (a) The
test loss of a residual CNN on CIFAR-5M is well described by a fit of the form
$\mathcal{L}\sim t^{-\alpha}+N^{-\beta}$ in the online training regime. (b)
The compute optimal strategy requires scaling up both model size and training
time simultaneously. (c) Transformer training on wikitext with 100M tokens
before data-repetition. Model performance is monotonic in width $N$. (d)
Wikitext with 5M subsampled tokens. Larger width $N$ is not always better as
wider models can overfit.

Large scale language and vision models have been shown to achieve better
performance as the number of parameters and number of training steps are
increased. Moreover, the scaling of various loss metrics (such as cross
entropy or MSE test loss) has been empirically observed to exhibit remarkably
regular, often power law behavior across several orders of magnitude (Hestness
et al., 2017; Kaplan et al., 2020). These findings are termed “neural scaling
laws”.

Neural scaling laws play a central role in modern deep learning practice, and
have substantial implications for the optimal trade-off between model size and
training time (Hoffmann et al., 2022), as well as architecture selection
(Alabdulmohsin et al., 2023). Understanding the origin of such scaling laws,
as well as their exponents, has the potential to offer insight into better
architectures, the design of better datasets (Sorscher et al., 2022), and the
failure modes and limitations of deep learning systems. Yet, many questions
about neural scaling laws remain open.

In this paper, we introduce and analyze a solvable model which captures many
important aspects of neural scaling laws. In particular, we are interested in
understanding the following empirically observed phenomena:

##### Test Loss Scales as a Power-law in Training Time and Model Size and
Compute.

In many domains of deep learning, the test loss of a model with $N$ trainable
parameters trained for $t$ iterations has been found to scale as
$\mathcal{L}(t,N,P)\approx\mathcal{L}_{0}+a_{t}t^{-r_{t}}+a_{N}N^{-r_{N}}$
(Kaplan et al., 2020; Hoffmann et al., 2022). These scaling law exponents
$r_{t},r_{N}$ generally depend on the dataset and architecture. We demonstrate
scaling laws on simple vision and language tasks in Figure 1. The compute is
proportional to the number of steps of gradient descent times the model size
$C\propto Nt$. Setting $N$ and $t$ optimally gives that test loss scales as a
power law in $C$. This is the compute optimal scaling law.

##### Compute-Optimal Training Time and Model Size Scaling Exponents Are
Different.

Hoffmann et al. (2022) observed that $r_{t}$ and $r_{N}$ are close but
slightly different, leading to asymmetric compute-optimal scaling of
parameters. For compute budget $C$, they scale model size $N\propto C^{c_{1}}$
and training time $t\propto C^{c_{2}}$ with $c_{2}>c_{1}$. This difference in
exponents led to a change in the scaling rule for large language models,
generating large performance gains.

##### Larger Models Train Faster.

Provided feature learning is held constant across model scales (i.e. adopting
mean-field or $\mu$P scaling), wider networks tend to train faster (Yang et
al., 2021) (Figure 1). If training proceeds in an online/one-pass setting
where datapoints are not repeated, then the wider models will also obtain
lower test loss at an equal number of iterations. This observation has been
found to hold both in overparameterized and underparameterized regimes
(Bordelon & Pehlevan, 2023; Vyas et al., 2023).

##### Models Accumulate Finite-Dataset and Finite-Width Corrections.

Early training can be well described by the learning curves for stochastic
gradient descent without reuse of samples (termed the online/ideal limiting
dynamics), however over time the effect of reusing data accumulates and leads
to worse test performance (Nakkiran et al., 2021b; Mignacco et al., 2020;
Ghosh et al., 2022). Similarly the gaps in model performance across various
model sizes also grow with training time (Yang et al., 2021; Vyas et al.,
2023). Figure 1 (d) shows overfitting and reversal of “wider is better”
phenomenon due to data reuse.

##### Scaling Exponents are Task-Dependent at Late Training Time, but not at
Early Time.

Prior works (Dyer & Gur-Ari, 2020; Atanasov et al., 2023; Roberts et al.,
2022; Bordelon & Pehlevan, 2023) predict early-time finite-width loss
corrections that go as $1/\textit{width}$ near the infinite width limit in
either lazy or feature-learning regimes. Bahri et al. (2021) et al provide
experiments demonstrating the $1/\textit{width}$ convergence. However, finite-
width models trained for a long time exhibit non-trivial exponents with
respect to model width (Kaplan et al., 2020; Vyas et al., 2023). See Figure 1
for examples of nontrivial scalings at late time on CIFAR-5M and Wikitext.

##### Ensembling is Not the Same as Going Wider.

Near the limit of infinite width, finite models can be thought of as noisy
approximations of the infinite-width model with noise that can be eliminated
through ensembling (Dyer & Gur-Ari, 2020; Geiger et al., 2020; Atanasov et
al., 2023). However recent experiments (Vyas et al., 2023) indicate that
ensembling is not enough to match performance of larger models.

These phenomena are not unique to deep networks, but can be observed in linear
models, or linearized neural networks operating in the lazy/kernel regime.
Though this regime does not capture feature learning, it has benefit of
analytical tractability. In this paper, we focus on such linearized models to
attempt to gain insight into the dynamics of training.

To attempt to explain these phenomena, we develop a mathematically tractable
model of neural scaling laws which allows one to simultaneously vary time,
model size, and dataset size. Our contributions are as follows:

  1. 1.

We analyze the learning dynamics of a structured and randomly projected linear
model trained with gradient descent or momentum. In an asymptotic limit of the
model, we obtain a dynamical mean field theory (DMFT) description of the
learning curve in terms of correlation functions, which measure the cross-time
correlation of training and test errors, and response functions which measure
sensitivity of the dynamics to small perturbations.

  2. 2.

We solve for the response functions exactly in Fourier domain. This solution
reveals faster training for larger models. The low frequency range of these
functions allow us to extract the long time limit of the loss.

  3. 3.

We show that the model and data corrections to the dynamics accumulate over
time. At early time, each of these corrections has a universal scaling,
consistent with prior works (Bahri et al., 2021).

  4. 4.

For power-law structured features we show that the model exhibits power law
scaling of test loss with time, model size and dataset size. While the data
and model exponents are the same, the time and model exponents are different
in general. We show that this gives rise to an asymmetric compute optimal
scaling strategy where training time increases faster than model size.

  5. 5.

Our theory explains why ensembling is not compute optimal as it gives less
benefit to performance than increase in model size.

  6. 6.

We observe in Section 5.1 that feature learning networks can obtain better
power law scalings, leading to a better compute optimal frontier. We
empirically study this phenomenon in Appendix K.

###  1.1 Related Works

The learning curves for linear models with structured (non-isotropic)
covariates, including infinite-width kernel regression, have been computed
using tools from statistical physics and random matrix theory (Bordelon et
al., 2020; Spigler et al., 2020; Canatar et al., 2021; Simon et al., 2021;
Bahri et al., 2021; Hastie et al., 2022). Mei & Montanari (2022) analyzed a
linear model with random projections of isotropic covariates. There, they
study the limiting effects of width and dataset size, and observe model-wise
and sample-wise double descent. In (Adlam & Pennington, 2020a) a related model
is used to study the finite-width neural tangent kernel (NTK) (Jacot et al.,
2018) of a given network. Further, (Adlam & Pennington, 2020b) extends this
analysis to understand the different sources of variance in the predictions of
random feature models and the effect of ensembling and bagging on the test
loss. Other works have extended this to models where an additional untrained
projection is applied to the structured covariates (Loureiro et al., 2021,
2022; Zavatone-Veth et al., 2022; Atanasov et al., 2023; Maloney et al., 2022;
Zavatone-Veth & Pehlevan, 2023; Ruben & Pehlevan, 2023; Simon et al., 2023).
Within this literature, which considered fully trained models, the works of
(Bordelon et al., 2020; Spigler et al., 2020) derived power-law decay rates
for power-law features which were termed resolution limited by (Bahri et al.,
2021) and recovered by (Maloney et al., 2022).

However, we also study the dependence on training time. The $t\to\infty$ limit
of our DMFT equations recovers the final losses computed in these prior works.
While these prior works find that the scaling exponents for model-size and
dataset-size are the same, we find that the test loss scales with a different
exponent with training time, leading to a different (model and task dependent)
compute optimal scaling strategy.

DMFT methods have been used to analyze the test loss dynamics for general
linear and spiked tensor models trained with high-dimensional random data
(Mannelli et al., 2019; Mignacco et al., 2020; Mignacco & Urbani, 2022) and
deep networks dynamics with random initialization (Bordelon & Pehlevan, 2022b;
Bordelon et al., 2023). High dimensional limits of SGD have been analyzed with
Volterra integral equations in the offline case (Paquette et al., 2021) or
with recursive matrix equations in the online case (Varre et al., 2021;
Bordelon & Pehlevan, 2022a). Random matrix approaches have also been used to
study test loss dynamics in linear regression with isotropic covariates by
(Advani et al., 2020). In this work, we consider averaging over both the
disorder in the sampled dataset and the random projection of the features
simultaneously.

Other models and hypotheses for scaling laws instead rely on a discrete
collection of subtasks or skills which are learned as compute grows (Caballero
et al., 2022; Arora & Goyal, 2023; Michaud et al., 2023). Our theory instead
focuses on spectral components of a data distribution.

##  2 Setup of the Model

We consider a “teacher-student” setting, where data sampled from a generative
teacher model is used to train a student random feature model. The teacher and
student models mismatch in a particular way that will be described below. This
mismatch is the key ingredient that leads to most of the phenomena that we
will discuss.

##### Teacher Model.

Take $\bm{x}\in\mathbb{R}^{D}$ to be drawn from a distribution $p(\bm{x})$
with a target function $y(\bm{x})$ expressible in terms of base features
$\bm{\psi}(\bm{x})\in\mathbb{R}^{M}$ up to noise:

$y(\bm{x})=\frac{1}{\sqrt{M}}\bm{w}^{\star}\cdot\bm{\psi}(\bm{x})+\sigma\epsilon(\bm{x}).$
(1)

Here $\psi(\bm{x})$ play the role of the infinite-width NTK eigenfunctions,
which form a complete basis for square-integrable functions $L^{2}[p]$. The
$\epsilon(\bm{x})$ function describes a component of $y$ with which is
uncorrelated with $\psi(\bm{x})$. We work in the eigenbasis of features as in
(Bordelon et al., 2020), so the covariance given by:

$\displaystyle\left<\psi_{k}(\bm{x})\psi_{\ell}(\bm{x})\right>_{\bm{x}\sim
p(\bm{x})}=\delta_{k\ell}\lambda_{k}.$ (2)

The power law structure in the $\lambda_{k}$ and $\bm{w}^{*}$ entries will
lead to power law scalings for the test loss and related quantities.

##### Student Model.

Our student model is motivated by a scenario where a randomly initialized
finite-width network is trained in the linearized or lazy regime (Chizat et
al., 2019; Jacot et al., 2018). Such training can be described through
learning linear combinations of the finite-width NTK features. These features
will span a lower-dimensional subspace of the space of square-integrable
functions, and relate to infinite-width NTK features in a complicated way.

To model this key aspect, the student model uses a projection of the
$\bm{\psi}(\bm{x})$ features, $\frac{1}{\sqrt{M}}\bm{A}\bm{\psi}(\bm{x})$
where $\bm{A}\in\mathbb{R}^{N\times M}$. These projected features represent
the finite-width (i.e. empirical) NTK eigenfunctions. This is motivated by the
fact that finite width kernel’s features can be linearly expanded in the basis
of infinite-width features, because infinite-width kernel eigenfunctions are
complete.

Our learned function then has the form:

$\displaystyle
f(\bm{x})=\frac{1}{\sqrt{NM}}\bm{w}\cdot\bm{A}\bm{\psi}(\bm{x}).$ (3)

Here, we will interpret $N$ as the model size with the $N\to\infty$ limit
recovering original kernel. Similar models were studied in (Maloney et al.,
2022; Atanasov et al., 2023).

We will focus on the setting where the elements of $\bm{A}$ are drawn iid from
a distribution of mean zero and variance one. See Appendix B for details on
the technical assumptions. The motivations for this choice are (1)
tractability and (2) it satisfies the constraint that as $N\to\infty$ the
student’s kernel approaches the infinite-width kernel $\bm{\psi}$. In more
realistic settings, such as when projecting the eigenfunctions of an infinite-
width NTK to a finite-width NTK, the form of the $\bm{A}$ matrix is generally
not known.

##### Training.

The model is trained on a random dataset
$\mathcal{D}=\\{\bm{x}_{\mu},y_{\mu}\\}_{\mu=1}^{P}$ of size $P$ with gradient
flow on MSE loss

$\displaystyle\frac{\partial}{\partial
t}\bm{w}(t)=\sum_{\mu=1}^{P}(y(\bm{x}_{\mu})-f(\bm{x}_{\mu},t))\bm{A}\bm{\psi}(\bm{x}_{\mu}).$
(4)

We explore extensions (momentum, discrete time, one-pass SGD in Appendix J).
We track the test and train loss

$\displaystyle\mathcal{L}(t)$
$\displaystyle=\mathbb{E}_{\bm{x}}[(f(\bm{x},t)-y(\bm{x}))^{2}],$ (5)
$\displaystyle\hat{\mathcal{L}}(t)$
$\displaystyle=\frac{1}{P}\sum_{\mu=1}^{P}(f_{\mu}(t)-y_{\mu})^{2}.$

In small size systems, these losses depend on the precise realization of the
data $\mathcal{D}$ and matrix $\bm{A}$. These two quantities can be viewed as
the disorder. For large systems, these losses approach a well-defined limit
independent of the specific realization of $\mathcal{D},\bm{A}$. We will use
this fact in the next section when analyzing the model.

##  3 DMFT for Scaling Laws

We next describe a theoretical approach for characterizing the learning curves
for this model. The full details of this approach is detailed in Appendices A,
B.

We derive a mean field theory for $M,N,P$ large. We analyze both the (1)
proportional regime where $N/M=\nu,P/M=\alpha$ and $M,N,P\to\infty$, and (2)
non-proportional regime where $M\to\infty$ first and $N,P\gg 1$. The theories
derived in these limits are structurally identical (App. F). 111While the
proportional limit is exact, the finite size $N,P$ theory will also contain
variability across realizations of disorder. When relevant, we show these in
experiments by plotting standard deviations over draws of data and projection
matrices $\bm{A}$. This variance decays as $\mathcal{O}(1/P+1/N)$.

Let $\bm{\Psi}\in\mathbb{R}^{P\times M}$ with
$\bm{\Psi}^{\mu}_{k}=\psi_{k}(\bm{x}^{\mu})$. Also define
$\bm{\Lambda}_{ij}=\lambda_{i}\delta_{ij}$.The discrepancy between the target
weights and the model’s effective weights is

$\displaystyle\bm{v}^{0}\equiv\bm{w}^{\star}-\frac{1}{\sqrt{N}}\bm{A}^{\top}\bm{w}(t).$
(6)

The test loss is then given by
$\mathcal{L}(t)=\frac{1}{M}\sum_{k}\lambda_{k}v^{0}_{k}(t)^{2}$. The
$\bm{v}^{0}$ vector has the following dynamics:

$\displaystyle\frac{d}{dt}\bm{v}^{0}(t)=-\left(\frac{1}{N}\bm{A}^{\top}\bm{A}\right)\left(\frac{1}{P}\bm{\Psi}^{\top}\bm{\Psi}\right)\bm{v}^{0}(t).$
(7)

Using DMFT, we characterize this limit by tracking $\bm{v}^{0}$ together with
the following random vectors:

$\displaystyle\bm{v}^{1}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{\Psi}\bm{v}^{0}(t),\quad\bm{v}^{2}(t)=\frac{1}{\sqrt{P}}\bm{\Psi}^{\top}\bm{v}^{1}(t),$
(8) $\displaystyle\bm{v}^{3}_{k}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{A}\bm{v}^{2}(t),\quad\bm{v}^{4}(t)=\frac{1}{\sqrt{N}}\bm{A}^{\top}\bm{v}^{3}(t).$

The key summary statistics (also called order parameters) are the correlation
functions:

$\displaystyle C_{0}(t,s)$
$\displaystyle=\frac{1}{M}\bm{v}^{0}(t)^{\top}\bm{\Lambda}\bm{v}^{0}(s),\
C_{1}(t,s)=\frac{1}{P}\bm{v}^{1}(t)\cdot\bm{v}^{1}(s),$ $\displaystyle
C_{2}(t,s)$ $\displaystyle=\frac{1}{M}\bm{v}^{2}(t)\cdot\bm{v}^{2}(t),\
C_{3}(t,s)=\frac{1}{N}\bm{v}^{3}(t)\cdot\bm{v}^{3}(s),$

as well as the response functions:

$\displaystyle
R_{1}(t,s)=\frac{1}{P}\mathrm{Tr}\left[\frac{\delta\bm{v}^{1}(t)}{\delta\bm{v}^{1}(s)}\right]\\!,R_{2,4}(t,s)=\frac{1}{M}\mathrm{Tr}\left[\frac{\delta\bm{v}^{2}(t)}{\delta\bm{v}^{4}(s)}\right],$
$\displaystyle
R_{3}(t,s)=\frac{1}{N}\mathrm{Tr}\left[\frac{\delta\bm{v}^{3}(t)}{\delta\bm{v}^{3}(s)}\right]\\!,R_{0,2}(t,s)=\frac{1}{M}\mathrm{Tr}\left[\bm{\Lambda}\frac{\delta\bm{v}^{0}(t)}{\delta\bm{v}^{2}(s)}\right].$

Here $\frac{\delta\bm{v}^{i}(t)}{\delta\bm{v}^{j}(s)}$ is the response of
$\bm{v}^{i}(t)$ to a kick in the dynamics of $\bm{v}^{j}$ at time $s$. See
appendix B.1.1 for details.

(a) $P=1000$ Test Loss Dynamics (b) Early Time Model Convergence (c) Late Time
Model Bottleneck (d) $N=1000$ Test Loss Dynamics (e) Early Time Data
Convergence (f) Late Time Data Bottleneck

Figure 2: Verification of the various bottleneck scalings for power-law
features with $a=1.5$ and $b=1.25$. Dashed black lines are DMFT solutions
while colors are simulations with standard deviation highlighted. (a) The loss
dynamics at large $\alpha$ will be bottlenecked by either time or finite
$\nu$. (b) Early in training, the loss converges to its limit as $N^{-1}$
(App. D). (c) At long times, the model’s asymptotic loss scales as
$N^{-(a-1)}$ (App. C). (d)-(f) The same results but for $N$ and $P$ switched.
The model exhibits $1/P$ corrections and early time and power law data
bottleneck scalings at late time.

The test loss $\mathcal{L}$ and train loss $\hat{\mathcal{L}}$ are related to
the time-time diagonal of $C_{0}(t,s)$ and $C_{1}(t,s)$ respectively

$\displaystyle\mathcal{L}(t)=C_{0}(t,t)+\sigma^{2}\ ,\
\mathcal{\hat{L}}(t)=C_{1}(t,t).$ (9)

These collective quantities concentrate over random draws of the disorder
(Sompolinsky & Zippelius, 1981). We show that these correlation and response
functions satisfy a closed set of integro-differential equations which depend
on $\alpha,\nu$ which we provide in the Appendices A.2.

Further, we show in Appendix A.3 that the response functions possess a time-
translation invariance property $R(t,s)=R(t-s)$. This enables exact analysis
in the Fourier domain
$R(\tau)=\int\frac{d\omega}{2\pi}e^{i\omega\tau}\mathcal{R}(\omega)$. These
response functions can then be used to solve for the correlation functions
$\\{C_{0},C_{1},C_{2},C_{3}\\}$.

To understand the convergence of the learned function $f$ along each
eigenfunction of the kernel, we introduce the transfer function222There are
dynamical analogues of the mode errors in (Bordelon et al., 2020; Canatar et
al., 2021) or learnabilities in (Simon et al., 2021). for mode $k$,
$H_{k}(t)\equiv\frac{\left<v_{k}^{0}(t)\right>}{w_{k}^{\star}}$. Our key
result is that the Fourier transform of $H_{k}$ can be simply expressed in
terms of the Fourier transforms of $R_{1},R_{3}$:

$\displaystyle\mathcal{H}_{k}(\omega)$
$\displaystyle=\frac{1}{i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}.$
(10)

These functions satisfy the self-consistent equations:

$\displaystyle\mathcal{R}_{1}(\omega)=1-\frac{1}{P}\sum_{k}\frac{\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}{i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)},$
(11)
$\displaystyle\mathcal{R}_{3}(\omega)=1-\frac{1}{N}\sum_{k}\frac{\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}{i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}.$

From these solved response functions $\mathcal{R}_{1},\mathcal{R}_{3}$, we can
compute local solutions to the correlation functions’ two-variable Fourier
transform $\mathcal{C}(\omega,\omega^{\prime})$ which are independent
equations for each pair of $\omega,\omega^{\prime}$. Information about the
early dynamics can be extracted from high frequencies $\omega\gg 1$ while
information about the late-time limit of the system can be extracted from
$\omega,\omega^{\prime}\to 0$ (App. C, D). For example, for the final test
loss,

$\displaystyle\lim_{t\to\infty}\mathcal{L}(t,\alpha,\nu)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\ \mathcal{C}_{0}(\omega,\omega^{\prime}).$ (12)

The full temporal trajectory can be obtained with an inverse Fourier transform
of $\mathcal{C}_{0}(\omega,\omega^{\prime})$. See Appendix A.4.

(a) $(a,b)=(2,1)\ ,\ L(C)\sim C^{-0.5}$ (b) $(a,b)=(2,1.5)\ ,\ L(C)\sim
C^{-0.4}$ (c) $(a,b)=(2,2)\ ,\ L(C)\sim C^{-0.33}$

Figure 3: Compute optimal scaling in our model is determined by tradeoff of
time and model-size bottlenecks. Solid colored lines are simulations with
power law features and in dashed black is the theoretical prediction of
compute optimal scaling. This Pareto frontier is defined as the minimum value
of $L$ at each compute $C$ over all possible choices of model size $N$.
Although the final losses do not depend on the spectral decay rate $b$ but
only on the task-power exponent $a$, the compute optimal scaling depends does
depend on $b$.

##  4 Results

Our results hold for any $\lambda_{k}$ and $w_{k}^{\star}$ and we provide some
simple analytically solvable examples in Appendix H. However, based on
empirical observations of NTK spectral decompositions on realistic datasets
(Bordelon & Pehlevan, 2021; Spigler et al., 2020; Bordelon & Pehlevan, 2022a;
Bahri et al., 2021; Maloney et al., 2022), here, we focus on the case of power
law features. In this setting, eigenvalues and target coefficients decay as a
power law in the index $k$

$\displaystyle(w^{\star}_{k})^{2}\lambda_{k}\sim k^{-a}\ ,\ \lambda_{k}\sim
k^{-b}.$ (13)

We will refer to $a$ as the task-power exponent and $b$ as the spectral decay
exponent333These power-law decay rates are also known as source and capacity
conditions in the kernel literature (Caponnetto & Vito, 2005; Cui et al.,
2021). See Figure 6 (a)-(b) for an example with a Residual CNN on CIFAR-5M.

##### Test loss power laws.

For power law features, the test loss will generally be bottlenecked by either
training time $t$ (steps of gradient descent), the size of the training set
$P$, or the size of the model $N$. We can derive bottleneck scalings from our
exact expressions for $\mathcal{L}(t,P,N)$ (Appendix I):

$\displaystyle\mathcal{L}(t,P,N)\approx\begin{cases}t^{-(a-1)/b}\ ,\
P,N\to\infty\ ,\ \text{(Time-Bottleneck)}\\\ P^{-(a-1)}\ ,\ t,N\to\infty\ ,\
\text{(Data-Bottleneck)}\\\ N^{-(a-1)}\ ,\ t,P\to\infty\ ,\ \text{(Model-
Bottleneck)}\end{cases}$ (14)

A consequence of this is an asymmetry in exponent between the model and data
bottlenecks compared to the time bottleneck. We verify this asymmetry in
Figure 2.

##### Bottlenecks as Rank-Constraints

All three of the bottleneck scalings arise due to rank constraints in the
effective dynamics. Heuristically, finite training time or the subsampling of
data/features leads to an approximate projection of the target function onto
the top $k_{\star}(t,P,N)$ eigenspace of the infinite-width kernel. The
components of the target function in the null-space of this projection are not
learned. This leads to an approximate test loss of the form

$\displaystyle\mathcal{L}\approx\sum_{k>k_{\star}}(w^{\star}_{k})^{2}\lambda_{k}\approx
k_{\star}^{-(a-1)}.$ (15)

For model and data bottlenecks we have that $k_{\star}\propto N$ and
$k_{\star}\propto P$ respectively (App. I). On the other hand, $k_{\star}$ for
the time bottleneck also depends on the structure of the features through the
exponent $b$. This is due to the fact that the $k$-th eigenfeature is learned
at a timescale $\tau_{k}\sim k^{b}$. Thus at time $t$, we have learned the
first $k_{\star}\approx t^{1/b}$ modes and the variance in the remaining modes
gives $\sim t^{-(a-1)/b}$. In the limit of $t\to\infty$ our data and model
bottleneck scalings agree with the resolution and variance-limited scalings
studied in (Bahri et al., 2021) as well as prior works on kernels and random
feature models (Bordelon et al., 2020; Maloney et al., 2022).

##### Asymmetric Compute Optimal Scaling Strategy

These results suggest that in our model, one should scale time and model size
differently at compute budget $C=Nt$,

$\displaystyle t\sim C^{\frac{b}{b+1}}\ ,\ \nu\sim C^{\frac{1}{b+1}}\ ,\
\implies\mathcal{L}^{\star}(C)\sim C^{-\frac{a-1}{1+b}}.$ (16)

We obtain the above scaling by approximating the loss as a sum of the three
terms in equation (14) and a constant, see Appendix M. This analysis suggests
that for features that have rapid decay in their eigenspectrum, it is
preferable to allocate greater resources toward training time rather than
model size as the compute budget increases. This is consistent with the
findings of (Hoffmann et al., 2022) for language models. In the limit as $b\to
1$, the time and parameter count should be scaled linearly together. We verify
this scaling rule and its $b$-dependence in Figure 3.

(a) Test Loss ($P=128$) (b) Train Loss ($P=128$) (c) Train-Test Gaps $N=512$

Figure 4: In a data limited regime, wider networks train faster but cannot
indefinitely improve generalization by making $N$ larger. (a) Test loss for
power-law features with $a=1.5$ and $b=1.25$ with $P=128$ and varying $N$. In
this regime, there are diminishing returns to making the model size larger.
(b) For $N<P$, the model is underparameterized and cannot achieve zero train
loss. For $N>P$, the train loss will eventually decay at exponential rate
which depends on $N$, despite the test loss saturating. (c) The train and test
losses gradually separate at a rate which depends on $P$.

##### Wider is Better Requires Sufficient Data

Larger models are not always better in terms of test loss for all time $t$, as
we showed in Figure 1 (c), especially if the dataset is limited. In Figure 4,
we illustrate that larger $N$ can improve convergence to a data-bottlenecked
loss for power law features. However, the loss may still be non-monotonic in
training time, motivating regularization or early stopping.

##### Gradual Buildup of Overfitting Effects

The exact gap between train and test losses can exactly be expressed in terms
of the DMFT order parameters:

$\displaystyle\mathcal{L}(t)-\mathcal{\hat{L}}(t)=-\frac{2}{P}\int_{0}^{t}dt^{\prime}\
R_{0,2}(t,t^{\prime})C_{1}(t,t^{\prime})$ (17)
$\displaystyle+\frac{1}{P^{2}}\int_{0}^{t}\int_{0}^{t}dt^{\prime}ds^{\prime}R_{0,2}(t,t^{\prime})R_{0,2}(t,s^{\prime})C_{1}(t^{\prime},s^{\prime}).$

We derive this relation in the Appendix E. At early time this gap goes as
$\mathcal{O}(1/P)$ (App. D, E). At late time, however, this picks up a
nontrivial task-dependent scaling with $P$ as we show in Figure 2 (e)-(f) and
App. C. In Figure 4 (c) we show this gradual accumulation of finite data on
the test-train loss gap. For larger datasets $P$ it takes longer training time
to begin overfitting (App. E).

##### Ensembling is Not Always Compute Optimal

(a) $E$ vs $N$ allocation of compute (b) Transfer Functions $H_{k}(\tau)$ for
$P=128$

Figure 5: Ensembling $E$ models of size $N$ improves performance by reducing
initialization variance by a factor of $E$ (see Appendix G.1) (a) However, at
fixed compute $C=NEt$, increasing the model size $N$ is preferable, since the
bias is also reduced. (b) The transfer functions ${H}_{k}(\tau)$ computed from
the DMFT determine the error as $E\to\infty$ depend on $N,P$ and saturate in
performance at long times, while the $N,P\to\infty$ curves decay
exponentially.

Ensembling a set of models means averaging their predictions over the same
datasets but with different intitialization seeds. This reduces test loss by
reducing the variance of the model output $f$ due to initialization. This
improvement can be predicted from an extension of our DMFT (App. G).
Analogously, bagging over $B$ datasets reduces variance due to sampling of
data.

(a) NTK Spectra for Varying Widths (b) Task-Power Decay (c) Lazy vs Feature
Compute Scaling

Figure 6: Our theory predicts time and compute scalings for linearized
networks on realistic datasets. (a) The initial NTK spectra and (b) task-power
distributions for ResNets of width $N$ on CIFAR-5M are well described by
powerlaws $\lambda_{k}\sim k^{-2.0}$ and $k^{-0.15}$ for large $k$. (c) The
predicted compute optimal scaling of for the ResNet obeys
$\mathcal{L}_{\star}(C)\sim C^{-0.05}$. However, for networks outside of the
kernel regime (dashed lines), feature learning can substantially alter the
observed scaling laws and improve the loss curves as a function of compute.

One might imagine that ensembling many finite sized models would allow one to
approach the performance of an infinite sized model ($N\to\infty$). If this
were possible, the compute optimal strategy could involve a tradeoff between
ensemble count and model size. However, recent experiments show that there is
a limited benefit from ensembling on large datasets when compared to
increasing model size (Vyas et al., 2023). We illustrate this in Figure 5 (a).
Our theory can explain these observations as it predicts the effect of
ensembling $E$ times on the learning dynamics as we show in App. G. The main
reason to prefer increasing $N$ rather than increasing $E$ is that larger $N$
has lower bias in the dynamics, whereas ensembling only reduces variance. The
bias of the model $\mathcal{B}$ has the form

$\displaystyle\mathcal{B}(t,N,P)=\sum_{k}\lambda_{k}(w_{k}^{\star})^{2}H_{k}(t,N,P)^{2},\vspace{-2pt}$
(18)

which depend on transfer function $H_{k}$ that we illustrate for power-law
features in Figure 5 (b). Since $H_{k}(\tau)$ depend on $N,P$, we see that
ensembling/bagging cannot recover the learning curve of the $N,P\to\infty$
system.

##  5 Tests on Realistic Networks

We now move beyond synthetic power-law datasets and consider realistic image
datasets and architectures. We take the CIFAR-5M dataset introduced in
(Nakkiran et al., 2021a) and consider the task of classfiying animate vs
inanimate objects. We plot the spectra of the finite-width NTK at
initialization across different widths for a Wide ResNet (Zagoruyko &
Komodakis, 2016) in Figure 6 a). Here the width parameter corresponds to the
number of channels in the hidden layers. Following (Canatar et al., 2021), we
define $C(k)$ as the fraction of the task captured by the top $k$ kernel
eigenmodes:

$C(k)\equiv\frac{\sum_{i\leq
k}\lambda_{i}(w_{i}^{*})^{2}}{\sum_{i}\lambda_{i}(w_{i}^{*})^{2}}.$ (19)

Then $1-C(k)$ is the portion of the task left unexplained. We plot this for
the initial NTKs across widths in Figure 6 b). We extract the spectral decay
exponent $b$ and the the task power exponent $a$ from these two curves.
Together, these give the learning scaling laws of the linearized neural
network model on this dataset. We plot the compute optimal scaling laws of
these linearized models in Figure 6 c). We also plot the predicted scaling law
$C^{-(a-1)/(1+b)}$ in blue and find excellent agreement.

###  5.1 The Role of Feature Learning

We also compare these scalings to those of the compute optimal learning curves
for feature-learning networks. We train several networks with different widths
and initialization seeds for 64 epochs through the dataset. We observe
substantially different compute-optimal scaling exponents in the dotted curves
of Figure 6 c). This means that although our random feature model does capture
the correct linearized scaling trends, which have all of the qualities
observed in realistic scaling laws, more is needed to capture the acceleration
of scaling induced by feature learning. Further analyses of the after-kernels
of feature learning networks are performed in Appendix K. We see that the
kernels continue to evolve substantially throughout training. This indicates
that a full explanation of the compute optimal scaling exponents will require
something resembling a mechanistic theory of kernel evolution (Long, 2021;
Fort et al., 2020; Atanasov et al., 2022; Bordelon & Pehlevan, 2022b).

##  6 Conclusion

We have presented a model that recovers a wide variety of phenomena observed
in more realistic deep learning settings. Our theory includes not just model
size and dataset size as parameters but also explicitly treats the temporal
dynamics of training. We observe different scaling exponents for performance
in terms of model size and number of time steps. Future work to incorporate
kernel evolution into this model could further shed insight into the improved
scaling laws in the feature-learning regime. Overall, our results provide a
theoretical interpretation of compute-optimal scaling as a competition between
the training dynamics of the infinite width/infinite data limit and finite
model-size bottleneck.

##  7 Broader Impact Statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

##  8 Acknowledgements

We are grateful to Yasaman Bahri, Stefano Mannelli, Francesca Mignacco, Jascha
Sohl-Dickstein, and Nikhil Vyas for useful conversations. We thank Clarissa
Lauditi and Jacob Zavatone-Veth for comments on the manuscript.

B.B. is supported by a Google PhD Fellowship. A.A. is supported by a Fannie
and John Hertz Fellowship. C.P. is supported by NSF grant DMS-2134157, NSF
CAREER Award IIS-2239780, and a Sloan Research Fellowship. This work has been
made possible in part by a gift from the Chan Zuckerberg Initiative Foundation
to establish the Kempner Institute for the Study of Natural and Artificial
Intelligence.

## References

  * Adlam & Pennington (2020a) Adlam, B. and Pennington, J.  The neural tangent kernel in high dimensions: Triple descent and a multi-scale theory of generalization.  In _International Conference on Machine Learning_ , pp. 74–84. PMLR, 2020a.
  * Adlam & Pennington (2020b) Adlam, B. and Pennington, J.  Understanding double descent requires a fine-grained bias-variance decomposition.  _Advances in neural information processing systems_ , 33:11022–11032, 2020b.
  * Advani et al. (2020) Advani, M. S., Saxe, A. M., and Sompolinsky, H.  High-dimensional dynamics of generalization error in neural networks.  _Neural Networks_ , 132:428–446, 2020.
  * Alabdulmohsin et al. (2023) Alabdulmohsin, I., Zhai, X., Kolesnikov, A., and Beyer, L.  Getting vit in shape: Scaling laws for compute-optimal model design.  _arXiv preprint arXiv:2305.13035_ , 2023.
  * Arora & Goyal (2023) Arora, S. and Goyal, A.  A theory for emergence of complex skills in language models.  _arXiv preprint arXiv:2307.15936_ , 2023.
  * Atanasov et al. (2022) Atanasov, A., Bordelon, B., and Pehlevan, C.  Neural networks as kernel learners: The silent alignment effect.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=1NvflqAdoom.
  * Atanasov et al. (2023) Atanasov, A., Bordelon, B., Sainathan, S., and Pehlevan, C.  The onset of variance-limited behavior for networks in the lazy and rich regimes.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=JLINxPOVTh7.
  * Bahri et al. (2021) Bahri, Y., Dyer, E., Kaplan, J., Lee, J., and Sharma, U.  Explaining neural scaling laws.  _arXiv preprint arXiv:2102.06701_ , 2021.
  * Bordelon & Pehlevan (2021) Bordelon, B. and Pehlevan, C.  Learning curves for sgd on structured features.  _arXiv preprint arXiv:2106.02713_ , 2021.
  * Bordelon & Pehlevan (2022a) Bordelon, B. and Pehlevan, C.  Learning curves for SGD on structured features.  In _International Conference on Learning Representations_ , 2022a.  URL https://openreview.net/forum?id=WPI2vbkAl3Q.
  * Bordelon & Pehlevan (2022b) Bordelon, B. and Pehlevan, C.  Self-consistent dynamical field theory of kernel evolution in wide neural networks.  _arXiv preprint arXiv:2205.09653_ , 2022b.
  * Bordelon & Pehlevan (2023) Bordelon, B. and Pehlevan, C.  Dynamics of finite width kernel and prediction fluctuations in mean field neural networks.  _arXiv preprint arXiv:2304.03408_ , 2023.
  * Bordelon et al. (2020) Bordelon, B., Canatar, A., and Pehlevan, C.  Spectrum dependent learning curves in kernel regression and wide neural networks.  In _International Conference on Machine Learning_ , pp. 1024–1034. PMLR, 2020.
  * Bordelon et al. (2023) Bordelon, B., Noci, L., Li, M. B., Hanin, B., and Pehlevan, C.  Depthwise hyperparameter transfer in residual networks: Dynamics and scaling limit, 2023.
  * Caballero et al. (2022) Caballero, E., Gupta, K., Rish, I., and Krueger, D.  Broken neural scaling laws.  _arXiv preprint arXiv:2210.14891_ , 2022.
  * Canatar et al. (2021) Canatar, A., Bordelon, B., and Pehlevan, C.  Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks.  _Nature communications_ , 12(1):2914, 2021.
  * Caponnetto & Vito (2005) Caponnetto, A. and Vito, E. D.  Fast rates for regularized least-squares algorithm.  2005\.
  * Cheng & Montanari (2022) Cheng, C. and Montanari, A.  Dimension free ridge regression.  _arXiv preprint arXiv:2210.08571_ , 2022.
  * Chizat et al. (2019) Chizat, L., Oyallon, E., and Bach, F.  On lazy training in differentiable programming.  _Advances in neural information processing systems_ , 32, 2019.
  * Cortes et al. (2012) Cortes, C., Mohri, M., and Rostamizadeh, A.  Algorithms for learning kernels based on centered alignment.  _The Journal of Machine Learning Research_ , 13(1):795–828, 2012.
  * Crisanti & Sompolinsky (2018) Crisanti, A. and Sompolinsky, H.  Path integral approach to random neural networks.  _Physical Review E_ , 98(6):062120, 2018.
  * Cui et al. (2021) Cui, H., Loureiro, B., Krzakala, F., and Zdeborová, L.  Generalization error rates in kernel regression: The crossover from the noiseless to noisy regime.  _Advances in Neural Information Processing Systems_ , 34:10131–10143, 2021.
  * Dyer & Gur-Ari (2020) Dyer, E. and Gur-Ari, G.  Asymptotics of wide networks from feynman diagrams.  In _International Conference on Learning Representations_ , 2020.
  * Fort et al. (2020) Fort, S., Dziugaite, G. K., Paul, M., Kharaghani, S., Roy, D. M., and Ganguli, S.  Deep learning versus kernel learning: an empirical study of loss landscape geometry and the time evolution of the neural tangent kernel.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Geiger et al. (2020) Geiger, M., Jacot, A., Spigler, S., Gabriel, F., Sagun, L., d’Ascoli, S., Biroli, G., Hongler, C., and Wyart, M.  Scaling description of generalization with number of parameters in deep learning.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(2):023401, 2020.
  * Ghosh et al. (2022) Ghosh, N., Mei, S., and Yu, B.  The three stages of learning dynamics in high-dimensional kernel methods.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=EQmAP4F859.
  * Hastie et al. (2022) Hastie, T., Montanari, A., Rosset, S., and Tibshirani, R. J.  Surprises in high-dimensional ridgeless least squares interpolation.  _Annals of statistics_ , 50(2):949, 2022.
  * Helias & Dahmen (2020) Helias, M. and Dahmen, D.  _Statistical field theory for neural networks_ , volume 970.  Springer, 2020.
  * Hestness et al. (2017) Hestness, J., Narang, S., Ardalani, N., Diamos, G., Jun, H., Kianinejad, H., Patwary, M. M. A., Yang, Y., and Zhou, Y.  Deep learning scaling is predictable, empirically.  _arXiv preprint arXiv:1712.00409_ , 2017.
  * Hoffmann et al. (2022) Hoffmann, J., Borgeaud, S., Mensch, A., Buchatskaya, E., Cai, T., Rutherford, E., Casas, D. d. L., Hendricks, L. A., Welbl, J., Clark, A., et al.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_ , 2022.
  * Jacot et al. (2018) Jacot, A., Gabriel, F., and Hongler, C.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Kaplan et al. (2020) Kaplan, J., McCandlish, S., Henighan, T., Brown, T. B., Chess, B., Child, R., Gray, S., Radford, A., Wu, J., and Amodei, D.  Scaling laws for neural language models.  _arXiv preprint arXiv:2001.08361_ , 2020.
  * Long (2021) Long, P. M.  Properties of the after kernel.  _arXiv preprint arXiv:2105.10585_ , 2021.
  * Loureiro et al. (2021) Loureiro, B., Gerbelot, C., Cui, H., Goldt, S., Krzakala, F., Mezard, M., and Zdeborová, L.  Learning curves of generic features maps for realistic datasets with a teacher-student model.  _Advances in Neural Information Processing Systems_ , 34:18137–18151, 2021.
  * Loureiro et al. (2022) Loureiro, B., Gerbelot, C., Refinetti, M., Sicuro, G., and Krzakala, F.  Fluctuations, bias, variance & ensemble of learners: Exact asymptotics for convex losses in high-dimension.  In _International Conference on Machine Learning_ , pp. 14283–14314. PMLR, 2022.
  * Maloney et al. (2022) Maloney, A., Roberts, D. A., and Sully, J.  A solvable model of neural scaling laws.  _arXiv preprint arXiv:2210.16859_ , 2022.
  * Mannelli et al. (2019) Mannelli, S. S., Krzakala, F., Urbani, P., and Zdeborova, L.  Passed & spurious: Descent algorithms and local minima in spiked matrix-tensor models.  In _international conference on machine learning_ , pp. 4333–4342. PMLR, 2019.
  * Martin et al. (1973) Martin, P. C., Siggia, E., and Rose, H.  Statistical dynamics of classical systems.  _Physical Review A_ , 8(1):423, 1973.
  * Mei & Montanari (2022) Mei, S. and Montanari, A.  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.
  * Michaud et al. (2023) Michaud, E. J., Liu, Z., Girit, U., and Tegmark, M.  The quantization model of neural scaling.  _arXiv preprint arXiv:2303.13506_ , 2023.
  * Mignacco & Urbani (2022) Mignacco, F. and Urbani, P.  The effective noise of stochastic gradient descent.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2022(8):083405, 2022.
  * Mignacco et al. (2020) Mignacco, F., Krzakala, F., Urbani, P., and Zdeborová, L.  Dynamical mean-field theory for stochastic gradient descent in gaussian mixture classification.  _Advances in Neural Information Processing Systems_ , 33:9540–9550, 2020.
  * Nakkiran et al. (2021a) Nakkiran, P., Neyshabur, B., and Sedghi, H.  The deep bootstrap framework: Good online learners are good offline generalizers.  In _International Conference on Learning Representations_ , 2021a.
  * Nakkiran et al. (2021b) Nakkiran, P., Neyshabur, B., and Sedghi, H.  The deep bootstrap framework: Good online learners are good offline generalizers.  In _International Conference on Learning Representations_ , 2021b.  URL https://openreview.net/forum?id=guetrIHLFGI.
  * Paquette et al. (2021) Paquette, C., Lee, K., Pedregosa, F., and Paquette, E.  Sgd in the large: Average-case analysis, asymptotics, and stepsize criticality.  In _Conference on Learning Theory_ , pp. 3548–3626. PMLR, 2021.
  * Roberts et al. (2022) Roberts, D. A., Yaida, S., and Hanin, B.  _The principles of deep learning theory_.  Cambridge University Press Cambridge, MA, USA, 2022.
  * Ruben & Pehlevan (2023) Ruben, B. S. and Pehlevan, C.  Learning curves for noisy heterogeneous feature-subsampled ridge ensembles.  _ArXiv_ , 2023.
  * Simon et al. (2021) Simon, J. B., Dickens, M., Karkada, D., and DeWeese, M. R.  The eigenlearning framework: A conservation law perspective on kernel regression and wide neural networks.  _arXiv preprint arXiv:2110.03922_ , 2021.
  * Simon et al. (2023) Simon, J. B., Karkada, D., Ghosh, N., and Belkin, M.  More is better in modern machine learning: when infinite overparameterization is optimal and overfitting is obligatory.  _arXiv preprint arXiv:2311.14646_ , 2023.
  * Sompolinsky & Zippelius (1981) Sompolinsky, H. and Zippelius, A.  Dynamic theory of the spin-glass phase.  _Physical Review Letters_ , 47(5):359, 1981.
  * Sorscher et al. (2022) Sorscher, B., Geirhos, R., Shekhar, S., Ganguli, S., and Morcos, A.  Beyond neural scaling laws: beating power law scaling via data pruning.  _Advances in Neural Information Processing Systems_ , 35:19523–19536, 2022.
  * Spigler et al. (2020) Spigler, S., Geiger, M., and Wyart, M.  Asymptotic learning curves of kernel methods: empirical data versus teacher–student paradigm.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(12):124001, 2020.
  * Varre et al. (2021) Varre, A. V., Pillaud-Vivien, L., and Flammarion, N.  Last iterate convergence of sgd for least-squares in the interpolation regime.  _Advances in Neural Information Processing Systems_ , 34:21581–21591, 2021.
  * Vyas et al. (2022) Vyas, N., Bansal, Y., and Nakkiran, P.  Limitations of the ntk for understanding generalization in deep learning.  _arXiv preprint arXiv:2206.10012_ , 2022.
  * Vyas et al. (2023) Vyas, N., Atanasov, A., Bordelon, B., Morwani, D., Sainathan, S., and Pehlevan, C.  Feature-learning networks are consistent across widths at realistic scales.  _arXiv preprint arXiv:2305.18411_ , 2023.
  * Yang et al. (2021) Yang, G., Hu, E. J., Babuschkin, I., Sidor, S., Liu, X., Farhi, D., Ryder, N., Pachocki, J., Chen, W., and Gao, J.  Tuning large neural networks via zero-shot hyperparameter transfer.  In Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , 2021.  URL https://openreview.net/forum?id=Bx6qKuBM2AD.
  * Zagoruyko & Komodakis (2016) Zagoruyko, S. and Komodakis, N.  Wide residual networks.  _arXiv preprint arXiv:1605.07146_ , 2016.
  * Zavatone-Veth & Pehlevan (2023) Zavatone-Veth, J. A. and Pehlevan, C.  Learning curves for deep structured gaussian feature models, 2023.
  * Zavatone-Veth et al. (2022) Zavatone-Veth, J. A., Tong, W. L., and Pehlevan, C.  Contrasting random and learned features in deep bayesian linear regression.  _Physical Review E_ , 105(6):064118, 2022.

##  Appendix A Derivation of Dynamical Model of Scaling Laws

We investigate the simplest possible model which can exhibit task-dependent
time, model size and finite data bottlenecks. We therefore choose to study a
linear model with projected features

$\displaystyle
f(\bm{x})=\frac{1}{\sqrt{M}}\bm{w}^{\top}\left(\frac{1}{\sqrt{N}}\bm{A}\bm{\psi}(\bm{x})\right)\
,\ y(\bm{x})=\frac{1}{\sqrt{M}}\bm{w}_{\star}\cdot\bm{\psi}(\bm{x}).$ (20)

The weights $\bm{w}$ are updated with gradient descent on a random training
dataset which has (possibly) noise corrupted target values
$y_{\mu}=y(\bm{x}_{\mu})+\sigma\epsilon_{\mu}$. This leads to the following
gradient flow dynamics

$\displaystyle\frac{\partial}{\partial
t}\bm{w}(t)=\frac{\sqrt{M}}{P\sqrt{N}}\sum_{\mu=1}^{P}(y_{\mu}-f_{\mu})\bm{A}\bm{\psi}_{\mu}=\frac{1}{\sqrt{N}}\bm{A}\left(\frac{1}{P}\sum_{\mu=1}^{P}\bm{\psi}_{\mu}\left[\bm{\psi}_{\mu}^{\top}\left(\bm{w}_{\star}-\frac{1}{\sqrt{N}}\bm{A}^{\top}\bm{w}\right)+\sqrt{M}\sigma\epsilon_{\mu}\right]\right).$
(21)

We introduce the variable
$\bm{v}^{0}=\bm{w}_{\star}-\frac{1}{\sqrt{N}}\bm{A}^{\top}\bm{w}$ to represent
the residual error of the learned weight vector. This residual error has the
following dynamics.

$\displaystyle\partial_{t}\bm{v}^{0}(t)=-\left(\frac{1}{N}\bm{A}^{\top}\bm{A}\right)\left[\left(\frac{1}{P}\bm{\Psi}^{\top}\bm{\Psi}\right)\bm{v}^{0}(t)+\frac{\sigma}{\alpha\sqrt{M}}\bm{\Psi}^{\top}\bm{\epsilon}\right].$
(22)

The entries of each matrix are treated as random with
$\Psi^{\mu}_{k}\sim\mathcal{N}(0,\lambda_{k})$ and
$A_{jk}\sim\mathcal{N}(0,1)$. To study the dynamical evolution of the test
error
$\mathcal{L}(t)=\frac{1}{M}\bm{v}^{0}(t)^{\top}\bm{\Lambda}\bm{v}^{0}(t)+\sigma^{2}$,
we introduce the sequence of vectors

$\displaystyle\bm{v}^{1}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{\Psi}\bm{v}^{0}(t)+\sigma\bm{\epsilon}\
,\ \bm{v}^{2}(t)=\frac{1}{\alpha\sqrt{M}}\bm{\Psi}^{\top}\bm{v}^{1}(t)$
$\displaystyle\bm{v}^{3}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{A}\bm{v}^{2}(t)\ ,\
\bm{v}^{4}(t)=\frac{1}{\nu\sqrt{M}}\bm{A}^{\top}\bm{v}^{3}(t).$ (23)

The train and test losses can be computed from the $\bm{v}^{0}$ and
$\bm{v}^{1}$ fields

$\displaystyle\hat{\mathcal{L}}(t)=\frac{1}{P}\sum_{\mu=1}^{P}v^{1}_{\mu}(t)^{2}\
,\
\mathcal{L}(t)=\frac{1}{M}\sum_{k=1}^{M}\lambda_{k}v^{0}_{k}(t)^{2}+\sigma^{2}.$
(24)

In the next section, we derive a statistical description of the dynamics in an
appropriate asymptotic limit using dynamical mean field theory methods.

###  A.1 DMFT Equations for the Asymptotic Limit

Standard field theoretic arguments such as the cavity or path integral methods
can be used to compute the effective statistical description of the dynamics
in the limit of large $M,N,P$ with fixed ratios $\alpha=P/M$ and
$\nu=\frac{N}{M}$ (see Appendix B). This computation gives us the following
statistical description of the dynamics.

$\displaystyle v^{1}(t)=u^{1}(t)+\frac{1}{\alpha}\int
dsR_{0,2}(t,s)v_{1}(s)+\sigma\epsilon\ ,\ u^{1}(t)\sim\mathcal{GP}(0,C_{0})\
,\ \epsilon\sim\mathcal{N}(0,1),$ (25) $\displaystyle
v^{2}_{k}(t)=u^{2}_{k}(t)+\lambda_{k}\int dsR_{1}(t,s)v^{0}_{k}(s)\ ,\
u^{2}_{k}(t)\sim\mathcal{GP}\left(0,\frac{1}{\alpha}\lambda_{k}C_{1}\right),$
$\displaystyle v^{3}(t)=u^{3}(t)+\frac{1}{\nu}\int dsR_{2,4}(t,s)v^{3}(s)\ ,\
u^{3}(t)\sim\mathcal{GP}\left(0,C_{2}\right),$ $\displaystyle
v^{4}_{k}(t)=u^{4}_{k}(t)+\int dsR_{3}(t,s)v^{2}_{k}(s)\ ,\
u^{4}_{k}(t)\sim\mathcal{GP}\left(0,\frac{1}{\nu}C_{3}\right),$
$\displaystyle\partial_{t}v^{0}_{k}(t)=-v^{4}_{k}(t).$

The correlation and response functions obey

$\displaystyle C_{0}(t,s)$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}\left<v^{0}_{k}(t)v^{0}_{k}(s)\right>\
,\ C_{1}(t,s)=\left<v^{1}(t)v^{1}(s)\right>\ ,\
C_{2}(t,s)=\frac{1}{M}\sum_{k=1}^{M}\left<v_{k}^{2}(t)v_{k}^{2}(s)\right>$
$\displaystyle R_{0,2}(t,s)$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}\left<\frac{\delta
v^{0}_{k}(t)}{\delta u^{2}_{k}(s)}\right>\ ,\
R_{2,4}(t,s)=\frac{1}{M}\sum_{k}\left<\frac{\delta v^{2}_{k}(t)}{\delta
u^{4}_{k}(s)}\right>$ $\displaystyle R_{1}(t,s)$
$\displaystyle=\left<\frac{\delta v^{1}(t)}{\delta u^{1}(s)}\right>\ ,\
R_{3}(t,s)=\left<\frac{\delta v^{3}(t)}{\delta u^{3}(s)}\right>$

These equations are exact in the joint proportional limit for any value of
$\alpha,\nu$.

###  A.2 Closing the Equations for the Order Parameters

Though we expressed the dynamics in terms of random fields, we stress in this
section that all of the dynamics for the correlation and response functions
close in terms of integro-dfferential equations. To shorten the expression, we
will provide the expression for $\beta=0$, but momementum can easily be added
back by making the substitution
$\partial_{t}\to\partial_{t}+\beta\partial^{2}_{t}$.

First, our closed integral equations for the response functions are

$\displaystyle R_{0,2,k}(t,s)=-\int
dt^{\prime}\Theta(t-t^{\prime})R_{3}(t^{\prime},s)-\lambda_{k}\int
dt^{\prime}dt^{\prime\prime}dt^{\prime\prime\prime}\Theta(t-t^{\prime})R_{3}(t^{\prime},t^{\prime\prime})R_{1}(t^{\prime\prime},t^{\prime\prime\prime})R_{0,2,k}(t^{\prime\prime\prime},s)$
$\displaystyle R_{1}(t,s)=\delta(t-s)+\frac{1}{\alpha}\int
dt^{\prime}R_{0,2}(t,t^{\prime})R_{1}(t^{\prime},s)$ $\displaystyle
R_{2,4,k}(t,s)=-\lambda_{k}\int
dt^{\prime}dt^{\prime\prime}R_{1}(t,t^{\prime})\Theta(t^{\prime}-t^{\prime\prime})-\lambda_{k}\int
dt^{\prime}dt^{\prime\prime}dt^{\prime\prime\prime}R_{1}(t,t^{\prime})\Theta(t^{\prime}-t^{\prime\prime})R_{3}(t^{\prime\prime},t^{\prime\prime\prime})R_{2,4,k}(t^{\prime\prime\prime},s)$
$\displaystyle R_{3}(t,s)=\delta(t-s)+\frac{1}{\nu}\int
dt^{\prime}R_{2,4}(t,t^{\prime})R_{3}(t^{\prime},s)$ $\displaystyle
R_{0,2}(t,s)=\frac{1}{M}\sum_{k}\lambda_{k}R_{0,2,k}(t,s)\ ,\
R_{2,4}(t,s)=\frac{1}{M}\sum_{k}R_{2,4,k}(t,s)$ (26)

We note that these equations imply causality in all of the response functions
since $R(t,s)=0$ for $t<s$. Once these equations are solved for the response
functions, we can determine the correlation functions, which satisfy

$\displaystyle\partial_{ts}^{2}C_{0,k}(t,s)=$ $\displaystyle-\lambda_{k}\int
dt^{\prime}dt^{\prime\prime}R_{3}(t,t^{\prime})R_{1}(t^{\prime},t^{\prime\prime})\partial_{s}C_{0,k}(t^{\prime\prime},s)$
$\displaystyle-\lambda_{k}\int
ds^{\prime}R_{3}(s,s^{\prime})R_{1}(s^{\prime},s^{\prime\prime})\partial_{t}C_{0,k}(t,s^{\prime\prime})$
$\displaystyle+\lambda_{k}^{2}\int
dt^{\prime}dt^{\prime\prime}ds^{\prime}ds^{\prime\prime}R_{3}(t,t^{\prime})R_{1}(t^{\prime},t^{\prime\prime})R_{3}(s,s^{\prime})R_{1}(s^{\prime},s^{\prime\prime})C_{0,k}(t^{\prime\prime},s^{\prime\prime})$
$\displaystyle-(w_{k}^{\star})^{2}\delta(t)\delta(s)-\frac{1}{\nu}C_{3}(t,s)-\frac{1}{\alpha}\int
dt^{\prime}ds^{\prime}R_{3}(t,t^{\prime})R_{3}(s,s^{\prime})C_{1}(t^{\prime},s^{\prime})$
$\displaystyle C_{1}(t,s)=$ $\displaystyle\int
dt^{\prime}R_{1}(t,t^{\prime})R_{1}(s,s^{\prime})C_{0}(t^{\prime},s^{\prime})$
$\displaystyle C_{2,k}(t,s)=$ $\displaystyle-\lambda_{k}\int
dt^{\prime}dt^{\prime\prime}R_{1}(t,t^{\prime})R_{3}(t^{\prime},t^{\prime\prime})C_{2,k}(t^{\prime\prime},s)-\lambda_{k}\int
ds^{\prime}ds^{\prime\prime}R_{1}(s,s^{\prime})R_{3}(s^{\prime},s^{\prime\prime})C_{2}(t,s^{\prime\prime})$
$\displaystyle+\lambda_{k}^{2}\int
dt^{\prime}dt^{\prime\prime}ds^{\prime}ds^{\prime\prime}R_{1}(t,t^{\prime})R_{3}(t^{\prime},t^{\prime\prime})R_{1}(s,s^{\prime})R_{3}(s^{\prime},s^{\prime\prime})C_{2,k}(t^{\prime\prime},s^{\prime\prime})$
$\displaystyle C_{3}(t,s)=$ $\displaystyle\int
dt^{\prime}ds^{\prime}R_{3}(t,t^{\prime})R_{3}(s,s^{\prime})C_{2}(t^{\prime},s^{\prime})$
(27)

Solving these closed equations provide the complete statistical
characterization of the limit. The test and train losses are given by the
time-time diagonal of $C_{0}(t,t),C_{1}(t,t)$.

###  A.3 Time-translation Invariant (TTI) Solution to Response Functions

From the structure of the above equations, the response functions are time-
translation invariant (TTI) since they are only functionals of TTI
$\delta(t-s)$ Dirac-Delta function and $\Theta(t-s)$ Heaviside step-function.
As a consequence, we write each of our response functions in terms of their
Fourier transforms

$\displaystyle R(t,s)=R(t-s)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\
e^{i\omega(t-s)}\mathcal{R}(\omega).$ (28)

Using the fact that

$\displaystyle\delta(\tau)=\int\frac{d\omega}{2\pi}e^{i\omega\tau}\ ,\
\Theta(\tau)=\lim_{\epsilon\to
0^{+}}e^{-\epsilon\tau}\Theta(\tau)=\lim_{\epsilon\to
0^{+}}\int\frac{d\omega}{2\pi}\frac{e^{i\omega\tau}}{\epsilon+i\omega}$ (29)

We will keep track of the regulator $\epsilon$ and consider $\epsilon\to
0^{+}$ at the end of the computation. The resulting DMFT equations for the
response functions have the following form in Fourier space

$\displaystyle\mathcal{R}_{1}(\omega)=1+\frac{1}{\alpha}\mathcal{R}_{2,4}(\omega)\mathcal{R}_{1}(\omega)$
$\displaystyle\mathcal{R}_{3}(\omega)=1+\frac{1}{\nu}\mathcal{R}_{2,4}(\omega)\mathcal{R}_{3}(\omega)$
$\displaystyle\mathcal{R}_{0,2}(\omega)=-\frac{1}{M}\sum_{k}\frac{\lambda_{k}}{\epsilon+i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}\mathcal{R}_{3}(\omega)$
$\displaystyle\mathcal{R}_{2,4}(\omega)=-\frac{1}{M}\sum_{k}\frac{\lambda_{k}}{\epsilon+i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}\mathcal{R}_{1}(\omega)$
(30)

where $\epsilon\to 0$ will be taken after. Combining these equations, we
arrive at the simple set of coupled equations

$\displaystyle\mathcal{R}_{1}(\omega)$
$\displaystyle=1-\frac{1}{P}\sum_{k}\frac{\lambda_{k}\mathcal{R}_{3}(\omega)\mathcal{R}_{1}(\omega)}{\epsilon+i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}$
$\displaystyle\mathcal{R}_{3}(\omega)$
$\displaystyle=1-\frac{1}{N}\sum_{k}\frac{\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}{\epsilon+i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}$
(31)

After solving these equations for all $\omega$, we can invert the dynamics of
$v^{0}_{k}(t)$ to obtain its Fourier transform

$\displaystyle v^{0}_{k}(\omega)$
$\displaystyle=\mathcal{H}_{k}(\omega)\left[w_{k}^{\star}-{u}^{4}_{k}(\omega)-\lambda_{k}\mathcal{R}_{3}(\omega){u}^{2}_{k}(\omega)\right]$
$\displaystyle\mathcal{H}_{k}(\omega)$
$\displaystyle\equiv\frac{1}{\epsilon+i\omega+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}$
(32)

where we defined the transfer functions $\mathcal{H}_{k}(\omega)$. From this
equation, we can compute $v^{0}_{k}(t)$ through inverse Fourier-transformation
and then compute the correlation function to calculate the test error. An
interesting observation is that the response functions
$\mathcal{R}_{1}(\omega),\mathcal{R}_{3}(\omega)$ alter the pole structure in
the transfer function, generating $\nu,\alpha$ dependent timescales of
convergence.

###  A.4 Fourier Representations for Correlation Functions

While the response functions are TTI, the correlation functions transparently
are not (if the time-time diagonal $C_{0}(t,t)$ did not evolve, then the loss
$\mathcal{L}(t)$ wouldn’t change!). We therefore define the need to define the
double Fourier transform $\mathcal{C}(\omega,\omega^{\prime})$ for each
correlation function $C(t,s)$

$\displaystyle\mathcal{C}(\omega,\omega^{\prime})=\int dtds\ e^{-i\omega
t-i\omega^{\prime}s}C(t,s)\ ,\
C(t,s)=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}e^{i\omega
t+i\omega^{\prime}s}\mathcal{C}(\omega,\omega^{\prime})$ (33)

Assuming that all response functions and transfer functions $\mathcal{H}_{k}$
have been solved for, the correlation functions satisfy the closed set of
linear equations.

$\displaystyle\mathcal{C}_{0}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}\mathcal{C}_{3}(\omega,\omega^{\prime})+\frac{1}{\alpha}\lambda_{k}\mathcal{R}_{3}(\omega)\mathcal{R}_{3}(\omega^{\prime})\mathcal{C}_{1}(\omega,\omega^{\prime})\right]$
$\displaystyle\mathcal{C}_{1}(\omega,\omega^{\prime})$
$\displaystyle=\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})\mathcal{C}_{0}(\omega,\omega^{\prime})$
$\displaystyle\mathcal{C}_{2}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[\frac{1}{\alpha}(i\omega)(i\omega^{\prime})\mathcal{C}_{1}(\omega,\omega^{\prime})+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})\left((w^{\star}_{k})^{2}+\frac{1}{\nu}\mathcal{C}_{3}(\omega,\omega^{\prime})\right)\right]$
$\displaystyle\mathcal{C}_{3}(\omega,\omega^{\prime})$
$\displaystyle=\mathcal{R}_{3}(\omega)\mathcal{R}_{3}(\omega^{\prime})\mathcal{C}_{2}(\omega,\omega^{\prime})$
(34)

These equations can be efficiently solved for all pairs of
$\omega,\omega^{\prime}$ after the response functions have been identified.
Then one can take an inverse Fourier transform in both indices.

##  Appendix B Field Theoretic Derivation of DMFT Equations

In this section, we derive the field theoretic description of our model. We
will derive this using both the Martin-Siggia-Rose (MSR) path integral method
(Martin et al., 1973) and the dynamical cavity method. For a recent review of
these topics in the context of neural networks, see (Helias & Dahmen, 2020).

###  B.1 Path Integral Derivation

With the MSR formalism, we evaluate the moment generating functional for the
field $\bm{v}^{0}(t)$:

$\displaystyle Z[\\{\bm{j}(t)\\}]$
$\displaystyle=\left<\int\mathcal{D}\bm{v}^{0}(t)\,\delta\left(\dot{\bm{v}}^{0}(t)+\frac{1}{NP}\bm{A}^{\top}\bm{A}\bm{\Psi}^{\top}\bm{\Psi}\bm{v}^{0}(t)\right)\exp\left(\int
dt\ \bm{j}(t)\cdot\bm{v}^{0}(t)\right)\right>_{\bm{A},\bm{\Psi}}.$ (35)

Note that at zero source, we have the important identity that

$Z[0]=1.$ (36)

We insert a Dirac delta functions to enforce the definitions of each of the
fields $\\{\bm{v}^{1},\bm{v}^{2},\bm{v}^{3},\bm{v}^{4}\\}$ as in equation 8.

$\displaystyle Z[\\{\bm{j}(t)\\}]=\int$
$\displaystyle\mathcal{D}[\bm{v}^{0},\dots,\bm{v}^{4},\hat{\bm{v}}^{1}\dots\hat{\bm{v}}^{4}]\,\delta(\dot{\bm{v}}^{0}+\bm{v}^{4})\,\exp\left(\int
dt\ \bm{j}(t)\cdot\bm{v}^{0}(t)\right)$ (37)
$\displaystyle\times\left<\exp\left[i\int
dt\left[\hat{\bm{v}}_{1}(t)\cdot\left(\bm{v}^{1}(t)-\frac{1}{\sqrt{M}}\bm{\Psi}\bm{v}^{0}(t)\right)+\hat{\bm{v}}_{2}(t)\cdot\left(\bm{v}^{2}(t)-\frac{1}{\sqrt{P}}\bm{\Psi}^{\top}\bm{v}^{1}(t)\right)\right]\right]\right>_{\bm{\Psi}}$
$\displaystyle\times\left<\exp\left[i\int
dt\left[\hat{\bm{v}}_{3}(t)\cdot\left(\bm{v}^{3}(t)-\frac{1}{\sqrt{M}}\bm{\bm{A}}\bm{v}^{2}(t)\right)+\hat{\bm{v}}_{4}(t)\cdot\left(\bm{v}^{4}(t)-\frac{1}{\sqrt{N}}\bm{\bm{A}}^{\top}\bm{v}^{3}(t)\right)\right]\right]\right>_{\bm{A}}.$

At this stage we can add sources $\tilde{\bm{j}}$ for each $\hat{\bm{v}_{i}}$
variable, yielding a $Z[\bm{j}(t),\tilde{\bm{j}}(t)]$. Interpreting each
source as modification of the respective evolution equation, we see that this
modified moment-generating function remains equal to unity at any value of
$\tilde{\bm{j}}$, $Z[0,\tilde{\bm{j}}(t)]=1$. As a consequence, all
correlation functions consisting only of $\hat{\bm{v}}^{i}$ variables vanish.
See (Crisanti & Sompolinsky, 2018) for further details and a worked example.

We now average over the sources of disorder. We assume that the entries of
$\bm{A}$ are i.i.d. with mean zero and variance 1. In the proportional limit,
we can replace the entries of $\bm{A}$ as a draw from a Gaussian
$\mathcal{N}(0,1)$ by appealing to Gaussian equivalence. We furhther justify
this in the cavity derivation in the next section. This allows us to evaluate
the averages over the matrix $\bm{A}$.

$\displaystyle\left<\exp\left(-\frac{i}{\sqrt{M}}\text{Tr}\bm{A}^{\top}\int
dt\left[\bm{\hat{}}{\bm{v}}^{3}(t)\bm{v}^{2}(t)^{\top}+\nu^{-1}\bm{v}^{3}(t)\bm{\hat{\bm{v}}}^{4}(t)^{\top}\right]\right)\right>_{\bm{A}}$
(38) $\displaystyle=\exp\left(-\frac{1}{2}\int
dtds[\hat{\bm{v}}^{3}(t)\cdot\hat{\bm{v}}^{3}(s)\underbrace{\frac{1}{M}\bm{v}^{2}(t)\cdot\bm{v}^{2}(s)}_{C_{2}(t,s)}+\nu^{-1}\hat{\bm{v}}^{4}(t)\cdot\hat{\bm{v}}^{4}(s)\underbrace{\frac{1}{N}\bm{v}^{3}(t)\cdot\bm{v}^{3}(s)}_{C_{3}(t,s)}]\right)$
$\displaystyle\quad\times\exp\left(-\int dtds\
\underbrace{\frac{1}{N}\hat{\bm{v}}^{3}(t)\cdot\bm{v}^{3}(s)}_{iR_{3}(s,t)}\,\bm{v}^{2}(t)\cdot\hat{\bm{v}}^{4}(s)\right).$

Similarly, we can calculate the averages over the data, which enters via the
design matrices $\bm{\Psi}$. Again in this proportional limit we can use a
Gaussian equivalence on $\bm{\Psi}$ to have it take the form
$\bm{\Psi}\sim\bm{\Phi}\bm{\Lambda}^{1/2}$ where $\bm{\Phi}$ has entries drawn
from a unit normal. Taking the average then gives us

$\displaystyle\left<\exp\left(-\frac{i}{\sqrt{M}}\text{Tr}\bm{\Psi}^{\top}\int
dt\left[\bm{\hat{}}{\bm{v}}^{1}(t)\bm{v}^{0}(t)^{\top}+\alpha^{-1}\bm{v}^{1}(t)\bm{\hat{\bm{v}}}^{2}(t)^{\top}\right]\right)\right>_{\bm{\Psi}}$
(39) $\displaystyle=\exp\left(-\frac{1}{2}\int
dtds[\hat{\bm{v}}^{1}(t)\cdot\hat{\bm{v}}^{1}(s)\underbrace{\frac{1}{M}\bm{v}^{0}(t)\cdot\bm{\Lambda}\bm{v}^{0}(s)}_{C_{0}(t,s)}+\alpha^{-1}\hat{\bm{v}}^{2}(t)\cdot\bm{\Lambda}\hat{\bm{v}}^{2}(s)\underbrace{\frac{1}{P}\bm{v}^{1}(t)\cdot\bm{v}^{1}(s)}_{C_{1}(t,s)}]\right),$
$\displaystyle\quad\times\exp\left(-\int dtds\
\underbrace{\frac{1}{P}\hat{\bm{v}}^{1}(t)\cdot\bm{v}^{1}(s)}_{iR_{1}(s,t)}\,\bm{v}^{0}(t)\cdot\bm{\Lambda}\cdot\hat{\bm{v}}^{2}(s)\right).$

We now insert delta functions for following bracketed terms:
$C_{0},C_{1},C_{2},C_{3}$ and $R_{1},R_{3}$ using the following identity (e.g.
for $C_{0}$ at times $s,t$):

$1=\int\frac{dC_{0}(s,t)d\hat{C}_{0}(s,t)}{2\pi
iM^{-1}}\exp\left[\frac{1}{2}M\int
dtds\,\hat{C}_{0}(t,s)\left(C_{0}(t,s)-\frac{1}{M}\bm{v}^{0}(t)\cdot\bm{\Lambda}\bm{v}^{0}(s)\right)\right].$
(40)

Here the $\hat{C}_{i},\hat{R}_{i}$ integrals are taken over the imaginary
axis. This yields a moment generating function (here we’ll take $\bm{j}=0$):

$\displaystyle
Z=\int\mathcal{D}[C_{1},\hat{C}_{1},\dots]\exp\left[MS[C_{0},C_{1},C_{2},C_{3},R_{1},R_{3},\hat{C}_{0},\hat{C}_{1},\hat{C}_{2},\hat{C}_{3},\hat{R}_{1},\hat{R}_{3}]\right].$
(41)

The constraint that $Z=1$ means that $S=0$ at the saddle point. $S$ here is
given by:

$\displaystyle S[\dots]$ $\displaystyle=\frac{1}{2}\int
dtds\Big{[}\hat{C}_{0}(t,s)C_{0}(t,s)+\alpha\hat{C}_{1}(t,s)C_{1}(t,s)+\hat{C}_{2}(t,s)C_{2}(t,s)+\nu\hat{C}_{3}(t,s)C_{3}(t,s)\Big{]}$
(42) $\displaystyle\quad+\int
dtds\Big{[}-R_{1}(t,s)\hat{R}_{1}(s,t)-R_{3}(t,s)\hat{R}_{3}(s,t)\Big{]}$
$\displaystyle\quad+\alpha\log\mathcal{Z}_{1}+\nu\log\mathcal{Z}_{3}+\frac{1}{M}\sum_{k}\log\mathcal{Z}_{0,2,4;k}.$

We have chosen to take $\hat{R}_{i}(s,t)$ to have a different sign and $s,t$
ordering convention than the $\hat{C}_{i}$ to simplify our notation later on.
We have also used that Equations (38), (39) factorize over their respective
indices, so each $\mathcal{Z}$ is a partition function over a single index.
The individual $\mathcal{Z}_{i}$ are given by:

$\displaystyle\mathcal{Z}_{1}=\int\mathcal{D}[v^{1},\hat{v}^{1}]\exp\left[i\int
dtds\left(\delta(t-s)-\alpha^{-1}\hat{R}_{1}(s,t)\right)v^{1}(t)\hat{v}^{1}(s)\right]$
(43) $\displaystyle\times\exp\left[-\frac{1}{2}\int
dtds(\hat{v}^{1}(t)\hat{v}^{1}(s)C_{0}(t,s)+v^{1}(t)v^{1}(s)\hat{C}_{1}(t,s))\right],$
$\displaystyle\mathcal{Z}_{3}=\int\mathcal{D}[v^{3},\hat{v}^{3}]\exp\left[i\int
dtds\left(\delta(t-s)-\nu^{-1}\hat{R}_{3}(s,t)\right)v^{3}(t)\hat{v}^{3}(s)\right]$
(44) $\displaystyle\times\exp\left[-\frac{1}{2}\int
dtds(\hat{v}^{3}(t)\hat{v}^{3}(s)C_{2}(t,s)+v^{3}(t)v^{3}(s)\hat{C}_{3}(t,s))\right],$
$\displaystyle\mathcal{Z}_{0,2,4;k}$
$\displaystyle=\int\mathcal{D}[v^{0,2,4},\hat{v}^{0,2,4}]\exp\left[-\frac{1}{2}\int
dtds\left(\alpha^{-1}\lambda_{k}\hat{v}^{2}_{k}(t)\hat{v}^{2}_{k}(s)C_{1}(t,s)+\nu^{-1}\hat{v}^{4}_{k}(t)\hat{v}^{4}_{k}(s)C_{3}(t,s)\right)\right]$
(45) $\displaystyle\qquad\times\exp\left[-\frac{1}{2}\int
dtds\left(v^{0}_{k}(t)v^{0}_{k}(s)\hat{C}_{0}(t,s)+v^{2}_{k}(t)v^{2}_{k}(s)\hat{C}_{2}(t,s)+v^{4}_{k}(t)v^{4}_{k}(s)\hat{C}_{4}(t,s)\right)\right]$
$\displaystyle\qquad\times\exp\left[-i\int
dtds\left(R_{3}(t,s)v^{2}_{k}(s)\hat{v}^{4}_{k}(t)+\lambda_{k}R_{1}(t,s)v^{0}_{k}(s)\hat{v}^{2}_{k}(t)\right)\right].$

In the large $M$ limit we evaluate this integral via saddle point. The saddle
point equations give:

$\displaystyle C_{i}(t,s)$ $\displaystyle=\langle
v_{i}(t)v_{i}(s)\rangle,\quad i=\\{0,1,2,3,4\\}$ (46)
$\displaystyle\hat{C}_{i}(t,s)$
$\displaystyle=\langle\hat{v}_{i}(t)\hat{v}_{i}(s)\rangle,\quad
i=\\{0,1,2,3,4\\}$ $\displaystyle R_{1}(t,s)$ $\displaystyle=-i\langle
v^{1}(t)\hat{v}^{1}(s)\rangle$ $\displaystyle R_{3}(t,s)$
$\displaystyle=-i\langle v^{3}(t)\hat{v}^{3}(s)\rangle$
$\displaystyle\hat{R}^{1}(t,s)$
$\displaystyle=-i\frac{1}{M}\sum_{k}\lambda_{k}\hat{v}_{k}^{2}(t)v_{k}^{0}(s)\equiv
R_{0,2}(t,s)$ $\displaystyle\hat{R}^{3}(t,s)$
$\displaystyle=-i\frac{1}{M}\sum_{k}\lambda_{k}\hat{v}_{k}^{4}(t)v_{k}^{2}(s)\equiv
R_{2,4}(t,s).$

Here $\langle\cdot\rangle$ denotes an average taken with respect to the
statistical ensemble given by the corresponding partition function
$\mathcal{Z}_{i}$. Following the discussion below Equation 37, we take
$\hat{C}_{i}(t,s)=0$, which will enforce
$\langle\hat{v}_{i}(t)\hat{v}_{i}(s)\rangle=0$ and lead to the correct
dynamical equations.

To evaluate the remaining, we can integrate out the $\hat{v}^{i}$ variables.
First let us look at $\mathcal{Z}_{1}$. Using the Hubbard-Stratonovich trick
we can write the action in terms linear in $\hat{v}^{1}$. This gives

$\displaystyle\mathcal{Z}_{1}$
$\displaystyle=\int\mathcal{D}[v^{1},\hat{v}^{1},u^{1}]\exp\left[i\int
dtds\,\hat{v}^{1}(t)\left[\delta(t-s)(v^{1}(s)-u(s))-\alpha^{-1}\hat{R}_{1}(t,s)v^{1}(s)\right]\right]$
(47) $\displaystyle\quad\times\exp\left[-\frac{1}{2}\int
dtds\,u(t)u(s)C^{-1}_{0}(t,s)+v^{1}(t)v^{1}(s)\hat{C}_{1}(t,s))\right]$

We now replace $\hat{C}_{1}$ by its saddle point value of $0$ and
$\hat{R}_{1}$ by $R_{0,2}$. Integrating over $\hat{v}$ gives a delta function:

$v^{1}(t)=u^{1}(t)+\frac{1}{\alpha}\int dsR_{0,2}(t,s)v_{1}(s),\quad
u^{1}(t)\sim\mathcal{GP}(0,C_{0}).$ (48)

Analogously for $v^{3}$ we get

$v^{3}(t)=u^{3}(t)+\frac{1}{\nu}\int dsR_{2,4}(t,s)v^{3}(s),\quad
u^{3}(t)\sim\mathcal{GP}(0,C_{2})$ (49)

For $\mathcal{Z}_{0,2,4;k}$ after replacing
$\hat{C}_{0},\hat{C}_{2},\hat{C}_{4}$ with their saddle point values we get:

$\displaystyle\mathcal{Z}_{0,2,4;k}$ $\displaystyle=\exp\left[-\frac{1}{2}\int
dtds\left(\alpha^{-1}\lambda_{k}\hat{v}^{2}_{k}(t)\hat{v}^{2}_{k}(s)C_{1}(t,s)+\nu^{-1}\hat{v}^{4}_{k}(t)\hat{v}^{4}_{k}(s)C_{3}(t,s)\right)\right]$
(50) $\displaystyle\qquad\times\exp\left[-\int
dtds\left(iR_{3}(t,s)v^{2}_{k}(s)\hat{v}^{4}_{k}(t)+i\lambda_{k}R_{1}(t,s)v^{0}_{k}(s)\hat{v}^{2}_{k}(t)\right)\right]$

Using the same Hubbard-Stratonovich trick on $\hat{v}^{2}_{k}$ gives:

$v^{2}_{k}(t)=u^{2}_{k}(t)+\lambda_{k}\int dsR_{1}(t,s)v_{k}^{0}(s),\quad
u^{2}_{k}(t)\sim\mathcal{GP}(0,\frac{1}{\alpha}\lambda_{k}C_{1}).$ (51)

On $\hat{v}^{4}_{k}$ we similarly get:

$v^{4}_{k}(t)=u^{4}_{k}(t)+\int dsR_{3}(t,s)\cdot v^{2}_{k}(s),\quad
u^{4}_{k}(t)\sim\mathcal{GP}(0,\frac{1}{\nu}C_{3})$ (52)

Lastly, the equations of motion for $v^{0}_{k}$ in terms of $v^{4}_{k}$ are
known:

$\partial_{t}v^{0}_{k}(t)=-v^{4}_{k}(t).$ (53)

One can easily add momentum by replacing $\partial_{t}v^{0}_{k}(t)$ with
$(\beta\partial^{2}_{t}+\partial_{t})v^{0}_{k}(t)$ without changing anything
else about the derivation.

####  B.1.1 Interpretation of the Response Functions

Following (Crisanti & Sompolinsky, 2018; Helias & Dahmen, 2020), we can
understand the $\langle\hat{v}^{a}(t)v^{b}(s)\rangle$ correlators by adding in
the single-site moment generating function (e.g. Equation 47) a source
$\tilde{j}^{b}(s)$ that couples to $\hat{v}^{b}$ at time $s$. As in the
discussion below equation 37, differentiating $\langle v^{a}(t)\rangle$ with
respect to that source corresponds to its response to a kick in the dynamics
of $v^{b}$ at time $s$. We denote this by:

$R_{i,j}(t,s)=\left<\frac{\delta v^{i}(t)}{\delta v^{j}(s)}\right>.$ (54)

###  B.2 Cavity Derivation

The cavity derivation relies on Taylor expanding the dynamics upon the
addition of a new sample or feature. We will work through each cavity step one
at a time by considering the influence of a single new base feature, new
sample, and new projected feature. In each step, the goal is to compute the
marginal statistics of the added variables. This requires tracking the linear
response to all other variables in the system.

##### Adding a Base Feature

Upon addition of a base feature with eigenvalue $\lambda_{0}$ so that there
are $M+1$ instead of $M$ features $\\{v^{0}_{k},v^{2}_{k},v^{4}_{k}\\}$ for
$k\in\\{0,1,...,M\\}$, we have a perturbation to both $v^{1}_{\mu}(t)$ and
$v^{3}_{n}(t)$. Denote the perturbed versions of the dynamics upon addition of
the $M+1$st feature as $\tilde{v}^{1}_{\mu}(t)$ and $\tilde{v}^{3}_{n}(t)$. At
large $M$ we can use linear-response theory to relate the dynamics at $M+1$
features to the dynamics of the original $M$ feature system

$\displaystyle\tilde{v}^{1}_{\mu}(t)\sim{v}^{1}_{\mu}(t)+\frac{1}{\sqrt{M}}\sum_{\nu=1}^{P}\int
ds\frac{\partial v^{1}_{\mu}(t)}{\partial
v^{1}_{\nu}(s)}\psi^{\nu}_{0}v^{0}_{0}(s)$
$\displaystyle\tilde{v}^{3}_{n}(t)\sim{v}^{1}_{\mu}(t)+\frac{1}{\sqrt{M}}\sum_{m=1}^{N}\int
ds\frac{\partial v^{3}_{n}(t)}{\partial v^{3}_{m}(s)}A_{m0}\ v^{2}_{0}(s)$
(55)

The next order corrections have a subleading influence on the dynamics. Now,
inserting these perturbed dynamics into the dynamics for the new $(M+1)$st set
of variables $\\{v^{2}_{0}(t),v^{4}_{0}(t)\\}$. For $v^{2}_{0}(t)$, we have

$\displaystyle
v^{2}_{0}(t)\sim\frac{1}{\alpha\sqrt{M}}\sum_{\mu=1}^{P}\psi^{\mu}_{0}v^{1}_{\mu}(t)+\frac{1}{\alpha
M}\sum_{\mu,\nu=1}^{P}\int ds\ \psi^{\mu}_{0}\frac{\partial
v^{1}_{\mu}(t)}{\partial v^{1}_{\nu}(s)}\psi^{\nu}_{0}v^{0}_{0}(s)$ (56)

There are now two key steps in simplifying the above expression in the
proportional limit:

  1. 1.

By the fact that the $v^{1}_{\mu}(t)$ dynamics are statistically independent
of the new feature $\psi^{\mu}_{0}$, we can invoke a central limit theorem for
the first term which is mean zero and variance $\mathcal{O}(1)$.

  2. 2.

Similarly, we can invoke a law of large numbers for the second term, which has
$\mathcal{O}(1)$ mean and variance on the order of $\mathcal{O}(M^{-1})$.
Therefore in the asymptotic limit it can be safely approximated by its mean.

We note in passing that neither of these steps require the $\psi^{\mu}_{0}$
variables to be Gaussian. Thus we obtain the following asymptotic statistical
description of the $v^{2}_{0}(t)$ random variable

$\displaystyle v^{2}_{0}(t)$ $\displaystyle\sim u^{2}_{0}(t)+\int
dsR_{1}(t,s)v^{0}_{0}(s)$ $\displaystyle u^{2}_{0}(t)$
$\displaystyle\sim\mathcal{GP}(0,\alpha^{-1}\lambda_{0}C_{1})\ ,\
R_{1}(t,s)\equiv\frac{1}{P}\sum_{\mu=1}^{P}\left<\frac{\partial
v^{1}_{\mu}(t)}{\partial v^{1}_{\mu}(s)}\right>.$ (57)

Following an identical argument for $v^{4}_{0}(t)$ we have

$\displaystyle v^{4}_{0}(t)$
$\displaystyle\sim\frac{1}{\nu\sqrt{M}}\sum_{n=1}^{M}A_{n0}v^{3}_{n}(t)+\frac{1}{\nu
M}\sum_{nm}\int ds\ A_{n0}\frac{\partial v^{3}_{n}(t)}{\partial
v^{3}_{m}(s)}A_{m0}v^{2}_{0}(s)$ $\displaystyle\sim u^{3}_{0}(t)+\int
dsR_{3}(t,s)v^{2}_{0}(s)$ $\displaystyle
u_{0}^{3}(t)\sim\mathcal{GP}(0,\nu^{-1}C_{3})\ ,\
R_{3}(t,s)=\frac{1}{N}\sum_{n=1}^{N}\left<\frac{\partial
v^{3}_{n}(t)}{\partial v^{3}_{n}(s)}\right>$ (58)

##### Adding a Sample

Next, we can consider the influence of adding a new data point. We will aim to
characterize a $P+1$ data point system in terms of the dynamics when $P$
points are present. Upon the addition of a new data point $\bm{\psi}^{0}$ the
field $v_{k}^{0}(t)$ will be perturbed to $\tilde{v}_{k}^{0}(t)$. Again
invoking linear response theory, we can expand the perturbed value around the
$P$-sample dynamics

$\displaystyle\tilde{v}^{0}_{k}(t)\sim
v^{0}_{k}(t)+\frac{1}{\alpha\sqrt{M}}\sum_{\ell=1}^{M}\int ds\ \frac{\partial
v^{0}_{k}(t)}{\partial v^{2}_{\ell}(s)}\psi^{0}_{\ell}v^{1}_{0}(s)$ (59)

Now, computing the dynamics of the new random variable $v^{1}_{0}(t)$

$\displaystyle v^{1}_{0}(t)$
$\displaystyle\sim\frac{1}{\sqrt{M}}\sum_{k=1}^{M}\psi^{0}_{k}v^{0}_{k}(t)+\frac{1}{\alpha
M}\int ds\sum_{k\ell}\psi^{0}_{k}\ \frac{\partial v^{0}_{k}(t)}{\partial
v^{2}_{\ell}(s)}\ \psi_{\ell}^{0}\ v^{1}_{0}(s)$ $\displaystyle\sim
u^{1}_{0}(t)+\frac{1}{\alpha}\int dsR_{0,2}(t,s)v^{1}_{0}(s)$ $\displaystyle
u^{1}_{0}(t)$
$\displaystyle\sim\mathcal{GP}\left(0,\frac{1}{M}\sum_{k}\lambda_{k}C^{0}_{k}\right)\
,\ R_{0,2}(t,s)=\frac{1}{M}\sum_{k=1}^{M}\lambda_{k}\left<\frac{\partial
v^{0}_{k}(t)}{\partial v^{2}_{k}(s)}\right>$ (60)

##### Adding a Projected Feature

Now, we finally consider the effect of introducing a single new projected
feature so that instead of $N$ we now have $N+1$ projected features. This
causes a perturbation to $\\{v^{2}_{k}(t)\\}$ which we

$\displaystyle\tilde{v}^{2}_{k}(t)\sim
v^{2}_{k}(t)+\frac{1}{\nu\sqrt{M}}\sum_{\ell=1}^{M}\int ds\ \frac{\partial
v^{2}_{k}(t)}{\partial v^{4}_{\ell}(s)}A_{0\ell}\ v^{3}_{0}(s)$ (61)

Now, we compute the dynamics for the added variable $v^{3}_{0}(t)$

$\displaystyle v^{3}_{0}(t)$
$\displaystyle\sim\frac{1}{\sqrt{M}}\sum_{k=1}^{M}A_{0k}v^{2}_{k}(t)+\frac{1}{\nu
M}\sum_{k\ell}\int dsA_{0k}\frac{\partial v^{2}_{k}(t)}{\partial
v^{4}_{\ell}(s)}A_{0\ell}\ v^{3}_{0}(s)$ $\displaystyle\sim
u^{3}_{0}(t)+\frac{1}{\nu}\int dsR_{2,4}(t,s)v^{3}_{0}(s)$ $\displaystyle
u^{3}_{0}(t)$ $\displaystyle\sim\mathcal{GP}(0,C_{2})\ ,\
R_{2,4}(t,s)=\frac{1}{M}\sum_{k=1}^{M}\left<\frac{\partial
v^{2}_{k}(t)}{\partial v^{4}_{k}(s)}\right>$ (62)

##### Putting it all together

Now, using the information gained in the previous sections, we can combine all
of the dynamics for each field into a closed set of stochastic processes. This
recovers the DMFT equations of Appendix A.2.

##  Appendix C Final Losses (the $t\to\infty$ Limit of DMFT)

In this section we work out exact expressions for the large time limit of
DMFT. By comparing with prior computations of the mean-field statics of this
problem computed in (Atanasov et al., 2023; Zavatone-Veth & Pehlevan, 2023;
Ruben & Pehlevan, 2023; Maloney et al., 2022; Simon et al., 2021), we show
that the large time and large $M$ limits commute, specifically that
$\lim_{M,N,P\to\infty}\lim_{t\to\infty}\mathcal{L}(M,N,P,t)=\lim_{t\to\infty}\lim_{M,N,P\to\infty}\mathcal{L}(M,N,P,t)$.
We invoke the final value theorem and use the response functions as before.

##### Final Value Theorem

We note that for functions which vanish at $t=-\infty$, that

$\displaystyle\lim_{\omega\to 0}i\omega\ \mathcal{H}(\omega)=-\lim_{\omega\to
0}\int_{-\infty}^{\infty}d\tau\left[\frac{\partial}{\partial\tau}e^{-i\omega\tau}\right]H(\tau)=\lim_{\omega\to
0}\int_{-\infty}^{\infty}d\tau\left[\frac{\partial}{\partial\tau}H(\tau)\right]e^{-i\omega\tau}=\lim_{\tau\to\infty}H(\tau)$
(63)

where we invoked integration by parts and used the assumption that
$\lim_{\tau\to-\infty}H(\tau)=0$, a condition that is satisfied for the
correlation and response functions in our theory. We can therefore use the
identity $\lim_{\tau\to\infty}H(\tau)=\lim_{\omega\to
0}i\omega\mathcal{H}(\omega)$ to extract the final values of our order
parameters.

$\displaystyle\lim_{\tau\to\infty}H_{k}(\tau)=\lim_{\omega\to
0}i\omega\mathcal{H}_{k}(\omega)=\lim_{\omega\to
0}\frac{1}{1+\lambda_{k}(i\omega)^{-1}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}.$
(64)

We also need to invoke a similar relationship for the final values of the
correlation functions

$\displaystyle\lim_{t,s\to\infty}C(t,s)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\ \mathcal{C}(\omega,\omega^{\prime})$
$\displaystyle\mathcal{C}(\omega,\omega^{\prime})=\int dt\int ds\ e^{-i(\omega
t+\omega^{\prime}s)}C(t,s)$ (65)

where $\mathcal{C}$ is the two-variable Fourier transform. The final value of
the test loss is
$\lim_{t\to\infty}\mathcal{L}(t)=\lim_{t,s\to\infty}C_{0}(t,s)$.

###  C.1 General Case (Finite $\nu,\alpha$)

Before working out the solution to the response functions, we note that the
following condition is always satisfied

$\displaystyle\nu(1-\mathcal{R}_{3}(\omega))=\alpha(1-\mathcal{R}_{1}(\omega)).$
(66)

For $\nu=\alpha$, this equation implies that
$\mathcal{R}_{1}=\mathcal{R}_{3}$. For $\nu\neq\alpha$, we can have either
$\mathcal{R}_{1}\to 0$ or $\mathcal{R}_{3}\to 0$ but not both. We consider
each of these cases below.

##### Over-parameterized Case $\nu>\alpha$:

In this case, the response function $\mathcal{R}_{1}\sim\mathcal{O}(i\omega)$
as $\omega\to 0$ and $\mathcal{R}_{3}\sim 1-\frac{\alpha}{\nu}$ as $\omega\to
0$. We thus define

$\displaystyle r\equiv\lim_{\omega\to
0}(i\omega)^{-1}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)$ (67)

Using the equation which defines $\mathcal{R}_{1}$, we find that the variable
$r$ satisfies the following relationship at $\omega\to 0$

$\displaystyle\alpha=\frac{1}{M}\sum_{k}\frac{\lambda_{k}r}{1+\lambda_{k}r}$
(68)

After solving this implicit equation, we can find the limiting value of
$i\omega\mathcal{H}_{k}(\omega)$ as

$\displaystyle
H_{k}^{\infty}=\lim_{\tau\to\infty}{H}_{k}(\tau)=\lim_{\omega\to
0}i\omega\mathcal{H}_{k}(\omega)=\frac{1}{1+\lambda_{k}r}$ (69)

Next, we can work out the scaling of the correlation functions in the limit of
low frequency. We define the following limiting quantities based on a scaling
analysis performed on our correlation functions for small $\omega$

$\displaystyle C_{0}^{\infty}$
$\displaystyle\equiv\lim_{t,s\to\infty}C_{0}(t,s)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\mathcal{C}_{0}(\omega,\omega^{\prime})$ (70)
$\displaystyle C_{1}^{\infty}$
$\displaystyle\equiv\int_{0}^{\infty}\int_{0}^{\infty}dt^{\prime}ds^{\prime}C_{1}(t^{\prime},s^{\prime})=\lim_{\omega,\omega^{\prime}\to
0}\mathcal{C}_{1}(\omega,\omega^{\prime})$ $\displaystyle C_{2}^{\infty}$
$\displaystyle\equiv\int_{0}^{\infty}\int_{0}^{\infty}dt^{\prime}ds^{\prime}C_{2}(t^{\prime},s^{\prime})=\lim_{\omega,\omega^{\prime}\to
0}\mathcal{C}_{2}(\omega,\omega^{\prime})$ $\displaystyle C_{3}^{\infty}$
$\displaystyle\equiv\int_{0}^{\infty}\int_{0}^{\infty}dt^{\prime}ds^{\prime}C_{3}(t^{\prime},s^{\prime})=\lim_{\omega,\omega^{\prime}\to
0}\mathcal{C}_{3}(\omega,\omega^{\prime})$

These limiting quantities satisfy the closed set of linear equations

$\displaystyle
C_{0}^{\infty}=\frac{1}{M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}C_{3}^{\infty}+\frac{1}{\alpha}\lambda_{k}\left(1-\frac{\alpha}{\nu}\right)^{2}C_{1}^{\infty}\right]$
(71) $\displaystyle
C_{1}^{\infty}=\frac{r^{2}}{(1-\frac{\alpha}{\nu})^{2}}C_{0}^{\infty}$
$\displaystyle C_{2}^{\infty}=\frac{1}{\alpha
M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}C_{1}^{\infty}+\frac{r^{2}}{(1-\frac{\alpha}{\nu})^{2}}\frac{1}{M}\sum_{k}\lambda_{k}^{2}(H_{k}^{\infty})^{2}[(w^{\star}_{k})^{2}+\nu^{-1}C_{3}^{\infty}]$
$\displaystyle
C_{3}^{\infty}=\left(1-\frac{\alpha}{\nu}\right)^{2}C_{2}^{\infty}$

These equations can be solved for
$\\{C_{0}^{\infty},C_{1}^{\infty},C_{2}^{\infty},C_{3}^{\infty}\\}$.
Simplifying the expressions to a two-variable system, we find

$\displaystyle C_{0}^{\infty}$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}\left(1-\frac{\alpha}{\nu}\right)^{2}C_{2}^{\infty}+\frac{1}{\alpha}\lambda_{k}r^{2}C_{0}^{\infty}\right]$
$\displaystyle C_{2}^{\infty}$
$\displaystyle=\frac{r^{2}}{\alpha(1-\alpha/\nu)^{2}M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}C_{0}^{\infty}+\frac{r^{2}}{(1-\frac{\alpha}{\nu})^{2}}\frac{1}{M}\sum_{k}\lambda_{k}^{2}(H_{k}^{\infty})^{2}\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}\left(1-\frac{\alpha}{\nu}\right)^{2}C_{2}^{\infty}\right]$

This expression recovers the ridgeless limit of the replica results of
(Atanasov et al., 2023; Zavatone-Veth & Pehlevan, 2023) and the random matrix
analysis of (Simon et al., 2023).

##### Under-parameterized Case $\nu<\alpha$:

Following the same procedure, we note that for $\nu<\alpha$ that
$\mathcal{R}_{3}\sim\mathcal{O}(i\omega)$ and $\mathcal{R}_{1}\sim
1-\frac{\nu}{\alpha}$. We thus find the following equation for
$r=\lim_{\omega\to\infty}(i\omega)^{-1}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)$.

$\displaystyle\nu=\frac{1}{M}\sum_{k}\frac{\lambda_{k}r}{\lambda_{k}r+1}$ (72)

where as before $H_{k}^{\infty}=\frac{1}{1+\lambda_{k}r}$. The analogous
scaling argument for small $\omega$ gives us the following set of well-defined
limiting quantities

$\displaystyle C_{0}^{\infty}$
$\displaystyle\equiv\lim_{t,s\to\infty}C_{0}(t,s)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\mathcal{C}(\omega,\omega^{\prime})$ (73)
$\displaystyle C_{1}^{\infty}$
$\displaystyle\equiv\lim_{t,s\to\infty}C_{1}(t,s)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\mathcal{C}_{1}(\omega,\omega^{\prime})$
$\displaystyle C_{2}^{\infty}$
$\displaystyle\equiv\lim_{t,s\to\infty}C_{2}(t,s)=\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\mathcal{C}_{2}(\omega,\omega^{\prime})$
$\displaystyle C_{3}^{\infty}$
$\displaystyle\equiv\int_{0}^{\infty}\int_{0}^{\infty}dtds\
C_{3}(t,s)=\lim_{\omega,\omega^{\prime}\to
0}\mathcal{C}_{3}(\omega,\omega^{\prime}).$

where these limiting correlation values satisfy

$\displaystyle C_{0}^{\infty}$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}C_{3}^{\infty}+\frac{1}{\alpha}\lambda_{k}\frac{r^{2}}{(1-\frac{\nu}{\alpha})^{2}}C_{1}^{\infty}\right]$
(74) $\displaystyle C_{1}^{\infty}$
$\displaystyle=\left(1-\frac{\nu}{\alpha}\right)^{2}C_{0}^{\infty}$
$\displaystyle C_{2}^{\infty}$ $\displaystyle=\frac{1}{\alpha
M}\sum_{k}\lambda_{k}(H_{k}^{\infty})^{2}C_{1}^{\infty}+\frac{1}{M}\left(1-\frac{\nu}{\alpha}\right)^{2}\sum_{k}\lambda_{k}^{2}(H_{k}^{\infty})^{2}\left[(w^{\star}_{k})^{2}+\frac{1}{\nu}C_{3}^{\infty}\right]$
$\displaystyle C_{3}^{\infty}$
$\displaystyle=\frac{r^{2}}{(1-\frac{\nu}{\alpha})^{2}}C_{2}^{\infty}.$

This is again a closed linear system of equations for the variables
$\\{C_{0}^{\infty},C_{1}^{\infty},C_{2}^{\infty},C_{3}^{\infty}\\}$. In the
next section, we recover the result for kernel regression where $\nu\to\infty$
and the learning curve for infinite data $\alpha\to\infty$ with respect to
model size $\nu$.

###  C.2 Learning Curves for Kernel Regression $\nu,t\to\infty$

In the $t\to\infty$ and $\nu\to\infty$ limit we recover the learning curve for
kernel regression with eigenvalues $\lambda_{k}$. To match the notation of
(Canatar et al., 2021), we define

$\displaystyle\lim_{\omega\to
0}(i\omega)^{-1}\mathcal{R}_{1}(\omega)\equiv\alpha\kappa^{-1}$ (75)

which generates the following self-consistent equation for $\kappa$

$\displaystyle
1=\frac{1}{M}\sum_{k}\frac{\lambda_{k}}{\lambda_{k}\alpha+\kappa}.$ (76)

Plugging this into the expression for the loss, we find

$\displaystyle(i\omega)(i\omega^{\prime})\mathcal{C}_{0}(\omega,\omega^{\prime})\sim\frac{1}{M}\sum_{k}\lambda_{k}\frac{\kappa^{2}}{(\kappa+\lambda_{k}\alpha)^{2}}[(w^{\star}_{k})^{2}+\alpha^{-1}\lambda_{k}\mathcal{C}_{1}(\omega,\omega^{\prime})]$
$\displaystyle\mathcal{C}_{1}(\omega,\omega^{\prime})=(i\omega)(i\omega^{\prime})\alpha^{2}\kappa^{-2}\mathcal{C}_{0}(\omega,\omega^{\prime})$
(77)

Letting $C_{\infty}\equiv\lim_{\omega,\omega^{\prime}\to
0}(i\omega)(i\omega^{\prime})\mathcal{C}(s,s^{\prime})$, we have

$\displaystyle C_{\infty}$
$\displaystyle=\frac{1}{M}\sum_{k}\lambda_{k}\frac{\kappa^{2}}{(\kappa+\lambda_{k}\alpha)^{2}}(w^{\star}_{k})^{2}+C_{\infty}\frac{\alpha}{M}\sum_{k}\frac{\lambda_{k}^{2}}{(\lambda_{k}\alpha+\kappa)^{2}}$
$\displaystyle=\frac{1}{1-\gamma}\sum_{k}\lambda_{k}(w^{\star}_{k})^{2}\frac{\kappa^{2}}{(\kappa+\lambda_{k}\alpha)^{2}}\
,\quad\gamma=\frac{\alpha}{M}\sum_{k}\frac{\lambda_{k}^{2}}{(\lambda_{k}\alpha+\kappa)^{2}}$
(78)

The variable $\kappa$ decreases from
$\left[\frac{1}{M}\sum_{k}\lambda_{k},0\right]$ as $\alpha\in[0,1]$. For
$\alpha>1$ we have $\kappa=0$. The quantity $\frac{1}{1-\gamma}$ comes from
overfitting due to variance from the randomly sampled dataset.

##  Appendix D Early Time Dynamics (High-Frequency Range)

In this section, we explore the early time dynamical effects of this model.
Similar to how the late time dynamical effects could be measured by examining
the low frequency $\omega\ll 1$ part of the response and correlation
functions, in this section, we analyze the high frequency components
$\omega\gg 1$. We start by noting the following expansions valid near
$\omega\to\infty$

$\displaystyle\mathcal{R}_{1}(\omega)$ $\displaystyle\sim
1-\frac{1}{\alpha(i\omega)}\left[\frac{1}{M}\sum_{k}\lambda_{k}\right]+\mathcal{O}(\omega^{-2})$
$\displaystyle\mathcal{R}_{3}(\omega)$ $\displaystyle\sim
1-\frac{1}{\nu(i\omega)}\left[\frac{1}{M}\sum_{k}\lambda_{k}\right]+\mathcal{O}(\omega^{-2})$
(79)

We let $c=\frac{1}{M}\sum_{k}\lambda_{k}$. These can be plugged into the
transfer function for mode $k$

$\displaystyle\mathcal{H}_{k}(\omega)\sim\frac{1}{i\omega+\lambda_{k}-c\alpha^{-1}(i\omega)^{-1}-c\nu^{-1}(i\omega)^{-1}}\sim\frac{1}{i\omega+\lambda_{k}}\left[1+c\alpha^{-1}(i\omega)^{-1}+c\nu^{-1}(i\omega)^{-1}\right]+\mathcal{O}\left(\omega^{-2}+(\omega^{\prime})^{-2}\right)$
(80)

Preforming an inverse Fourier transform, we find the following early time
asymptotics

$\displaystyle H_{k}(t)\sim
e^{-\lambda_{k}t}+\frac{c}{\lambda_{k}}\left(\frac{1}{\alpha}+\frac{1}{\nu}\right)\left[1-e^{-\lambda_{k}t}\right]$
(81)

We see from this expression that the early time corrections always scale as
$1/\alpha$ or $1/\nu$ and that these corrections build up over time. A similar
expansion can be performed for $\mathcal{C}(\omega,\omega^{\prime})$ with
$\omega,\omega^{\prime}\gg 1$ which also gives leading corrections which scale
as $1/\alpha$ and $1/\nu$.

##  Appendix E Buildup of Overfitting Effects

In this section, we derive a formula for the gap between test loss
$\mathcal{L}(t)$ and train loss $\hat{\mathcal{L}}(t)$. We start from the
following formula

$\displaystyle
v_{1}(t)=u_{1}(t)+\frac{1}{\alpha}\int_{0}^{t}dsR_{0,2}(t,s)v_{1}(s)$ (82)

Moving the $v_{1}(t)$ term to the other side, and using the fact that
$\left<u_{1}(t)u_{1}(s)\right>=C_{1}(t,s)$, we find the following relationship
between train and test loss

$\displaystyle\mathcal{L}(t)$
$\displaystyle=\left<u_{1}(t)u_{1}(t)\right>=\left<v_{1}(t)v_{1}(t)\right>-\frac{2}{\alpha}\int_{0}^{t}dt^{\prime}R_{0,2}(t,t^{\prime})\left<v_{1}(t)v_{1}(t^{\prime})\right>$
$\displaystyle+\frac{1}{\alpha^{2}}\int_{0}^{t}dt^{\prime}\int_{0}^{t}ds^{\prime}R_{0,2}(t,t^{\prime})R_{0,2}(t,s^{\prime})\left<v_{1}(t^{\prime})v_{1}(s^{\prime})\right>$
$\displaystyle=\hat{\mathcal{L}}(t)-\frac{2}{\alpha}\int_{0}^{t}dt^{\prime}R_{0,2}(t,t^{\prime})C_{1}(t,t^{\prime})+\frac{1}{\alpha^{2}}\int_{0}^{t}dt^{\prime}\int_{0}^{t}ds^{\prime}R_{0,2}(t,t^{\prime})R_{0,2}(t,s^{\prime})C_{1}(t^{\prime},s^{\prime}).$
(83)

To get a sense of these expressions at early and late timescales, we
investigate the Fourier transforms at high $\omega\gg 1$ and low $\omega\ll 1$
frequencies respectively.

###  E.1 High Frequency Range / Early Time

The relationship between Fourier transforms at high frequencies $\omega\gg 1$
is

$\displaystyle\mathcal{C}_{0}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})}\
\mathcal{C}_{1}(\omega,\omega^{\prime})\sim\mathcal{C}_{1}(\omega,\omega^{\prime})+\frac{c}{\alpha(i\omega^{\prime})}\mathcal{C}_{1}(\omega,\omega^{\prime})+\frac{c}{\alpha(i\omega^{\prime})}\mathcal{C}_{1}(\omega,\omega^{\prime})+\mathcal{O}((i\omega)^{-2}+(i\omega^{\prime})^{-2})$
(84)

where $c=\frac{1}{M}\sum_{k}\lambda_{k}$. Taking a Fourier transform back to
real time gives us the following early time differential equation for the
test-loss train loss gap

$\displaystyle\partial^{2}_{ts}[C_{0}(t,s)-C_{1}(t,s)]=\partial^{2}_{ts}C_{1}(t,s)+\frac{c}{\alpha}(\partial_{t}+\partial_{s})C_{1}(t,s).$
(85)

The above equation should hold for early times. We note that
$C_{0}(t,t)-C_{1}(t,t)=\mathcal{L}(t)-\mathcal{\hat{L}}(t)$ exactly recovers
the test-train gap.

###  E.2 Low Frequency Range/Late Time

At late time/low frequency, as we showed in Appendix C, the behavior of the
$C_{1}$ correlation function depends on whether the model is over-
parameterized or under-parameterized. In the overparameterized case, the
asymptotic train loss is zero while the asymptotic test loss is nonzero. In
the underparameterized case, we have a limiting value for both the test and
train loss which can be computed from the expressions in Appendix C.

##  Appendix F Non-Proportional (Dimension-Free) Limit

We can imagine a situation where the original features are already infinite
dimensional ($M\to\infty$ is taken first). This would correspond more
naturally to the connection between infinite dimensional RKHS’s induced by
neural networks at infinite width (Bordelon et al., 2020; Canatar et al.,
2021; Cheng & Montanari, 2022). Further, we will assume a trace class kernel
$K(\bm{x},\bm{x}^{\prime})=\bm{\psi}(\bm{x})\cdot\bm{\psi}(\bm{x}^{\prime})$
for the base features $\bm{\psi}$ which diagonalizes over the data
distribution $p(\bm{x})$ as

$\displaystyle\int
K(\bm{x},\bm{x}^{\prime})\phi_{k}(\bm{x}^{\prime})p(\bm{x}^{\prime})d\bm{x}^{\prime}=\lambda_{k}\phi_{k}(\bm{x})\quad,\quad\sum_{k=1}^{\infty}\lambda_{k}<\infty.$
(86)

As before, we are concerned with the test and train losses

$\displaystyle\mathcal{L}(t)=\sum_{k=1}^{\infty}\lambda_{k}v^{0}_{k}(t)^{2}\
,\ \mathcal{\hat{L}}(t)=\frac{1}{P}\sum_{\mu=1}^{P}v^{1}_{\mu}(t)^{2}.$ (87)

The appropriate scaling of our four fields of interest in this setting are

$\displaystyle\bm{v}^{1}(t)=\bm{\Psi}\bm{v}^{0}(t)\ ,\
\bm{v}^{2}=\frac{1}{P}\bm{\Psi}^{\top}\bm{v}^{1}(t)$
$\displaystyle\bm{v}^{3}(t)=\bm{A}\bm{v}^{2}(t)\ ,\
\bm{v}^{4}=\frac{1}{N}\bm{A}^{\top}\bm{v}^{3}(t).$ (88)

Following the cavity argument given in the previous section, we can
approximate the the correlation and response functions as concentrating to
arrive at the following field description of the training dynamics

$\displaystyle\partial_{t}v^{0}_{k}(t)=-v^{4}_{k}(t)$ (89) $\displaystyle
v^{1}(t)=u^{1}(t)+\frac{1}{P}\int ds\ R_{0,2}(t,s)v^{1}(s)\ ,\
u^{1}(t)\sim\mathcal{GP}(0,C_{0})$ $\displaystyle
v^{2}_{k}(t)=u^{2}_{k}(t)+\int ds\ R_{1}(t,s)v^{0}_{k}(s)\ ,\
u^{2}_{k}(t)\sim\mathcal{GP}\left(0,\frac{1}{P}\lambda_{k}C^{1}\right)$
$\displaystyle v^{3}(t)=u^{3}(t)+\frac{1}{N}\int ds\ R_{2,4}(t,s)v^{3}(s)\ ,\
u^{3}(t)\sim\mathcal{GP}(0,C_{2})$ $\displaystyle
v^{4}_{k}(t)=u^{4}_{k}(t)+\int ds\ R_{3}(t,s)v^{2}_{k}(s)\ ,\
u^{4}_{k}(t)\sim\mathcal{N}\left(0,\frac{1}{N}C_{3}\right).$

which are exactly the same equations as in the proportional limit except with
the substitution $\nu\to N$ and $\alpha\to P$. The correlation and response
functions have the form

$\displaystyle
C_{0}(t,s)=\sum_{k=1}^{\infty}\lambda_{k}\left<v^{0}_{k}(t)v^{0}_{k}(s)\right>\
,\ C_{1}(t,s)=\left<v^{1}(t)v^{1}(s)\right>$ (90) $\displaystyle
C_{2}(t,s)=\sum_{k=1}^{\infty}v^{2}_{k}(t)v^{2}_{k}(s)\ ,\
C_{3}(t,s)=\left<v^{3}(t)v^{3}(s)\right>$ $\displaystyle
R_{0,2}(t,s)=\sum_{k=1}^{\infty}\lambda_{k}\left<\frac{\delta
v^{0}_{k}(t)}{\delta u^{2}_{k}(s)}\right>\ ,\ R_{1}(t,s)=\left<\frac{\delta
v^{1}(t)}{\delta u^{1}(s)}\right>$ $\displaystyle
R_{2,4}(t,s)=\sum_{k=1}^{\infty}\left<\frac{\delta v^{2}_{k}(t)}{\delta
u^{4}_{k}(s)}\right>\ ,\ R_{3}(t,s)=\left<\frac{\delta v^{3}(t)}{\delta
u^{3}(s)}\right>$

which will all be $\mathcal{O}(1)$ under this scaling.

##  Appendix G Effect of Ensembling and Bagging on Dynamics

###  G.1 Definition of Bias and Variance

We adopt the language of the fine-grained bias-variance decomposition in
(Adlam & Pennington, 2020b). There, a given learned function generally depends
on both the dataset $\mathcal{D}$ and initialization seed $\theta_{0}$. We
write this as $f_{\mathcal{D},\mathcal{\theta}_{0}}$. The role of random
initialization is played by the $\bm{A}$ matrix in our setting. For a given
function, its variance over datasets and its variance over initializations are
respectively given by

$\displaystyle\mathrm{Var}_{\mathcal{D}}f$
$\displaystyle\equiv\mathbb{E}_{\mathcal{D}}(f_{\mathcal{D},\mathcal{\theta}_{0}}-\mathbb{E}_{\mathcal{D}}[f_{\mathcal{D},\mathcal{\theta}_{0}}])^{2}$
(91) $\displaystyle\mathrm{Var}_{\theta_{0}}f$
$\displaystyle\equiv\mathbb{E}_{\theta_{0}}(f_{\mathcal{D},\mathcal{\theta}_{0}}-\mathbb{E}_{\theta_{0}}[f_{\mathcal{D},\mathcal{\theta}_{0}}])^{2}$
(92)

Here $\mathbb{E}_{\mathcal{D}}[f_{\mathcal{D},\mathcal{\theta}_{0}}]$ and
$\mathbb{E}_{\theta_{0}}[f_{\mathcal{D},\mathcal{\theta}_{0}}]$ can be viewed
as infinitely bagged or infinitely ensembled predictors respectively. The bias
of a function over datasets or initializations is given by the test error of
$\mathbb{E}_{\mathcal{D}}[f_{\mathcal{D},\mathcal{\theta}_{0}}]$,
$\mathbb{E}_{\theta_{0}}[f_{\mathcal{D},\mathcal{\theta}_{0}}]$ respectively.
The irreducible bias is given by
$\mathbb{E}_{\mathcal{D},\theta_{0}}[f_{\mathcal{D},\mathcal{\theta}_{0}}]$.

###  G.2 Derivation

In this section, we consider the effect of ensembling over $E$ random initial
conditions and bagging over $B$ random datasets. We let $\bm{v}^{0}_{e,b}(t)$
represent the weight discrepancy for model $e$ on dataset $b$. Here $e$ runs
from $1$ to $E$ and $b$ runs from $1$ to $B$. The $(e,b)$th vector has
dynamics:

$\displaystyle\frac{d}{dt}\bm{v}^{0}_{e,b}(t)=-\left(\frac{1}{N}\bm{A}_{e}^{\top}\bm{A}_{e}\right)\left(\frac{1}{N}\bm{\Psi}_{b}^{\top}\bm{\Psi}_{b}\right)\bm{v}^{0}_{e,b}(t).$
(93)

Ensembling and bagging would correspond to averaging these $\bm{v}^{0}$s over
these $EB$ systems

$\displaystyle\bar{\bm{v}}^{0}(t)=\frac{1}{EB}\sum_{e=1}^{E}\sum_{b=1}^{B}\bm{v}^{0}_{e,b}(t).$
(94)

The key vectors to track for this computation are

$\displaystyle\bm{v}^{1}_{e,b}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{\Psi}_{b}\bm{v}^{0}_{e,b}(t)+\sigma\bm{\epsilon}_{b}\
,\
\bm{v}^{2}_{e,b}(t)=\frac{1}{\alpha\sqrt{M}}\bm{\Psi}_{b}^{\top}\bm{v}^{1}_{e,b}(t)$
$\displaystyle\bm{v}^{3}_{e,b}(t)$
$\displaystyle=\frac{1}{\sqrt{M}}\bm{A}_{e}\bm{v}^{2}_{e,b}(t)\ ,\
\bm{v}^{4}_{e,b}(t)=\frac{1}{\nu\sqrt{M}}\bm{A}_{e}^{\top}\bm{v}^{3}_{e,b}(t).$
(95)

We can further show that the $\bm{v}^{0}_{e,b}$ and
$\bm{v}^{0}_{e^{\prime},b^{\prime}}$ have response functions that decouple
across $e,b$. Intuitively, giving the dynamical system $e,b$ a kick should not
alter the trajectory of the separate $e^{\prime},b^{\prime}$ dynamical system,
even if they share disorder $\\{\bm{\Psi},\bm{A}\\}$. The DMFT description of
the proportional limit yields the following integral equations for the $v$
fields:

$\displaystyle\partial_{t}v_{e,b,k}^{0}(t)=-v^{4}_{e,b,k}(t)$ $\displaystyle
v^{1}_{e,b}(t)=u^{1}_{e,b}(t)+\frac{1}{\alpha}\int
dsR_{0,2}(t,s)v^{1}_{e,b}(s)$ $\displaystyle
v^{2}_{e,b,k}(t)=u^{2}_{e,b,k}(t)+\lambda_{k}\int
dsR_{1}(t,s)v^{0}_{e,b,k}(s)$ $\displaystyle
v^{3}_{e,b}(t)=u^{3}_{e,b}(t)+\frac{1}{\nu}\int dsR_{2,4}(t,s)v^{1}_{e,b}(s)$
$\displaystyle v^{4}_{e,b,k}(t)=u^{4}_{e,b,k}(t)+\int
dsR_{3}(t,s)v^{2}_{e,b,k}(s).$ (96)

Here, the response functions $R$ are to be computed within a single system. In
what follows, we will use $\langle\cdot\rangle$ to denote averages over the
disorder, and explicitly write out any averages over the ensemble members and
datasets.

The Gaussian variables in the DMFT have the following covariance

$\displaystyle\left<u^{1}_{e,b}(t)u^{1}_{e^{\prime},b^{\prime}}(s)\right>=\delta_{b,b^{\prime}}C^{0}_{e,e^{\prime}}(t,s)$
$\displaystyle\left<u^{2}_{e,b,k}(t)u^{2}_{e^{\prime},b^{\prime},k}(s)\right>=\delta_{b,b^{\prime}}\frac{\lambda_{k}}{\alpha}C_{1,e,e^{\prime}}(t,s)$
$\displaystyle\left<u^{3}_{e,b}(t)u^{3}_{e^{\prime},b^{\prime}}(s)\right>=\delta_{e,e^{\prime}}C_{2,b,b^{\prime}}(t,s)$
$\displaystyle\left<u^{4}_{e,b,k}(t)u^{4}_{e^{\prime},b^{\prime},k}(s)\right>=\delta_{e,e^{\prime}}\frac{1}{\nu}C_{3,b,b^{\prime}}(t,s)$
(97)

The covariances above
$C_{0,e,e^{\prime}},C_{1,e,e^{\prime}},C_{2,b,b^{\prime}},C_{3,b,b^{\prime}}$
allow for different ensemble or dataset index but not both. We will use
$C_{0},C_{1},C_{2},C_{3}$ etc to represent the correlation functions within a
single system. For instance,
$C^{0}_{e,e^{\prime}}(t,s)=\frac{1}{M}\sum_{k}\lambda_{k}\left<v^{0}_{e,b}(t)v^{0}_{e^{\prime},b}(s)\right>$
while
$C^{0}=\frac{1}{M}\sum_{k}\lambda_{k}\left<v^{0}_{e,b}(t)v^{0}_{e,b}(s)\right>$.
The correlation function of interest is thus

$\displaystyle\mathcal{C}_{0,k,e,e^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle=\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[(w_{k}^{\star})^{2}+\frac{1}{\nu}\delta_{e,e^{\prime}}\mathcal{C}_{3}(\omega,\omega^{\prime})+\frac{\lambda_{k}}{\alpha}\mathcal{C}_{1,e,e^{\prime}}(\omega,\omega^{\prime})\right]$
$\displaystyle\mathcal{C}_{1,e,e^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle=\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})\mathcal{C}^{0}_{e,e^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle\mathcal{C}_{2,b,b^{\prime},k}(\omega,\omega^{\prime})$
$\displaystyle=(i\omega)(i\omega^{\prime})\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\frac{\lambda_{k}}{\alpha}\delta_{b,b^{\prime}}\mathcal{C}_{1}(\omega,\omega^{\prime})+\lambda_{k}^{2}\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})\left[(w_{k}^{\star})^{2}+\frac{1}{\nu}\mathcal{C}_{3,b,b^{\prime}}(\omega,\omega^{\prime})\right]$
$\displaystyle\mathcal{C}_{3,b,b^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle=\mathcal{R}_{3}(\omega)\mathcal{R}_{3}(\omega)\mathcal{C}_{2,b,b^{\prime}}(\omega,\omega^{\prime})$
(98)

We can combine the first two equations and the second two equations to
identify the structure of the cross-ensemble and cross-dataset (across-system)
correlations in terms of the marginal (within-system) correlation statistics

$\displaystyle\mathcal{C}_{0,e,e^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{1-\gamma_{0}(\omega,\omega^{\prime})}\frac{1}{M}\sum_{k}\lambda_{k}\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[(w_{k}^{\star})^{2}+\frac{1}{\nu}\delta_{e,e^{\prime}}\mathcal{C}_{3}(\omega,\omega^{\prime})\right]$
$\displaystyle\gamma_{0}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{\alpha
M}\sum_{k}\lambda_{k}^{2}\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})$
$\displaystyle\mathcal{C}_{2,b,b^{\prime}}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{1-\gamma_{2}(\omega,\omega^{\prime})}\frac{1}{M}\sum_{k}\lambda_{k}\left[\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})\lambda_{k}(w_{k}^{\star})^{2}+\frac{1}{\alpha}(i\omega)(i\omega^{\prime})\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\delta_{b,b^{\prime}}\mathcal{C}_{1}(\omega,\omega^{\prime})\right]$
$\displaystyle\gamma_{2}(\omega,\omega^{\prime})$
$\displaystyle=\frac{1}{\nu}\mathcal{R}_{3}(\omega)\mathcal{R}_{3}(\omega)\frac{1}{M}\sum_{k}\lambda_{k}^{2}$
(99)

These equations give the necessary cross-ensemble and cross-dataset
correlations. Now we can consider the effect of ensembling and bagging on the
dynamics. To do so, consider the Fourier transform of the bagged-ensembled
error $\bar{v}^{0}_{k}(t)=\frac{1}{EB}\sum_{eb}v^{0}_{k,e,b}(t)$, which has
the Fourier transform

$\displaystyle{\bar{v}}^{0}_{k}(\omega)=\mathcal{H}_{k}(\omega)\left[w^{\star}_{k}-\frac{1}{EB}\sum_{e,b}\left(u^{4}_{e,b,k}(\omega)+\mathcal{R}_{3}(\omega)u^{2}_{e,b,k}(\omega)\right)\right]$
(100)

Computing the correlation function for this bagged-ensembled field random
variable, we find

$\displaystyle\left<{\bar{v}}^{0}_{k}(\omega){\bar{v}}^{0}_{k}(\omega^{\prime})\right>$
$\displaystyle=\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[(w^{\star}_{k})^{2}+\frac{1}{\nu
E^{2}B^{2}}\sum_{e,e^{\prime},b,b^{\prime}}\delta_{e,e^{\prime}}\mathcal{C}_{3,b,b^{\prime}}(\omega,\omega^{\prime})+\frac{\lambda_{k}}{\alpha
E^{2}B^{2}}\sum_{ee^{\prime}bb^{\prime}}\delta_{b,b^{\prime}}\mathcal{C}_{1,e,e^{\prime}}(\omega,\omega^{\prime})\right]$
$\displaystyle=\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\left[(w^{\star}_{k})^{2}+\frac{1}{\nu
EB^{2}}\sum_{b,b^{\prime}}\mathcal{C}_{3,b,b^{\prime}}(\omega,\omega^{\prime})+\frac{\lambda_{k}}{\alpha
E^{2}B}\sum_{ee^{\prime}}\mathcal{C}_{1,e,e^{\prime}}(\omega,\omega^{\prime})\right]$
$\displaystyle=\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})(w^{\star}_{k})^{2}$
$\displaystyle+\frac{1}{\nu
E}\frac{\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\mathcal{R}_{3}(\omega)\mathcal{R}_{3}(\omega^{\prime})\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})}{1-\gamma_{2}(\omega,\omega^{\prime})}\
\left[\frac{1}{M}\sum_{\ell}\lambda_{\ell}^{2}(w_{\ell}^{\star})^{2}\right]$
$\displaystyle+\frac{1}{\nu\alpha
EB}\frac{\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\mathcal{C}_{1}(\omega,\omega^{\prime})(i\omega)(i\omega^{\prime})}{(1-\gamma_{2}(\omega,\omega^{\prime}))}\frac{1}{M}\sum_{\ell}\lambda_{\ell}\mathcal{H}_{\ell}(\omega)\mathcal{H}_{\ell}(\omega^{\prime})$
$\displaystyle+\frac{\lambda_{k}}{\alpha
B}\frac{\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\mathcal{R}_{1}(\omega)\mathcal{R}_{1}(\omega^{\prime})}{1-\gamma_{0}(\omega,\omega^{\prime})}\frac{1}{M}\sum_{\ell}\lambda_{\ell}\mathcal{H}_{\ell}(\omega)\mathcal{H}_{\ell}(\omega^{\prime})(w^{\star}_{\ell})^{2}$
$\displaystyle+\frac{\lambda_{k}}{\alpha\nu
EB}\frac{\mathcal{H}_{k}(\omega)\mathcal{H}_{k}(\omega^{\prime})\mathcal{C}_{3}(\omega,\omega^{\prime})}{1-\gamma_{0}(\omega,\omega^{\prime})}\frac{1}{M}\sum_{\ell}\lambda_{\ell}\mathcal{H}_{\ell}(\omega)\mathcal{H}_{\ell}(\omega^{\prime})$
(101)

The first term is the irreducible bias for mode $k$ which is the loss for mode
$k$ when the learned function is averaged over all possible datasets and all
possible projections. We see that the second term scales as $\frac{1}{\nu E}$
which will persist even if $B\alpha\to\infty$. Similarly, there is a term that
is order $\frac{1}{\alpha B}$ which will persist even if $\nu E\to\infty$.
Lastly, there are two terms which depend on both $B,E$. This is similar to the
variance that is explained by the interaction of the dataset and the random
projection (Adlam & Pennington, 2020b). The test loss is then a Fourier
transform of the above function

$\displaystyle\bar{\mathcal{L}}(t)=\frac{1}{M}\sum_{k}\lambda_{k}\left<\bar{v}^{0}_{k}(t)^{2}\right>.$
(102)

If $E,B\to\infty$, then we obtain the stated irreducible bias of the main
paper

$\displaystyle\lim_{E,B\to\infty}\bar{\mathcal{L}}(t)=\frac{1}{M}\sum_{k}\lambda_{k}(w^{\star}_{k})^{2}H_{k}(t)^{2}.$
(103)

This is the error of the mean output function over all possible datasets and
random projections of a certain size.

###  G.3 Ensembling is Not Always Compute Optimal

For a compute budget $C=NEt$, we find that ensembling does not provide as much
benefit as increasing the size of the model. From the results in the last
section, we note that ensembling reduces the variance. For this section, we
consider the $P\to\infty$ limit. We let $\mathcal{B}(N,t)$ represent the bias
and $\mathcal{V}(N,t)$ represent the variance within a single ensemble. The
loss at fixed compute then takes the form

$\displaystyle\mathcal{L}(\nu,C,t)=\mathcal{B}(\nu,t)+\frac{1}{\nu
E}\mathcal{V}(\nu,t).$ (104)

For any $\nu$ which satisfies the condition that

$\displaystyle\frac{\partial}{\partial\nu}\mathcal{B}(\nu,t)\leq 0\ ,\
\frac{\partial}{\partial\nu}\mathcal{V}(\nu,t)\leq 0$ (105)

we have that ensembling is strictly dominated by increasing $\nu$.

##  Appendix H White Bandlimited Model

Figure 7: The white bandlimited model ($\lambda_{k}=1$) with $\alpha=0.8$ and
varying model size $\nu$ with no explicit noise $\sigma=0$ exhibits double
descent at late time. Optimal early stopping, like optimal regularization,
recovers monotonic scaling with $\nu$.

To gain intuition for the model, we can first analyze the case where
$\lambda_{k}=1$, which has a simpler DMFT description since each of the $M$
features are statistically identical. We illustrate the dependence of the loss
on model size $\nu$ and training time $t$ for $\alpha<1$ in Figure 7. We note
that the loss can be non-monotonic in $\nu$ at late training times, but that
monotonicity is maintained for optimal early stopping, similar to results on
optimal regularization in linear models (Advani et al., 2020) and random
feature models (Mei & Montanari, 2022; Simon et al., 2023).

###  H.1 Derivation

In the case of all $\lambda_{k}=1$ we have the following definitions

$\displaystyle\mathcal{R}_{1}(\omega)$
$\displaystyle=1-\frac{1}{\alpha}\frac{\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}{i\omega+\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}$
(106) $\displaystyle\mathcal{R}_{3}(\omega)$
$\displaystyle=1-\frac{1}{\nu}\frac{\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}{i\omega+\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}$

Writing $\mathcal{R}_{1}=1-\frac{\nu}{\alpha}(\mathcal{R}_{3}-1)$ allows us to
solve for $\mathcal{R}_{3}$ exactly:

$\displaystyle\mathcal{R}_{3}\left(i\omega+\mathcal{R}_{3}(1+\frac{\nu}{\alpha}(\mathcal{R}_{3}-1))\right)=i\omega+\mathcal{R}_{3}(1+\frac{\nu}{\alpha}(\mathcal{R}_{3}-1))-\frac{1}{\nu}\left(\mathcal{R}_{3}(1+\frac{\nu}{\alpha}(\mathcal{R}_{3}-1))\right).$
(107)

This is a cubic equation that can be solved for $\mathcal{R}_{3}$ as a
function of $\omega$. In the limit of $\alpha\to\infty$ this simplifies to:

$\displaystyle\mathcal{R}_{3}i\omega+\mathcal{R}_{3}^{2}$
$\displaystyle=i\omega+\left(1-\frac{1}{\nu}\right)\mathcal{R}_{3}$ (108)
$\displaystyle\Rightarrow\mathcal{R}_{3}$
$\displaystyle=\frac{1}{2}[(1-\nu^{-1}-i\omega)+\sqrt{(1-\nu^{-1}-i\omega)^{2}+4i\omega}].$

###  H.2 Timescale Corrections in The Small $\nu$ Regime

By expanding the above in the limit of small $\nu$ we get that
$\mathcal{R}_{3}$ goes as

$\displaystyle\mathcal{R}_{3}\sim\frac{i\omega}{\nu^{-1}-1+i\omega}\ ,\ \nu\to
0$ (109)

From this approximate response function, we find that the transfer function
takes the form

$\displaystyle H(\tau)$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{e^{i\omega\tau}}{i\omega+\frac{i\omega}{\nu^{-1}-1+i\omega}}=\int\frac{d\omega}{2\pi}\frac{(\nu^{-1}-1+i\omega)e^{i\omega\tau}}{i\omega\left[\nu^{-1}+i\omega\right]}$
$\displaystyle=(1-\nu)+\nu e^{-\tau/\nu},$ (110)

where in the last line, we used the residue theorem. We note that in this
perturbative approximation that this transfer function is always greater than
the transfer function at $\nu\to\infty$ which is $e^{-\tau}$. Thus finite
$\nu$ leads to higher bias in this regime. We define bias and variance
precisely in Appendix G.1.

###  H.3 Timescale corrections in fully expressive regime $\nu>1$

For $\nu\gg 1$, we can approximate $R_{3}(\omega)\sim
1-\nu^{-1}(1+i\omega)^{-1}$, we have

$\displaystyle H(\tau)$
$\displaystyle\sim\int\frac{d\omega}{2\pi}\frac{e^{i\omega\tau}}{i\omega+1-\nu^{-1}(i\omega+1)^{-1}}=\int\frac{d\omega}{2\pi}\frac{e^{i\omega\tau}(1+i\omega)}{(i\omega+1-\nu^{-1/2})(i\omega+1+\nu^{-1/2})}$
$\displaystyle=\frac{1}{2}e^{-\tau(1+\sqrt{\nu})}+\frac{1}{2}e^{-\tau(1-\sqrt{\nu})}=e^{-\tau}\cosh\left(\tau/\sqrt{\nu}\right)$
(111)

where we used the residue theorem after closing the contour in the upper half-
plane. In Figure 8, we show that this perturbative approximation does capture
a slowdown in the dynamics for large but finite $\nu$.

Figure 8: Slower timescales in the $\nu>1$ regime for white bandlimited
features.

##  Appendix I Power-Law Bottleneck Scalings

In this section we calculate the scaling of the loss with the various limiting
resources (time, model size, and data) when using power law features. Since
the power-law features give a trace class kernel (i.e.
$\sum_{k=1}^{\infty}\lambda_{k}<\infty$), we use the non-proportional limit
formalism in Appendix F, which gives an expression for $\mathcal{L}(t,N,P)$
with $M$ already considered infinite. While the resulting expressions are not
a formal proportional thermodynamic limit and finite $N,P$ corrections exist
in the form of fluctuations from one random realization of the system to
another. These corrections decay rapidly enough at finite $N,P$ for this mean
field theory to be accurate and descriptive in realistic systems (Bordelon et
al., 2020; Simon et al., 2023; Cheng & Montanari, 2022). We plot this
variability of random finite size experiments as highlighted standard
deviations in the main text figures.

###  I.1 Time Bottleneck

The time bottleneck is defined as the limiting dynamics in the absence of any
model or data finite size effects. To eliminate those effects, we simply study
the $\alpha,\nu\to\infty$ limit

$\displaystyle\mathcal{L}_{\infty}(t)=\lim_{P,N\to\infty}\mathcal{L}(t,P,N).$
(112)

In this limit, the response functions simplify to
$\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)\to 1$ so that

$\displaystyle\mathcal{H}_{k}(\omega)=\frac{1}{i\omega+\lambda_{k}}\implies
H_{k}(\tau)=e^{-\lambda_{k}\tau}\Theta(\tau).$ (113)

Further, in this limit, we have that
$C_{0}(t,s)=\frac{1}{M}\sum_{k}\lambda_{k}H_{k}(t)H_{k}(s)(w^{\star}_{k})^{2}$
since all the variance terms (which depend on $\nu^{-1},\alpha^{-1}$) drop
out. Thus we have the following loss at time $t$,

$\displaystyle\mathcal{L}(t)=\sum_{k}\lambda_{k}(w_{k}^{\star})^{2}e^{-2\lambda_{k}t}\sim\int_{1}^{\infty}dkk^{-a}\exp\left(-2k^{-b}t\right)\sim
t^{-(a-1)/b}.$ (114)

where the final scaling with time can be obtained through either change of
variables or steepest descent methods (Bordelon & Pehlevan, 2022a).

###  I.2 Model Bottleneck

In this section we take $\alpha,t\to\infty$. This leaves us with the following
equation for $r\equiv\lim_{\omega\to 0}(i\omega)^{-1}\mathcal{R}_{3}(\omega)$.

$\displaystyle N$
$\displaystyle=\sum_{k}\frac{\lambda_{k}r}{\lambda_{k}r+1}\approx\int_{1}^{\infty}\frac{dk}{k^{b}/r+1}\approx
r^{1/b}\implies r\approx N^{b}.$ (115)

Now, the large time limit of the transfer functions $H_{k}(\tau)$ can be
obtained from the final-value theorem

$\displaystyle\lim_{t\to\infty}H_{k}(\tau)=\lim_{\omega\to
0}\frac{i\omega}{i\omega+\lambda_{k}ri\omega}=\frac{1}{1+\lambda_{k}r}\approx\Theta(k-N).$
(116)

Now, integrating over the eigenvalue density to get the total loss gives

$\displaystyle\mathcal{L}(t)\sim\int_{1}^{\infty}dk\ k^{-a}\
\Theta(k-N)=\int_{N}^{\infty}k^{-a}=N^{-(a-1)}.$ (117)

Thus we expect a powerlaw scaling of the form $\mathcal{L}\sim N^{-(a-1)}$ in
this regime.

###  I.3 Data Bottleneck

In this section we take $\nu,t\to\infty$. This leaves us with the following
equation for $r\equiv\lim_{\omega\to 0}(i\omega)^{-1}\mathcal{R}_{1}(\omega)$.

$\displaystyle P$
$\displaystyle=\sum_{k}\frac{\lambda_{k}r}{\lambda_{k}r+1}\approx\int_{1}^{\infty}\frac{dk}{k^{b}/r+1}\approx
r^{1/b}\implies r\approx P^{b}.$ (118)

Now, the large time limit of the transfer functions $H_{k}(\tau)$ can again be
obtained from the final-value theorem

$\displaystyle\lim_{t\to\infty}H_{k}(\tau)=\lim_{\omega\to
0}\frac{i\omega}{i\omega+\lambda_{k}ri\omega}=\frac{1}{1+\lambda_{k}r}\approx\Theta(k-P).$
(119)

Now, integrating over the eigenvalue density to get the total loss gives

$\displaystyle\mathcal{L}(t)\sim\int_{1}^{\infty}dk\ k^{-a}\
\Theta(k-P)=\int_{P}^{\infty}k^{-a}=P^{-(a-1)}.$ (120)

The loss will therefore scale as $P^{-(a-1)}$ in this data-bottleneck regime.

##  Appendix J Optimization Extensions

###  J.1 Discrete Time

In this section, we point out that DMFT can also completely describe discrete
time training as well. In this section we consider discrete time gradient
descent with learning rate $\eta$

$\displaystyle\bm{v}^{0}(t+1)=\bm{v}^{0}(t)-\eta\bm{v}^{4}(t)$
$\displaystyle\bm{v}^{4}(t)=\frac{1}{\nu\sqrt{M}}\bm{A}^{\top}\bm{v}^{3}(t)\
,\ \bm{v}^{3}(t)=\frac{1}{\sqrt{M}}\bm{A}\bm{v}^{2}(t)$
$\displaystyle\bm{v}^{2}(t)=\frac{1}{\alpha\sqrt{M}}\bm{\Psi}^{\top}\bm{v}^{1}(t)\
,\ \bm{v}^{1}(t)=\frac{1}{\sqrt{M}}\bm{\Psi}\bm{v}^{0}(t)$ (121)

Following either the MSR or cavity derivation, we obtain an analgous set of
limiting DMFT equations defined for integer times $t\in\mathbb{Z}$,

$\displaystyle v^{0}_{k}(t+1)=v^{0}_{k}(t)-\eta
v^{4}_{k}(t)+\delta(t+1)w^{\star}_{k}$ $\displaystyle
v^{1}(t)=u^{1}(t)+\alpha^{-1}\sum_{s}R_{0,2}(t,s)v^{1}(s)$ $\displaystyle
v^{2}_{k}(t)=u^{2}_{k}(t)+\sum_{s}R_{1}(t,s)v^{0}_{k}(s)$ $\displaystyle
v^{3}(t)=u^{3}(t)+\sum_{s}R_{2,4}(t,s)v^{3}(s)$ $\displaystyle
v^{4}_{k}(t)=u^{4}_{k}(t)+\sum_{s}R_{3}(t,s)v^{2}_{k}(s)$ (122)

The delta function in this context is defined as

$\displaystyle\delta(t+1)=\begin{cases}1&t=-1\\\ 0&\text{else}\end{cases}$
(123)

ensures that the initial condition $v_{k}^{0}(0)=w^{\star}_{k}$ is satisfied.
These iteration equations can be closed for the response functions and
correlation functions and solved over $T\times T$ matrices.

Alternatively, we can also solve this problem in an analogous frequency space.
Analogous to the Fourier transform method, the equations in discrete time can
be closed in terms of the $Z$-transform

$\displaystyle v(z)=\sum_{t=-\infty}^{\infty}z^{-t}v(t)$ (124)

Applying this transform gives us the following expression for the $v^{0}_{k}$
fields.

$\displaystyle v^{0}_{k}(z)=\frac{zw^{\star}_{k}-\eta
u^{4}_{k}(z)-\eta\mathcal{R}_{3}(z)u^{2}_{k}(z)}{z-1+\eta\lambda_{k}\mathcal{R}_{1}(z)\mathcal{R}_{3}(z)}\equiv\mathcal{H}_{k}(z)\left[zw^{\star}_{k}-\eta
u^{4}_{k}(z)-\eta\mathcal{R}_{3}(z)u^{2}_{k}(z)\right]$ (125)

Similar to the Fourier case, the final losses can be extracted as the $z\to 1$
limit of these objects.

###  J.2 Momentum

As mentioned in appendix B, it is straightforward to extend the DMFT treatment
beyond just gradient descent dynamics to include a momentum term with momentum
$\beta$.

We first consider this replacement in continuous time. This requires applying
the following replacement:

$\partial_{t}v^{0}_{k}(t)=-v^{4}_{k}(t)\to(\beta\partial_{t}^{2}+\partial_{t})v^{0}_{k}(t)=-v^{4}_{k}(t).$
(126)

This slightly modifies the expressions for the response functions. For
example, in Fourier space the response functions become:

$\displaystyle\mathcal{R}_{0,2}(\omega)=-\frac{1}{M}\sum_{k}\frac{\lambda_{k}}{\epsilon+i\omega+\beta(i\omega)^{2}+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}\mathcal{R}_{3}(\omega)$
(127)
$\displaystyle\mathcal{R}_{2,4}(\omega)=-\frac{1}{M}\sum_{k}\frac{\lambda_{k}}{\epsilon+i\omega+\beta(i\omega)^{2}+\lambda_{k}\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)}\mathcal{R}_{1}(\omega).$

In discrete time, momentum updates can be expressed as

$\displaystyle\bm{v}^{0}(t+1)=\bm{v}^{0}(t)-\eta\bm{b}(t)$
$\displaystyle\bm{b}(t)=\bm{v}^{4}(t)+\mu\bm{b}(t-1)$ (128)

where $\bm{b}(t)$ is the filtered version of the loss gradient (the
$\bm{v}^{4}(t)$ field) with momentum coefficient $\mu$ and $\eta$ is the
learning rate. The dependence on the $\bm{b}(t)$ field can be eliminated by
turning this into a second order difference equation

$\displaystyle\bm{v}^{0}(t+1)-\bm{v}^{0}(t)-\mu\left(\bm{v}^{0}(t)-\bm{v}^{0}(t-1)\right)=-\eta\bm{v}^{4}(t).$
(129)

Again, the final result can be expressed in terms of the $Z$-transformed
transfer functions $\mathcal{H}_{k}(z)$ which have the form

$\displaystyle\mathcal{H}_{k}(z)=\frac{1}{z-1-\mu+\mu
z^{-1}+\eta\mathcal{R}_{1}(z)\mathcal{R}_{3}(z)}.$ (130)

###  J.3 One Pass SGD

In this section we derive online SGD with projected features. At each step a
random batch of $P=\alpha M$ samples are collected (independent of previous
samples), giving a matrix $\bm{\Psi}(t)$ of sampled features. The update at
step $t$ is

$\displaystyle\bm{v}^{0}(t+1)=\bm{v}^{0}(t)+\eta\left(\frac{1}{N}\bm{A}^{\top}\bm{A}\right)\left(\frac{1}{P}\bm{\Psi}(t)^{\top}\bm{\Psi}(t)\right)\bm{v}^{0}(t).$
(131)

The DMFT limit gives the following statistical description of the fields,
which decouple over time for the $v^{1}(t),v^{2}_{k}(t)$ but remain coupled
across time for $v^{3}(t),v^{4}_{k}(t)$

$\displaystyle v^{1}(t)$ $\displaystyle=u^{1}(t)\ ,\
u^{1}(t)\sim\mathcal{N}(0,C_{0}(t,t)\delta(t-s)),$ (132) $\displaystyle
v^{2}_{k}(t)$ $\displaystyle=u^{2}_{k}(t)+\lambda_{k}v^{0}_{k}(t),$
$\displaystyle v^{3}(t)$
$\displaystyle=u^{3}(t)+\frac{1}{\nu}\sum_{s}R_{2,4}(t,s)v^{3}(s),$
$\displaystyle v^{4}_{k}(t)$
$\displaystyle=u^{4}_{k}(t)+\sum_{s}R_{3}(t,s)v^{2}_{k}(s),$ $\displaystyle
v^{0}(t+1)$ $\displaystyle=v^{0}_{k}(t)-\eta v^{4}_{k}(t).$

This system cannot exhibit overfitting effects as we have the statistical
equivalence between the covariance of $v^{1}$ and the test loss:

$\displaystyle\hat{\mathcal{L}}(t)=\left<v^{1}(t)^{2}\right>=\left<u^{1}(t)^{2}\right>=C_{0}(t,t)=\mathcal{L}(t)$
(133)

We note that this is very different than the case where data is reused at
every step, which led to a growing gap between train and test loss as we
derive in Appendix E.

##  Appendix K Kernel Analysis of Feature Learning Networks

(a) Final NTK Spectra for Varying Widths (b) Final NTK Task-Power Decay (c)
Losses for Feature Learning Networks

Figure 9: a) The observed power law spectrum on a held out test set of the
after-kernel for a width $N=128$ ResNet trained on CIFAR-5m. Early on in
training, the spectrum flattens quite rapidly. At later times, the spectral
decay remains relatively constant. b) The fraction of the task unexplained, as
defined in Equation 19. Throughout training, the top eigenmode of the after-
kernel explains more and more of the task. c) The test loss of the network. We
see that the observed scaling of this quantity is faster than that predicted
from analyzing the after-kernel. d) The kernel-target alignment of the after
kernel improves throughout training time. The error bars here denote different
ensemble members. Their relatively small size implies that the kernel
trajectory is relatively deterministic over different initialization seeds. e)
The norm of the after-kernel throughout training is relatively constant for
this task.

In Section 5.1, we observed that feature learning networks can achieve better
loss and compute-optimal scaling. In such settings, it may be useful to
observe the after kernel, namely the NTK at the end of training. This object
can often shed insight into the structure of the learned network function
(Atanasov et al., 2022; Long, 2021) and its generalization. In some cases, it
has been observed that the final kernel stabilizes during the course of
training (Fort et al., 2020), potentially allowing one to potentially deduce
scaling laws from the spectrum and task-model alignment of this after-kernel,
though other papers have observed contrary results (Vyas et al., 2022).

Motivated by this, we study the NTKs of the finite-width networks trained for
64 epochs with the animate-inanimate CIFAR-5m discrimination task. We observe
in Figure 9 a) that the spectrum becomes flatter, with a decay exponent of
close to $1.4$ down from $2.0$ for the initial kernel.

The fraction of the task power unexplained is also observed to have a lower
exponent in Figure 9 b), however there is also the presence of a low rank
spike indicative of the kernel aligning to this discrimination tasks.

From these scalings we can obtain the $a$ and $b$ exponents and get a
prediction for the scaling of the test loss. We plot this in grey in Figure 9
c). The observed scaling (in black) is much better than that predicted by the
after-kernel. This is an indication the the after kernel continues evolving in
this task, improving the scaling exponent of the test loss.

The kernel-target alignment (Cortes et al., 2012), as measured by

$A=\frac{\bm{y}^{\top}\bm{K}\bm{y}}{\bm{y}^{\top}\bm{y}|\bm{K}|_{F}},$ (134)

is plotted in 9 d). Here $\bm{y}$ is the target labels on a held-out test set,
and $\bm{K}$ is the gram matrix of the after-kernel on this test set. We
indeed observe a consistent increase in this quantity across time. This gives
an indication that understanding the evolution of the after-kernel will be
useful

##  Appendix L Numerical Recipes

###  L.1 Iteration of DMFT Equations on Time $\times$ Time matrices

The simplest way to solve the DMFT equations is to iterate them from a
reasonable initial condition (Mignacco et al., 2020; Bordelon & Pehlevan,
2022b). We solve in discrete time for $T\times T$ matrices
$\\{\bm{R}_{0,2},\bm{R}_{1},\bm{R}_{2,4},\bm{R}_{3},\bm{C}_{0},\bm{C}_{1},\bm{C}_{2},\bm{C}_{3}\\}$
which have entries $[\bm{R}]_{t,s}=R(t,s),[\bm{C}]_{t,s}=C(t,s)$, etc. We let
$\bm{\Theta}(t,s)=\eta\Theta(t-s)$ where $\eta$ is the learning rate.

  1. 1.

Solve for the response functions by updating the closed equations as matrices
by iterating the equations.

$\displaystyle\bm{R}_{0,2,k}\leftarrow-\left[\bm{\Theta}^{-1}+\lambda_{k}\bm{R}_{3}\bm{R}_{1}\right]^{-1}\bm{R}_{3},$
(135)
$\displaystyle\bm{R}_{0,2}\leftarrow\frac{1}{M}\sum_{k}\lambda_{k}\bm{R}_{0,2,k},$
$\displaystyle\bm{R}_{1}\leftarrow\left[\bm{I}-\alpha^{-1}\bm{R}_{0,2}\right]^{-1},$
$\displaystyle\bm{R}_{2,4,k}\leftarrow-\lambda_{k}\left[\bm{I}+\lambda_{k}\bm{R}_{1}\bm{\Theta}\bm{R}_{3}\right]^{-1}\bm{R}_{1}\bm{\Theta}\
,\ \bm{R}_{2,4}=\frac{1}{M}\sum_{k}\bm{R}_{2,4,k},$
$\displaystyle\bm{R}_{3}\leftarrow\left[\bm{I}-\nu^{-1}\bm{R}_{2,4}\right]^{-1}.$

  2. 2.

Once these response functions have converged, we can iterate the equations for
the correlation functions

$\displaystyle\bm{C}_{0,k}\leftarrow\left[\bm{I}+\lambda_{k}\bm{\Theta}\bm{R}_{3}\bm{R}_{1}\right]^{-1}\left[(w^{\star}_{k})^{2}\bm{1}\bm{1}^{\top}+\bm{\Theta}\left(\nu^{-1}\bm{C}_{3}+\frac{\lambda_{k}}{\alpha}\bm{R}_{3}\bm{C}_{1}\bm{R}_{3}^{\top}\right)\bm{\Theta}^{\top}\right]\left[\bm{I}+\lambda_{k}\bm{\Theta}\bm{R}_{3}\bm{R}_{1}\right]^{-1\top},$
(136)
$\displaystyle\bm{C}_{0}\leftarrow\frac{1}{M}\sum_{k=1}^{M}\lambda_{k}\bm{C}_{0,k},$
$\displaystyle\bm{C}_{1}\leftarrow\bm{R}_{1}\bm{C}_{0}\bm{R}_{1}^{\top},$
$\displaystyle\bm{C}_{2,k}\leftarrow\left[\bm{I}+\lambda_{k}\bm{R}_{1}\bm{\Theta}\bm{R}_{3}\right]^{-1}\left(\frac{\lambda_{k}}{\alpha}\bm{C}_{1}+\bm{R}_{1}\left[(w^{\star}_{k})^{2}\lambda_{k}^{2}\bm{1}\bm{1}^{\top}+\frac{\lambda_{k}^{2}}{\nu}\bm{\Theta}\bm{C}_{3}\bm{\Theta}^{\top}\right]\bm{R}_{1}^{\top}\right)\left[\bm{I}+\lambda_{k}\bm{R}_{1}\bm{\Theta}\bm{R}_{3}\right]^{-1\top},$
$\displaystyle\bm{C}_{2}\leftarrow\frac{1}{M}\sum_{k}\bm{C}_{2,k}.$

After iterating these equations, one has the discrete time solution to the
DMFT order parameters and any other observable can then be calculated.

###  L.2 Fourier Transform Method

To accurately compute the Fourier transforms in the model/data bottleneck
regime ($\alpha<1$ or $\nu<1$) we have that
$\mathcal{R}_{1}(\omega)\mathcal{R}_{3}(\omega)\sim i\omega r$ as $\omega\to
0$ so we must resort to analyzing the principal part and the delta-function
contribution to the integral. Construct a shifted and non-divergent version of
the function $\mathcal{H}(\omega)$.

$\displaystyle\mathcal{H}(\omega)$
$\displaystyle=\mathcal{\tilde{H}}(\omega)+\frac{1}{\epsilon+i\omega(1+r)}$