# Depthwise Hyperparameter Transfer in Residual Networks: Dynamics and Scaling
Limit

Blake Bordelon∗ ¶ , Lorenzo Noci ♯ , Mufan (Bill) Li § ,  Boris Hanin† § &
Cengiz Pehlevan ¶
¶ Harvard University
♯ ETH Zürich
§ Princeton University
Equal first authors.Equal senior authors.

###### Abstract

The cost of hyperparameter tuning in deep learning has been rising with model
sizes, prompting practitioners to find new tuning methods using a proxy of
smaller networks. One such proposal uses $\mu$P parameterized networks, where
the optimal hyperparameters for small width networks _transfer_ to networks
with arbitrarily large width. However, in this scheme, hyperparameters do not
transfer across depths. As a remedy, we study residual networks with a
residual branch scale of $1/\sqrt{\text{depth}}$ in combination with the
$\mu$P parameterization. We provide experiments demonstrating that residual
architectures including convolutional ResNets and Vision Transformers trained
with this parameterization exhibit transfer of optimal hyperparameters across
width and depth on CIFAR-10 and ImageNet. Furthermore, our empirical findings
are supported and motivated by theory. Using recent developments in the
dynamical mean field theory (DMFT) description of neural network learning
dynamics, we show that this parameterization of ResNets admits a well-defined
feature learning joint infinite-width and infinite-depth limit and show
convergence of finite-size network dynamics towards this limit.

##  1 Introduction

(a) ResNet, $\mu$P (b) This work: $\frac{1}{\sqrt{L}}$-ResNet + $\mu$P

Figure 1: Learning rates transfer across both depth and width in our proposed
parameterization but not in $\mu$P or standard parameterization (Fig. B.2).
Loss is plotted after $20$ epochs on CIFAR-10. All the missing datapoints
indicate that the corresponding run diverged. In particular, all depth $30$
runs diverged when using $\mu$P without $1/\sqrt{\text{depth}}$ scaling.

Increasing the number of parameters in a neural network has led to consistent
and often dramatic improvements in model quality (Kaplan et al., 2020;
Hoffmann et al., 2022; Zhai et al., 2022; Klug & Heckel, 2023; OpenAI, 2023).
To realize these gains, however, it is typically necessary to conduct by
trial-and-error a grid search for optimal choices of hyperparameters, such as
learning rates. Doing so directly in SOTA models, which may have hundreds of
billions (Brown et al., 2020) or even trillions of parameters (Fedus et al.,
2022), often incurs prohibitive computational costs.

To combat this, an influential recent line of work by Yang & Hu (2021)
proposed the so-called $\mu$P parameterization, which seeks to develop
principles by which optimal hyperparameters from small networks can be reused
— or _transferred_ — to larger networks (Yang et al., 2021). The $\mu$P
prescription focuses on transfer from narrower to wider models, but does not
always transfer across depth (see Figure 1(a)). This leads us to pose the
following problem:

Question: Can we transfer hyperparameters simultaneously across depth and
width?

In this work, we provide an affirmative answer for a particular flexible class
of residual architectures (including ResNets and transformers) with residual
branches scaled like $1/\sqrt{\text{depth}}$ (see Equation 1 and Table 1 for
the exact parameterization). More specifically, our empirical contributions
can be summarized as follows:

  * •

We show that a simple variant of the $\mu$P parameterization allows for
transfer of learning rates simultaneously across depth and width in residual
networks with $1/\sqrt{\text{depth}}$ scaling on the residual branches (see
Equation 1 and Figure 1(b)). Moreover, we observe transfer not only for
training with a fixed learning rate, but also with complex learning rate
schedules (Figure 3(a-d)).

  * •

We show that the lack of learning rate transfer across depth in $\mu$P ResNets
can be partly overcome by adding normalization layers (Figure 2). However, we
find on CIFAR-10, Tiny ImageNet, ImageNet, that the transfer can be
significantly improved with our $1/\sqrt{\text{depth}}$ residual branches
(Figures 3 and B.5).

  * •

While most of our experiments concern residual convolutional networks, such as
wide ResNets, we also present evidence that our prescription leads to learning
rate transfer in Vision Transformers (ViTs), both with and _without_
normalization (Figure 3).

  * •

We demonstrate transfer for learning rates, as well as other hyperparameters
such as momentum coefficients and regularization strengthsin convolutional
residual networks (Figure 4).

(a) ResNet + $\mu$P, BatchNorm (b) $\frac{1}{\sqrt{L}}$ ResNet + $\mu$P,
BatchNorm

Figure 2: Effect of normalization layers. The above examples are taken after
$20$ epochs on CIFAR-10. Normalization layers can slightly improve consistency
across depths in standard ResNets, but consistency is much more reliable with
the $1/\sqrt{L}$ scaling.

The reason hyperparameters can transfer is due to the existence of a scaling
limit where the network predictions and internal features of the model move at
a rate independent of the size of the network. The key underlying idea is that
two neural networks will exhibit similar training dynamics (and hence similar
optimal hyperparameters) if they are both finite size approximations to, or
discretizations of, such a consistent scaling limit. In fact, the success of
the mean field or $\mu$P scaling (Yang & Hu, 2021; Bordelon & Pehlevan, 2022b)
is precisely due to the width-independent scale of feature updates. In the
same spirit, we propose our parameterization so that the rate of feature
updates and prediction updates are also independent of the depth of the
network. Our main theoretical contributions can be summarized as follows:

  * •

Extending the infinite-width limit results of Bordelon & Pehlevan (2022b), we
show the infinite-width-and-depth limit of a residual network with
$1/\sqrt{\text{depth}}$ residual branch scaling and $\mu$P initialization can
be characterized by a dynamical mean field theory (DMFT) (Section 4, 4.2)
where the layer index takes on a continuous “layer time.” In general, the
network’s internal features evolve during training in this limit. We show that
a lazy limit of the DMFT dynamics, where features are frozen, can also be
analyzed as a special case. See Section 4.3.1 and Proposition 2.

  * •

We provide exact solutions of the DMFT dynamics in the rich (i.e. non-kernel)
regime for a simple deep linear network. See Section 4.3.2.

The rest of the article is organized as follows. In Section 2, we define our
parameterization and scaling rules. In Section 3 demonstrate our empirical
results on hyperparameter transfer. Then in Section 4, we proceed to
characterize the limiting process at infinite width and depth and the
approximation error at finite $N,L$. Next we provide a brief overview of
related works in Section 5, before concluding with a discussion of limitations
and future directions in Section 6.

##  2 Framework for Depthwise Hyperparameter Transfer

We consider the following type of residual network (and suitable
architectural/optimization generalizations, see Appendix J) of width $N$ and
depth $L$ with inputs $\bm{x}\in\mathbb{R}^{D}$ which are mapped to
$N$-dimensional preactivations $\bm{h}^{\ell}\in\mathbb{R}^{N}$ and outputs
$f(\bm{x})\in\mathbb{R}$:

$\displaystyle
f(\bm{x})=\frac{\beta_{L}}{\gamma}\bm{w}^{L}\cdot\phi(\bm{h}^{L}(\bm{x}))\
,\quad\bm{h}^{\ell+1}(\bm{x})$
$\displaystyle=\bm{h}^{\ell}(\bm{x})+\beta_{\ell}\bm{W}^{\ell}\phi(\bm{h}^{\ell}(\bm{x}))\
,\quad\bm{h}^{1}(\bm{x})=\beta_{0}\bm{W}^{0}\bm{x},$ (1)

where $\gamma$ is a scaling factor, $\phi(\cdot)$ the activation function, and
the weights are initialized as
$\bm{W}^{\ell}_{ij}\sim\mathcal{N}(0,\sigma_{\ell}^{2})$ with a corresponding
learning rate $\eta(t)$. See Table 1.

| SP (PyTorch) |  $\mu$P / Mean Field |  $\mu$P$+1/\sqrt{L}$-Residual (Ours)
---|---|---|---
Branch Scale $\beta_{\ell}$ | $1$ | $\begin{cases}N^{-1/2},&\ell>0\\\ D^{-1/2},&\ell=0\end{cases}$ | $\begin{cases}N^{-1/2},&\ell=L\\\ (LN)^{-1/2},&0<\ell<L\\\ D^{-1/2},&\ell=0\end{cases}$
Output Scale $\gamma$ | $1$ | $\gamma_{0}N^{1/2}$ | $\gamma_{0}N^{1/2}$
LR Schedule $\eta(t)$ | $\eta_{0}(t)$ | $\eta_{0}(t)\gamma_{0}^{2}N$ | $\eta_{0}(t)\gamma_{0}^{2}N$
Weight Variance $\sigma_{\ell}^{2}$ | $\begin{cases}N^{-1},&\ell>0\\\ D^{-1},&\ell=0\end{cases}$ | $1$ | $1$
Table 1: In the framework of Equation 1, different choices of coefficients
leads to standard parameterization (SP, PyTorch default (Paszke et al.,
2019)), $\mu$P / Mean Field parameterization, and our
$\mu$P$+1/\sqrt{L}$-residual branch scaling which transfers over width and
depth.

Compared to what is commonly done in practical residual architectures, such as
ResNets and Transformers, our parameterization uses the $\sqrt{L}$ factor
proposed in Hayou et al. (2021); Fischer et al. (2023), and the coefficient
$\gamma$ from $\mu$P, but do not use normalization layers (Table 1). The
parameters $\bm{\theta}=\\{\bm{W}^{0},...,\bm{w}^{L}\\}$ are updated with a
gradient based learning rule. For concreteness consider SGD which we use in
many experiments (we also allow for gradient flow and other optimizers in
Appendix J):

$\displaystyle\bm{\theta}(t+1)=\bm{\theta}(t)-\eta(t)\
\mathbb{E}_{\bm{x}\sim\mathfrak{B}_{t}}\nabla_{\bm{\theta}}\
\mathcal{L}[f(\bm{\theta},\bm{x})],$ (2)

where $\mathfrak{B}_{t}$ is the minibatch at iteration $t$, $\eta(t)$ is the
learning rate schedule, and the loss function $\mathcal{L}$ depends on
parameters through the predictions $f$ of the model.

The choice of how $(\gamma,\eta(t))$ should be scaled as width $N$ and depth
$L$ go to infinity determines the stability and feature learning properties of
the network. If $\gamma,\eta\sim\mathcal{O}(1)$111We use $\mathcal{O}(1)$ to
denote a constant independent of width $N$ and depth $L$. (_kernel regime_),
then the network evolution is stable under gradient descent, but the network
does not learn features at infinite width (Jacot et al., 2018; Lee et al.,
2019; Yang & Hu, 2021). In this work, we study the richer feature learning
(_mean field_) parameterization $\gamma=\gamma_{0}\sqrt{N}$,
$\eta(t)=\eta_{0}(t)\gamma_{0}^{2}\ N$, where $\gamma_{0},\eta_{0}(t)$ are
independent of $N$. In this regime, the scaling factor of $1/\gamma$ in the
last layer enables feature learning at infinite width, and the constant
$\gamma_{0}$ controls the rate of feature learning. At fixed depth $L$, this
mean field parameterization is capable of generating consistent feature
learning dynamics across network widths $N$ (Yang & Hu, 2021; Bordelon &
Pehlevan, 2022b).

One contribution of this work is identifying how $\eta_{0},\gamma_{0}$ should
be scaled with depth $L$ in order to obtain a stable feature-learning large
$L$ limit. We argue that after the $\frac{1}{\sqrt{L}}$ factors have been
introduced the correct scaling is simply

$\displaystyle\eta=\eta_{0}\gamma_{0}^{2}N\ ,\ \gamma=\gamma_{0}\sqrt{N}\ ,\
\eta_{0},\gamma_{0}\sim\mathcal{O}(1).$ (3)

We show both theoretically and experimentally that this parameterization
exhibits consistent learning dynamics and hyperparameter optima across large
network depths $L$. To justify this choice, we will characterize the large
width and depth $N,L\to\infty$ limit of this class of models.

##  3 Experimental Verification of Hyperparameter Transfer

We report here a variety of experiments on hyperparameter transfer in residual
networks. We begin with a basic modification of equation 1 by replacing fully
connected layers with convolutional layers. Apart from the proposed
$1/\sqrt{L}$ scaling, notice that this convolutional residual model is similar
to widely adopted implementations of ResNets (He et al., 2016), with the main
difference being the absence of normalization layers, and a single convolution
operator per residual block (details in Appendix A). We use this residual
convolutional architecture in figures 1, 4, B.2 and B.1. We also experiment
with a practical ResNet (in figures 2, B.4, B.5), which has two convolutional
layers per block, and normalization layers (see Appendix A.2). Finally, we
also perform selected experiments on Vision Transformers (Dosovitskiy et al.,
2020), to demonstrate the broad applicability of our theoretical insights
(Figure 3). Below we list the the main takeaways from our experiments.

(a) $\frac{1}{\sqrt{L}}$ ViT + $\mu$P (b) $\frac{1}{\sqrt{L}}$ ViT + $\mu$P,
LayerNorm (c) ViT + $\mu$P, LayerNorm (d) $\frac{1}{\sqrt{L}}$ ViT + $\mu$P,
Cosine Annealing (e) $\frac{1}{\sqrt{L}}$ ViT ImageNet Dynamics (f)
$\frac{1}{\sqrt{L}}$ ViT + $\mu$P, ImageNet

Figure 3: ViTs trained with Adam also exhibit learning rate transfer, with or
without LayerNorm and $1/\sqrt{L}$-scaling. The above examples are taken after
20 epochs on Tiny ImageNet for (a)-(d) and ImageNet after $10$ epochs for
(e)-(f). In (a)-(c), the learning rate is linearly increased to the target in
the first $1000$ steps. In (d), after the same warm-up steps, the learning
rate is decreased to zero with a cosine schedule. Notice that in addition to
the transfer across all the settings, the $1/\sqrt{L}$ models show a stronger
non-saturating benefit of deeper models. (e)-(f) ViTs on first few epochs of
ImageNet. Also, notice that (a), (d), (e) and (f) refer to models _without
normalization layers_.

$\mu$P does not Transfer at Large Depths. While $\mu$P exhibits transfer over
network widths, it does not immediately give transfer over network depths. In
Figure 1(a), we show the train loss after $20$ epochs on CIFAR-10. Notice how
the network transfers over widths but not depths for $\mu$P. For completeness,
in Fig. B.2 we show that also SP does not transfer.

$1/\sqrt{L}$-scaled Residual Networks Transfer. We repeat the same experiment,
but now scale the residual branch by ${1}/{\sqrt{L}}$. The networks transfer
over both width and depth (Figure 1 (b)).

The Role of Normalization Layers. Normalization layers, such as LayerNorm and
BatchNorm, are commonly used instead of the $1/\sqrt{L}$-scaling in residual
architectures. In Fig. 2, we repeat the learning rate transfer experiment with
BatchNorm placed after the convolutional layer. In Fig. B.4 we repeat the same
experiment with LayerNorm. In both cases, while it is known that the presence
of normalization layers certainly helps the trainability of the model at
larger depth (Daneshmand et al., 2021; Joudaki et al., 2023), we observe that
the $1/\sqrt{L}$-scaled version exhibits a more consistent transfer across
both width and depth.

Empirical Verification on Vision Transformers with Adam. Despite our theory
being primarily on (S)GD-based training, we empirically test whether the
transfer results can be extended to Vision Transformers (ViTs) trained with
$\mu$P-parameterized Adam optimizer (Yang et al., 2021). In particular, we
consider Pre-LN Transformers (Xiong et al., 2020), and adapt it to our setting
by scaling all the residual branches by $1/\sqrt{L}$ (as in Noci et al.
(2022)), and make the architecture compatible with the $\mu$P framework of
Yang et al. (2021) (see Appendix A.3 for details). We test the variants both
with and without LayerNorm. Results for 20 epochs training with the Tiny
ImageNet dataset are shown in Fig. 3 (a-c), where we also see a good transfer
across all the settings. The observation of a good transfer with depth for
Pre-LN Transformers was empirically observed for translation tasks in Yang et
al. (2021). Here, we additionally notice that while the performance benefits
of Pre-LN Transformers saturates at large depth, the models with
$1/\sqrt{L}$-scaling exhibits a more consistent improvement with depth.
Finally, in Fig. 3(e-f) we successfully test scaling the learning rate
transfer experiment of $1/\sqrt{L}$-scaled ViTs (without LayerNorm) to
ImageNet (see Appendix B for ImageNet experiments with ResNets).

Other Hyperparameters and Learning Rate Schedules also Transfer. Given the
standard practice of using learning rate schedulers in optimizing Transformers
(Dosovitskiy et al., 2020; Touvron et al., ), we show that linear warm-up (Fig
3(a-c)), optionally followed by a cosine decay (Fig 3(d)) also transfer
consistently. Finally, in Figure 4, we plot the learning dynamics for
momentum, weight decay, and the feature learning rate $\gamma_{0}$, which
interpolates between the kernel and feature learning regimes (Bordelon &
Pehlevan, 2022b). Despite the usual correct transfer, we highlight the
_consistency_ of the dynamics, in the sense that learning curves with the same
hyperparameter but different depths are remarkably similar, especially at the
beginning of training.

(a) Momentum (b) Feature Learning Rate $\gamma_{0}$

Figure 4: Other hyperparameters also transfer. This example shows training
dynamics during the first $20$ epochs on CIFAR-10 (architecture details in
Appendix A). The dynamics at two different depths are provided for (a)
momentum (b) feature learning rate $\gamma_{0}$.

##  4 Convergence to the Large Width and Depth Limit

We argued that width/depth-invariant feature and prediction updates are
crucial for hyperparameter transfer. To formalize this, let $f_{N,L}(t)$ be
the network output after $t$ steps of optimization, where we explicitly write
the dependence on width $N$ and depth $L$. For hyperparameter transfer we want
$f_{N,L}(t)$ to be close to $f_{N^{\prime},L^{\prime}}(t)$. One strategy to
achieve this is to adopt a parameterization which makes finite networks as
close as possible to a limit $f_{\infty,\infty}(t)$. Following this line of
thought, we desire that a parameterization:

  1. 1.

Have a well-defined large depth and width limit, i.e.
$\lim_{N,L\to\infty}f_{N,L}(t)$.

  2. 2.

Have internal features / kernels evolve asymptotically _independent of width
$N$ or depth $L$_.

  3. 3.

Have minimal finite size approximation error.

These desiderata222In addition, an ideal parameterization would also have the
property that scaling up the width and depth of the model would monotonically
improve performance, which empirically seems to hold in our parameterization.
motivated the $\frac{1}{\sqrt{L}}$-$\mu$P ResNet in Equation 1 and Table 1.
The existence of the limit for the forward pass at initialization $t=0$ was
shown in Hayou & Yang (2023); Cirone et al. (2023). Here, we study the
learning dynamics in this model with $\mu$P parameterization, establish that
the limit exists and show that feature updates are width and depth
independent. Further, we characterize the limit throughout training in
Appendices C, D, E. We attempt to visualize convergence to the limit in Figure
5 and Appendix Figure B.6. See Appendix G and H for theoretical analysis of
finite size approximation error and Figure B.6 for additional tests of the
sources of finite width/depth approximation error.

(a) $N=512$ SGD Dynamics (b) Convergence in $L,N$

Figure 5: Approximation of the joint $N\to\infty$, $L\to\infty$ limit requires
both sufficiently large $N$ and $L$. CNNs are compared after $250$ steps on
CIFAR-10 with batchsize $32$ with $\gamma_{0}=1.0$. Since we cannot compute
the infinite width and depth predictor, we use a proxy
$f_{\approx\infty,\infty}$: the ensemble averaged (over random inits)
predictor of networks with $N=512$ and $L=128$. Error bars are standard
deviation computed over $15$ random initializations. (a) SGD dynamics at large
width are strikingly consistent across depths. (b) The convergence of
$f_{N,L}$ to a large width and depth proxy $f_{\approx\infty,\infty}$ is
bottlenecked by width $N$ for small $N$ while for large $N$, it decays like
$\mathcal{O}(L^{-1})$.

###  4.1 Primer on the DMFT in the $N\to\infty$ Limit

To derive our theoretical results, we use a technique known as dynamical mean
field theory (DMFT) Bordelon & Pehlevan (2022b). DMFT describes the dynamics
of neural networks at (or near) infinite width (preserving feature learning).
In this setting one tracks a set of kernel order parameters

$\displaystyle\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\frac{1}{N}\phi(\bm{h}^{\ell}(\bm{x},t))\cdot\phi(\bm{h}^{\ell}(\bm{x}^{\prime},s))\
,\quad
G^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\frac{1}{N}\bm{g}^{\ell}(\bm{x},t)\cdot\bm{g}^{\ell}(\bm{x}^{\prime},s),$
(4)

where $\bm{g}^{\ell}(\bm{x},t)\equiv N\gamma_{0}\frac{\partial
f(\bm{x},t)}{\partial\bm{h}^{\ell}(\bm{x},t)}$ are the back-propagation
signals. DMFT reveals the following facts about the infinite width limit:

  1. 1.

The dynamical kernels $\\{\Phi^{\ell}$, $G^{\ell}\\}$ and network outputs $f$
concentrate over random initializations of network weights throughout
training.

  2. 2.

The preactivations and gradient fields of each neuron
$\\{h_{i}^{\ell},g_{i}^{\ell}\\}$ become i.i.d. random variables drawn from a
single-site (i.e. neuron marginal) density which depends on the kernels.

  3. 3.

The kernels can be computed as averages over this single-site density.

The last two facts provide a self-consistency interpretation: given the
kernels, one can compute or sample the distribution of preactivations
$\\{h^{\ell},g^{\ell}\\}$. Given the density of $\\{h^{\ell},g^{\ell}\\}$, one
can compute the kernels. Prior work on this limit has demonstrated that these
equations can be approximately solved at finite $L$ with a Monte Carlo
procedure (but with computational complexity cubic in the training iterations)
(Bordelon & Pehlevan, 2022b).

###  4.2 The Large-Width-and-Depth Limit DMFT Equations

From the finite-$L$ DMFT equations (App. C), we calculate the large depth
$L\to\infty$ limit, generating a continuum process over the layer time
$\tau=\frac{\ell}{L}\in[0,1]$, a concept introduced in Li et al. (2022).
333Alternatively, one can also take the large-$L$ limit to get a limiting
distribution over order parameters (kernels, output logits, etc) at fixed $N$,
which yields the same limiting equations when $N\to\infty$ (see App. E).

######  Proposition 1 (Informal)

Consider network training dynamics of a ResNet as in Equation 1 with our $\mu
P$ $\frac{1}{\sqrt{L}}$ scaling and the appropriate choice of learning rate
and $\gamma_{0}$ in the infinite width and depth limit. The preactivation
$h(\tau;\bm{x},t)$ drawn from the marginal density of neurons in a layer at
layer-time $\tau$, data point $\bm{x}$ and training time $t$ obeys a
stochastic integral equation

$\displaystyle h(\tau;\bm{x};t)$
$\displaystyle=h(0;\bm{x};t)+\int_{0}^{\tau}du(\tau^{\prime};\bm{x};t)+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime}\int_{0}^{t}ds\int
d\bm{x}^{\prime}C_{h}(\tau^{\prime};\bm{x},\bm{x}^{\prime};t,s)g(\tau^{\prime};\bm{x}^{\prime};s)$
(5)

where $du(\tau;\bm{x};t)$ is zero mean Brownian motion. The covariance of $du$
is the feature kernel $\Phi$ and the deterministic operator $C_{h}$ can be
computed from the deterministic limiting kernels $\Phi,G$. The gradients
$g(\tau)$ satisfy an analogous integral equation (See App. D.7 for formulas
for $C_{h}$ and $g$). The stochastic variables $h,g$ and the kernels $\Phi,G$
satisfy a closed system of equations.

A full proposition can be found in Appendix D.7. We derive this limit in the
Appendix D, E. These equations are challenging to solve in the general case,
as the $h$ and $g$ distributions are non-Gaussian. Its main utility in the
present work is demonstrating that there is a well-defined feature-learning
infinite width and depth limit and that, in the limit, the predictions, and
kernels evolve by $\mathcal{O}(1)$.

###  4.3 Special Exactly Solveable Cases of the $N,L\to\infty$ DMFT equations

Though the resulting DMFT equations for the kernels are intractable in the
general case, we can provide analytical solutions in special cases to verify
our approach. We first write down the training dynamics in the lazy limit by
characterizing the neural tangent kernel (NTK) at initialization. Next, we
discuss deep linear ResNets where the feature learning dynamics gives ODEs
that can be closed while preserving Gaussianity of preactivations.

####  4.3.1 The Lazy Learning Limit

A special case of the above dynamics potentially of interest is the limit of
$\gamma_{0}\to 0$ where the feature kernels are frozen. In this limit the NTK
$K$ is static and can be computed at initialization and used to calculate the
full trajectory of the network dynamics (see Appendix F.1).

######  Proposition 2 (Informal)

Consider a ResNet with the $\frac{1}{\sqrt{L}}$-$\mu$P parameterization in the
$N,L\to\infty$ limit. In the $\gamma_{0}\to 0$ limit, the feature kernel
$\Phi(\tau)$, gradient kernel $G(\tau)$ and NTK $K$ do not evolve during
training and satisfy the following differential equations over layer time
$\tau$

$\displaystyle\partial_{\tau}H(\tau;\bm{x},\bm{x}^{\prime})$
$\displaystyle=\Phi(\tau;\bm{x},\bm{x}^{\prime})\ ,\
\Phi(\tau;\bm{x},\bm{x}^{\prime})=\left<\phi(h)\phi(h^{\prime})\right>_{h,h^{\prime}\sim\mathcal{N}(0,\bm{H}_{\bm{x},\bm{x}^{\prime}}(\tau))}$
$\displaystyle\partial_{\tau}G(\tau;\bm{x},\bm{x}^{\prime})$
$\displaystyle=-G(\tau;\bm{x},\bm{x}^{\prime})\left<\dot{\phi}(h)\dot{\phi}(h^{\prime})\right>_{h,h^{\prime}\sim\mathcal{N}(0,\bm{H}_{\bm{x},\bm{x}^{\prime}}(\tau))}$
$\displaystyle\bm{H}_{\bm{x},\bm{x}^{\prime}}(\tau)$
$\displaystyle=\begin{bmatrix}H(\tau;\bm{x},\bm{x})&H(\tau;\bm{x},\bm{x}^{\prime})\\\
H(\tau;\bm{x}^{\prime},\bm{x})&H(\tau;\bm{x}^{\prime},\bm{x}^{\prime})\end{bmatrix}\
,\quad K(\bm{x},\bm{x}^{\prime})=\int_{0}^{1}d\tau\
G(\tau,\bm{x},\bm{x}^{\prime})\Phi(\tau,\bm{x},\bm{x}^{\prime}).$ (6)

The dynamics of the neural network outputs $f(\bm{x},t)$ can be computed from
$K(\bm{x},\bm{x}^{\prime})$.

The forward ODE for $\Phi$ at initialization was studied by Hayou & Yang
(2023); Cirone et al. (2023). Figure 6 shows this limiting NTK for a ReLU
network. The finite depth kernels converge to the limiting large $L$ kernel
(blue line in (a)) at a rate $\mathcal{O}(L^{-2})$ in square error.

(a) ReLU NTKs at Initialization (b) Initial NTK Convergence in $L$

Figure 6: Convergence of the depth $L$ NTK to the limiting kernel at
initialization. (a) The infinite width NTK for a ReLU residual multilayer
perpectron (MLP) for two inputs $\bm{x},\bm{x}^{\prime}$ separated by angle
$\theta\in[-\pi,\pi]$ over varying depth $L$. The kernel converges as
$L\to\infty$. (b) The convergence rate of the initialization averaged NTK
$\left<K\right>$ with network widths $N$ and network depths $L$ at
initialization. For sufficiently large $N$, the initial NTK will converge at
rate $\mathcal{O}(L^{-2})$ in square error.

####  4.3.2 Deep Linear ResNets

Next we explore the dynamics of feature learning in a solveable infinite depth
ResNet with linear activations. In this case, the preactivations $h$ and
gradient fields $g$ remain Gaussian over the full time-course of learning,
rendering the DMFT equations tractable without any need for Monte Carlo
methods. As an example, in Figure 7 we show the kernel updates in a deep
linear ResNet after a single training step on a single training example
$(\bm{x},y)$. The theoretical solution to the kernel ODE at initialization
$H_{0}(\tau)=e^{\tau}H_{0}(0)$ and after one step $H_{1}(\tau)$ is (see
Appendix F.2.2 for derivation)

$\displaystyle
H_{1}(\tau)=e^{\tau}H_{1}(0)+2\eta_{0}^{2}\gamma_{0}^{2}y^{2}e^{\tau}\
\left[\frac{1}{3}\tau^{3}+\frac{1}{2}\tau^{2}+(1+e^{-1})(\tau^{2}e^{\tau}-\tau
e^{\tau}+e^{\tau}-1)\right]$ (7)

For networks with depths on the order of $\sim 100$ the infinite depth
solution gives an accurate description of the profile of the kernel updates.
More general equations for infinite depth linear networks are provided in
Appendix F.2.

(a) Kernels at initialization (b) Kernel Updates after 1 Step

Figure 7: We plot the feature kernels $H(\tau)=H^{\ell}|_{\ell=\tau L}$ in
width $N=1000$ residual linear networks before and after taking a gradient
step in our parameterization. The kernel profiles at initialization
$H_{0}(\tau)$ and the kernel changes $H_{1}(\tau)-H_{0}(\tau)$ converge for
large $L$. We compare to the theoretical result expected from the $L\to\infty$
limit of the DMFT equations (dashed black).

##  5 Related Works

The $1/\sqrt{L}$-scaled residual model (He et al., 2016) adopted in this work
belongs to a rapidly expanding research domain. This approach has several
advantageous characteristics in signal propagation, in terms of dynamical
isometry (Tarnowski et al., 2019), as well as stable backward and forward
kernels (Balduzzi et al., 2017; Hayou et al., 2021; Fischer et al., 2023). As
a result, this has paved the way for the development of reliable
initialization methods (Taki, 2017; Zhang et al., 2019), including
Transformers without normalization layers (Huang et al., 2020; Noci et al.,
2022).

The study of neural network scaling limits started with MLPs in the infinite-
width kernel regime (Jacot et al., 2018; Neal, 2012; Matthews et al., 2018;
Lee et al., 2017), where the preactivations are also Gaussian with
deterministic covariance. However, the joint depth and width limit has been
shown to be non-Gaussian (Hanin & Nica, 2020; Noci et al., 2021) and have a
stochastic covariance (Li et al., 2021). In particular, the feature covariance
of particular classes of networks have a limit described by stochastic
differential equations (SDEs) (Li et al., 2022; Noci et al., 2023), where the
layer time is the depth-to-width ratio in this case. On the other hand, Hayou
& Yang (2023) study the forward signal propagation in ${1}/{\sqrt{L}}$-scaled
ResNets at initialization, where the preactivations converge to a Gaussian
distribution with a _deterministic_ covariance described by an ODE.

Successful hyperparameter transfer across network widths in $\mu$P was shown
in Yang et al. (2021). Width-consistent early time prediction and
representation dynamics were also observed in an empirical study of Vyas et
al. (2023). The infinite width limits of feature learning networks have been
characterized with a mean-field PDE in two layer networks (Chizat & Bach,
2018; Mei et al., 2019; Chizat & Bach, 2020), Tensor programs (Yang, 2020;
Yang & Hu, 2021), and DMFT techniques (Bordelon & Pehlevan, 2022b; a; 2023),
which have been proven useful to understand the gradient dynamics in other
settings as well (Mignacco et al., 2020; Gerbelot et al., 2022). This present
work extends the DMFT technique to ResNets to analyze their feature learning
behavior at large depth $L$.

Understanding feature learning in non-residual feedforward networks close to
the asymptotic regime has also been studied with finite width and depth
corrections from the corresponding infinite limit, both in the lazy (Chizat et
al., 2019) NTK (Hanin, 2018; Hanin & Nica, 2019) and NNGP limits (Antognini,
2019; Yaida, 2020; Li & Sompolinsky, 2021; Zavatone-Veth et al., 2021; Segadlo
et al., 2022; Hanin, 2022), and more recently in the rich feature learning
regime (Bordelon & Pehlevan, 2023). At any finite width and depth the analysis
is even more challenging and has been studied in special cases using
combinatorial approaches and hypergeometric functions (Noci et al., 2021;
Zavatone-Veth & Pehlevan, 2021; Hanin & Zlokapa, 2023).

##  6 Discussion

In this work, we have proposed a simple parameterization of residual networks.
In this parameterization, we have seen empirically that hyperparameter
transfer is remarkably consistent across both width and depth for a variety of
architectures, hyperparameters, and datasets. The experimental evidence is
supported by theory, which shows that dynamics of all hidden layers are both
non-trivial and independent of width and depth while not vanishing in the
limit. We believe that these results have the potential to reduce
computational costs of hyperparameter tuning, allowing the practitioner to
train a large depth and width model only once with near optimal
hyperparameters.

Limitations. While we demonstrate the existence of an infinite width and depth
feature learning limit of neural network dynamics using statistical field
theory methods and empirical tests, our results are derived at the level of
rigor of physics, rather than formal proof. We also do not consider joint
scaling limits in which the dataset size and number of optimization steps can
scale with width $N$ or depth $L$. Instead, we treat time and sample size as
$\mathcal{O}(1)$ quantities. Due to computational limits, most of our
experiments are limited to $10$-$20$ epochs on each dataset.

Future Directions. This paper explores one possible scaling limit of infinite
width and depth networks. Other parameterizations and learning rate scaling
rules may also be worth exploring in future work (Jelassi et al., 2023; Li et
al., 2022; Noci et al., 2023). It would be especially interesting to see if an
stochastic limit for the kernels can be derived for the training dynamics in
ResNets without the $1/\sqrt{L}$ scaling factor. This parameterization makes
it especially convenient to study depth-width trade-offs while scaling neural
networks since hyperparameter dependencies are consistent across widths/depths
(see Figure B.1). It would also be interesting to study our limiting dynamics
in further detail to characterize what kinds of features are learned. Finally,
it remains an open problem theoretically why hyperparameters transfer across
widths and depths even when the actual predictions/losses differ across width
and depth.

## Acknowledgements

BB is supported by a Google PhD research fellowship. BB and CP are also
supported by NSF Award DMS-2134157. 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. BH and
MBL are supported by NSF grant DMS-2133806. BH is further supported by NSF
CAREER grant DMS-2143754, NSF grant DMS-1855684 and an ONR MURI on Foundations
of Deep Learning. We would like to thank Jacob Zavatone-Veth, Mary Letey, Alex
Atanasov for comments on early versions of this manuscript.

## References

  * Antognini (2019) Joseph M Antognini.  Finite size corrections for neural network gaussian processes.  _arXiv preprint arXiv:1908.10030_ , 2019.
  * Balduzzi et al. (2017) David Balduzzi, Marcus Frean, Lennox Leary, JP Lewis, Kurt Wan-Duo Ma, and Brian McWilliams.  The shattered gradients problem: If resnets are the answer, then what is the question?  In _International Conference on Machine Learning_ , pp. 342–350. PMLR, 2017.
  * Bordelon & Pehlevan (2022a) Blake Bordelon and Cengiz Pehlevan.  The influence of learning rule on representation dynamics in wide neural networks.  _arXiv preprint arXiv:2210.02157_ , 2022a.
  * Bordelon & Pehlevan (2022b) Blake Bordelon and Cengiz Pehlevan.  Self-consistent dynamical field theory of kernel evolution in wide neural networks.  _Advances in Neural Information Processing Systems_ , 35:32240–32256, 2022b.
  * Bordelon & Pehlevan (2023) Blake Bordelon and Cengiz Pehlevan.  Dynamics of finite width kernel and prediction fluctuations in mean field neural networks.  _arXiv preprint arXiv:2304.03408_ , 2023.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Chizat & Bach (2018) Lenaic Chizat and Francis Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  _Advances in neural information processing systems_ , 31, 2018.
  * Chizat & Bach (2020) Lenaic Chizat and Francis Bach.  Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss.  In _Conference on Learning Theory_ , pp. 1305–1338. PMLR, 2020\.
  * Chizat et al. (2019) Lenaic Chizat, Edouard Oyallon, and Francis Bach.  On lazy training in differentiable programming.  _Advances in neural information processing systems_ , 32, 2019.
  * Cho & Saul (2009) Youngmin Cho and Lawrence Saul.  Kernel methods for deep learning.  _Advances in neural information processing systems_ , 22, 2009.
  * Cirone et al. (2023) Nicola Muca Cirone, Maud Lemercier, and Cristopher Salvi.  Neural signature kernels as infinite-width-depth-limits of controlled resnets.  _arXiv preprint arXiv:2303.17671_ , 2023.
  * Daneshmand et al. (2021) Hadi Daneshmand, Amir Joudaki, and Francis Bach.  Batch normalization orthogonalizes representations in deep random networks.  _Advances in Neural Information Processing Systems_ , 34:4896–4906, 2021.
  * Dosovitskiy et al. (2020) Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Fedus et al. (2022) William Fedus, Barret Zoph, and Noam Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity.  _The Journal of Machine Learning Research_ , 23(1):5232–5270, 2022.
  * Fischer et al. (2023) Kirsten Fischer, David Dahmen, and Moritz Helias.  Optimal signal propagation in resnets through residual scaling.  _arXiv preprint arXiv:2305.07715_ , 2023.
  * Gerbelot et al. (2022) Cedric Gerbelot, Emanuele Troiani, Francesca Mignacco, Florent Krzakala, and Lenka Zdeborova.  Rigorous dynamical mean field theory for stochastic gradient descent methods.  _arXiv preprint arXiv:2210.06591_ , 2022.
  * Hanin (2018) Boris Hanin.  Which neural net architectures give rise to exploding and vanishing gradients?  _Advances in neural information processing systems_ , 31, 2018.
  * Hanin (2022) Boris Hanin.  Correlation functions in random fully connected neural networks at finite width.  _arXiv preprint arXiv:2204.01058_ , 2022.
  * Hanin & Nica (2019) Boris Hanin and Mihai Nica.  Finite depth and width corrections to the neural tangent kernel.  _arXiv preprint arXiv:1909.05989_ , 2019.
  * Hanin & Nica (2020) Boris Hanin and Mihai Nica.  Products of many large random matrices and gradients in deep neural networks.  _Communications in Mathematical Physics_ , 376(1):287–322, 2020.
  * Hanin & Zlokapa (2023) Boris Hanin and Alexander Zlokapa.  Bayesian interpolation with deep linear networks.  _Proceedings of the National Academy of Sciences_ , 120(23):e2301345120, 2023.
  * Hayou (2023) Soufiane Hayou.  On the infinite-depth limit of finite-width neural networks.  _Transactions on Machine Learning Research_ , 2023.  ISSN 2835-8856.  URL https://openreview.net/forum?id=RbLsYz1Az9.
  * Hayou & Yang (2023) Soufiane Hayou and Greg Yang.  Width and depth limits commute in residual networks.  _arXiv preprint arXiv:2302.00453_ , 2023.
  * Hayou et al. (2021) Soufiane Hayou, Eugenio Clerico, Bobby He, George Deligiannidis, Arnaud Doucet, and Judith Rousseau.  Stable resnet.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1324–1332. PMLR, 2021.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Hoffmann et al. (2022) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_ , 2022.
  * Huang et al. (2020) Xiao Shi Huang, Felipe Perez, Jimmy Ba, and Maksims Volkovs.  Improving transformer optimization through better initialization.  In _International Conference on Machine Learning_ , pp. 4475–4483. PMLR, 2020.
  * Itô (1944) Kiyosi Itô.  109\. stochastic integral.  _Proceedings of the Imperial Academy_ , 20(8):519–524, 1944.
  * Jacot et al. (2018) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Jelassi et al. (2023) Samy Jelassi, Boris Hanin, Ziwei Ji, Sashank J Reddi, Srinadh Bhojanapalli, and Sanjiv Kumar.  Depth dependence of $\mu$p learning rates in relu mlps.  _arXiv preprint arXiv:2305.07810_ , 2023.
  * Joudaki et al. (2023) Amir Joudaki, Hadi Daneshmand, and Francis Bach.  On the impact of activation and normalization in obtaining isometric embeddings at initialization.  _arXiv preprint arXiv:2305.18399_ , 2023.
  * Kaplan et al. (2020) Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models.  _arXiv preprint arXiv:2001.08361_ , 2020.
  * Klug & Heckel (2023) Tobit Klug and Reinhard Heckel.  Scaling laws for deep learning based image reconstruction.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=op-ceGueqc4.
  * Lee et al. (2017) Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  _arXiv preprint arXiv:1711.00165_ , 2017.
  * Lee et al. (2019) Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  _Advances in neural information processing systems_ , 32, 2019.
  * Li et al. (2021) Mufan Li, Mihai Nica, and Dan Roy.  The future is log-gaussian: Resnets and their infinite-depth-and-width limit at initialization.  _Advances in Neural Information Processing Systems_ , 34:7852–7864, 2021.
  * Li et al. (2022) Mufan Li, Mihai Nica, and Dan Roy.  The neural covariance sde: Shaped infinite depth-and-width networks at initialization.  _Advances in Neural Information Processing Systems_ , 35:10795–10808, 2022.
  * Li & Sompolinsky (2021) Qianyi Li and Haim Sompolinsky.  Statistical mechanics of deep linear neural networks: The backpropagating kernel renormalization.  _Physical Review X_ , 11(3):031059, 2021.
  * Martin et al. (1973) Paul Cecil Martin, ED Siggia, and HA Rose.  Statistical dynamics of classical systems.  _Physical Review A_ , 8(1):423, 1973.
  * Matthews et al. (2018) Alexander G de G Matthews, Mark Rowland, Jiri Hron, Richard E Turner, and Zoubin Ghahramani.  Gaussian process behaviour in wide deep neural networks.  _arXiv preprint arXiv:1804.11271_ , 2018.
  * Mei et al. (2019) Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Mean-field theory of two-layers neural networks: dimension-free bounds and kernel limit.  In _Conference on Learning Theory_ , pp. 2388–2464. PMLR, 2019\.
  * Mignacco et al. (2020) Francesca Mignacco, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborová.  Dynamical mean-field theory for stochastic gradient descent in gaussian mixture classification.  _Advances in Neural Information Processing Systems_ , 33:9540–9550, 2020.
  * Neal (2012) Radford M Neal.  _Bayesian learning for neural networks_ , volume 118.  Springer Science & Business Media, 2012.
  * Noci et al. (2021) Lorenzo Noci, Gregor Bachmann, Kevin Roth, Sebastian Nowozin, and Thomas Hofmann.  Precise characterization of the prior predictive distribution of deep relu networks.  _Advances in Neural Information Processing Systems_ , 34:20851–20862, 2021.
  * Noci et al. (2022) Lorenzo Noci, Sotiris Anagnostidis, Luca Biggio, Antonio Orvieto, Sidak Pal Singh, and Aurelien Lucchi.  Signal propagation in transformers: Theoretical perspectives and the role of rank collapse.  _Advances in Neural Information Processing Systems_ , 35:27198–27211, 2022.
  * Noci et al. (2023) Lorenzo Noci, Chuning Li, Mufan Bill Li, Bobby He, Thomas Hofmann, Chris Maddison, and Daniel M Roy.  The shaped transformer: Attention models in the infinite depth-and-width limit.  _arXiv preprint arXiv:2306.17759_ , 2023.
  * OpenAI (2023) OpenAI.  Gpt-4 technical report, 2023.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Rudin (1953) Walter Rudin.  _Principles of mathematical analysis_.  1953\.
  * Segadlo et al. (2022) Kai Segadlo, Bastian Epping, Alexander van Meegen, David Dahmen, Michael Krämer, and Moritz Helias.  Unified field theoretical approach to deep and recurrent neuronal networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2022(10):103401, 2022.
  * Taki (2017) Masato Taki.  Deep residual networks and weight initialization.  _arXiv preprint arXiv:1709.02956_ , 2017.
  * Tarnowski et al. (2019) Wojciech Tarnowski, Piotr Warchoł, Stanisław Jastrzebski, Jacek Tabor, and Maciej Nowak.  Dynamical isometry is achieved in residual networks in a universal way for any activation function.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 2221–2230. PMLR, 2019.
  * (53) Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al.  Llama: open and efficient foundation language models, 2023.  _URL https://arxiv. org/abs/2302.13971_.
  * Vyas et al. (2023) Nikhil Vyas, Alexander Atanasov, Blake Bordelon, Depen Morwani, Sabarish Sainathan, and Cengiz Pehlevan.  Feature-learning networks are consistent across widths at realistic scales, 2023.
  * Xiong et al. (2020) Ruibin Xiong, Yunchang Yang, Di He, Kai Zheng, Shuxin Zheng, Chen Xing, Huishuai Zhang, Yanyan Lan, Liwei Wang, and Tieyan Liu.  On layer normalization in the transformer architecture.  In _International Conference on Machine Learning_ , pp. 10524–10533. PMLR, 2020.
  * Yaida (2020) Sho Yaida.  Non-gaussian processes and neural networks at finite widths.  In _Mathematical and Scientific Machine Learning_ , pp. 165–192. PMLR, 2020.
  * Yang et al. (2021) Ge Yang, Edward Hu, Igor Babuschkin, Szymon Sidor, Xiaodong Liu, David Farhi, Nick Ryder, Jakub Pachocki, Weizhu Chen, and Jianfeng Gao.  Tuning large neural networks via zero-shot hyperparameter transfer.  _Advances in Neural Information Processing Systems_ , 34:17084–17097, 2021.
  * Yang (2020) Greg Yang.  Tensor programs ii: Neural tangent kernel for any architecture.  _arXiv preprint arXiv:2006.14548_ , 2020.
  * Yang & Hu (2021) Greg Yang and Edward J Hu.  Tensor programs iv: Feature learning in infinite-width neural networks.  In _International Conference on Machine Learning_ , pp. 11727–11737. PMLR, 2021.
  * Zavatone-Veth & Pehlevan (2021) Jacob Zavatone-Veth and Cengiz Pehlevan.  Exact marginal prior distributions of finite bayesian neural networks.  _Advances in Neural Information Processing Systems_ , 34:3364–3375, 2021.
  * Zavatone-Veth et al. (2021) Jacob Zavatone-Veth, Abdulkadir Canatar, Ben Ruben, and Cengiz Pehlevan.  Asymptotics of representation learning in finite bayesian neural networks.  _Advances in neural information processing systems_ , 34:24765–24777, 2021.
  * Zhai et al. (2022) Xiaohua Zhai, Alexander Kolesnikov, Neil Houlsby, and Lucas Beyer.  Scaling vision transformers.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 12104–12113, 2022.
  * Zhang et al. (2019) Hongyi Zhang, Yann N Dauphin, and Tengyu Ma.  Fixup initialization: Residual learning without normalization.  _arXiv preprint arXiv:1901.09321_ , 2019.

## Appendix

##  Appendix A Experimental Details

###  A.1 Simplified ResNet

We design a residual architecture following the footprint of Eq. 1, and
replacing the weight matrices with $1$-strided convolutional layers with
$3\times 3$ kernel size. We max-pool the feature map at three selected layers
equally distributed across depth. For instance, for a depth=$12$ model, we
perform average pooling after the third, 6-th and 9-th layer. After each
average pooling, we apply an extra convolutional layer to double the number of
channels, which play the role of width in a convolutional network. Hence, also
the width doubled three times across the depth (one time per average pooling).
The final feature map is then flattened and passed through a final readout
layer. The reported width in all the figures is the last width before the
flattening operations. In the experiment with the $1/\sqrt{L}$-scaling, we
divide by the total number of layers instead of the number of residual blocks
(i.e. $L=3+\text{\\# blocks}$). This extra factor does not affect the scaling
limit.

###  A.2 ResNet

Compared to the simplified ResNet architecture of the previous section, the
ResNet family as proposed in He et al. (2016) are augmented with two
convolutional layers per residual block, each followed by batch normalization
normalization. The subsampling is performed by $2$-strided convolutional
layers, for a total of $4$ times across the depth. For further details, we
refer to He et al. (2016). In some experiments (e.g. Fig. B.4), we replace
batch norm with LayerNorm applied across each two dimensional channel.

###  A.3 Vision Transformer

We use the standard ViT architecture as in Dosovitskiy et al. (2020), i.e. we
tokenize the input images by splitting it into patches of fixed size, and
embed it to $N$ dimension using a linear layer. Then we apply a number of the
usual Transformer blocks, each one composed of two residual blocks, a Softmax-
based self attention block, and an MLP block with GeLU activation function. If
used, we place LayerNorm at the beginning of each residual branch (Pre-LN
configuration). Finally, we apply the following modification to make it
compatible with $\mu$P (Yang et al., 2021):

  * •

We re-scale the Softmax logits (i.e. the product of queries and keys) of each
attention layer by $1/N_{q}$. The usual convention would be to use instead
$1/\sqrt{N_{q}}$.

  * •

Initialize the queries to zero.

  * •

Initialize all the other weights to have standard deviation of $1/\sqrt{N}$,
where $N$ here is the model dimension.

###  A.4 Details

We always use data augmentation, in particular the following random
transformations:

  * •

CIFAR-10: horizontal flips, 10-degree rotations, affine transformations, color
jittering.

  * •

Tiny ImageNet and ImageNet: crops, horizontal flips.

Across all the experiments, we fix the batch size to 64. Unless stated
otherwise, all the hyperparameters that are not tuned are set to a default
value: $\beta,\gamma_{0},\sigma_{\ell}=1$. For ViTs, we use a patch size We
use the $\mu$P implementation of the $\mu$-Readout layer, and optimizers (SGD,
and Adam under $\mu$P parametrizstion) as in the released $\mu$P package.

##  Appendix B Additional Experiments

###  B.1 Depth Width Trade-Off

While tuning the architecture at small widths/depths, it is interesting to
study along which direction — width or depth — the architecture should be
scaled. To test this, we train the residual convolutional model (Sec. A.1) at
relatively large scale (up to almost a billion parameters), for 20 epochs on
CIFAR-10 with fixed learning rate of $0.046$ (using data augmentation, see
Sec. A). We report the results in Fig. B.1, where we plot the train loss as a
function of parameters. In particular, in Fig. B.1(a) we highlight models of
same widths, while in Fig. B.1(a) models of equal depth. We find that the
Pareto optimal curve coincides by the depth $12$ model for the range of
widths/depths used for the experiment. As shown in Fig. B.1, we then fit power
laws to the curves of constant depth, which seems to predict an optimal depth
30 at very large widths. Notice that we fix the number of epochs. We
hypothesize the optimal curve to shift to different trade-offs as the number
of epochs increases. We also expect different trade-offs for different
datasets and models. However, extensively studying width-depth trade-off is
beyond the scope of this work, thus we leave it as a potential direction for
future work.

(a) Envelope of Depth-Width Trade-off (b) Optimal Intermediate Depths

Figure B.1:  Since all models have hyperparameter transfer across depths and
widths, one can explore a trade-off between depth and width at the optimal
learning rate. We illustrate this idea with a simple experiment on CIFAR-10.
(a) We plot the final train loss as a function of parameters for different
widths. (b) intermediate depth of $\sim 9-12$ is preferred at fixed parameter
count.

(a) ResNet, SP (b) $\frac{1}{\sqrt{L}}$ ResNet SP

Figure B.2: Learning rates transfer across both depth and width in our
proposed parameterization but not in standard parameterization (SP) or $\mu$P.
Loss is plotted after $20$ epochs on CIFAR-10. All the missing datapoints
indicate that the corresponding run diverged. In particular, all depth $30$
runs diverged when using $\mu$P without $1/\sqrt{\text{depth}}$ scaling.

###  B.2 How Long Should you Train for Accurate Transfer?

An important question is after how many steps the optimal hyperparameters are
fixed and can hence be selected without training for longer. In Fig. B.3(a),
we plot the training loss profiles for two different depths and 4 different
learning rates for a residual convolutional network of width $512$ trained on
CIFAR-10. Notice how after the first few epochs it is already possible to
decide the optimal learning rate. We confirm this for a ViT trained on Tiny
ImageNet in Fig. B.3(b-c), where the optimal learning rate at two selected
epochs (3 and 9) is the same.

(a) LR dynamics (b) Epoch 3 (c) Epoch 9

Figure B.3: (a) dynamics of train loss for different depths and learning
rates. The network is the residual convolutional network trained on CIFAR-10.
(b-c) learning rate tranfer plots for the same ViT at two different
checkpoints (epoch $3$ and $5$). Notice how the optimal learning rate is
determined after the first few epochs of training.

###  B.3 Other Experiments

In Fig. B.2, we show how the learning rate does not transfer under SP for the
residual model considered in this work, both with and without the
$1/\sqrt{L}$-scaling. In Fig. B.4, we train models from the ResNets family,
replacing BatchNorm with LayerNorm, concluding that the version with
$1/\sqrt{L}$ scaling exhibits a better transfer. In Fig. B.5 we show our
results on learning rate transfers for ImageNet, and in Fig. B.7, we show how
also weight decay transfer under the proposed $1/\sqrt{L}$-scaling and has
very consistent dynamics. The model used is the simplified residual
convolutional model, trained for 20 epochs on CIFAR-10.

(a) ResNet + $\mu$P, LayerNorm (b) $\frac{1}{\sqrt{L}}$ ResNet + $\mu$P,
LayerNorm

Figure B.4: Effect of normalization layers. The above examples are taken after
$20$ epochs on CIFAR-10. Normalization layers can slightly improve consistency
across depths in standard ResNets, but consistency is much more reliable with
the $1/\sqrt{L}$ scaling.

(a) Dynamics (b) Learning Rate Transfer

Figure B.5: Learning rate + depth sweep on ImageNet for $5$ epochs with depth
extrapolations of ResNet-18 _without any normalization layers_. Training is
performed with SGD with momentum $0.9$ and batchsize $128$. (a) The loss
curves for different learning rates (each learning rate is a distinct
linestyle) and for different depths (colors). Dynamics are strikingly
consistent across depths for each learning rate. (b) The optimal learning rate
is preserved across depths.

(a) Fixed $L$ Convergence in $N$ (b) Average Square Error (c) Error of
Ensembled Logits

Figure B.6: Further numerical support of our finite $N,L$ error decomposition
in Appendix H through ensembling logits over $E=15$ random seeds in the same
experimental setting as Figure 5. (a) At each depth $L$, the finite width $N$
network converges in square error as $\mathcal{O}_{N,L}(N^{-1})$, independent
of $L$. (b) The full expected square error from Figure 5 shows a convergence
rate of $\mathcal{O}(L^{-1})$ for sufficiently large $N$, consistent with
theory. (c) The ensemble averaged logits $\left<f_{L,N}\right>$ converge to
the infinite depth proxy at a faster rate, which is consistent with the
theoretical $\mathcal{O}(L^{-2})$ for sufficiently large $N$.

(a) Weight Decay Dynamics (b) Final Loss vs Depth

Figure B.7: Weight decay dynamics also transfer across depths.

##  Appendix C Derivation of the $N\to\infty$ limit at fixed $L$

In this section, we characterize the feature learning dynamics of the
$N\to\infty$ limit of fixed depth $L$ residual networks. Our derivation is an
extension of the approach presented in (Bordelon & Pehlevan, 2022b). Before
taking the limit, we first identify our order parameters and weight space
dynamics. Then we take the $N\to\infty$ limit using a technique known as
Dynamical Mean Field Theory (DMFT).

###  C.1 Finite $N,L$ Function and Weight Dynamics

We start by computing the neural tangent kernel (NTK). For the sake of
brevity, we will first write the equations where the read-in $\bm{W}^{0}$
weights are held fixed and will add their contribution to the dynamics in
Appendix J.3. Excluding these two weight matrices from the present discussion
does not alter any of the scalings with width $N$ or depth $L$. The dynamics
of the prediction $f$ under gradient flow take the form

$\displaystyle\frac{d}{dt}f(\bm{x};t)=\eta_{0}\
\mathbb{E}_{\bm{x}^{\prime}}K(\bm{x},\bm{x}^{\prime};t)\Delta(\bm{x}^{\prime};t)$
(8)

where $\Delta(\bm{x};t)=-\frac{\partial\mathcal{L}(f(\bm{x},t))}{\partial
f(\bm{x},t)}$ is an error signal and $K$ is the NTK of this model, which takes
the form

$\displaystyle
K(\bm{x},\bm{x}^{\prime})=\gamma_{0}^{2}N\sum_{\ell=1}^{L}\frac{\partial
f(\bm{x})}{\partial\bm{W}^{\ell}}\cdot\frac{\partial
f(\bm{x}^{\prime})}{\partial\bm{W}^{\ell}}.$ (9)

Explicitly calculating the gradients of the output with respect to the weights
$\bm{W}^{\ell}$ gives

$\displaystyle N\gamma_{0}\frac{\partial
f(\bm{x})}{\partial\bm{W}^{\ell}}=\frac{1}{\sqrt{NL}}\bm{g}^{\ell+1}(\bm{x})\phi(\bm{h}^{\ell}(\bm{x}))^{\top}\
,\ g_{i}^{\ell}\equiv N\gamma_{0}\frac{\partial f(\bm{x})}{\partial
h_{i}^{\ell}(\bm{x})}\sim\mathcal{O}_{N,L}(1)$ (10)

Using this above formula, we find that the NTK is $\mathcal{O}(1)$ with
respect to $N,L$

$\displaystyle
K(\bm{x},\bm{x}^{\prime})=\frac{1}{L}\sum_{\ell=1}^{L}G^{\ell+1}(\bm{x},\bm{x}^{\prime})\Phi^{\ell}(\bm{x},\bm{x}^{\prime})$
(11)

where we introduced the $\mathcal{O}(1)$ feature kernel $\Phi^{\ell}$ and
gradient kernels $G^{\ell}$ which have the form

$\displaystyle\Phi^{\ell}(\bm{x},\bm{x}^{\prime})=\frac{1}{N}\phi(\bm{h}^{\ell}(\bm{x}))\cdot\phi(\bm{h}^{\ell}(\bm{x}^{\prime}))\
,\
G^{\ell}(\bm{x},\bm{x}^{\prime})=\frac{1}{N}\bm{g}^{\ell}(\bm{x})\cdot\bm{g}^{\ell}(\bm{x}^{\prime})$
(12)

The base cases for the feature and gradient kernels are given by

$\displaystyle\Phi^{0}(\bm{x},\bm{x}^{\prime})=\frac{1}{D}\bm{x}\cdot\bm{x}^{\prime}\
,\ G^{L+1}(\bm{x},\bm{x}^{\prime})=1.$ (13)

We next start by computing the dynamics of the weights under gradient flow

$\displaystyle\frac{d}{dt}\bm{W}^{\ell}=\frac{\eta_{0}\gamma_{0}}{\sqrt{LN}}\
\mathbb{E}_{\bm{x}}\
\Delta(\bm{x};t)\bm{g}^{\ell+1}(\bm{x},t)\phi(\bm{h}^{\ell}(\bm{x},t))^{\top}.$
(14)

Integrating this equation from time $0$ to time $t$, we get the following

$\displaystyle\bm{h}^{\ell+1}(\bm{x},t)=$
$\displaystyle\bm{h}^{\ell}(\bm{x},t)+\frac{1}{\sqrt{NL}}\bm{W}^{\ell}(0)\phi(\bm{h}^{\ell})$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\int_{0}^{t}ds\mathbb{E}_{\bm{x}^{\prime}}\Delta(\bm{x}^{\prime},s)\bm{g}^{\ell+1}(\bm{x}^{\prime},s)\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
(15)

where
$\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\frac{1}{N}\phi(\bm{h}^{\ell}(\bm{x},t))\cdot\phi(\bm{h}^{\ell}(\bm{x}^{\prime};s))$
is the double time-index generalization of the feature kernel. Similarly, we
obtain a backward pass relation of the form

$\displaystyle\bm{g}^{\ell}(\bm{x},t)=$
$\displaystyle\bm{g}^{\ell+1}(\bm{x},t)+\frac{1}{\sqrt{NL}}\dot{\phi}(\bm{h}^{\ell}(\bm{x},t))\odot\bm{W}^{\ell}(0)^{\top}\bm{g}^{\ell+1}$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\dot{\phi}(\bm{h}^{\ell}(\bm{x},t))\odot\int_{0}^{t}ds\mathbb{E}_{\bm{x}^{\prime}}\Delta(\bm{x}^{\prime},s)G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)\phi(\bm{h}^{\ell}(\bm{x}^{\prime},s))$
(16)

These equations define update rules for the preactivations $\bm{h}^{\ell}$ and
the gradients $\bm{g}^{\ell}$ in terms of the feature and gradient kernels
$\\{\Phi^{\ell},G^{\ell}\\}$. However the above equations explicitly depend on
the random initial weights $\bm{W}^{\ell}(0)$ through which time-varying
signals propagate. Our ultimate goal is to gain insight into the effective
dynamics over random initializations in the large $N$ limit. To do that, we
will invoke dynamical mean field theory. In the next section, we characterize
how these terms behave at large $N$ using a saddle point argument.

###  C.2 Intuitive Explanation of Large Width Feature Learning Limit

In the limit of $N\to\infty$ with $\gamma_{0}$ held fixed, the above equations
can be simplified. The two key effects that occur at large width are

  * •

Kernels (feature kernels, gradient kernels, NTK) and output logits concentrate
by a law of large numbers.

  * •

At each layer, each of the hidden neuron’s activities becomes independent
draws from a single site stochastic process which is characterized by the
kernels.

Once the single site density is known, the kernels can be computed as averages
over the single site density. The next section derives this result formally
from a Martin-Siggia-Rose Path Integral derivation of the Dynamical Mean Field
Theory (Martin et al., 1973).

###  C.3 Path-integral and Saddle Point Equations in $N\to\infty$ Limit

This computation sketches the DMFT derivation of the limiting large $N$
process. Examples of detailed computations of the DMFT action $\mathcal{S}$
can be found in the Appendix of Bordelon & Pehlevan (2022b; a). Our notation
is chosen to match the derivations found in their Appendix.

To characterize the effect of the random initial weights, we will attempt to
characterize the distribution of
$\bm{h}^{\ell}(\bm{x},t),\bm{g}^{\ell}(\bm{x},t)$ by tracking the moment
generating functional of these fields

$\displaystyle
Z[\\{\bm{j}^{h},\bm{j}_{g}\\}]=\mathbb{E}_{\\{\bm{W}^{\ell}(0)\\}}\exp\left(\sum_{\ell=1}^{L}\int
d\bm{x}\int dt\
[\bm{j}_{h}^{\ell}(\bm{x},t)\cdot\bm{h}^{\ell}(\bm{x},t)+\bm{j}_{g}(\bm{x},t)\cdot\bm{g}^{\ell}(\bm{x},t)]\right)$
(17)

Moments of the fields can be obtained with differentiation near zero source,
for example

$\displaystyle\lim_{\\{\bm{j}_{h}^{\ell},\bm{j}_{g}^{\ell}\\}\to
0}\frac{\delta}{\delta
j^{\ell}_{h,i_{1}}(\bm{x}_{1},t_{1})}...\frac{\delta}{\delta
j^{\ell}_{h,i_{k}}(\bm{x}_{k},t_{k})}Z[\\{\bm{j}^{h},\bm{j}_{g}\\}]=\left<h_{i_{1}}^{\ell_{1}}(\bm{x}_{1},t_{1})...h_{i_{k}}^{\ell_{k}}(\bm{x}_{k},t_{k})\right>.$
(18)

Now, following Bordelon & Pehlevan (2022b), we introduce the following useful
auxiliary fields

$\displaystyle\bm{\chi}^{\ell+1}(\bm{x},t)=\frac{1}{\sqrt{N}}\bm{W}^{\ell}(0)\phi(\bm{h}^{\ell}(\bm{x},t))\
,\
\bm{\xi}^{\ell}(\bm{x},t)=\frac{1}{\sqrt{N}}\bm{W}^{\ell}(0)^{\top}\bm{g}^{\ell+1}(\bm{x},t)$
(19)

The choice to introduce these auxiliary fields is apparent once one realizes
that the $h$ and $g$ fields can be regarded as functions of $\chi,\xi$ and the
kernels $\Phi,G$. We enforce the constraints defining these fields with Dirac-
delta functions by repeatedly multiplying by unity

$\displaystyle 1$ $\displaystyle=\int d\bm{\chi}^{\ell+1}(\bm{x},t)\
\delta\left(\bm{\chi}^{\ell+1}(\bm{x},t)-\frac{1}{\sqrt{N}}\bm{W}^{\ell}(0)\phi(\bm{h}^{\ell}(\bm{x},t))\right)$
$\displaystyle=\int\frac{d\hat{\bm{\chi}}^{\ell+1}(\bm{x},t)d\bm{\chi}^{\ell+1}(\bm{x},t)}{(2\pi)^{N}}\exp\left(i\hat{\bm{\chi}}^{\ell+1}(\bm{x},t)\cdot\left[\bm{\chi}^{\ell+1}(\bm{x},t)-\frac{1}{\sqrt{N}}\bm{W}^{\ell}(0)\phi(\bm{h}^{\ell}(\bm{x},t))\right]\right)$
(20)

for all $\ell,\bm{x},t$. And similarly for the $\bm{g}^{\ell}$ fields. After
this insertion, the averages over $\\{\bm{W}^{\ell}(0)\\}$ can be performed
with Gaussian integration. Concretely, we obtain the following average over
Gaussian initialization $\bm{W}^{\ell}(0)$

$\displaystyle\ln\mathbb{E}\exp\left(-\frac{i}{\sqrt{N}}\text{Tr}\
\bm{W}^{\ell}(0)^{\top}\int d\bm{x}\int
dt\left[\bm{\hat{\chi}}^{\ell+1}(\bm{x},t)\phi(\bm{h}^{\ell}(\bm{x},t))^{\top}+\bm{g}^{\ell+1}(\bm{x},t)\bm{\xi}^{\ell}(\bm{x},t)^{\top}\right]\right)$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\int dtdsd\bm{x}d\bm{x}^{\prime}\
\bm{\hat{\chi}}^{\ell+1}(\bm{x},t)\cdot\bm{\hat{\chi}}^{\ell+1}(\bm{x}^{\prime},s)\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle-\frac{1}{2}\int dtdsd\bm{x}d\bm{x}^{\prime}\
\bm{\hat{\xi}}^{\ell}(\bm{x},t)\cdot\bm{\hat{\xi}}^{\ell}(\bm{x}^{\prime},s)G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle-i\int dtdsd\bm{x}d\bm{x}^{\prime}\
\bm{\hat{\chi}}^{\ell+1}(\bm{x},t)\cdot\bm{g}^{\ell+1}(\bm{x}^{\prime},s)\
A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$ (21)

where we introduced the following order parameters

$\displaystyle\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=\frac{1}{N}\phi(\bm{h}^{\ell}(\bm{x},t))\cdot\phi(\bm{h}^{\ell}(\bm{x}^{\prime},s))\
,\
G^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\frac{1}{N}\bm{g}^{\ell}(\bm{x},t)\cdot\bm{g}^{\ell}(\bm{x}^{\prime},s)$
(22) $\displaystyle A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=-\frac{i}{N}\phi(\bm{h}^{\ell}(\bm{x},t))\cdot\bm{\hat{\xi}}^{\ell}(\bm{x}^{\prime},s)$
(23)

As we did with the $\bm{\chi},\bm{\xi}$ fields, we can enforce the definition
of the above order parameters with integral representations of Dirac-delta
functions. For example, the $\Phi$ kernels can be enforced by inserting

$\displaystyle
1=\int\frac{d\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)d\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)}{2\pi
iN^{-1}}$
$\displaystyle\exp\left(N\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\right)$
$\displaystyle\exp\left(-\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\phi(\bm{h}^{\ell}(\bm{x},t))\cdot\phi(\bm{h}^{\ell}(\bm{x}^{\prime},t^{\prime}))\right)$
(24)

where the $\hat{\Phi}$ integral is performed along the imaginary axis
$(-i\infty,i\infty)$ in the complex plane. The same trick is performed for
$G^{\ell}$ and $A^{\ell}$ which have conjugate variables
$\hat{G}^{\ell},B^{\ell}$ respectively. After inserting these new variables,
we find that the moment generating function can be expressed as an integral
over a collection of order parameters $\bm{q}$,

$\displaystyle\bm{q}=\text{Vec}\\{\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s),\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t.s),G^{\ell}(\bm{x},\bm{x}^{\prime},t,s),\hat{G}^{\ell}(\bm{x},\bm{x}^{\prime},t,s),A^{\ell}(\bm{x},\bm{x}^{\prime},t,s)B^{\ell}(\bm{x},\bm{x}^{\prime},t,s)\\},$
(25)

which is vectorized over all layers $\ell$ time points $t,s$ and samples
$\bm{x},\bm{x}^{\prime}$. The moment generating function takes the form

$\displaystyle Z[\\{\bm{j}_{h},\bm{j}_{g}\\}]=\int
d\bm{q}\exp\left(N\mathcal{S}[\bm{q},\bm{j}_{h},\bm{j}_{g}]\right)$ (26)

where $\mathcal{S}$ is the $\mathcal{O}_{N}(1)$ DMFT action which has the
following form

$\displaystyle\mathcal{S}=$ $\displaystyle\sum_{\ell}\int
d\bm{x}d\bm{x}^{\prime}dtds\left[\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t.s)+G^{\ell}(\bm{x},\bm{x}^{\prime},t,s)\hat{G}^{\ell}(\bm{x},\bm{x}^{\prime},t,s)\right]$
$\displaystyle-\sum_{\ell}\int d\bm{x}d\bm{x}^{\prime}dtds\
A^{\ell}(\bm{x},\bm{x}^{\prime},t,s)B^{\ell}(\bm{x}^{\prime},\bm{x},s,t)+\frac{1}{N}\sum_{i=1}^{N}\mathcal{Z}_{i}[\bm{q},\bm{j}_{h},\bm{j}_{g}]$
(27)

The functions $\mathcal{Z}_{i}$ are single-site moment generating functions
(MGFs) which express the distribution of $h_{i},g_{i}$ for a given value of
the order parameters $\bm{q}$ and sources $\bm{j}_{h},\bm{j}_{g}$.

$\displaystyle\mathcal{Z}_{i}=\int\prod_{\ell,\bm{x},t}d\hat{\chi}^{\ell}(\bm{x},t)d\hat{\chi}^{\ell}(\bm{x},t)d\hat{\chi}^{\ell}(\bm{x},t)d\hat{\chi}^{\ell}(\bm{x},t)\exp\left(-\sum_{\ell}\mathcal{H}_{i}^{\ell}[\hat{\chi},\hat{\xi},\chi,\xi,j_{h,i},j_{g,i}]\right)$
(28)

where $\mathcal{H}_{i}^{\ell}$ are the single-site Hamiltonians

$\displaystyle\mathcal{H}_{i}^{\ell}=$ $\displaystyle\frac{1}{2}\int
d\bm{x}d\bm{x}^{\prime}dtds\left[\Phi^{\ell-1}(\bm{x},\bm{x}^{\prime};t,s)\hat{\chi}^{\ell}(\bm{x},t)\hat{\chi}^{\ell}(\bm{x}^{\prime},s)+G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)\hat{\xi}^{\ell}(\bm{x},t)\hat{\xi}^{\ell}(\bm{x}^{\prime},s)\right]$
$\displaystyle+\int
d\bm{x}d\bm{x}^{\prime}dtds\left[\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\phi(h^{\ell}(\bm{x},t))\phi(h^{\ell}(\bm{x}^{\prime},t^{\prime}))+\hat{G}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)g^{\ell}(\bm{x},t)g^{\ell}(\bm{x}^{\prime},s)\right]$
$\displaystyle+i\int
d\bm{x}d\bm{x}^{\prime}dtds\left[A^{\ell-1}(\bm{x},\bm{x}^{\prime};t,s)\hat{\chi}^{\ell}(\bm{x},t)g^{\ell}(\bm{x}^{\prime},s)+B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\hat{\xi}^{\ell}(\bm{x},t)\phi(h^{\ell}(\bm{x}^{\prime},s))\right]$
$\displaystyle-\int
d\bm{x}dt\left[i\hat{\chi}^{\ell}(\bm{x},t)\chi^{\ell}(\bm{x},t)+i\hat{\xi}^{\ell}(\bm{x},t)\xi^{\ell}(\bm{x},t)+j_{h,i}^{\ell}(\bm{x},t)h^{\ell}(\bm{x},t)+j_{g,i}^{\ell}(\bm{x},t)g^{\ell}(\bm{x},t)\right]$
(29)

In the above formulas, the $h,g$ fields should be regarded as functions of the
$\chi,\xi$ fields.

Since the moment generating function has the form $Z=\int
d\bm{q}\exp\left(N\mathcal{S}[\bm{q}]\right)$, the $N\to\infty$ limit can thus
be characterized by the steepest-descent method (saddle-point integration over
$\bm{q}$). The result is that all of the order parameters take on definite
values (concentrate) at infinite width. The saddle point equations
$\frac{\partial\mathcal{S}}{\partial\bm{q}}=0$ state that each of the kernels
takes on definite values

$\displaystyle\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left<\phi(h^{\ell}(\bm{x},t))\phi(h^{\ell}(\bm{x}^{\prime},s))\right>_{i}$
$\displaystyle G^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left<g^{\ell}(\bm{x},t)g^{\ell}(\bm{x}^{\prime},s)\right>_{i}$
$\displaystyle A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=-\frac{i}{N}\sum_{i=1}^{N}\left<\phi(h^{\ell}(\bm{x},t))\hat{\xi}^{\ell}(\bm{x}^{\prime},s)\right>_{i}$
$\displaystyle B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=-\frac{i}{N}\sum_{i=1}^{N}\left<g^{\ell+1}(\bm{x},t)\hat{\chi}^{\ell+1}(\bm{x}^{\prime},s)\right>_{i}$
$\displaystyle\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=-\frac{1}{2N}\sum_{i=1}^{N}\left<\hat{\chi}^{\ell}(\bm{x},t)\hat{\chi}^{\ell}(\bm{x}^{\prime},s)\right>_{i}$
$\displaystyle\hat{G}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=-\frac{1}{2N}\sum_{i=1}^{N}\left<\hat{\xi}^{\ell}(\bm{x},t)\hat{\xi}^{\ell}(\bm{x}^{\prime},s)\right>_{i}$
(30)

where $\left<\right>_{i}$ denotes averaging over a single site stochastic
process defined by the moment generating function $\mathcal{Z}_{i}$. At zero
source $\bm{j}\to 0$ all of the single site MGF functions $\mathcal{Z}_{i}$
are identical so all averages $\left<\cdot\right>_{i}$ are also identical. At
zero source, the kernels have the simple expressions such as
$\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<\phi(h^{\ell}(\bm{x},t))\phi(h^{\ell}(\bm{x}^{\prime},s))\right>$.

####  C.3.1 Hubbard-Stratonovich transform

Now that we know the kernels $\Phi^{\ell},G^{\ell}$ take on definite values at
$N\to\infty$, we can perform a Hubbard-Stratonovich transformation to
introduce Gaussian sources

$\displaystyle\exp\left(-\frac{1}{2}\int d\bm{x}d\bm{x}^{\prime}dtds\
\hat{\chi}^{\ell}(\bm{x},t)\hat{\chi}^{\ell}(\bm{x}^{\prime},t^{\prime})\Phi^{\ell-1}(\bm{x},\bm{x}^{\prime};t,s)\right)$
$\displaystyle=\left<\exp\left(-i\int d\bm{x}dt\
\hat{\chi}^{\ell}(\bm{x},t)u^{\ell}(\bm{x},t)\right)\right>_{u^{\ell}\sim\mathcal{GP}(0,\Phi^{\ell-1})}$
(31) $\displaystyle\exp\left(-\frac{1}{2}\int d\bm{x}d\bm{x}^{\prime}dtds\
\hat{\xi}^{\ell}(\bm{x},t)\hat{\xi}^{\ell}(\bm{x}^{\prime},t^{\prime})G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)\right)$
$\displaystyle=\left<\exp\left(-i\int d\bm{x}dt\
\hat{\xi}^{\ell}(\bm{x},t)r^{\ell}(\bm{x},t)\right)\right>_{r^{\ell}\sim\mathcal{GP}(0,G^{\ell+1})}$
(32)

After introducing these Gaussian sources $u^{\ell},r^{\ell}$, we have
linearized the single site Hamiltonians $\mathcal{H}_{i}^{\ell}$ with respect
to $\hat{\chi}$ and $\hat{\xi}$. Integration over these variables gives us the
following Dirac-Delta function constraints

$\displaystyle\chi^{\ell}(\bm{x},t)$ $\displaystyle=u^{\ell}(\bm{x},t)+\int
d\bm{x}ds\ A^{\ell-1}(\bm{x},\bm{x}^{\prime};t,s)g^{\ell}(\bm{x}^{\prime},s)$
$\displaystyle\xi^{\ell}(\bm{x},t)$ $\displaystyle=r^{\ell}(\bm{x},t)+\int
d\bm{x}ds\
B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\phi(h^{\ell}(\bm{x}^{\prime},s))$ (33)

Using similar manipulations, we can simplify the expressions for
$A^{\ell},B^{\ell},\hat{\Phi},\hat{G}$ at zero source $\bm{j}=0$

$\displaystyle A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=\left<\frac{\delta\phi(h^{\ell}(\bm{x},t))}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>$ $\displaystyle
B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$ $\displaystyle=\left<\frac{\delta
g^{\ell+1}(\bm{x},t)}{\delta u^{\ell}(\bm{x}^{\prime},s)}\right>$
$\displaystyle\hat{\Phi}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$ $\displaystyle=0\
,\ \hat{G}^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=0$ (34)

Plugging the dynamics for $\chi,\xi$ into the formulas for
$h^{\ell},g^{\ell}$, we obtain the following final equations for the
stochastic process of $h^{\ell},g^{\ell}$.

$\displaystyle h^{\ell+1}(\bm{x},t)$
$\displaystyle=h^{\ell}(\bm{x},t)+\frac{1}{\sqrt{L}}u^{\ell+1}(\bm{x},t)+\frac{1}{\sqrt{L}}\int
d\bm{x}^{\prime}ds\
A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)g^{\ell+1}(\bm{x}^{\prime},s)$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\mathbb{E}_{\bm{x}^{\prime}}\
\int_{0}^{t}ds\
\Delta(\bm{x}^{\prime},s)\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\
g^{\ell+1}(\bm{x}^{\prime},s)\ ,\ u^{\ell}\sim\mathcal{GP}(0,\Phi^{\ell-1})$
$\displaystyle g^{\ell}(\bm{x},t)$
$\displaystyle=g^{\ell+1}(\bm{x},t)+\frac{1}{\sqrt{L}}\dot{\phi}(h^{\ell}(\bm{x},t))z^{\ell}(\bm{x},t)\
,\ $ $\displaystyle z^{\ell}(\bm{x},t)$ $\displaystyle=r^{\ell}(\bm{x},t)+\int
d\bm{x}^{\prime}ds\
B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\phi(h^{\ell}(\bm{x}^{\prime},s))\ ,\
r^{\ell}\sim\mathcal{GP}(0,G^{\ell+1})$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{\sqrt{L}}\mathbb{E}_{\bm{x}^{\prime}}\int_{0}^{t}ds\
\Delta(\bm{x}^{\prime},s)G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)\phi(h^{\ell}(\bm{x}^{\prime},s))$
(35a)

where the order parameters satisfy the the following saddle point equations

$\displaystyle\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<\phi(h^{\ell}(\bm{x},t))\phi(h^{\ell}(\bm{x}^{\prime},s))\right>\
,\
G^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<g^{\ell}(\bm{x},t)g^{\ell}(\bm{x}^{\prime},s)\right>$
$\displaystyle
A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<\frac{\delta\phi(h^{\ell}(\bm{x},t))}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>\ ,\
B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<\frac{\delta
g^{\ell+1}(\bm{x},t)}{\delta u^{\ell+1}(\bm{x}^{\prime},s)}\right>$
$\displaystyle\frac{df(\bm{x},t)}{dt}=\eta_{0}\mathbb{E}_{\bm{x}^{\prime}}K(\bm{x},\bm{x}^{\prime},t)\Delta(\bm{x}^{\prime},t)\
,\ \Delta(\bm{x},t)=-\frac{\partial\mathcal{L}}{\partial f(\bm{x},t)}$
$\displaystyle
K(\bm{x},\bm{x}^{\prime},t)=\frac{1}{L}\sum_{\ell=1}^{L}G^{\ell+1}(\bm{x},\bm{x}^{\prime};t,s)\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)$
(36a)

Up to the residual skip connections and the factors of $\frac{1}{\sqrt{L}}$,
this agrees with the $N\to\infty$ limit DMFT derived in Bordelon & Pehlevan
(2022b). At fixed $L$ this above equation can be solved with a self-consistent
Monte Carlo sampling procedure. For deep linear networks, the equations close
since all $h,g$ fields are Gaussian. In the next section, we analyze the
behavior of the above equations at large depth $L$ so that we can take a large
depth limit.

##  Appendix D The $L\to\infty$ limit from the DMFT Saddle Point Equations

In this section, we consider how the infinite width dynamics behaves at large
depth. We will show that the limiting dynamics takes the form of a system of
ODEs for the kernels in terms of layer time $\tau=\ell/L$. First, we will show
that the response functions need to be rescaled with respect to $L$. Then we
will take a continuum limit over layers to arrive at the final ODEs for the
$\Phi$ and $G$ correlation functions.

###  D.1 Unfolding the Layer Recursion

A useful first step which we will use in the following analysis is “unfolding”
the recurrence to get explicit formulas for $h^{\ell}$ in terms of all of the
Gaussian sources $u$ and all of the gradient fields $g^{\ell}$.

$\displaystyle h^{\ell}(\bm{x},t)=$ $\displaystyle
h^{1}(\bm{x},t)+\frac{1}{\sqrt{L}}\sum_{k=0}^{\ell}u^{k}(\bm{x},t)+\frac{1}{\sqrt{L}}\sum_{k=1}^{\ell-1}\int
d\bm{x}^{\prime}ds\
A^{k}(\bm{x},\bm{x}^{\prime};t,s)g^{k+1}(\bm{x}^{\prime},s)$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\mathbb{E}_{\bm{x}^{\prime}}\int_{0}^{t}ds\
\Delta(\bm{x}^{\prime},s)\sum_{k=0}^{\ell-1}\Phi^{k}(\bm{x},\bm{x}^{\prime},t,s)g^{k+1}(\bm{x}^{\prime},s)$
(37)

and similarly for the backward pass, we have

$\displaystyle g^{\ell}(\bm{x},t)=$ $\displaystyle
g^{L}(\bm{x},t)+\frac{1}{\sqrt{L}}\sum_{k=\ell}^{L}\dot{\phi}(h^{k}(\bm{x},t))r^{k}(\bm{x},t)$
$\displaystyle+\frac{1}{\sqrt{L}}\sum_{k=\ell}^{L}\dot{\phi}(h^{k}(\bm{x},t))\int
d\bm{x}^{\prime}ds\
B^{k}(\bm{x},\bm{x}^{\prime};t,s)\phi(h^{k}(\bm{x}^{\prime},s))$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=\ell}^{L}\dot{\phi}(h^{k}(\bm{x},t))\mathbb{E}_{\bm{x}^{\prime}}\int_{0}^{t}ds\
\Delta(\bm{x}^{\prime},s)\sum_{k=1}^{\ell-1}G^{k}(\bm{x},\bm{x}^{\prime},t,s)\phi(h^{k}(\bm{x}^{\prime},s))$
(38)

A very useful fact that can be immediately gleaned from the above equations is
that

$\displaystyle\forall k\leq\ell\ ,\ \frac{\delta h^{\ell}(\bm{x},t)}{\delta
u^{k}(\bm{x},t)}\sim\mathcal{O}\left(\frac{1}{\sqrt{L}}\right)$
$\displaystyle\forall k\geq\ell\ ,\ \frac{\delta g^{\ell}(\bm{x},t)}{\delta
r^{k}(\bm{x},t)}\sim\mathcal{O}\left(\frac{1}{\sqrt{L}}\right)$ (39)

These facts will be utilized to obtain the scale of the response functions in
the large depth limit.

###  D.2 Characterizing Response Functions at Large Depth

First, we will study the response functions $A^{\ell},B^{\ell}$. We note that
from the saddle point equations,

$\displaystyle
A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=\left<\frac{\delta\phi(h^{\ell}(\bm{x},t^{\prime}))}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>.$ (40)

Since $\frac{\delta\phi(h)}{\delta r}=\dot{\phi}(h)\frac{\delta h}{\delta r}$,
it suffices to consider the scale of the derivative $\frac{\delta
h^{\ell}(\bm{x},t^{\prime})}{\delta r(\bm{x}^{\prime},s)}$

$\displaystyle A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)=$ $\displaystyle\
\frac{1}{\sqrt{L}}\sum_{k=1}^{\ell}\int
d\bm{x}^{\prime\prime}dt^{\prime}A^{k}(\bm{x},\bm{x}^{\prime\prime};t,t^{\prime})\left<\dot{\phi}(h^{\ell}(\bm{x},t))\frac{\delta
g^{k+1}(\bm{x}^{\prime\prime},t^{\prime})}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=0}^{\ell-1}\mathbb{E}_{\bm{x}^{\prime}}\int_{0}^{t}dt^{\prime}\Delta(\bm{x}^{\prime},s)\Phi^{k}(\bm{x},\bm{x}^{\prime};t,s)\
\left<\dot{\phi}(h^{\ell}(\bm{x},t))\frac{\delta
g^{k}(\bm{x}^{\prime\prime},t^{\prime})}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>$ (41)

Next, we use the fact that $\frac{\delta
g^{k}(\bm{x}^{\prime\prime},t^{\prime})}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\sim\mathcal{O}(L^{-1/2})$. Therefore, the second
sum (the one which depends on $\eta_{0}\gamma_{0}L^{-1}$) scales as
$\mathcal{O}(\eta_{0}\gamma_{0}L^{-1/2})$. Solving the $A^{\ell}$ recursion
iteratively from the first layer, we arrive at the conclusion that

$\displaystyle
A^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\sim\mathcal{O}(\eta_{0}\gamma_{0}L^{-1/2}).$
(42)

Indeed, as a sanity check, plugging this in gives a consistent solution to the
above equation. Repeating an identical argument for the backward pass reveals
an identical scaling for $B^{\ell}$ with depth $L$

$\displaystyle
B^{\ell}(\bm{x},\bm{x}^{\prime};t,s)\sim\mathcal{O}(\eta_{0}\gamma_{0}L^{-1/2})$
(43)

We therefore need to introduce a rescaled version of the response functions

$\displaystyle
A^{\ell}(\bm{x},\bm{x}^{\prime},t,s)\leftarrow\frac{\sqrt{L}}{\eta_{0}\gamma_{0}}\
A^{\ell}(\bm{x},\bm{x}^{\prime},t,s)$ (44)
$\displaystyle{B}^{\ell}(\bm{x},\bm{x}^{\prime},t,s)\leftarrow\frac{\sqrt{L}}{\eta_{0}\gamma_{0}}\
B^{\ell}(\bm{x},\bm{x}^{\prime},t,s)$ (45)

After this rescaling, we have

$\displaystyle
A^{\ell}(\bm{x},\bm{x}^{\prime},t,s)=\frac{\sqrt{L}}{\eta_{0}\gamma_{0}}\left<\frac{\delta\phi(h^{\ell}(\bm{x},t))}{\delta
r^{\ell}(\bm{x}^{\prime},s)}\right>$ (46)

To calculate this rescaled response function, we can derive a closed set of
equations for all of the cross layer correlators

$\displaystyle\frac{\delta h^{\ell}(\bm{x},t)}{\delta
r^{\ell^{\prime}}(\bm{x}^{\prime},s)}$
$\displaystyle=\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=1}^{\ell}\int
d\bm{x}^{\prime\prime}ds\
C_{h}^{k}(\bm{x},\bm{x}^{\prime\prime},t,s)\frac{\delta
g^{k}(\bm{x}^{\prime\prime},s)}{\delta r^{\ell^{\prime}}(\bm{x}^{\prime},s)}$
$\displaystyle\frac{\delta g^{\ell}(\bm{x},t)}{\delta
r^{\ell^{\prime}}(\bm{x}^{\prime},s)}$
$\displaystyle=\frac{1}{\sqrt{L}}\Theta(\ell^{\prime}-\ell)\delta(\bm{x}-\bm{x}^{\prime})\delta(t-s)\dot{\phi}(h^{\ell^{\prime}}(\bm{x},t))$
$\displaystyle+\frac{1}{\sqrt{L}}\sum_{k=\ell}^{L}r^{k}(\bm{x},t)\frac{\delta\dot{\phi}(h^{k}(\bm{x},t))}{\delta
r^{\ell^{\prime}}(\bm{x}^{\prime},s)}$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=\ell}^{L}\dot{\phi}(h^{k}(\bm{x},t))\int
d\bm{x}^{\prime\prime}dsC_{g}^{k}(\bm{x},\bm{x}^{\prime\prime},t,s)\frac{\delta\phi(h^{k}(\bm{x}^{\prime\prime},s))}{\delta
r^{\ell^{\prime}}(\bm{x}^{\prime},s)}$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=\ell}^{L}\frac{\delta\dot{\phi}(h^{k}(\bm{x},t))}{\delta
r^{\ell^{\prime}}(\bm{x}^{\prime},s)}\int
d\bm{x}^{\prime\prime}dsC_{g}^{k}(\bm{x},\bm{x}^{\prime\prime},t,s)\phi(h^{k}(\bm{x}^{\prime\prime},s))$
(47)

where we introduced the shorthand

$\displaystyle C_{h}^{\ell}(\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle=A^{\ell-1}(\bm{x},\bm{x}^{\prime},t,s)+\Phi^{\ell-1}(\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime})$
$\displaystyle C_{h}^{\ell}(\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle=B^{\ell}(\bm{x},\bm{x}^{\prime},t,s)+G^{\ell+1}(\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime})$
(48)

Similar equations exist for all derivatives with respect to
$u^{\ell}(\bm{x},t)$. Once the above equations are solved, then one can
compute the response functions $A^{\ell},B^{\ell}$ by averaging the necessary
field derivatives.

###  D.3 Continuum Layer-time Limit at Infinite Depth

In this section, we argue the infinite depth limit can be viewed as a
continuum limiting process defined by ODEs for the kernels and SDEs for the
preactivation and gradient fields for a layertime $\tau=\ell/L$. We start by
examining the asymptotic structure of the preactivation $h^{\ell}$ dynamics

$\displaystyle h^{\ell}(\bm{x};t)$ $\displaystyle\sim
h^{1}(\bm{x};t)+\frac{1}{\sqrt{L}}\sum_{k=1}^{\ell}u^{k}(\bm{x};t)$
$\displaystyle+\frac{\eta_{0}\gamma_{0}}{L}\sum_{k=0}^{\ell-1}\int_{0}^{t}ds\
\mathbb{E}_{\bm{x}^{\prime}}C_{h}^{k}(\bm{x},\bm{x}^{\prime};t,s)\
g^{k+1}(\bm{x}^{\prime},s)\Delta(\bm{x}^{\prime};s)$ (49)

The Gaussian source variables $u^{k}$ are all uncorrelated across layers. So
the first sum is also a mean zero Gaussian with variance

$\displaystyle\left<\left(\frac{1}{\sqrt{L}}\sum_{k=1}^{\ell}u^{k}(\bm{x},t)\right)\left(\frac{1}{\sqrt{L}}\sum_{k=1}^{\ell}u^{k}(\bm{x}^{\prime},s)\right)\right>=\frac{1}{L}\sum_{k=1}^{\ell}\Phi^{k-1}(\bm{x},\bm{x}^{\prime};t,s)\sim\mathcal{O}(1)$
(50)

At large depth $L$, this term will behave as integrated Brownian motion since
each layer’s contribution is independent and order
$\sqrt{d\tau}\sim\frac{1}{\sqrt{L}}$. Next, we reason about the scale of
feature learning update

$\displaystyle\frac{\eta_{0}^{2}\gamma_{0}^{2}}{L^{2}}$
$\displaystyle\sum_{k=0}^{\ell-1}\sum_{k^{\prime}=0}^{\ell-1}\int_{0}^{t}ds^{\prime}\int_{0}^{t}ds\
\mathbb{E}_{\bm{x}^{\prime},\bm{x}^{\prime\prime}}\
\Delta(\bm{x}^{\prime};s)\Delta(\bm{x}^{\prime\prime};s^{\prime})$
$\displaystyle\times C_{h}^{k}(\bm{x},\bm{x}^{\prime};t,s)\
C_{h}^{k^{\prime}}(\bm{x},\bm{x}^{\prime};t,s)\
\left<g^{k+1}(\bm{x}^{\prime},s)g^{k^{\prime}+1}(\bm{x}^{\prime\prime},s^{\prime})\right>\sim\mathcal{O}(1)$
(51)

The above double sum is $\mathcal{O}(1)$ because we have
$\left<g^{k}g^{k^{\prime}}\right>\sim\mathcal{O}(1)$ for all $k,k^{\prime}$.
This fact arises from the skip connections, which generate long-range
correlations across layers. This behavior is very different from non-residual
networks where $g^{\ell}$ and $g^{\ell^{\prime}}$ are independent for all
$\ell\neq\ell^{\prime}$ (Bordelon & Pehlevan, 2022b).

To take the continuum limit, we define functions of a continuous layer time
variable $\tau\in[0,1]$,

$\displaystyle h(\tau,\bm{x},t)$ $\displaystyle=h^{\ell}(\bm{x},t)|_{\ell=\tau
L}\ ,\ g(\tau,\bm{x},t)=g^{\ell}(\bm{x},t)|_{\ell=\tau L}$
$\displaystyle\Phi(\tau;\bm{x},\bm{x}^{\prime};t,s)$
$\displaystyle=\Phi^{\ell}(\bm{x},\bm{x}^{\prime};t,s)|_{\ell=\tau L}\ ,\
G(\tau;\bm{x},\bm{x}^{\prime};t,s)=G^{\ell}(\bm{x},\bm{x}^{\prime};t,s)|_{\ell=\tau
L}.$ (52)

The key idea when taking the continuum limit is the following relations based
on a notion an infinitesimal layer increment $d\tau\sim\frac{1}{L}$. Let
$F^{\ell}$ be a deterministic function of the layer index and let $F(\tau)$ be
its continuous analogue which satisfies $F(\tau)=F^{\ell}|_{\ell=\tau L}$.
Then,

$\displaystyle\lim_{L\to\infty}\frac{1}{L}\sum_{\ell=1}^{\tau
L}F^{\ell}\to\int_{0}^{\tau}d\tau^{\prime}F(\tau)$
$\displaystyle\lim_{L\to\infty}L\left[F^{\tau L+1}-F^{\tau
L}\right]\to\frac{\partial}{\partial\tau}F(\tau).$ (53)

We will use these relations when deriving continuum limits of the DMFT
process. We note that standard error analysis for Euler discretizations can be
used to argue that finite depth averages converge to the $\tau$ integral at
rate

$\displaystyle\left|\frac{1}{L}\sum_{\ell=1}^{\tau
L}F^{\ell}-\int_{0}^{\tau}d\tau^{\prime}F(\tau^{\prime})\right|\sim\mathcal{O}(L^{-1})$
(54)

This is merely the approximation error of a Riemann sum which is proportional
to the step size (Rudin, 1953), which in this case is $1/L$. This idea can be
used to generate a convergence rate of the NTK and the network output
predictions at large depth as we discuss in H.

###  D.4 Preactivation Dynamics

As discussed in the main text, the large $L$ limit for the dynamics of
$h^{\ell}$ can be viewed as a stochastic integral over $\tau$ and over
Brownian motion for the $u^{k}$ sum

$\displaystyle h(\tau;\bm{x},t)$ $\displaystyle=\ h(0;\bm{x},t)+\
\int_{0}^{\tau}du(\tau^{\prime},\bm{x},t)$
$\displaystyle+\eta_{0}\gamma_{0}\int
d\bm{x}^{\prime}\int_{0}^{t}ds\int_{0}^{\tau}d\tau^{\prime}\left[\Phi(\tau^{\prime};\bm{x},\bm{x}^{\prime};t,s)\Delta(\bm{x}^{\prime};s)p(\bm{x}^{\prime})+A(\tau^{\prime},\bm{x},\bm{x}^{\prime};t,s)\right]g(\tau^{\prime},\bm{x}^{\prime},s)$
(55)

where $du$ is Brownian motion term with covariance structure

$\displaystyle\left<du(\tau,\bm{x},t)du(\tau^{\prime},\bm{x}^{\prime},s)\right>=\delta(\tau-\tau^{\prime})\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)d\tau.$
(56)

Similarly, the gradient fields have an analogous structure

$\displaystyle g(\tau;\bm{x},t)=$ $\displaystyle\
g(\tau;\bm{x},t)+\int_{\tau}^{1}dr(\tau^{\prime},\bm{x},t)\dot{\phi}(h(\tau^{\prime},\bm{x},t))$
$\displaystyle+\eta_{0}\gamma_{0}\mathbb{E}_{\bm{x}^{\prime}}\int_{0}^{t}ds\
\int_{\tau}^{1}d\tau^{\prime}\dot{\phi}(h(\tau^{\prime},\bm{x},t))G(\tau^{\prime};\bm{x},\bm{x}^{\prime};t,s)\Delta(\bm{x}^{\prime};s)\phi(h(\tau^{\prime};\bm{x}^{\prime},s))$
$\displaystyle+\eta_{0}\gamma_{0}\int d\bm{x}^{\prime}\int_{0}^{t}ds\
\int_{\tau}^{1}d\tau^{\prime}\dot{\phi}(h(\tau^{\prime},\bm{x},t))B(\tau^{\prime},\bm{x},\bm{x}^{\prime};t,s)\phi(h(\tau^{\prime};\bm{x}^{\prime},s)).$
(57)

where $dr$ is Brownian motion with covariance structure

$\displaystyle\left<dr(\tau,\bm{x},t)dr(\tau^{\prime},\bm{x}^{\prime},s)\right>=\delta(\tau-\tau^{\prime})G(\tau,\bm{x},\bm{x}^{\prime},t,s)d\tau.$
(58)

We note that, due to the nonlinearities, the $h,g$ fields are non-Gaussian at
finite $\gamma_{0}$ as in the case of mean field networks Bordelon & Pehlevan
(2022b). We note that an alternative to the above integral formulation is a
stochastic differential equations

$\displaystyle dh(\tau,\bm{x},t)$ $\displaystyle=du(\tau,\bm{x},t)$
$\displaystyle+d\tau\ \eta_{0}\gamma_{0}\int d\bm{x}^{\prime}\int_{0}^{t}ds\
[\Phi(\tau;\bm{x},\bm{x}^{\prime};t,s)\Delta(\bm{x}^{\prime};s)p(\bm{x}^{\prime})+A(\tau;\bm{x},\bm{x}^{\prime};t,s)]g(\tau,\bm{x}^{\prime},s)$
(59)

The $\gamma_{0}\to 0$ limit of the $h$ equation at $t,s=0$ coincides with the
SDE derived at initialization by Hayou (2023); Hayou & Yang (2023). However,
we now see that feature learning term (the term which depends on $\gamma_{0}$)
adds a drift term to the layerwise dynamics compared to the diffusion term
from the initial weights.

###  D.5 Kernel Layerwise ODEs in $\tau$ in Feature Learning Regime

We can start from the stochastic difference equation

$\displaystyle dh(\tau,\bm{x},t)$
$\displaystyle=du(\tau,\bm{x}^{\prime},t)+d\tau\ \eta_{0}\gamma_{0}\int
d\bm{x}^{\prime}\int_{0}^{t}ds\
C_{h}(\tau,\bm{x},\bm{x}^{\prime},t,s)g(\tau,\bm{x}^{\prime},s)$
$\displaystyle C_{h}(\tau,\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle=A(\tau,\bm{x},\bm{x}^{\prime},t,s)+\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime}).$
(60)

and derive an ODE for the feature kernel at layer-time $\tau$. Applying Ito’s
lemma (Itô, 1944), we obtain the following feature kernel dynamics

$\displaystyle\frac{\partial}{\partial\tau}\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)=$
$\displaystyle\frac{1}{2}\int
d\tilde{\bm{x}}d\tilde{\bm{x}}^{\prime}dt^{\prime}ds^{\prime}\Phi(\tau,\tilde{\bm{x}},\tilde{\bm{x}}^{\prime},t^{\prime},s^{\prime})\left<\frac{\partial^{2}}{\partial
h(\tau,\tilde{\bm{x}},t^{\prime})\partial
h(\tau,\tilde{\bm{x}}^{\prime},s^{\prime})}[\phi(h(\tau,\bm{x},t))\phi(h(\tau,\bm{x}^{\prime},s))]\right>$
$\displaystyle+\eta_{0}\gamma_{0}\int
d\bm{x}^{\prime\prime}dt^{\prime}C_{h}(\tau,\bm{x},\bm{x}^{\prime\prime},t^{\prime})\left<g(\tau,\bm{x}^{\prime\prime},t^{\prime})\dot{\phi}(h(\tau,\bm{x},t))\phi(h(\tau,\bm{x}^{\prime},s))\right>$
$\displaystyle+\eta_{0}\gamma_{0}\int
d\bm{x}^{\prime\prime}ds^{\prime}C_{h}(\tau,\bm{x}^{\prime},\bm{x}^{\prime\prime},s^{\prime})\left<g(\tau,\bm{x}^{\prime\prime},s^{\prime})\phi(h(\tau,\bm{x},t))\dot{\phi}(h(\tau,\bm{x}^{\prime},s))\right>$
(61)

The first term on the right hand side comes from signal propagation, which
persists even in the lazy limit. The second set of terms come from feature
learning corrections. This equation should be integrated from $\tau=0$ to
$\tau=1$. Similarly, there is a backward pass ODE for the $G$ kernel

$\displaystyle\frac{\partial}{\partial\tau}G(\tau,\bm{x},\bm{x}^{\prime},t,s)=-\left<\dot{\phi}(h(\tau,\bm{x},t))\dot{\phi}(h(\tau,\bm{x}^{\prime},s))\right>G(\tau,\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle-\eta_{0}\gamma_{0}\int
d\bm{x}^{\prime\prime}\int_{0}^{t^{\prime}}dt^{\prime}C_{g}(\tau,\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\left<\dot{\phi}(h(\tau,\bm{x},t))\phi(h(\tau,\bm{x}^{\prime},s))g(\tau,\bm{x}^{\prime\prime},t^{\prime})\right>$
$\displaystyle-\eta_{0}\gamma_{0}\int
d\bm{x}^{\prime\prime}\int_{0}^{s^{\prime}}ds^{\prime}C_{g}(\tau,\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\left<\phi(h(\tau,\bm{x},t))\dot{\phi}(h(\tau,\bm{x}^{\prime},s))g(\tau,\bm{x}^{\prime\prime},s^{\prime})\right>$
$\displaystyle
C_{g}(\tau,\bm{x},\bm{x}^{\prime},t,s)=B(\tau,\bm{x},\bm{x}^{\prime},t,s)+G(\tau,\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime})$
(62)

which should be integrated from $\tau=1$ to $\tau=0$.

We note again, that convergence to the ODE will occur at rate
$\mathcal{O}(L^{-1})$ in absolute error.

###  D.6 Response Function Equations in the Continuum Limit

The response functions can also be determined in the $L\to\infty$ limit using
the following relations

$\displaystyle
A(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{1}{\eta_{0}\gamma_{0}}\left<\frac{\delta\phi(h(\tau,\bm{x},t))}{\delta\
dr(\tau,\bm{x}^{\prime},s)}\right>\ ,\
B(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{1}{\eta_{0}\gamma_{0}}\left<\frac{\delta
g(\tau,\bm{x},t)}{\delta\ du(\tau,\bm{x}^{\prime},s)}\right>$ (63)

The reason that the variation is computed with respect to the differential
Brownian motion $du,dr$ is to match the scale of the increment
$\frac{1}{\sqrt{L}}$, which appears in the finite $L$ response functions
$A^{\ell}=\frac{\sqrt{L}}{\eta_{0}\gamma_{0}}\left<\frac{\delta\phi(h^{\ell})}{\delta
r^{\ell}}\right>$. Using our continuum limit equations, we obtain the
following field derivatives

$\displaystyle\frac{\delta h(\tau,\bm{x},t)}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}=$ $\displaystyle\
\Theta(\tau-\tau^{\prime})\delta(\bm{x}-\bm{x}^{\prime})\delta(t-s)$
$\displaystyle+\eta_{0}\gamma_{0}\int d\tau^{\prime\prime}\int
d\bm{x}^{\prime\prime}\int
dt^{\prime}C_{h}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\frac{\delta
g(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
g(\tau,\bm{x},t)}{\delta du(\tau^{\prime},\bm{x}^{\prime},s)}=$
$\displaystyle\ \eta_{0}\gamma_{0}\int
d\tau^{\prime\prime}\dot{\phi}(h(\tau^{\prime\prime},\bm{x},t))\int
d\bm{x}^{\prime\prime}\int
dt^{\prime}C_{g}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\frac{\delta\phi(h(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime}))}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
h(\tau,\bm{x},t)}{\delta dr(\tau^{\prime},\bm{x}^{\prime},s)}=$
$\displaystyle\ \eta_{0}\gamma_{0}\int d\tau^{\prime\prime}\int
d\bm{x}^{\prime\prime}\int
dt^{\prime}C_{h}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\frac{\delta
g(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
dr(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
g(\tau,\bm{x},t)}{\delta dr(\tau^{\prime},\bm{x}^{\prime},s)}=$
$\displaystyle\
\Theta(\tau^{\prime}-\tau)\delta(\bm{x}-\bm{x}^{\prime})\delta(t-s)\
\dot{\phi}(h(\tau,\bm{x},t))$ $\displaystyle+\eta_{0}\gamma_{0}\int
d\tau^{\prime\prime}\dot{\phi}(h(\tau^{\prime\prime},\bm{x},t))\int
d\bm{x}^{\prime\prime}\int
dt^{\prime}C_{h}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\frac{\delta\phi(h(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime}))}{\delta
dr(\tau^{\prime},\bm{x}^{\prime},s)}$ (64)

Once these recursions are solved, the response functions can be computed by
averaging the necessary field derivatives at each layer time $\tau$

$\displaystyle
A(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<\frac{\delta\phi(h(\tau,\bm{x},t))}{\delta
dr(\tau,\bm{x}^{\prime},s)}\right>\ ,\
B(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<\frac{\delta g(\tau,\bm{x},t)}{\delta
du(\tau,\bm{x}^{\prime},s)}\right>.$ (65)

###  D.7 Final Result: Infinite Width and Depth Limit

In this section we combine the results obtained in the preceeding sections to
state our main theoretical result which is the infinite width and depth limit
of the learning dynamics in our ResNet. Below we state the result for gradient
flow, but minor modifications can also characterize discrete time SGD,
momentum, weight decay etc as we discuss in Appendices J.

###### Proposition 3

Consider the $N,L\to\infty$ limit of the ResNet $\mu$P-$\frac{1}{\sqrt{L}}$
model presented in the main text trained with gradient flow on data
distribution $p(\bm{x})$. The preactivation $h(\tau,\bm{x},t)$ and gradient
$g(\tau,\bm{x},t)$ stochastic fields obey the equations

$\displaystyle h(\tau,\bm{x},t)$
$\displaystyle=h(0,\bm{x},t)+\int_{0}^{\tau}du(\tau^{\prime},\bm{x},t)+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime}\int
d\bm{x}^{\prime}\int_{0}^{t}ds\
C_{h}(\tau^{\prime},\bm{x},\bm{x}^{\prime},t,s)g(\tau^{\prime},\bm{x}^{\prime},s)$
$\displaystyle g(\tau,\bm{x},t)$
$\displaystyle=g(1,\bm{x},t)+\int_{\tau}^{1}\dot{\phi}(h(\tau^{\prime},\bm{x},t))dr(\tau^{\prime},\bm{x},t)$
$\displaystyle+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime}\dot{\phi}(h(\tau^{\prime},\bm{x},t))\int
d\bm{x}^{\prime}\int_{0}^{t}ds\
C_{g}(\tau^{\prime},\bm{x},\bm{x}^{\prime},t,s)\phi(h(\tau^{\prime},\bm{x}^{\prime},s))$
(66)

where the zero mean Brownian motions $du(\tau,\bm{x},t),dr(\tau,\bm{x},t)$
have the following covariance structure

$\displaystyle\left<du(\tau,\bm{x},t)du(\tau^{\prime},\bm{x}^{\prime},s)\right>=\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)\delta(\tau-\tau^{\prime})d\tau$
$\displaystyle\left<dr(\tau,\bm{x},t)dr(\tau^{\prime},\bm{x}^{\prime},s)\right>=G(\tau,\bm{x},\bm{x}^{\prime},t,s)\delta(\tau-\tau^{\prime})d\tau,$
(67)

the deterministic kernels $C_{h}$ and $C_{g}$ have the form

$\displaystyle C_{h}(\tau,\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle=\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime})+A(\tau,\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle C_{g}(\tau,\bm{x},\bm{x}^{\prime},t,s)$
$\displaystyle=G(\tau,\bm{x},\bm{x}^{\prime},t,s)\Delta(\bm{x}^{\prime},s)p(\bm{x}^{\prime})+B(\tau,\bm{x},\bm{x}^{\prime},t,s)$
(68)

where the kernels $\Phi,G$ and response functions $A,B$ are computed as
averages over the $h,g$ stochastic process

$\displaystyle\Phi(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<\phi(h(\tau,\bm{x},t))\phi(h(\tau,\bm{x}^{\prime},s))\right>\
,\
G(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<g(\tau,\bm{x},t)g(\tau,\bm{x}^{\prime},s)\right>$
$\displaystyle
A(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{1}{\eta_{0}\gamma_{0}}\left<\frac{\delta\phi(h(\tau,\bm{x},t))}{\delta
dr(\tau,\bm{x}^{\prime},s)}\right>\ ,\
B(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{1}{\eta_{0}\gamma_{0}}\left<\frac{\delta
g(\tau,\bm{x},t)}{\delta du(\tau,\bm{x}^{\prime},s)}\right>.$ (69)

From the feature and gradient kernels, we can compute the dynamical NTK and
the dynamics of network predictions

$\displaystyle K(\bm{x},\bm{x}^{\prime},t)$ $\displaystyle=\int_{0}^{1}d\tau\
G(\tau,\bm{x},\bm{x}^{\prime},t,t)\Phi(\tau,\bm{x},\bm{x}^{\prime},t,t)$
$\displaystyle\frac{d}{dt}f(\bm{x},t)$ $\displaystyle=\eta_{0}\int
d\bm{x}^{\prime}p(\bm{x}^{\prime})K(\bm{x},\bm{x}^{\prime},t)\Delta(\bm{x}^{\prime},t)\
,\ \Delta(\bm{x},t)=-\frac{\partial\mathcal{L}}{\partial f(\bm{x},t)}.$ (70)

This is a closed system of equations relating preactiation and gradient
marginals, the kernel dynamics and the network prediction dynamics. In this
limit, the kernels, response functions and network outputs are all
deterministic.

##  Appendix E Infinite Depth then Infinite Width Limit

In this section, we compute the limiting process in the other direction, by
directly working with the distribution of order parameters $\bm{q}$ at finite
$N$. We show that at large depth $L$ the DMFT action converges to a limiting
object $\mathcal{S}_{\infty}(\bm{q})$. We the show that the large $N$ limit of
the distribution recovers the same saddle point derived in Appendix D for the
large $N$ then large $L$ limit. Thus, this section constitutes a theoretical
argument that the width and depth limits commute for the whole course of
training dynamics. The derivations are performed at the level of rigor of
physics, and we leave open for future work a rigorous result about
commutativity of limits throughout training.

###  E.1 Large $L$ Limit at Fixed $N$

We derived the finite $N,L$ distribution of order parameters in Appendix
Section C. The result was that the log-density of the order parameters
$\bm{q}$ over random initialization of network weights was given by the DMFT
action $\ln p(\bm{q})\propto N\mathcal{S}_{L}(\bm{q})$. Suppressing
unnecessary (for the purposes of this section) sample and time indices, the
action at depth $\mathcal{S}_{L}$ had the form

$\displaystyle\mathcal{S}_{L}=$
$\displaystyle\sum_{\ell}\left[\Phi^{\ell}\hat{\Phi}^{\ell}+G^{\ell}\hat{G}^{\ell}-A^{\ell}B^{\ell}\right]+\mathcal{Z}[\bm{q}]$
(71)

where the single site density $\mathcal{Z}$ has the form

$\displaystyle\mathcal{Z}=$
$\displaystyle\int\prod_{\ell}dh^{\ell}d\hat{h}^{\ell}dg^{\ell}d\hat{g}^{\ell}$
$\displaystyle\mathbb{E}_{\\{u^{\ell},r^{\ell}\\}}\exp\left(i\sum_{\ell}\hat{h}^{\ell}\left(h^{\ell+1}-h^{\ell}-\frac{1}{\sqrt{L}}u^{\ell}-\frac{1}{\sqrt{L}}A^{\ell}g^{\ell+1}-\frac{\eta\gamma}{L}\Phi^{\ell}\Delta
g^{\ell+1}\right)\right)$
$\displaystyle\exp\left(i\sum_{\ell}\hat{g}^{\ell}\left(g^{\ell}-g^{\ell+1}-\frac{1}{\sqrt{L}}r^{\ell}\dot{\phi}(h^{\ell})-\frac{1}{\sqrt{L}}\dot{\phi}(h^{\ell})B^{\ell}\phi(h^{\ell})-\frac{\eta_{0}\gamma_{0}}{L}\dot{\phi}(h^{\ell})G^{\ell+1}\Delta\phi(h^{\ell})\right)\right)$
$\displaystyle\exp\left(-\sum_{\ell}\phi(h^{\ell})\hat{\Phi}^{\ell}\phi(h^{\ell})-\sum_{\ell}g^{\ell}\hat{G}^{\ell}g^{\ell}\right)$
(72)

where the averages are over $u^{\ell},r^{\ell}$ which are Gaussians with
covariances $\Phi^{\ell-1}$ and $G^{\ell+1}$ respectively. In the above
equation we work with $h,g$ directly rather than working with the intermediate
fields $\chi,\xi$ which only depend on the initialization.

To obtain a proper large depth limit, we rescale the conjugate order
parameters as

$\displaystyle\hat{\Phi}^{\ell}\to\frac{1}{L}\hat{\Phi}^{\ell}\ ,\
\hat{G}^{\ell}\to\frac{1}{L}\hat{G}^{\ell}\ ,\
A^{\ell}\to\frac{\eta_{0}\gamma_{0}}{\sqrt{L}}A^{\ell}\ ,\
B^{\ell}\to\frac{\eta_{0}\gamma_{0}}{\sqrt{L}}B^{\ell}$ (73)

We note that such a rescaling will not alter the distribution of real network
observables $\\{\Phi,G,K,f,...\\}$ (Bordelon & Pehlevan, 2022b; 2023). After
this rescaling, we take $L\to\infty$ which gives a continuum limit of the
resulting layer sums, we obtain the following simplified infinite depth action
in terms of the kernels and features at layer time $\tau\in[0,1]$

$\displaystyle\mathcal{S}_{\infty}$ $\displaystyle=\int
d\tau[\Phi(\tau)\hat{\Phi}(\tau)+G(\tau)\hat{G}(\tau)-\eta_{0}^{2}\gamma_{0}^{2}A(\tau)B(\tau)]+\ln\mathcal{Z}$
$\displaystyle\mathcal{Z}$
$\displaystyle=\int\mathcal{D}\hat{h}\mathcal{D}\hat{g}\mathcal{D}h\mathcal{D}g\
\mathbb{E}_{u(\tau),r(\tau)}\exp\left(i\int\hat{h}(\tau)[dh(\tau)-du(\tau)-\eta_{0}\gamma_{0}C_{h}(\tau)g(\tau)d\tau]\right)$
$\displaystyle\exp\left(i\int\hat{g}(\tau)[-dg(\tau)-dr(\tau)\dot{\phi}(h(\tau))-\eta_{0}\gamma_{0}\dot{\phi}(h(\tau))C_{g}(\tau)\phi(h(\tau))d\tau]\right)$
$\displaystyle\exp\left(-\int
d\tau[\phi(h(\tau))\hat{\Phi}(\tau)\phi(h(\tau))+g(\tau)\hat{G}(\tau)g(\tau)]\right).$
(74)

where $C_{h},C_{g}$ are given in the main text and the $du(\tau),dr(\tau)$ are
Brownian motion processes over layer time $\tau$. We therefore find that the
distribution of order parameters at finite width $N$ and infinite depth
$L\to\infty$ is

$\displaystyle\lim_{L\to\infty}p_{L}(\bm{q})=p_{\infty}(\bm{q})\propto\exp\left(N\mathcal{S}_{\infty}(\bm{q})\right)$
(75)

This constitutes taking the infinite depth limit $L\to\infty$ at fixed $N$. At
this stage, the order parameters do not obey a deterministic set of equations,
but are samples from the above density. In the next section, we study infinite
width of this density and find that we recover the same result that we derived
previously in the other direction.

###  E.2 Saddle Point Equations on Infinite Depth Action

This section considers the large $N$ limit of the infinite depth stochastic
process. In particular, we evaluate the saddle point of $\mathcal{S}_{\infty}$
and show that it recovers the same limit as we obtained by first taking a
saddle point and then taking the large $L$ limit. This shows that the limits
commute.

$\displaystyle\frac{\delta\mathcal{S}_{\infty}}{\delta\hat{\Phi}(\tau)}$
$\displaystyle=\Phi(\tau)-\left<\phi(h(\tau))\phi(h(\tau))\right>=0$
$\displaystyle\frac{\delta\mathcal{S}_{\infty}}{\delta\hat{G}(\tau)}$
$\displaystyle=G(\tau)-\left<g(\tau)g(\tau)\right>=0$
$\displaystyle\frac{\delta\mathcal{S}_{\infty}}{\delta A(\tau)}$
$\displaystyle=-\eta_{0}^{2}\gamma_{0}^{2}B(\tau)-i\eta_{0}\gamma_{0}\left<\hat{h}(\tau)g(\tau)\right>=0$
$\displaystyle\frac{\delta\mathcal{S}_{\infty}}{\delta B(\tau)}$
$\displaystyle=-\eta_{0}^{2}\gamma_{0}^{2}A(\tau)-i\eta_{0}\gamma_{0}\left<\hat{g}(\tau)\dot{\phi}(h(\tau))\phi(h(\tau))\right>=0$
(76)

For the last two expressions, we recognizing that

$\displaystyle\frac{\delta}{\delta
du(\tau^{\prime})}\exp\left(-i\int\hat{h}(\tau)du(\tau)\right)=-i\hat{h}(\tau)\exp\left(-i\int\hat{h}(\tau)du(\tau)\right)$
$\displaystyle\frac{\delta}{\delta
dr(\tau^{\prime})}\exp\left(-i\int\hat{g}(\tau)\dot{\phi}(h(\tau))dr(\tau)\right)=-i\hat{g}(\tau)\dot{\phi}(h(\tau))\exp\left(-i\int\hat{g}(\tau)\dot{\phi}(h(\tau))dr(\tau)\right)$
(77)

We can thus express the saddle point equations in terms of the left hand side
derivatives with respect to Gaussian sources. Performing integration by parts,
we arrive at the usual response function formulas

$\displaystyle
A(\tau)=\eta_{0}^{-1}\gamma_{0}^{-1}\left<\frac{\delta\phi(h(\tau))}{\delta
dr(\tau))}\right>\ ,\ B(\tau)=\eta_{0}^{-1}\gamma_{0}^{-1}\left<\frac{\delta
g(\tau)}{\delta du(\tau)}\right>$ (78)

Lastly, we can integrate over $\hat{h}(\tau)$ and $\hat{g}(\tau)$ which give
Dirac-Delta functions that define the $h$ and $g$ stochastic processes

$\displaystyle dh(\tau)$
$\displaystyle=du(\tau)+\eta_{0}\gamma_{0}C_{h}(\tau)g(\tau)\ ,\
dg(\tau)=-dr(\tau)\dot{\phi}(h(\tau))-\eta_{0}\gamma_{0}\dot{\phi}(h(\tau))C_{g}(\tau)\phi(h(\tau))d\tau$
(79)

The rest of the indices (over samples and training times) can be incorporated
into the above argument and the result is easily verified to match that
derived in Appendix D from the large width then large depth limit. This
exercise thus verifies that the large width and large depth limit continue to
commute throughout training, extending the commutativity at initialization
proven by Hayou & Yang (2023).

##  Appendix F Exactly Solveable Cases

In this section, we give exact solutions to the dynamics which do not require
any non-Gaussian integrals. The two cases where this can be achieved are

  1. 1.

the lazy limit $\gamma_{0}\to 0$ of network training where the features are
frozen

  2. 2.

networks with linear activations, which preserve Gaussianity of $h,g$ at any
$\gamma_{0}$.

We discuss these two limits in great detail below.

###  F.1 Lazy Limit

In the lazy limit $\gamma_{0}\to 0$, the flow for the preactivations
$h(\tau,\bm{x},t)=h(\tau,\bm{x})$ and the gradient fields
$g^{\ell}(\tau,\bm{x},t)=g^{\ell}(\tau,\bm{x})$ are constant over training
time. Similarly, the kernels $\Phi^{\ell},G^{\ell}$ are also constant
throughout the dynamics. We thus have the following forward and backward
dynamics

$\displaystyle h(\tau,\bm{x})$
$\displaystyle=h(0,\bm{x})+\int_{0}^{\tau}du(\tau^{\prime},\bm{x})$ (80)
$\displaystyle g(\tau,\bm{x})$
$\displaystyle=g(1,\bm{x})+\int_{\tau}^{1}\dot{\phi}(h(\tau^{\prime},\bm{x}))dr(\tau^{\prime},\bm{x})$
(81)

We note that the above equations indicate that $h(\tau,\bm{x})$ are all
jointly Gaussian with zero mean. Since $h$ is Gaussian, it makes sense to
characterize its covariance
$H(\tau,\bm{x},\bm{x}^{\prime})=\left<h(\tau,\bm{x})h(\tau,\bm{x}^{\prime})\right>$
at layer-time $\tau$. This object satisfies the differential equation

$\displaystyle\frac{\partial}{\partial\tau}H(\tau,\bm{x},\bm{x}^{\prime})$
$\displaystyle=\Phi(\tau,\bm{x},\bm{x}^{\prime})$
$\displaystyle\Phi(\tau,\bm{x},\bm{x}^{\prime})$
$\displaystyle=\left<\phi(h(\tau,\bm{x}))\phi(h(\tau,\bm{x}^{\prime}))\right>$
$\displaystyle h(\tau,\bm{x},\bm{x}^{\prime})$
$\displaystyle\sim\mathcal{N}(0,\bm{H}(\tau))$ (82)

Since $h(\tau,\bm{x})$ are zero-mean jointly Gaussian random variables, the
average to compute $\Phi$ can be computed in closed form for polynomials and
several practically used nonlinearities (ReLU, ERF, etc) (Cho & Saul, 2009).
Now, to compute the backward pass, we must

$\displaystyle\frac{\partial}{\partial\tau}G(\tau,\bm{x},\bm{x}^{\prime})$
$\displaystyle=-\left<\dot{\phi}(h(\tau,\bm{x}))\dot{\phi}(h(\tau^{\prime},\bm{x}^{\prime}))\right>G(\tau,\bm{x},\bm{x}^{\prime})$
(83)

Again, the bivariate Gaussian integral
$\left<\dot{\phi}(h(\tau,\bm{x}))\dot{\phi}(h(\tau^{\prime},\bm{x}^{\prime}))\right>$
has a closed form expression for many regularly used activation functions
$\phi$. Given that we have solved for $\Phi(\tau)$ and $G(\tau)$, the dynamics
of the predictor is

$\displaystyle\frac{d}{dt}f(\bm{x},t)$
$\displaystyle=\eta_{0}\mathbb{E}_{\bm{x}^{\prime}}K(\bm{x},\bm{x}^{\prime})\Delta(\bm{x}^{\prime},t)\
,\ \Delta(\bm{x},t)=-\frac{\partial\mathcal{L}}{\partial f(\bm{x},t)}$
$\displaystyle K(\bm{x},\bm{x}^{\prime})$ $\displaystyle=\int_{0}^{1}d\tau\
G(\tau,\bm{x},\bm{x}^{\prime})\Phi(\tau,\bm{x},\bm{x}^{\prime}).$ (84)

This is an exact description of the network predictor dynamics in the
$N\to\infty,L\to\infty$ and $\gamma_{0}\to 0$ limit. However, the
$\gamma_{0}\to 0$ limit lacks a fundamental phenomenon of deep neural
networks, namely feature learning. To gain analytical insight into the feature
learning dynamics, we will study linear networks in the next section.

###  F.2 Infinite Depth ResNets with Linear Activations

In the case where the activations are linear, we can also close the equations
at the level of the kernels (since activations remain Gaussian even after
feature learning). To provide an exactly solveable model of the dynamics, we
consider computing the kernels after one step of feature learning.

####  F.2.1 Full DMFT equations in Deep Linear Networks

The DMFT equations for the preactivation processes have the form

$\displaystyle h(\tau,\bm{x},t)$
$\displaystyle=h(0,\bm{x},t)+\int_{0}^{\tau}du(\tau^{\prime},\bm{x},t)+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime}\int
d\bm{x}^{\prime}\int ds\
C_{h}(\tau^{\prime},\bm{x},\bm{x}^{\prime},t,s)g(\tau^{\prime},\bm{x}^{\prime},s)$
$\displaystyle g(\tau,\bm{x},t)$
$\displaystyle=g(1,\bm{x},t)+\int_{\tau}^{1}dr(\tau^{\prime},\bm{x},t)+\eta_{0}\gamma_{0}\int_{\tau}^{1}d\tau^{\prime}\int
d\bm{x}^{\prime}\int ds\
C_{g}(\tau^{\prime},\bm{x},\bm{x}^{\prime},t,s)h(\tau^{\prime},\bm{x}^{\prime},s)$
(85)

where $C_{h},C_{g}$ have the usual expressions. The field derivatives (to get
the response functions) have the form

$\displaystyle\frac{\delta h(\tau,\bm{x},t)}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}$
$\displaystyle=\Theta(\tau-\tau^{\prime})\delta(\bm{x}-\bm{x}^{\prime})\delta(t-s)$
$\displaystyle+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime\prime}\int
dx^{\prime\prime}\int
dt^{\prime}C_{h}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime\prime},t,t^{\prime})\frac{\delta
g(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
g(\tau,\bm{x},t)}{\delta du(\tau^{\prime},\bm{x}^{\prime},s)}$
$\displaystyle=\eta_{0}\gamma_{0}\int_{\tau}^{1}d\tau^{\prime\prime}\int
dx^{\prime\prime}\int
dt^{\prime}C_{g}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime},t,t^{\prime})\frac{\delta
h(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
du(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
h(\tau,\bm{x},t)}{\delta dr(\tau^{\prime},\bm{x}^{\prime},s)}$
$\displaystyle=\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime\prime}\int
dx^{\prime\prime}\int
dt^{\prime}C_{h}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime},t,t^{\prime})\frac{\delta
g(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
dr(\tau^{\prime},\bm{x}^{\prime},s)}$ $\displaystyle\frac{\delta
g(\tau,\bm{x},t)}{\delta dr(\tau^{\prime},\bm{x}^{\prime},s)}$
$\displaystyle=\Theta(\tau^{\prime}-\tau)\delta(\bm{x}-\bm{x}^{\prime})\delta(t-s)$
$\displaystyle+\eta_{0}\gamma_{0}\int_{\tau}^{1}d\tau^{\prime\prime}\int
dx^{\prime\prime}\int
dt^{\prime}C_{g}(\tau^{\prime\prime},\bm{x},\bm{x}^{\prime},t,t^{\prime})\frac{\delta
h(\tau^{\prime\prime},\bm{x}^{\prime\prime},t^{\prime})}{\delta
dr(\tau^{\prime},\bm{x}^{\prime},s)}$ (86)

Since the activations are linear, these derivatives are no longer stochastic
and we have the equal layer time $\tau=\tau^{\prime}$ derivatives as the
response functions

$\displaystyle A(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{\delta
h(\tau,\bm{x},t)}{\delta dr(\tau,\bm{x}^{\prime},s)}\ ,\
B(\tau,\bm{x},\bm{x}^{\prime},t,s)=\frac{\delta g(\tau,\bm{x},t)}{\delta
du(\tau,\bm{x}^{\prime},s)}$ (87)

Next, we compute the kernels as simple Gaussian averages

$\displaystyle
H(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<h(\tau,\bm{x},t)h(\tau,\bm{x},t)\right>\
,\
G(\tau,\bm{x},\bm{x}^{\prime},t,s)=\left<g(\tau,\bm{x},t)g(\tau,\bm{x},t)\right>$
(88)

In the next section, we show a simple example where these equations become
simple linear differential equations.

####  F.2.2 One step One Sample Feature Kernel Equations

Before training, the initial conditions for the kernels $H_{0}(\tau)$ and
$G_{0}(\tau)$ are

$\displaystyle\partial_{\tau}H_{0}(\tau)=H_{0}(\tau)\ ,\
\partial_{\tau}G_{0}(\tau)=-G_{0}(\tau)$ $\displaystyle\implies
H_{0}(\tau)=e^{\tau}\ ,\ G_{0}(\tau)=e^{-\tau}+1-e^{-1}$ (89)

After a step of gradient descent, the preactivation at layer time $\tau$ has
the form

$\displaystyle
h(\tau)=h(0)+\int_{0}^{\tau}du(\tau^{\prime})+\eta_{0}\gamma_{0}\int_{0}^{\tau}d\tau^{\prime}[A(\tau^{\prime})+H_{0}(\tau^{\prime})y]g_{0}(\tau^{\prime})$
(90)

We let $H(\tau),G(\tau)$ represent the kernels after one step. The response
function after one step has satisfies the differential equation

$\displaystyle\partial_{\tau}{A}(\tau)=A(\tau)+H_{0}(\tau)y\ ,\
A(0)=0\implies{A}(\tau)=\tau e^{\tau}\ y$ (91)

The kernel $H(\tau)$ after one step therefore satisfies

$\displaystyle H(\tau)$
$\displaystyle=H(0)+\int_{0}^{\tau}d\tau^{\prime}H(\tau^{\prime})$
$\displaystyle+\eta_{0}^{2}\gamma_{0}^{2}\int_{0}^{\tau}\int_{0}^{\tau}d\tau^{\prime}d\tau^{\prime\prime}G_{0}(\tau^{\prime},\tau^{\prime\prime})[A(\tau^{\prime})+H_{0}(\tau^{\prime})y][A(\tau^{\prime\prime})+H_{0}(\tau^{\prime\prime})y]$
$\displaystyle=H(0)+\int_{0}^{\tau}d\tau^{\prime}H(\tau^{\prime})+\eta_{0}^{2}\gamma_{0}^{2}y^{2}\int_{0}^{\tau}\int_{0}^{\tau}d\tau^{\prime}d\tau^{\prime\prime}G_{0}(\tau^{\prime},\tau^{\prime\prime})e^{\tau^{\prime}+\tau^{\prime\prime}}[\tau^{\prime}+1][\tau^{\prime\prime}+1].$
(92)

where
$G_{0}(\tau,\tau^{\prime})=\left<g_{0}(\tau)g_{0}(\tau^{\prime})\right>$.
Differentiating both sides with respect to $\tau$ gives the integro-
differential equation

$\displaystyle\partial_{\tau}H(\tau)=H(\tau)+2\eta_{0}^{2}\gamma_{0}^{2}y^{2}e^{\tau}[\tau+1]\int_{0}^{\tau}d\tau^{\prime}G_{0}(\tau,\tau^{\prime})e^{\tau^{\prime}}[\tau^{\prime}+1]$
(93)

Note that $G_{0}(\tau,\tau^{\prime})=G_{0}(\tau)=e^{-\tau}+1-e^{-1}$ for
$\tau^{\prime}<\tau$. So we can pull this out of the integral and we are left
with

$\displaystyle\int_{0}^{\tau}d\tau^{\prime}e^{\tau^{\prime}}[\tau^{\prime}+1]=\tau
e^{\tau}$ (94)

We are finally left with the ODE

$\displaystyle\partial_{\tau}H(\tau)=H(\tau)+2\eta_{0}^{2}\gamma_{0}^{2}y^{2}\
e^{2\tau}[e^{-\tau}+1-e^{-1}][\tau^{2}+\tau]$ (95)

Integrating gives

$\displaystyle
H(\tau)=e^{\tau}H(0)+2\eta_{0}^{2}\gamma_{0}^{2}e^{\tau}\left[\frac{1}{3}\tau^{3}+\frac{1}{2}\tau^{2}+(1+e^{-1})(\tau^{2}e^{\tau}-\tau
e^{\tau}+e^{\tau}-1)\right]$ (96)

This solution gives the profile of the feature kernels after a single step of
gradient descent. This is verified numerically in Figure 7. At small $\tau$
this goes as
$H(\tau)-H_{0}(\tau)\sim\eta_{0}^{2}\gamma_{0}^{2}[1+(1+e^{-1})]\tau^{2}$.

##  Appendix G Finite Width Convergence: Fluctuations do not scale with $L$

At fixed $L$, one can compute systematic asymptotic corrections in
$\mathcal{O}_{N}(N^{-1})$ to the statistical properties of $\bm{q}$. This
follows closely from the finite width analysis around mean field dynamics
performed by Bordelon & Pehlevan (2023). The idea of their expansion is to
extract from $\mathcal{S}$ information about finite size deviations from the
infinite width DMFT process. The key fact which appears in the next to leading
order terms is that at finite width $N$, the order parameters $\bm{q}$ become
random with covariance

$\displaystyle\text{Cov}(\bm{q})\sim\frac{1}{N}\bm{\Sigma}^{L}+\mathcal{O}_{N}(N^{-2})\
,\ \bm{\Sigma}^{L}=\left[-\nabla^{2}\mathcal{S}(\bm{q}_{\infty})\right]^{-1}$
(97)

The matrix $\frac{1}{N}\bm{\Sigma}^{L}$ gives the leading order (in powers of
$1/N$) covariance of the order parameters over random initialization at width
$N$ depth $L$. The goal of this section is to show that the entries of this
covariance matrix $\bm{\Sigma}^{L}$ that correspond to the NTK and
(consequently the predictor $f$) do not diverge with $L$ for the architecture
of this paper (residual networks with $1/\sqrt{L}$ branch scaling), but rather
approach a well defined limit at rate $\mathcal{O}(L^{-1})$. We again stress
that this is unique to this architecture. Non-residual deep networks networks
have propagator entries which scale linearly with $L$ since variances
accumulate rather than average out over layers (Bordelon & Pehlevan, 2022a;
2023). We note that the predictor $f$ can be included in the set of order
parameters $\bm{q}$ at depth $L$, so that the finite width deviation scales as

$\displaystyle\left<(f_{N,L}-f_{\infty,L})^{2}\right>\sim\frac{1}{N}\Sigma_{f}^{L}$
(98)

where averages are taken over random initialization.

###  G.1 Components of the Hessian

Since we are primarily interested in the scaling of the finite width effects
with respect to $N$ and $L$, we will suppress the sample and time indices of
the kernels and instead focus on the dependence of the fluctuations on layer
index $\ell$. This will not change the analysis of the scaling but will save
us from tracking a tremendous amount of indices. We will let $\Phi^{\ell}$ to
represent the feature kernel at layer $\ell$ and $G^{\ell}$ the gradient
kernel at layer $\ell$ with no time or sample indices. Under this simplified
notation, the Hessian has entries of the form

$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{\Phi}^{\ell}\partial\hat{\Phi}^{\ell^{\prime}}}$
$\displaystyle=\left<\phi(h^{\ell})\phi(h^{\ell})\phi(h^{\ell^{\prime}})\phi(h^{\ell^{\prime}})\right>-\Phi^{\ell}\Phi^{\ell^{\prime}}\equiv\kappa^{\ell,\ell^{\prime}}_{\Phi}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{\Phi}^{\ell}\partial\Phi^{\ell^{\prime}}}$
$\displaystyle=\delta_{\ell,\ell^{\prime}}-\frac{\partial}{\partial\Phi^{\ell^{\prime}}}\left<\phi(h^{\ell})\phi(h^{\ell})\right>\equiv\delta_{\ell,\ell^{\prime}}-\frac{1}{L}D^{\ell,\ell^{\prime}}_{\Phi}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\Phi^{\ell}\partial\Phi^{\ell^{\prime}}}$
$\displaystyle=0$ (99)

First, we note that since the $h^{\ell}$ fields have long-range correlations,
the $\kappa^{\ell,\ell^{\prime}}$ is non-vanishing for all
$\ell,\ell^{\prime}$. We next note that
$D_{\Phi}^{\ell,\ell^{\prime}}\sim\mathcal{O}(1)$ since $\frac{\partial
h^{\ell}}{\partial\Phi^{\ell^{\prime}}}\sim\mathcal{O}(L^{-1})$ and by Price’s
theorem,

$\displaystyle
D^{\ell,\ell^{\prime}}_{\Phi}=\frac{1}{2}\left<\frac{\partial^{2}}{\partial
u^{\ell^{\prime}}\partial
u^{\ell^{\prime}}}[\phi(h^{\ell})\phi(h^{\ell})]\right>+2L\left<\dot{\phi}(h^{\ell})\phi(h^{\ell})\frac{\partial
h^{\ell}}{\partial\Phi^{\ell^{\prime}}}\right>\sim\mathcal{O}(1)$ (100)

Similarly, we also need to calculate the terms involving $G^{\ell}$

$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{G}^{\ell}\partial\hat{G}^{\ell^{\prime}}}$
$\displaystyle=\left<g^{\ell}g^{\ell}g^{\ell^{\prime}}g^{\ell^{\prime}}\right>-G^{\ell}G^{\ell^{\prime}}\equiv\kappa^{\ell,\ell^{\prime}}_{G}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{G}^{\ell}\partial
G^{\ell^{\prime}}}$
$\displaystyle=\delta_{\ell,\ell^{\prime}}-\frac{\partial}{\partial
G^{\ell^{\prime}}}\left<g^{\ell}g^{\ell}\right>\equiv\delta_{\ell,\ell^{\prime}}-\frac{1}{L}D^{\ell,\ell^{\prime}}_{G}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial G^{\ell}\partial
G^{\ell^{\prime}}}$ $\displaystyle=0$ (101)

Lastly, we must consider the cross terms involving both $\\{\Phi^{\ell}\\}$
and $\\{G^{\ell}\\}$

$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{\Phi}^{\ell}\partial\hat{G}^{\ell^{\prime}}}$
$\displaystyle=\left<\phi(h^{\ell})\phi(h^{\ell})g^{\ell^{\prime}}g^{\ell^{\prime}}\right>-\Phi^{\ell}G^{\ell^{\prime}}\equiv\kappa^{\ell,\ell^{\prime}}_{\Phi,G}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{\Phi}^{\ell}\partial
G^{\ell^{\prime}}}$ $\displaystyle=-\frac{\partial}{\partial
G^{\ell^{\prime}}}\left<\phi(h^{\ell})\phi(h^{\ell})\right>\equiv-\frac{1}{L}D^{\ell,\ell^{\prime}}_{\Phi,G}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\hat{G}^{\ell}\partial\Phi^{\ell^{\prime}}}$
$\displaystyle=-\frac{\partial}{\partial\Phi^{\ell^{\prime}}}\left<g^{\ell}g^{\ell}\right>\equiv-\frac{1}{L}D^{\ell,\ell^{\prime}}_{G,\Phi}$
$\displaystyle\frac{\partial^{2}\mathcal{S}}{\partial\Phi^{\ell}\partial
G^{\ell^{\prime}}}$ $\displaystyle=0$ (102)

The $\kappa$ tensor represents an instantaneous source of variance for the
kernels, while the $D$ tensor represents the sensitivity of the kernel at
layer $\ell$ to a perturbation in the kernel at layer $\ell^{\prime}$. All
entries of each $D$ matrix are $\mathcal{O}(1)$. Now, to compute the
propagator, we construct block matrices

$\displaystyle\bm{\kappa}=\begin{bmatrix}\bm{\kappa}_{\Phi}&\bm{\kappa}_{\Phi,G}\\\
\bm{\kappa}_{G,\Phi}&\bm{\kappa}_{G}\end{bmatrix}\ ,\
\bm{D}=\begin{bmatrix}\bm{D}_{\Phi}&\bm{D}_{\Phi,G}\\\
\bm{D}_{G,\Phi}&\bm{D}_{G}\end{bmatrix}$ (103)

Now, to obtain the propagator $\bm{\Sigma}$, we must compute the inverse of
the Hessian $\nabla^{2}\mathcal{S}$

$\displaystyle\bm{\Sigma}=-\begin{bmatrix}\bm{0}&(\bm{I}-\frac{1}{L}\bm{D})^{\top}\\\
\bm{I}-\frac{1}{L}\bm{D}&\bm{\kappa}\end{bmatrix}^{-1}=\begin{bmatrix}\left(\bm{I}-\frac{1}{L}\bm{D}\right)^{-1}\bm{\kappa}\left[\left(\bm{I}-\frac{1}{L}\bm{D}\right)^{-1}\right]^{\top}&-\left(\bm{I}-\frac{1}{L}\bm{D}\right)^{-1}\\\
-\left[\left(\bm{I}-\frac{1}{L}\bm{D}\right)^{-1}\right]^{\top}&\bm{0}\end{bmatrix}$
(104)

The observables of interest are the top-left block which gives the joint
covariance over all $\\{\Phi^{\ell},G^{\ell}\\}$ variables. We introduce the
following $2\times 2$ block matrices

$\displaystyle\bm{\kappa}^{\ell,\ell^{\prime}}=\begin{bmatrix}\kappa_{\Phi}^{\ell,\ell^{\prime}}&\kappa_{\Phi
G}^{\ell,\ell^{\prime}}\\\
\kappa_{G\Phi}^{\ell,\ell^{\prime}}&\kappa_{G}^{\ell,\ell^{\prime}}\end{bmatrix}\
,\
\bm{\Sigma}^{\ell,\ell^{\prime}}=\begin{bmatrix}\Sigma_{\Phi}^{\ell,\ell^{\prime}}&\Sigma_{\Phi
G}^{\ell,\ell^{\prime}}\\\
\Sigma_{G\Phi}^{\ell,\ell^{\prime}}&\Sigma_{G}^{\ell,\ell^{\prime}}\end{bmatrix}\
,\
\bm{D}^{\ell,\ell^{\prime}}=\begin{bmatrix}D_{\Phi}^{\ell,\ell^{\prime}}&D_{\Phi
G}^{\ell,\ell^{\prime}}\\\
D_{G\Phi}^{\ell,\ell^{\prime}}&D_{G}^{\ell,\ell^{\prime}}\end{bmatrix}$ (105)

The above equation gives us we get the following $2\times 2$ matrix equations

$\displaystyle\bm{\Sigma}^{\ell,\ell^{\prime}}-\frac{1}{L}\sum_{k}\bm{D}^{\ell,k}\bm{\Sigma}^{k,\ell^{\prime}}-\frac{1}{L}\sum_{k}\bm{\Sigma}^{\ell,k}[\bm{D}^{\ell^{\prime},k}]^{\top}+\frac{1}{L^{2}}\sum_{kk^{\prime}}\bm{D}^{\ell,k}\bm{\Sigma}^{k,k^{\prime}}[\bm{D}^{k^{\prime},\ell^{\prime}}]^{\top}=\bm{\kappa}^{\ell,\ell^{\prime}}$
(106)

We see that the solutions $\bm{\Sigma}^{\ell,\ell^{\prime}}$ to the above
equation will be $\mathcal{O}_{L}(1)$. This is precisely due to the factors of
$1/L$ which appear in the above sums. One way to see this is to consider a
large $L\gg 1$ limit where the above layer sums converge to integrals over
$\tau$

$\displaystyle\bm{\Sigma}(\tau,\tau^{\prime})-\int
d\tau^{\prime\prime}\bm{D}(\tau,\tau^{\prime\prime})\bm{\Sigma}(\tau^{\prime\prime},\tau^{\prime})-\int
d\tau^{\prime\prime}\bm{\Sigma}(\tau^{\prime\prime},\tau^{\prime})\bm{D}(\tau,\tau^{\prime\prime})^{\top}$
$\displaystyle+\int
d\tau^{\prime\prime}d\tau^{\prime\prime\prime}\bm{D}(\tau,\tau^{\prime\prime})\bm{\Sigma}(\tau^{\prime\prime},\tau^{\prime\prime\prime})\bm{D}(\tau^{\prime},\tau^{\prime\prime\prime})^{\top}=\bm{\kappa}(\tau,\tau^{\prime})$
(107)

So we find that the covariance of kernels at each pair of layers
$\bm{\Sigma}(\tau,\tau^{\prime})$ is $\mathcal{O}_{L}(1)$.

###  G.2 Variance of the NTK and Predictor

Using the above fact that
$\bm{\Sigma}^{\ell,\ell^{\prime}}\sim\mathcal{O}_{L}(1)$, we can reason about
the variance of the neural tangent kernel $K$

$\displaystyle\text{Var}(K)$
$\displaystyle=\frac{1}{NL^{2}}\sum_{\ell\ell^{\prime}}\bm{\Sigma}_{\Phi}^{\ell,\ell^{\prime}}G^{\ell+1}G^{\ell^{\prime}+1}+\frac{1}{NL^{2}}\sum_{\ell\ell^{\prime}}\bm{\Sigma}_{G}^{\ell,\ell^{\prime}}\Phi^{\ell-1}\Phi^{\ell^{\prime}-1}$
$\displaystyle+\frac{2}{NL^{2}}\sum_{\ell\ell^{\prime}}\bm{\Sigma}_{\Phi
G}^{\ell,\ell^{\prime}}G^{\ell+1}\Phi^{\ell^{\prime}-1}\sim\mathcal{O}_{L,N}(N^{-1}).$
(108)

This demonstrates that the NTK will have initialization variance that scales
only with the width $N$ but not the depth $L$. Using the fact that
$\frac{\partial f}{\partial t}=K\Delta$, we see that $f$ will inherit
fluctuations on the same scale of $\mathcal{O}_{N,L}(N^{-1/2})$. The depth $L$
NTK and predictor covariance converge to the infinite depth propagator with
rate $\mathcal{O}(L^{-1})$.

##  Appendix H Finite Width And Depth Error Analysis

In this section, we combine the results of the previous sections. We consider
the following error decomposition where we introduce the norm
$||z||=\sqrt{\left<z^{2}\right>}$ and $\left<\right>$ denotes an average over
random initializations of the network. Applying the triangle inequality on
this norm, we have

$\displaystyle||f_{N,L}-f_{\infty,\infty}||$
$\displaystyle\leq||f_{N,L}-\left<f_{N,L}\right>||\quad\text{(Finite Width
Variance)}$
$\displaystyle+||\left<f_{N,L}\right>-f_{\infty,L}||\quad\text{(Mean Predictor
Finite Width Error)}$
$\displaystyle+||f_{\infty,L}-f_{\infty,\infty}||\quad\text{(Finite Depth
Error at Infinite Width)}$ (109)

From the preeceeding section G, the first term has asymptotic behavior

$\displaystyle\left<\left(f_{N,L}-\left<f_{N,L}\right>\right)^{2}\right>\sim\frac{1}{N}\Sigma_{f}^{L}+\mathcal{O}(N^{-2})=\frac{1}{N}\Sigma_{f}^{\infty}+\mathcal{O}\left(\frac{1}{NL}+\frac{1}{N^{2}}\right)$
(110)

where the last term comes from the fact that $\Sigma^{L}$ approximates the
infinite depth propagator $\Sigma^{\infty}$ with error $1/L$. Next, we can
apply results developed from Bordelon & Pehlevan (2023) which show that the
mean predictor satisfies
$\left<f_{N,L}\right>=f_{\infty,L}+\mathcal{O}(N^{-1})$ to find that mean
finite width error scales as

$\displaystyle|\left<f_{N,L}\right>-f_{\infty,L}|^{2}\sim\mathcal{O}(N^{-2})$
(111)

Lastly, the finite depth error in the infinite width DMFT arises from the
$1/L$ discretization effect. Thus the square error goes as

$\displaystyle|f_{\infty,L}-f_{\infty,\infty}|^{2}\sim\mathcal{O}(L^{-2}).$
(112)

which follows from the discretization error of the limiting large depth
process for the NTK and thus the predictor. Altogether, we have

$\displaystyle||f_{N,L}-f_{\infty,\infty}||\sim\sqrt{\frac{1}{N}\Sigma^{\infty}_{f}+\mathcal{O}\left(\frac{1}{NL}\right)}+\mathcal{O}\left(\frac{1}{N}\right)+\mathcal{O}\left(\frac{1}{L}\right)$
(113)

We see that at large width and depth this is dominated by the first term which
comes from the finite width fluctuations. However, if the model is ensembled,