# Random Fully Connected Neural Networks as Perturbatively Solvable
Hierarchies

Boris Hanin111BH gratefully acknowledges support by the NSF through
DMS-2143754, DMS-1855684, and DMS-2133806 as well support from an ONR MURI on
Foundations of Deep Learning.
Department of Operations Research
and Financial Engineering
Princeton University
<EMAIL_ADDRESS>
###### Abstract

This article considers fully connected neural networks with Gaussian random
weights and biases as well as $L$ hidden layers, each of width proportional to
a large parameter $n$. For polynomially bounded non-linearities we give sharp
estimates in powers of $1/n$ for the joint cumulants of the network output and
its derivatives. Moreover, we show that network cumulants form a
perturbatively solvable hierarchy in powers of $1/n$ in that $k$-th order
cumulants in one layer have recursions that depend to leading order in $1/n$
only on $j$-th order cumulants at the previous layer with $j\leq k$. By
solving a variety of such recursions, however, we find that the depth-to-width
ratio $L/n$ plays the role of an effective network depth, controlling both the
scale of fluctuations at individual neurons and the size of inter-neuron
correlations. Thus, while the cumulant recursions we derive form a hierarchy
in powers of $1/n$, contributions of order $1/n^{k}$ often grow like $L^{k}$
and are hence non-negligible at positive $L/n$. We use this to study a
somewhat simplified version of the exploding and vanishing gradient problem,
proving that this particular variant occurs if and only if $L/n$ is large.
Several key ideas in this article were first developed at a physics level of
rigor in a recent monograph of Daniel A. Roberts, Sho Yaida, and the author.
This article not only makes these ideas mathematically precise but also
significantly extends them, opening the way to obtaining corrections to all
orders in $1/n$.

##  1 Introduction

We live in an era of big data and cheap computation. This has led to
remarkable progress in domains ranging from self-driving cars [35] to
automatic drug discovery [33] and machine translation [8]. Underpinning many
of these exciting practical developments is a class of computational models
called neural networks. While they were originally developed in the 1940’s and
1950’s [26, 54], the complexity of state-of-the-art neural nets is
unprecedented. And yet, despite their empirical utility, a theoretical
understanding of how they work and how to make them better is nascent. In
fact, it is sometimes said that neural networks are simply too complicated to
allow for a rigorous understanding of their key features.

This article adds to the growing body of literature to the contrary. Namely,
in the simplest important setting of fully connected networks, we develop a
flexible set of probabilistic tools for studying correlation functions of
random fully connected neural networks at finite width.

A random fully connected neural network, defined precisely in §2.1, is a
random field whose law is determined by a fixed and typically non-linear
function $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ as well as two structural
parameters: a depth $L\in\mathbb{N}$ and a width $n\in\mathbb{N}$. For a fixed
depth $L$, in the limit of infinite width $n\rightarrow\infty$ random neural
networks converge in distribution to Gaussian processes (see Theorem 2.2). To
study finite width effects we therefore consider the higher cumulants for the
distribution of the outputs. Our approach is recursive in the network depth
$L$ and perturbatively in the reciprocal $1/n$ of the network width. Our main
results are:

  * •

We give in Theorem 3.1 sharp estimates in powers of $1/n$ for the the size of
all joint cumulants of networks outputs and their derivatives. These can be
viewed as quantitative results refinements of the central limit theorem at
fixed $L$ and large $n$ (cf Theorems 2.2 and Theorem 3.2).

  * •

We derive exact recursions with respect to depth that describe cumulants at
one layer in terms of cumulants at the previous layer (see e.g. Corollary 3.4
and Proposition 11.2). These recursions take the form of perturbatively
solvable hierarchies in the sense that the recursion for the $k$-th cumulant
at layer $\ell+1$ involves, at leading order in $1/n$ only the $j$-th
cumulants at layer $\ell$ with $j\leq k$ (see Corollaries 3.3 and 3.4). This
is distinct from but similar in spirit to the remarkable article [29], which
derived a dynamical perturbative hierarchy for the so-called neural tangent
kernel. As in our work, the perturbative parameter is $1/n$. However, the
emphasis in [29] is primarily on training dynamics while ours is on
understanding the effect of depth.

  * •

Solving some special cases of the cumulant recursions to leader order in $1/n$
reveals that it is the depth-to-width ratio $L/n$, which we term the effective
network depth, rather than the apparent depth $L$, that is a more informative
measure of neural network depth and complexity. Mathematically, this suggests
a non-trivial double scaling limit for random neural networks in which

$n,L\rightarrow\infty\qquad\text{and}\qquad L/n\rightarrow\xi\in[0,\infty).$

For non-linear networks this scaling limit has only started to be considered
[17, 21, 53, 40]. Even in the very special case of product of $L$ iid random
$n\times n$ matrices (sometimes called deep linear networks) the simultaneous
large $n,L$ regime has revealed a range of interesting and not fully
understood properties (see e.g. references in [1, 2] as well as [16, 21, 22,
42]). In contrast to the $\xi=0$ regime typically considered in previous work
on neural networks (c.f. e.g. [12, 13, 32, 41]), we show that when $\xi>0$ our
double scaling limit is capable of exhibiting non-Gaussian and non-linear
effects. We find the following effects to leading order in $1/n$:

    * –

We prove in Corollary 4.1 that for any non-linearity $\sigma$ from the
$K_{*}=0$ universality class (defined in §4.2.2) and for $k=2,3,4$, the
$2k^{th}$ cumulants of the output of a random neural network with non-
linearity $\sigma$ grow like $(L/n)^{k/2-1}$. This implies, for instance, that
both the correlations between neurons and the fluctuations of a single neuron
grow like the effective depth $L/n$ at large $L,n$ (see Remark 4.2 and just
after Theorem 4.4). Since the components of the output of a depth $0$ neural
network are simply iid Gaussians, we see that the effective depth $L/n$ can be
thought of as a measure of how close a random fully connected network is to a
Gaussian process. Related questions were considered in [4, 15].

    * –

We show in Theorems 4.4 and 4.5 that, for random neural networks initialized
as in practice (see §4.1), the variance of the gradient of the network output
with respect to either its input or a trainable parameter in its first layer
grows like $L/n$ at large $L$. As explained in §4.4, this gives the first
mathematical characterization for fairly general fully connected networks of
the so-called exploding and vanishing gradient problem (EVGP). Previously,
this problem was solved rigorously in the important special case of random
ReLU networks [17, 21] and solved at a physics level of rigor by a somewhat
different but related set of ideas in the monograph [53] of Roberts, Yaida,
and the author.

The remainder of the introduction is structured as follows. We begin by giving
in §2.1 the precise definition of random fully connected neural networks. We
then formulate and motivate the main question taken up in this article in
§2.2.

##  2 Background and Motivation

###  2.1 What is a (Random) Fully Connected Neural Network?

Neural networks are parameterized families of functions. The simplest kind of
networks are called fully connected. Each such network is specified by its
architecture, which consists of an input dimension $n_{0}\in\mathbb{Z}_{+}$, a
network depth $L\in\mathbb{Z}_{+}$, hidden layer widths
$n_{1},\ldots,n_{L}\in\mathbb{Z}_{+}$, an output dimension
$n_{L+1}\in\mathbb{Z}_{+}$, and a non-linearity
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$. The functions computed by a fully
connected network with a given architecture are all maps from
$\mathbb{R}^{n_{0}}$ to $\mathbb{R}^{n_{L+1}}$ that associate to each network
input $x_{\alpha}\in\mathbb{R}^{n_{0}}$ an output
$z_{\alpha}^{{}_{(L+1)}}\in\mathbb{R}^{n_{L+1}}$ through a sequence of
intermediate representations
$z_{\alpha}^{{}_{(\ell)}}\in\mathbb{R}^{n_{\ell}}$ as follows:

$z_{\alpha}^{(\ell+1)}:=\begin{cases}b^{(\ell+1)}+W^{(\ell+1)}\sigma(z_{\alpha}^{(\ell)}),&\quad\ell\geq
1\\\ b^{(1)}+W^{(1)}x_{\alpha},&\quad\ell=0\end{cases},\qquad
W^{(\ell+1)}\in\mathbb{R}^{n_{\ell+1}\times
n_{\ell}},\,b^{(\ell+1)}\in\mathbb{R}^{n_{\ell+1}}.$ (2.1)

In this recursion, the univariate function $\sigma$ applied to a vector
$z_{\alpha}^{{}_{(\ell)}}\in\mathbb{R}^{n_{\ell}}$ is short-hand for applying
it separately to each component. The entries of the matrices $W^{(\ell)}$ and
the components of the vectors $b^{(\ell)}$ are called the weights and biases
in layer $\ell$, respectively. One typically refers to

$z_{\alpha}^{(\ell)}=\left(z_{1;\alpha}^{(\ell)},\ldots,z_{n_{\ell};\alpha}^{(\ell)}\right)\in\mathbb{R}^{n_{\ell}}$

as the vector of pre-activations at layer $\ell$ corresponding to the input
$x_{\alpha}$. The most popular choices of $\sigma$ in practice include
$\mathrm{ReLU}(t):=\max\left\\{0,t\right\\}$ as well the hyperbolic tangent
and their variations. We will analyze these cases in detail later (see §4.2),
but for our general results make only the following mild assumption

######  Assumption 2.1.

There exists $r\geq 1$ so that the $r$-th derivative of $\sigma$ exists almost
everywhere and grows at most polynomially:

$\exists k\geq 1\text{ s.t.
}\mathrm{sup}_{x\in\mathbb{R}}\left|(1+\left|x\right|)^{-k}\frac{d^{r}}{dx^{r}}\sigma(x)\right|<\infty.$

The primary objects of study in this article are random fully connected neural
networks, obtained by choosing network weights and biases to be independent
centered Gaussians:

$W_{ij}^{(\ell)}\sim\mathcal{N}(0,C_{W}/n_{\ell-1}),\qquad
b_{i}^{(\ell)}\sim\mathcal{N}(0,C_{b})\qquad\text{independent}.$ (2.2)

Here $C_{b}\geq 0,C_{W}>0$ are fixed constants. The $1/n_{\ell-1}$ scaling in
the weight variance ensures that the moments of the outputs
$z_{\alpha}^{{}_{(L+1)}}$ remain uniformly bounded as
$n_{1},\ldots,n_{L},\rightarrow\infty$ (see e.g. Theorem 2.2 and (3.6)). While
we carry out a substantial portion of our analysis for any choice of
$C_{b},C_{W}$, as we will see in §4.1, there often exist distinguished
$\sigma$-dependent settings of $C_{b},C_{W}$ for which random neural networks
at infinite width $n_{1},\ldots,n_{L}\rightarrow\infty$ are well-behaved at
large depth $L$.

###  2.2 Statement and Motivation for Questions Addressed in this Article

The main problem we take up in the present article is to characterize the
finite dimensional distributions of a random neural network $x_{\alpha}\mapsto
z_{\alpha}^{{}_{(L+1)}}$ (and its derivatives with respect to $x_{\alpha}$) in
the regime where the input dimension $n_{0}$ is arbitrary with the inputs
$x_{\alpha}$ satisfying

$\frac{1}{n_{0}}\left|\left|x_{\alpha}\right|\right|_{2}^{2}<\infty,$

the output dimension $n_{L+1}$ is fixed, and the hidden layer widths
$n_{\ell}$ are large but finite:

$\exists c,C>0\text{ s.t. }\qquad cn\leq n_{1},\ldots,n_{L}\leq
Cn,\qquad\qquad n\gg 1.$ (2.3)

Our approach will be to describe the random field $x_{\alpha}\mapsto
z_{\alpha}^{{}_{(L+1)}}$ perturbatively in $1/n$ and recursively in $L$.
Before proceeding to the technical statements of our results, we pause to
address below the following motivational questions:

  * •

Why study random neural networks?

  * •

Why treat $1/n$ as a perturbative parameter?

  * •

Why study networks at finite width?

  * •

What role does $\sigma$ play?

The primary use of a neural network with a fixed architecture is to find a
setting of its parameters $W^{{(\ell)}},b^{{(\ell)}}$ giving rise to a mapping
$x_{\alpha}\mapsto z_{\alpha}^{{}_{(L+1)}}$ as in (2.1) that matches a dataset
$\left\\{(x_{i},f(x_{i}))\right\\}$ of values for some otherwise unknown
function $f:\mathbb{R}^{n_{0}}\rightarrow\mathbb{R}^{n_{L+1}}$. The
optimization procedure for finding these weights and biases is almost always
some variant of gradient descent starting with $W^{{(\ell)}},b^{{(\ell)}}$
drawn at random from the distribution (2.2). Thus, studying the properties of
random neural networks gives direct insights into the starting conditions for
neural network optimization. For instance, understanding the behavior of
neural networks at the start of training gives a principled way to set
optimization hyperparameters (e.g. the variances $C_{b},C_{W}$ and the step
size used for gradient descent). We refer the interested reader our discussion
of criticality in §4 as well as to §3.2 in [53] the articles articles [25, 43,
52, 67] for more on this point.

Next, to address why $1/n$ can reasonably be treated perturbatively, we recall
that networks used in practice often have a very large number of parameters.
This is reflected in the fact that both the layer widths $n_{\ell}$ and the
network depth $L$ are often big. Thus, it is sensible to first understand
various limits of neural networks in which the number of parameters tends to
infinity. The most well-studied regime of this type (though not the only
option cf eg [45, 55, 57, 58, 64]) is the infinite width limit, also known as
the NTK regime. By definition, this regime is accessed by fixing the depth
$L$, the input and output dimensions $n_{0},n_{L+1}$, and the non-linearity
$\sigma$, the initialization scheme (2.2) and considering the limit when
$n_{1},\ldots,n_{L}\rightarrow\infty$. From the random matrix theory point of
view this is the free probability regime. In view of the relation (2.3) the
NTK regime is obtained by taking $n\rightarrow\infty$ at fixed $L$ and has two
salient features:

  * •

At the start of training, neural networks converge to a Gaussian process (see
Theorem 2.2 below) [36, 44, 48, 49, 62, 63, 65].

  * •

For the purposes of optimization of a squared loss, the network can be
replaced by its linearization at the start of training (see [9, 12, 32, 37,
41]).

Taken together, these two points show that at least at infinite width and
finite depth, it is the structure of the network at initialization that
determines not only the start of training but really the entire training
trajectory. However, the infinite width limit is too rigid to capture the
ability of real work networks to learn data-dependent features (see e.g. [20,
29, 53, 64]). Only finite width networks can capture these effects! Since the
starting point for our analysis of random neural networks at finite width is
their infinite width behavior, we take this opportunity to record the
following result about random neural networks in the NTK regime.

######  Theorem 2.2 (Random Networks at fixed $L$ and Infinite Width are
Gaussian Processes).

Fix a non-negative integer $r\geq 0$, and suppose
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is $r$-times differentiable and that
its $r$-th derivative is polynomially bounded:

$\exists k\geq 1\text{ s.t.
}\sup_{x\in\mathbb{R}}\left|(1+\left|x\right|)^{-k}\frac{d^{r}}{dx^{r}}\sigma(x)\right|<\infty.$

Then the finite-dimensional distributions of the stochastic process
$x_{\alpha}\mapsto z_{\alpha}^{(L+1)}$ and its derivatives of order up to $r$
converge to those of a centered Gaussian process with $n_{L+1}$ iid
components. The limiting variance of each component

$K_{\alpha\beta}^{(L+1)}:=\lim_{n_{1},\ldots,n_{L}\rightarrow\infty}\mathrm{Cov}\left(z_{i;\alpha}^{(L+1)},z_{i;\beta}^{(L+1)}\right),\qquad
x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{0}},$

satisfies the recursion

$\displaystyle
K_{\alpha\beta}^{(\ell+1)}=\begin{cases}C_{b}+C_{W}\left\langle\sigma(z_{\alpha})\sigma(z_{\beta})\right\rangle_{K^{(\ell)}},&\quad\ell\geq
1\\\
C_{b}+\frac{C_{W}}{n_{0}}\sum_{j=1}^{n_{0}}x_{j;\alpha}x_{j;\beta},&\quad\ell=0\end{cases}.$
(2.4)

In the statement of Theorem 2.2 and henceforth we reserve the symbol
$\left\langle f(z_{\alpha},z_{\beta})\right\rangle_{\kappa}$ to denote the
expectation of $f(z_{\alpha},z_{\beta})$ with respect to the Gaussian
distribution

$\left(z_{\alpha},z_{\beta}\right)\sim\mathcal{N}\left(0,\left(\begin{array}[]{cc}\kappa_{\alpha\alpha}&\kappa_{\alpha\beta}\\\
\kappa_{\alpha\beta}&\kappa_{\beta\beta}\end{array}\right)\right),$

where $\kappa_{\alpha\beta}=\kappa(x_{\alpha},x_{\beta})$ is a given
covariance function. The conclusion in Theorem 2.2 is not new, having been
obtained many times and under a variety of different assumptions (including
for more general architectures) [19, 36, 44, 51, 63]. We refer the interested
reader to [19] for a discussion of prior work and note only that convergence
of the derivatives of the field $z_{\alpha}^{{}_{(L+1)}}$ to its Gaussian
limit does not seem to have been previously considered. We give a short proof
that includes convergence of derivatives along the lines of the arguments in
[19, 36] in Appendix §A.

Let us remark that the relation between $K^{{(\ell+1)}}$ to $K^{(\ell)}$
supplied by the recursion (2.4) is in general non-linear due to the presence
of $\sigma$ and depends crucially on the weight and bias variances
$C_{W},C_{b}$. Given $\sigma$, finding a choice of $C_{b},C_{W}$ so that this
recursion is well-behaved (e.g. not exponentially growing or decaying) at
large $\ell$ is an important practical matter [30, 36, 50]. We explain this in
§4.1, where we also point out that the different possible large $\ell$ limits
suggest the existence of universality classes of random neural networks.

Finally, as we alluded to above, the tremendous simplification of random
neural networks in the NTK regime comes at a steep cost in terms of the
descriptive power of the resulting model in the sense that such models cannot
learn data-dependent features. We again refer to reader to articles such as
[10, 45, 55, 57, 58, 64] in which a different initialization scheme leads to
infinite width limits capable of feature learning. Since we do not treat
feature learning in this article, we will not elaborate further on this point,
referring the interested reader to [9, 20, 29, 53, 60, 64] instead.

Finally, before formulating the new results in this article, we point the
interested to a final tranche of physics-style references such as [3, 11, 14,
34, 38, 46, 47, 56], which analyze a variety of questions related to either
very wide or infinitely wide networks.

##  3 Results

In this section we give precise formulations of our results on random neural
networks. We start in §3.1 by stating a structural result, Theorem 3.1, on the
size of the joint cumulants of a random neural networks and its derivatives in
powers of $1/n$. We then extend this order of magnitude estimate to obtain in
Theorem 3.2 a general prescription for obtaining series expansions in powers
of $1/n$ for expectation values of functions of a random depth $\ell+1$ neural
network in terms of those of a random depth $\ell$ network. An application of
Theorem 3.2 yields explicit recursions for the $4^{th},6^{th},8^{th}$ joint
cumulants for components $z_{i;\alpha}^{{}_{(\ell)}}$ of the vector of pre-
activations at layer $\ell$ corresponding to a fixed network input
$x_{\alpha}.$ These recursions are recorded in Corollary 3.4. See §5 for an
overview of the proof of Theorems 3.1 and 3.2 as well as Corollary 3.4.

After stating these results, we describe in §4 the important notion of
criticality in random neural networks, which corresponds to choosing values of
the weight and bias variance parameters $C_{b},C_{W}$ in (2.2) depending on
$\sigma$ in such a way that the infinite width covariance
$K_{\alpha\beta}^{{}_{(\ell)}}$ (and in fact the higher cumulants at finite
width as well) is well-behaved at large $\ell$. The variety of resulting large
$\ell$ behaviors determines universality classes of random neural networks, as
we explain in §4.2.

We then turn in §4.2.1 and §4.2.2 to solving to the $2k$-th cumulant
recursions for $k=1,2,3,4$ from Corollary 3.4 in random networks tuned to
criticality with non-linearities either from the universality class of ReLU or
of tanh. We will see that, at large network depth $L$ and leading order in
$1/n$, these cumulants depend only on the effective network depth $L/n$. We
formalize our belief that this is a universal phenomenon in Conjecture 4.3. We
then solve in §4.4 a simplified version of the so-called exploding and
vanishing gradient problem (EVGP), which concerns the empirical variance of
parameter gradients in a random neural network. More precisely, for non-
linearities in the $K_{*}=0$ universality class in prove (see Theorem 4.5)
that to leading order in $1/n$ the empirical variance over weights in the
first layer of network gradients grows linearity in the effective network
depth $L/n$. These results depend on Theorem 4.4, which gives for non-
linearities in the $K_{*}=0$ universality class asymptotics at large $\ell$
for joint fourth cumulants between the values of partial derivatives of the
output of a random neural network.

###  3.1 Cumulants of Random Neural Networks at Finite Width

Since in the infinite $n$ limit, the field $z_{\alpha}^{{}_{(L+1)}}$ is
Gaussian (see Theorem 2.2), it is natural to measure perturbations around this
regime by considering the behavior of the cumulants of
$z_{\alpha}^{{}_{(L+1)}}$ and its derivatives. Let us therefore agree that,
given random variables $X_{1},\ldots,X_{k}$ with finite moments defined on the
same probability space, we will denote their mixed cumulant by

$\kappa\left(X_{1},\ldots,X_{k}\right):=i^{k}\frac{\partial^{k}}{\partial
t_{1}\cdots\partial
t_{k}}\bigg{|}_{t=0}\log\mathbb{E}\left[\exp\left[-i(t_{1}X_{1}+\cdots+t_{k}X_{k})\right]\right].$
(3.1)

Thus, for example, $\kappa(X_{1})=\mathbb{E}\left[X_{1}\right]$ and
$\kappa(X_{1},X_{2})=\mathrm{Cov}(X_{1},X_{2}).$ We refer the reader to §6.1
for background on cumulants. Our first result, Theorem 3.1, gives estimates on
the order in $1/n$ of the cumulants of $z_{\alpha}^{{}_{(L+1)}}$ and its
derivatives. To state it, we fix a finite collection

$\left\\{x_{\alpha},\,\alpha\in\mathcal{A}\right\\}\subseteq\mathbb{R}^{n_{0}},$

of $\left|\mathcal{A}\right|$ distinct network inputs. Moreover, we fix a
collection of $p$ directional derivatives:

$D=\left(d_{1},\ldots,d_{p}\right),\qquad
d_{j}:=\nabla_{v_{j}}=\sum_{i=1}^{n_{0}}v_{ij}\partial_{x_{i}}.$ (3.2)

and for any multi-index $J=\left(j_{1},\ldots,j_{p}\right)\in\mathbb{N}^{p}$
set

$D_{\alpha}^{J}:=d_{1}^{j_{1}}\cdots d_{m}^{j_{p}}\bigg{|}_{x=x_{\alpha}}$

for the corresponding differential operator of order
$\left|J\right|:=j_{1}+\cdots+j_{p}.$

######  Theorem 3.1 (order of magnitude for cumulants of random neural
networks).

Fix $r,L\geq 1$ and suppose that $\sigma:\mathbb{R}\rightarrow\mathbb{R}$
satisfies Assumption (2.1) with this value of $r$. Suppose further that one of
the following two conditions holds:

  * •

$\sigma$ is smooth

  * •

the limiting covariance matrix

$\left(\lim_{n\rightarrow\infty}\mathrm{Cov}\left(D_{\alpha_{1}}^{J_{1}}z_{1;\alpha_{1}}^{(\ell)},D_{\alpha_{2}}^{J_{2}}z_{1;\alpha_{2}}^{(\ell)}\right)\right)_{\begin{subarray}{c}\left|J_{1}\right|,\left|J_{2}\right|\leq
r\\\ \alpha_{1},\alpha_{2}\in\mathcal{A}\end{subarray}}$ (3.3)

of derivatives of order at most $r$ in the directional derivatives
$d_{1},\ldots,d_{p}$ of the scalar field $z_{1;\alpha}^{(\ell)}$ is strictly
positive definite in the infinite width limit for all $\ell\leq L$.

Then, for each $k,\ell\geq 1$, as $n\rightarrow\infty$

$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell)},\ldots,D_{\alpha_{k}}^{J_{k}}z_{i_{k};\alpha_{k}}^{(\ell)}\right)=\begin{cases}0,&\quad
k\text{ odd}\\\ O(n^{-\frac{k}{2}+1}),&\quad k\text{ even}\end{cases},$ (3.4)

where the implicit constant in the error term depends on $k$, the inputs
$x_{\alpha_{1}},\ldots,x_{\alpha_{k}},$ the multi-indices
$J_{1},\ldots,J_{k}$, the weight and bias variances $C_{b},C_{W},$ the non-
linearity $\sigma$, and the layer index $\ell$.

We prove Theorem 3.1 in §7. At a physics level of rigor, Theorem 3.1 with
$k=4$ and no derivatives was already derived in the breakthrough work of Yaida
[61]. In fact, Yaida’s original article went much further: it obtained a
recursive formula with respect to $\ell$ for the fourth cumulant
$\kappa(z_{i_{1};\alpha_{1}}^{{}_{(\ell)}},\ldots,z_{i_{4};\alpha_{4}}^{{}_{(\ell)}})$
at layer $\ell$ in terms of the second and fourth cumulants at layer $\ell-1$.
This is analogous to the recursion (2.4) for the infinite width covariance
$K_{\alpha_{1}\alpha_{2}}^{{}_{(\ell)}}$. This theme was then picked up and
significantly expanded upon in the physics monograph [53], which computes,
among other things, at order $1/n$ the leading corrections to the field
$z_{\alpha}^{{}_{(\ell)}}$ and its derivatives with respect to both
$x_{\alpha}$ and model parameters. We will reproduce some of these recursions
and obtain new ones of a similar flavor below.

Compared to this prior work the main novelty of Theorem 3.1 is two-fold.
First, it gives sharp estimates in powers of $1/n$ for cumulants of all orders
(the sharpness can already be seen when $\ell=2$). Second, it treats in a
uniform way the cumulants for not only the values but also all derivatives of
$z_{\alpha}^{{}_{(\ell)}}$.

In order to put Theorem 3.1 into context, we take this opportunity to make an
important remark. Namely, it is no accident that the order of magnitude
$O(n^{-k+1})$ for the $2k$-th cumulant in Theorem 3.1 is the same as that of
the $k$-th cumulant of an average of $n$ iid random variables. To see why, let
us denote by $\mathcal{F}^{(\ell)}$ the sigma algebra generated by the weights
and biases in layers up to and including $\ell$. Since we’ve assumed weights
and biases to be Gaussian and independent for different neurons, we find that
conditional on $\mathcal{F}^{(\ell)}$ the vectors

$z_{i;\mathcal{A}}^{(\ell+1)}:=\left(z_{i;\alpha}^{(\ell+1)},\,\alpha\in\mathcal{A}\right)$

are independent centered Gaussians with the covariance

$\Sigma_{\alpha\alpha^{\prime}}^{(\ell)}:=\mathrm{Cov}\left(z_{i;\alpha}^{(\ell+1)},z_{i;\alpha^{\prime}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right)=C_{b}+\frac{C_{W}}{n_{\ell}}\sum_{j=1}^{n_{\ell}}\sigma\left(z_{j;\alpha}^{(\ell)}\right)\sigma\left(z_{j;\alpha^{\prime}}^{(\ell)}\right).$
(3.5)

Put another way, we have the following equality in distribution:

$\left(z_{i;\alpha}^{(\ell+1)}\right)_{i=1}^{n_{\ell+1}}~{}\stackrel{{\scriptstyle
d}}{{=}}~{}\left(\left(\Sigma^{(\ell)}\right)^{1/2}Z_{i}\right)_{i=1}^{n_{\ell+1}},$
(3.6)

where $Z_{i}$ are i.i.d. standard Gaussians and are independent of the random
covariance matrix

$\Sigma^{(\ell)}:=\left(\Sigma_{\alpha\alpha^{\prime}}^{(\ell)}\right)_{\begin{subarray}{c}\alpha,\alpha^{\prime}\in\mathcal{A}\end{subarray}}.$

Inspecting (3.5) shows this conditional covariance has the structure of an
average of $n_{\ell}$ identically distributed random matrices, each depending
only on the pre-activations $z_{i;\mathcal{A}}^{{}_{(\ell)}}$ of a single
neuron at layer $\ell$. We will refer to such random variables as collective
observables. The key point is that while the $z_{i;\mathcal{A}}^{{}_{(\ell)}}$
are not independent for different $i$ at finite width, we will show that they
are sufficiently weakly correlated that the cumulants of
$\Sigma_{\alpha\alpha^{\prime}}^{{}_{(\ell)}}$ have the same order in $n$ as
if they were exactly independent (see Theorem 7.3 and Lemma 7.5). This
reasoning applies equally well to derivatives of
$z_{i;\alpha}^{{}_{(\ell+1)}}$ and explains the order of magnitude of the
estimates in Theorem 3.1.

Collective observables such as $\Sigma_{\alpha\alpha^{\prime}}^{{}_{(\ell)}}$
are a convenient book-keeping device for studying the fully distribution of

$D^{\leq
r}z_{i;\mathcal{A}}^{(\ell+1)}:=\left(D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)},\quad\alpha\in\mathcal{A},\,J\in\mathbb{N}^{p},\,\left|J\right|\leq
r\right).$

Indeed, the conditional Gaussian structure (3.6) means the cumulants of
$D^{\leq r}z_{i;\mathcal{A}}^{{}_{(\ell+1)}}$ are easily expressible in terms
of those of
$D_{\alpha}^{J}D_{\alpha^{\prime}}^{J^{\prime}}\Sigma_{\alpha\alpha^{\prime}}^{{}_{(\ell)}}$.
For example, using Wick’s Theorem and the multi-linearity of cumulants reveals
that the fourth cumulant

$\displaystyle\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},D_{\alpha_{2}}^{J_{2}}z_{i_{2};\alpha_{2}}^{(\ell+1)},D_{\alpha_{3}}^{J_{3}}z_{i_{3};\alpha_{3}}^{(\ell+1)},D_{\alpha_{4}}^{J_{4}}z_{i_{4};\alpha_{4}}^{(\ell+1)}\right)$

equals

$\delta_{i_{1}i_{2}}\delta_{i_{3}i_{4}}\kappa_{(\alpha_{1}\alpha_{2})(\alpha_{3}\alpha_{4})}^{(J_{1}J_{2})(J_{3}J_{4}),(\ell)}+\delta_{i_{1}i_{3}}\delta_{i_{2}i_{4}}\kappa_{(\alpha_{1}\alpha_{3})(\alpha_{2}\alpha_{4})}^{(J_{1}J_{3})(J_{2}J_{4}),(\ell)}+\delta_{i_{1}i_{4}}\delta_{i_{2}i_{3}}\kappa_{(\alpha_{1}\alpha_{4})(\alpha_{2}\alpha_{3})}^{(J_{1}J_{4})(J_{2}J_{3}),(\ell)},$

where we’ve abbreviated

$\kappa_{(\alpha_{1}\alpha_{1}^{\prime}),\ldots,(\alpha_{k}\alpha_{k}^{\prime})}^{{(J_{1}J_{1}^{\prime}),\ldots,(J_{k}J_{k}^{\prime}),(\ell)}}:=\kappa\left(D_{\alpha_{1}}^{J_{1}}D_{\alpha_{1}^{\prime}}^{J_{1}^{\prime}}\Sigma_{\alpha_{1}\alpha_{1}^{\prime}}^{(\ell)},\ldots,D_{\alpha_{k}}^{J_{k}}D_{\alpha_{k}^{\prime}}^{J_{k}^{\prime}}\Sigma_{\alpha_{k}\alpha_{k}^{\prime}}^{(\ell)}\right).$
(3.7)

The general pattern is that odd mixed cumulants in the $z_{\alpha}^{(\ell)}$
and its derivatives vanish by symmetry and $2k$-th cumulants are finite sums
of cumulants
$\kappa_{(\alpha_{1}\alpha_{1}^{\prime}),\ldots,(\alpha_{k}\alpha_{k}^{\prime})}^{{}_{(J_{1}J_{1}^{\prime}),\ldots,(J_{k}J_{k}^{\prime}),(\ell)}}$.
In Corollary 3.4 and Proposition 10.2 below we will obtain a range of
recursions for these cumulants and thereby implicitly for the cumulants of the
original field $z_{\alpha}^{{}_{(\ell)}}$ as well.

To obtain such recursions note that, at first glance, the estimate (3.4) only
gives the order of magnitude in powers of $1/n$ of the cumulants of
$\kappa(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell)},\ldots,D_{\alpha_{k}}^{J_{k}}z_{i_{k};\alpha_{k}}^{{}_{(\ell)}})$
but does not seem to provide information about their structural dependence on
the remaining model parameters $C_{b},C_{W},\sigma,\ell$. However, such
information can be obtained by combining (3.4) with an additional argument. To
state the result, let us agree to write

$\kappa_{\alpha_{1}\alpha_{2}}^{(\ell)}:=\mathrm{Cov}\left(z_{i;\alpha_{1}}^{(\ell)},z_{i;\alpha_{2}}^{(\ell)}\right),\qquad\qquad\kappa_{\alpha_{1}\alpha_{2}}^{J_{1}J_{2},(\ell)}:=D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\kappa_{\alpha_{1}\alpha_{2}}^{(\ell)}.$
(3.8)

Physicists might refer to $\kappa_{\alpha_{1}\alpha_{2}}^{{}_{(\ell)}}$ as a
dressed two point function. Moreover, let us denote by
$\left\langle\cdot\right\rangle_{\kappa^{(\ell)}}$ the expectation with
respect to a collection of centered jointly Gaussian random vectors

$D_{\mathcal{A}}^{\leq
r}z_{i}=\left(D_{\alpha}^{J}z_{i;\alpha},\,\alpha\in\mathcal{A},\,\left|J\right|\leq
r\right)$

with the same covariance

$\mathrm{Cov}\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}},\,D_{\alpha_{2}}^{J_{2}}z_{i_{2};\alpha_{2}}\right):=\mathrm{Cov}\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell)},\,D_{\alpha_{2}}^{J_{2}}z_{i_{2};\alpha_{2}}^{(\ell)}\right)=\delta_{i_{1}i_{2}}\kappa_{\alpha_{1}\alpha_{2}}^{J_{1}J_{2},(\ell)}$

as the true vectors of derivatives $D_{\mathcal{A}}^{\leq
r}z_{i;\mathcal{A}}^{{}_{(\ell)}}$ in each component separately but zero
covariance for different $i$.

######  Theorem 3.2 (perturbative expansions of expectations of observables at
finite width).

Fix $r\geq 0$ and suppose that $f$ is both a continuous function and a
tempered distribution. Then for any $q_{*}\geq 0$ we have

$\displaystyle\mathbb{E}\left[f\left(D^{\leq
r}z_{1,\mathcal{A}}^{(\ell+1)},\ldots,D^{\leq
r}z_{m,\mathcal{A}}^{(\ell+1)}\right)\right]$
$\displaystyle\quad=\sum_{q=0}^{2q_{*}}\frac{(-1)^{q}}{2^{q}q!}\bigg{\langle}\mathbb{E}\bigg{[}\bigg{(}\sum_{\begin{subarray}{c}\left|J\right|,\left|J^{\prime}\right|\leq
r\\\
\alpha,\alpha^{\prime}\in\mathcal{A}\end{subarray}}\Delta_{\alpha\alpha^{\prime}}^{JJ^{\prime},(\ell)}\sum_{j=1}^{m}\partial_{D_{\alpha}^{J}z_{j;\alpha}}\partial_{D_{\alpha^{\prime}}^{J^{\prime}}z_{j;\alpha^{\prime}}}\bigg{)}^{q}\bigg{]}f\left(D_{\mathcal{A}}^{\leq
r}z_{1},\ldots,D_{\mathcal{A}}^{\leq
r}z_{m}\right)\bigg{\rangle}_{\kappa^{(\ell+1)}}$
$\displaystyle\quad+O(n^{-q_{*}-1}),$ (3.9)

where the sum is over multi-indices $J,J^{\prime}\in\mathbb{N}^{p}$ of order
at most $r$, we’ve set

$\Delta_{\alpha\alpha^{\prime}}^{JJ^{\prime},(\ell)}:=D_{\alpha}^{J}D_{\alpha^{\prime}}^{J^{\prime}}\Sigma_{\alpha\alpha^{\prime}}^{(\ell)}-\mathbb{E}\left[D_{\alpha}^{J}D_{\alpha^{\prime}}^{J^{\prime}}\Sigma_{\alpha\alpha^{\prime}}^{(\ell)}\right],$
(3.10)

and the derivatives $\partial_{D_{\alpha}^{J}z_{j;\alpha}}$ are interpreted in
the weak sense if $f$ is not differentiable.

We prove Theorem 3.2 in §9 and give the main idea of the proof in §5. By
substituting various polynomials in $z_{\alpha}^{{}_{(\ell+1)}}$ and its
derivatives for $f$ into the perturbative expansion (3.9), it is now
straightforward in principle to deduce recursions for the cumulants
$\kappa_{(\alpha_{1}\alpha_{1}^{\prime}),\ldots,(\alpha_{k}\alpha_{k}^{\prime})}^{{}_{(J_{1}J_{1}^{\prime}),\ldots,(J_{k}J_{k}^{\prime}),(\ell+1)}}$
at layer $\ell+1$ in terms of objects of the same type at layer $\ell$. In
particular, we have the following

######  Corollary 3.3 (cumulants in random neural networks form a hierarchy to
leading order in $1/n$).

With the assumptions of Theorem 3.1,the mixed cumulant

$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},\ldots,D_{\alpha_{2k}}^{J_{2k}}z_{i_{2k};\alpha_{2k}}^{(\ell+1)}\right)$

equals

$\sum_{\begin{subarray}{c}j\leq k\\\ J_{i}^{\prime},\,i=1,\ldots,2j\\\
\left|J_{1}^{\prime}\right|+\cdots+\left|J_{2j}^{\prime}\right|\\\
\quad\leq\left|J_{1}\right|+\cdots+\left|J_{2k}\right|\end{subarray}}C(J_{i}^{\prime},K_{\alpha_{i}\alpha_{i^{\prime}}}^{(\ell)}i,i^{\prime}=1,\ldots,2j)~{}\kappa\left(D_{\alpha_{1}}^{J_{1}^{\prime}}z_{1;\alpha_{1}}^{(\ell)},\ldots,D_{\alpha_{2j}}^{J_{2j}^{\prime}}z_{2j;\alpha_{2j}}^{(\ell)}\right)+O\left(n^{-k}\right),$
(3.11)

where the sum if over multi-indices $J_{i}^{\prime}$, the constants
$C(J_{i}^{\prime},K_{\alpha_{i}\alpha_{j}}^{(\ell)}i,j=1,\ldots,2k)$ depend
only on the multi-indices $J_{i}^{\prime}$ and the infinite width covariance
$K^{(\ell)}$, while the implicit constant in the error term depends on $k$,
the inputs $x_{\alpha_{1}},\ldots,x_{\alpha_{k}},$ the multi-indices
$J_{1},\ldots,J_{k}$, the weight and bias variances $C_{b},C_{W},$ the non-
linearity $\sigma$, and the layer index $\ell$.

We do not know how to efficiently compute the coefficients $C$ in the
recursion (3.11) for arbitrary $k$. Instead, we will compute them by hand for
small values of $k$ in the two settings:

  * •

We fix a single input $x_{\alpha}\in\mathbb{R}^{n_{0}}$ and obtain a recursion
for the cumulants:

$\kappa_{2k;\alpha}^{(\ell)}:=\kappa_{(\alpha\alpha)\cdots(\alpha\alpha)}^{(00)\cdots(00),(\ell)}=\frac{1}{(2k-1)!!}\kappa\bigg{(}\underbrace{z_{i;\alpha}^{(\ell+1)},\ldots,z_{i;\alpha}^{(\ell+1)}}_{2k\text{
times}}\bigg{)}$ (3.12)

when $k=2,3,4$ (the pairs $(\alpha\alpha)$ and $(00)$ both appear $k$ times).
See Corollary 3.4.

  * •

We again fix a single input $x_{\alpha}\in\mathbb{R}^{n_{0}}$ and consider any
two partial derivatives

$d_{j}:=\sum_{i=1}^{n_{0}}c_{i,j}\partial_{x_{j}}\bigg{|}_{x=x_{\alpha}},\qquad
j=1,2.$

We obtain recursions for

$\displaystyle\kappa_{(i_{1}i_{1}^{\prime})(i_{2}i_{2}^{\prime});\alpha}^{(\ell)}:=\kappa\left(d_{i_{1}}d_{i_{1}^{\prime}}\Delta_{\alpha\alpha}^{(\ell)},d_{i_{2}}d_{i_{2}^{\prime}}\Delta_{\alpha\alpha}^{(\ell)}\right).$
(3.13)

These cumulants are the building blocks for understanding the fourth mixed
cumulants of
$z_{i;\alpha}^{{}_{(\ell+1)}},d_{1}z_{i;\alpha}^{{}_{(\ell+1)}},d_{2}z_{i;\alpha}^{{}_{(\ell+1)}}$.
For example,

$\kappa\left(d_{1}z_{i;\alpha}^{(\ell+1)},d_{1}z_{i;\alpha}^{(\ell+1)},d_{2}z_{i;\alpha}^{(\ell+1)},d_{2}z_{i;\alpha}^{(\ell+1)}\right)=\kappa_{(11)(22);\alpha}^{(\ell)}+2\kappa_{(12)(12);\alpha}^{(\ell)}.$

This will be done in the course of proving Theorem 4.5. We refer the
interested reader to Proposition 10.2 for the precise recursions and to
Proposition 11.2 their solutions for non-linearities from the $K_{*}=0$
universality class (including $\tanh$).

In order to facilitate a compact form for the recursions described in first
bullet point, let us write

$T_{i,j;\alpha}^{(\ell)}:=C_{W}^{j}\left\langle\partial_{z}^{i}\left\\{\left(\sigma^{2}(z)-\left\langle\sigma^{2}(z)\right\rangle_{K_{\alpha\alpha}^{(\ell)}}\right)^{j}\right\\}\right\rangle_{K_{\alpha\alpha}^{(\ell)}},$
(3.14)

with the derivatives interpreted in the weak sense if $\sigma$ is not
sufficiently smooth and where we remind the reader of our standing notation

$\left\langle
f(z)\right\rangle_{K}=\int_{-\infty}^{\infty}f(zK^{1/2})e^{-\frac{z^{2}}{2}}\frac{dz}{\sqrt{2\pi}},\qquad
K\geq 0.$

######  Corollary 3.4.

Fix $r\geq 1$ and suppose that $\sigma:\mathbb{R}\rightarrow\mathbb{R}$
satisfies Assumption (2.1) with this value of $r$. Consider a depth $L$ random
neural network with input dimension $n_{0},$ hidden layer widths
$n_{1},\ldots,n_{L}$, output dimension $n_{L+1}$ and non-linearity $\sigma$.
Fix $x_{\alpha}\in\mathbb{R}^{n_{0}}$ and define

$\chi_{||;\alpha}^{(\ell)}:=\frac{1}{2}T_{2,1;\alpha}^{(\ell)}=\frac{C_{W}}{2}\left\langle\partial_{z}^{2}\sigma(z)^{2}\right\rangle_{K_{\alpha\alpha}^{(\ell)}},$

where the second derivative is interpreted in the weak sense if $\sigma$ is
not twice differentiable. For each $\ell=1,\ldots,L$, in the notation of
(3.12), the fourth cumulant satisfies

$\displaystyle\kappa_{4;\alpha}^{(\ell+1)}$
$\displaystyle=\frac{T_{0,2;\alpha}^{(\ell)}}{n_{\ell}}+\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}\kappa_{4,\alpha}^{(\ell)}+O(n^{-2}).$
(3.15)

Further, the $6$-th cumulant satisfies

$\displaystyle\kappa_{6;\alpha}^{(\ell+1)}$
$\displaystyle=\frac{T_{0,3;\alpha}}{n_{\ell}^{2}}+\frac{3T_{2,2;\alpha}^{(\ell)}}{2n_{\ell}}\chi_{||;\alpha}^{(\ell)}\kappa_{4;\alpha}^{(\ell)}-\frac{3T_{4,1;\alpha}^{(\ell)}}{8}\left(\chi_{||;\alpha}^{(\ell)}\kappa_{4;\alpha}^{(\ell)}\right)^{2}+\left(\chi_{||;\alpha}^{(\ell)}\right)^{3}\kappa_{6;\alpha}^{(\ell)}+O(n^{-3}).$
(3.16)

Finally, the $8$-th cumulant satisfies:

$\displaystyle\kappa_{8;\alpha}^{(\ell+1)}$
$\displaystyle=\frac{1}{n_{\ell}^{3}}\left(T_{0,4;\alpha}^{(\ell)}-3\left(T_{0,2;\alpha}^{(\ell)}\right)^{2}\right)$
$\displaystyle+\frac{1}{n_{\ell}^{2}}\left[2T_{2,3;\alpha}^{(\ell)}\chi_{||;\alpha}^{(\ell)}-12T_{0,2;\alpha}^{(\ell)}\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}+\frac{3}{2}\left(T_{2,2;\alpha}^{(\ell)}\right)^{2}-\frac{3}{2}T_{4,1;\alpha}^{(\ell)}T_{0,2;\alpha}^{(\ell)}\right]\kappa_{4;\alpha}^{(\ell)}$
$\displaystyle-\frac{1}{n_{\ell}}\left[2T_{2,2;\alpha}^{(\ell)}T_{4,1;\alpha}^{(\ell)}\chi_{||;\alpha}^{(\ell)}-\frac{1}{2}T_{4,2;\alpha}^{(\ell)}\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}+\left(\chi_{||;\alpha}^{(\ell)}\right)^{4}\right]\left(\kappa_{4;\alpha}^{(\ell)}\right)^{2}$
$\displaystyle+\frac{1}{n_{\ell}}\left[5T_{0,2;\alpha}^{(\ell)}T_{4,1;\alpha}^{(\ell)}\chi_{||;\alpha}^{(\ell)}+12T_{2,2;\alpha}^{(\ell)}\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}\right]\kappa_{6;\alpha}^{(\ell)}$
$\displaystyle+\frac{3}{32}\left(T_{4,1;\alpha}^{(\ell)}\right)^{2}\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}\left(\kappa_{4;\alpha}^{(\ell)}\right)^{3}-\frac{1}{2}\left(\chi_{||;\alpha}^{(\ell)}\right)^{3}T_{4,1;\alpha}^{(\ell)}\kappa_{4;\alpha}^{(\ell)}\kappa_{6;\alpha}^{(\ell)}$
$\displaystyle+\left(\chi_{||;\alpha}^{(\ell)}\right)^{4}\kappa_{8;\alpha}^{(\ell)}+O(n^{-4}).$
(3.17)

The initial condition for the recursions (3.15)-(3.17) is that
$\kappa_{2k;\alpha}^{{}_{(1)}}=0$ for all $k\geq 2$.

######  Remark 3.5.

Note that for $k=2,3,4$, we therefore see that to leading order in $1/n$ the
recursions for $\kappa_{2k;\alpha}^{{}_{(\ell+1)}}$ depends only on
$\kappa_{2j;\alpha}^{{}_{(\ell)}}$ for $j\leq k$. This allows us to interpret
(3.15) - (3.17) as a forming the start of hierarchy in powers of $1/n$
describing the cumulants of the output of a random neural network.

Let us briefly compare Corollary 3.4 to results in prior work:

  * •

In the special case when $\sigma$ is $1-$homogeneous (i.e. is linear, ReLU,
leaky ReLU, etc, see (4.5)), the full distribution of a neuron pre-activation
$z_{i;\alpha}^{{}_{(\ell)}}$ can be worked out in closed form. Namely, as we
explain in §4.2.1 and Appendix C, is has the same distribution as a Gaussian
with an independent random variance given by a product of independent weighted
chi-squared random variables. This was first pointed out in [17, 21] and
described in the language of special functions (namely Meijer G functions) in
[66]. For such non-linearities obtaining the recursions (3.15)-(3.17) is not
new.

  * •

The breakthrough work of Yaida [61] was the first to obtain, at a physics
level of rigor, the recursion (3.15) and probe its solutions at large $\ell$.

  * •

The ideas of Yaida in [61] then seeded the development in the monograph [53]
of Roberts, Yaida, and the author a much richer analysis, producing at a
physics level of rigor many recursions similar in flavor to (3.15)-(3.17) that
describe the the behavior of objects such as network derivatives at the start
of training, the NTK at the start of training, and even the change in the NTK
and the resulting output of a trained network. Many of these results go far
beyond what we are currently capable of doing mathematically.

  * •

The analysis in [53] never required studying cumulants
$\kappa_{2k;\alpha}^{{}_{(\ell)}}$ for $k\geq 3$, and while the techniques
developed there can certainly be used to obtain the recursions (3.16) and
(3.17) we take a rather different approach in this article that produces such
recursions more directly.

The functional $\chi_{||;\alpha}^{{}_{(\ell)}}$ plays a fundamental role in
the recursive description of random neural networks supplied by Corollary 3.4,
whose proof is in §9. In the following section we explain a principled
procedure, called tuning to criticality, that reveals its origin (as well as
that of a similar object we denote $\chi_{\perp;\alpha}^{{}_{(\ell)}}$) and
explains how to choose $C_{b},C_{W}$ so that these functionals are
approximately equal to $1$ at large $\ell$. As we will see, such a choice will
ensure that the recursions in Corollary 3.4 and their infinite width
counterpart (2.4) have well-behaved solutions at large $\ell$. We will then
return in §4.2.1 and §4.2.2 to solving the recursions from Corollary 3.4 in
random networks tuned to criticality (see (4.7) and Corollary 4.1).

##  4 Criticality and Universality in Random Neural Networks at Large Depth

In the previous section we presented two kinds of results about the structure
of random neural networks at large but finite width. The first, Theorem 3.1,
concerned the order of magnitude for cumulants of the output of such a random
network and its derivatives. The second, Theorem 3.2 and Corollary 3.4,
spelled out recursions with respect to the layer index $\ell$ that describe,
to leading order in $1/n$, network cumulants at layer $\ell+1$ in terms of
those at layer $\ell$. Our purpose in forthcoming sections is to analyze these
recursions at large $\ell$ and to apply this analysis to obtain results about
the structure of gradients in deep fully connected networks. Before doing
this, we must take a step back and ask: for which $\sigma,C_{b},C_{W}$ are the
recursions (2.4) describing the infinite width covariance $K^{(\ell)}$ well-
behaved at large $\ell$?

In section §4.1, we recall a more or less canonical answer to this question
whose roots are in the early articles [51, 52] and that was recently spelled
out in the generality presented here in [53]. This procedure, called tuning to
criticality, prescribes combinations of $\sigma,C_{b},C_{W}$ for which
$K^{(\ell)}$ is indeed well-behaved at large $\ell$. As we shall see below,
the term criticality is meant to be evocative of it’s use in the analysis of
2d systems in statistical mechanics in that tuning to criticality consists of
choosing $C_{b},C_{W}$ so that the infinite width covariance function
$K^{{(\ell)}}$ is as close to constant as a function $\ell$ as possible.

At a high level, there are two reasons to ask that $K^{(\ell)}$ be slowly
varying as a function of $\ell$. First, it arguably only makes sense to study
perturbative corrections in $1/n$ recursively in $\ell$ if the limiting
$n\rightarrow\infty$ covariance structure does not change too rapidly between
consecutive layers. Second, and perhaps more importantly, as explained and
thoroughly validated in [50, 52], deep fully connected networks (without
residual connections [24], batch normalization [31], etc) are numerically
stable enough for gradient-based optimization to succeed only if they are
tuned to criticality.

We discuss in §4.2 how considerations underlying criticality naturally give
rise to a notion of universality classes for random neural networks. Even the
correct definition of universality is still not fully understood. Unlike in
random matrix theory, universality for random neural networks depends not on
the statistics of the individual weights and biases (though this is also an
interesting direction to consider e.g. [19]) but rather on the effect of the
non-linearity $\sigma$ on the behavior of the infinite width covariance
$K^{(\ell)}$ at large values of the depth $\ell$.

Before giving the details, we take this opportunity to emphasize, as we have
elsewhere, that the definitions of criticality and universality, the approach
to solving the recursions for $\kappa_{2k;\alpha}^{{}_{(\ell)}}$ from
Corollary 3.4, and the resulting lessons learned about the role of the
effective network depth $L/n$ closely follow the ideas developed in the
monograph [53]. Though we pursue them in a somewhat different way, the author
would nonetheless like to acknowledge that his co-authors Dan Roberts and Sho
Yaida in the book deserve significant credit.

###  4.1 Tuning to Criticality

As originally explained in [51, 52] and recently spelled out in a definitive
way in [53], tuning a neural network to criticality means seeking choices of
$(C_{b},C_{W})$ that lead to critical fixed points of the form
$(K_{*},K_{*},K_{*})$ for the recursion (2.4), viewed as a dynamical system
describing
$(K_{\alpha\alpha}^{{}_{(\ell)}},K_{\beta\beta}^{{}_{(\ell)}},K_{\alpha\beta}^{{}_{(\ell)}})$
with time parameter $\ell$. Specifically, criticality requires

$\displaystyle\exists K_{*}\geq 0\qquad\text{s.t.}\qquad$ $\displaystyle
K_{*}=C_{b}+C_{W}\left\langle\sigma^{2}(z)\right\rangle_{K_{*}}$ ($*$)
$\displaystyle\forall\ell\geq 1\qquad$ $\displaystyle\frac{\partial
K_{\alpha\alpha}^{(\ell)}}{\partial
K_{\alpha\alpha}^{(1)}}\bigg{|}_{K_{\alpha\alpha}^{(1)}=K_{*}}=1$ ($||$)
$\displaystyle\forall\ell\geq 1\qquad$ $\displaystyle\frac{\partial
K_{\alpha\beta}^{(\ell)}}{\partial
K_{\alpha\beta}^{(1)}}\bigg{|}_{K_{\alpha\alpha}^{(1)}=K_{\alpha\alpha}^{(1)}=K_{\alpha\beta}^{(1)}=K_{*}}=1,$
($\perp$)

where

$K_{\alpha\alpha}^{(1)}=C_{b}+C_{W}K_{\alpha\beta}^{(0)},\qquad
K_{\alpha\beta}^{(0)}:=\frac{1}{n_{0}}\sum_{j=1}^{n_{0}}x_{j;\alpha}x_{j;\beta},\qquad
x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{0}}.$

The intuitive meaning of these conditions is as follows. Due to Theorem 2.2,
the first guarantees the existence of a fixed point $K_{*}$ for of the
recursion

$K_{\alpha\alpha}^{(\ell+1)}=C_{b}+C_{W}\left\langle\sigma^{2}(z)\right\rangle_{K_{\alpha\alpha}^{(\ell)}}$
(4.1)

of the infinite width variance. In particular, ($*$ ‣ 4.1) implies

$K_{\alpha\alpha}^{(1)}=C_{b}+\frac{C_{W}}{n_{0}}\left|\left|x_{\alpha}\right|\right|^{2}=K_{*}\qquad\Longrightarrow\qquad
K_{\alpha\alpha}^{(\ell)}=\lim_{n\rightarrow\infty}\mathrm{Var}\left[z_{i;\alpha}^{(\ell)}\right]=K_{*}\quad\forall\ell\geq
1.$

Thus, if a network is tuned to criticality, there is a critical radius

$K_{\text{crit}}^{2}:=n_{0}C_{W}^{-1}(K_{*}-C_{b})$

such that for inputs $x_{\alpha}$ on the sphere of radius $K_{\text{crit}}$
the variance of $z_{i;\alpha}^{{}_{(\ell)}}$ is independent of $\ell$ in the
infinite width limit. In non-critical networks, we expect this variance to
either grow or decay exponentially in $\ell$, leading to numerical
instabilities. The second condition ($||$ ‣ 4.1) considers the infinite width
limit of the variance of $z_{i;\alpha}^{{}_{(\ell)}}$ for an input
$x_{\alpha}$ for which $K_{\alpha\alpha}^{{}_{(1)}}$ is close to $K_{*}$.
Specifically, condition ($||$ ‣ 4.1) requires for all $\ell\geq 1$ that

$K_{\alpha\alpha}^{(1)}=\mathrm{Var}[z_{i;\alpha}^{(1)}]=K_{*}+\delta
K\qquad\Longrightarrow\qquad
K_{\alpha\alpha}^{(\ell)}=\lim_{n\rightarrow\infty}\mathrm{Var}\left[z_{i;\alpha}^{(\ell)}\right]=K_{*}+\delta
K+o(\delta K).$

This guarantees that the fixed point $K_{*}$ of the recursion (4.1) is
critical and hence that for inputs near the sphere of radius $K_{\text{crit}}$
the variance of the resulting pre-activations $z_{i;\alpha}^{{}_{(\ell)}}$ is
approximately constant in $\ell$ at large $n$. The final condition ($\perp$ ‣
4.1) considers instead the covariance between two inputs on the sphere of
radius $K_{\text{crit}}$. Namely, given two nearby network inputs
$x_{\alpha},x_{\beta}\in\mathbb{R}^{n_{0}}$ with

$K_{\alpha\alpha}^{(1)}=K_{\beta\beta}^{(1)}=K_{*},\qquad
K_{\alpha\beta}^{(1)}=C_{b}+\frac{C_{W}}{n_{0}}\sum_{j=1}^{n_{0}}x_{j;\alpha}x_{j;\beta}=K_{*}-\delta
K,$

the third condition asks that

$K_{\alpha\beta}^{(\ell)}=\lim_{n\rightarrow\infty}\mathrm{Cov}\left(z_{i;\alpha}^{(\ell)},z_{i;\beta}^{(\ell)}\right)=K_{*}-\delta
K+o(\delta K),\quad\forall\ell.$

This ensures that the covariance between pre-activations
$z_{i;\alpha}^{{}_{(\ell)}}$ and $z_{i;\beta}^{{}_{(\ell)}}$ corresponding to
two nearby inputs on the $K_{\text{crit}}$-sphere are approximately
independent of $\ell$ at large $n$. A simple computation directly from the
recursion (2.4) shows that

$\displaystyle\chi_{||}(K)$ $\displaystyle:=\frac{\partial
K_{\alpha\alpha}^{(\ell+1)}}{\partial
K_{\alpha\alpha}^{(\ell)}}\bigg{|}_{K_{\alpha\alpha}^{(\ell)}=K}=\frac{C_{W}}{2}\left\langle\partial_{z}^{2}(\sigma^{2}(z))\right\rangle_{K}$
(4.2) $\displaystyle\chi_{\perp}(K)$ $\displaystyle:=\frac{\partial
K_{\alpha\beta}^{(\ell+1)}}{\partial
K_{\alpha\beta}^{(\ell)}}\bigg{|}_{K_{\alpha\alpha}^{(\ell)}=K_{\beta\beta}^{(\ell)}=K_{\alpha\beta}^{(\ell)}=K}=C_{W}\left\langle(\partial_{z}\sigma(z))^{2}\right\rangle_{K}.$
(4.3)

Hence, all together, tuning to criticality requires

$\boxed{K_{*}\geq 0\text{ s.t.
}K_{*}=C_{b}+C_{W}\left\langle\sigma^{2}(z)\right\rangle_{K_{*}}\qquad\text{and}\qquad\chi_{||}(K_{*})=\chi_{\perp}(K_{*})=1.}$
(4.4)

###  4.2 Universality Classes of Random Neural Networks: Two Examples In
Search of a General Definition

We now turn to discussing the notion of universality classes for random neural
networks. To start, recall from Theorem 2.2 that the behavior at large depth
$\ell$ of random fully connected neural networks at infinite width is
completely specified by the asymptotics of the limiting covariance function
$K^{{(\ell)}}$. Observe, moreover, that the coefficients in the recursions for
$k=2,3,4$ of the cumulants $\kappa_{2k;\alpha}^{{}_{(\ell+1)}}$ from Corollary
3.4, which by Theorem 3.1 determine the behavior of random neural networks at
finite width to the first four orders in $1/n$, are completely determined by
$\sigma$, the infinite width covariance $K^{(\ell)}$, and cumulants
$\kappa_{2j;\alpha}^{{}_{(\ell)}},\,j\leq k$. It is therefore in terms of the
large $\ell$ behavior of $K^{{(\ell)}}$ that we should hope to define
universality classes of random neural networks at large depth. At present it
is not clear what the correct general definition of such a universality class
should be. We content ourselves instead with studying two important classes of
examples.

####  4.2.1 The Universality Class of ReLU

The most popular non-linearities used in practice are positively homogeneous
of degree $1$, i.e. have the form

$\sigma(t)=(a_{-}{\bf 1}_{\left\\{t<0\right\\}}+a_{+}{\bf
1}_{\left\\{t>0\right\\}})t,\qquad a_{-},a_{+}\in\mathbb{R},\quad a_{-}\neq
a_{+},\quad a_{-}^{2}+a_{+}^{2}\neq 0.$ (4.5)

Such non-linearities include the ReLU ($a_{-}=0,a_{+}=1$) and the leaky ReLU
($a_{-}\in(0,1),a_{+}=1$). A direct computation, left to the reader, shows
that criticality is achieved if and only if

$K_{*}\geq 0\text{ is arbitrary}\quad\text{and}\quad
C_{b}=0,\,C_{W}=\frac{2}{a_{+}^{2}+a_{-}^{2}}.$

Thus, the first property of the ReLU universality class is that setting
$(C_{b},C_{W})=(0,2/(a_{+}^{2}+a_{-}^{2}))$ allows all non-negative $K_{*}$ to
satisfy ($*$). In fact, at criticality, a simple symmetrization argument shows
that the variance of neuron pre-activations is preserved exactly even at
finite width

$\mathrm{Var}\left[z_{i;\alpha}^{(\ell)}\right]=\mathrm{Var}\left[z_{i;\alpha}^{(1)}\right]=\frac{C_{W}}{n_{0}}\left|\left|x_{\alpha}\right|\right|^{2}\qquad\forall\ell,n_{0},\ldots,n_{\ell}\geq
1,\,x_{\alpha}\in\mathbb{R}^{n_{0}}$ (4.6)

and, relatedly, that we have

$\chi_{||;\alpha}^{(\ell)}:=\chi_{||}(K_{\alpha\alpha}^{(\ell)})=1=\chi_{\perp;}(K_{\alpha\alpha}^{(\ell)})=:\chi_{\perp;\alpha}^{(\ell)},\qquad\forall\ell\geq
1,\,x_{\alpha}\in\mathbb{R}^{n_{0}}.$

The remarkable property (4.6) is much stronger than the criticality condition
$(*)$, which requires only that this condition holds for some value of
$n_{0}^{-1}\left|\left|x_{\alpha}\right|\right|^{2}$ and only in the limit
when $n\rightarrow\infty$. It implies that the cumulant recursions from
Corollary 3.4 for $1-$homogeneous non-linearities have constant coefficients
and are therefore particularly simple to solve. For instance, we find at
leading order in $1/n$

$\displaystyle\kappa_{4;\alpha}^{(\ell+1)}$
$\displaystyle=\frac{C_{W}^{2}}{n_{\ell}}\left[\left\langle\sigma(z)^{4}\right\rangle_{K_{\alpha\alpha}^{(\ell)}}-\left\langle\sigma(z)^{2}\right\rangle_{K_{\alpha\alpha}^{(\ell)}}^{2}\right]+\left(\chi_{||;\alpha}^{(\ell)}\right)^{2}\kappa_{4;\alpha}^{(\ell)}$
$\displaystyle=\left(\frac{2}{(a_{+}^{2}+a_{-}^{2})n_{0}}\left|\left|x_{\alpha}\right|\right|^{2}\right)^{2}\left(6\frac{a_{+}^{4}+a_{-}^{4}}{(a_{+}^{2}+a_{-}^{2})^{2}}-1\right)\sum_{\ell^{\prime}=1}^{\ell}\frac{1}{n_{\ell^{\prime}}},$

which shows that while $\kappa_{4,\alpha}^{{}_{(\ell)}}$ is suppressed by one
power of $1/n$ relative to the infinite width variance
$K_{\alpha\alpha}^{{}_{(\ell)}}$, it also grows one order faster in $\ell$.
This illustrates an important and general theme: depth amplifies finite width
effects. It is the effective depth $\ell/n$ of neurons at layer $\ell$ that
measures the distance to the infinite width Gaussian regime.

Moreover, in the special setting of $1$-homogeneous non-linearities there is a
simple method for obtaining the full distribution of the pre-activation vector
$z_{\alpha}^{{}_{(\ell)}}$ at a single input at any finite values of
$n_{0},\ldots,n_{\ell}$. This was first pointed out in [17, 21, 66] and is
briefly reviewed in Appendix C. A key takeaway is that if we take the hidden
layer widths $n_{1}=\cdots=n_{L}=n$, then we have following convergence in
distribution to product of independent normal and log-normal random variables:

$\lim_{\begin{subarray}{c}n,L\rightarrow\infty\\\
L/n\rightarrow\xi\in[0,\infty)\end{subarray}}z_{i;\alpha}^{(L)}\quad\stackrel{{\scriptstyle
d}}{{=}}\quad\left(\frac{2\left|\left|x_{\alpha}\right|\right|^{2}}{(a_{+}^{2}+a_{-}^{2})n_{0}}\right)^{1/2}Z_{1}\exp\left[-\mu(\xi,a_{+},a_{-})+\sigma(\xi,a_{+},a_{-})Z_{2}\right],$
(4.7)

where

$\mu(\xi,a_{+},a_{-})=\sigma^{2}(\xi,a_{+},a_{-}):=\frac{\xi}{4}\left(6\frac{a_{+}^{4}+a_{-}^{4}}{(a_{+}^{2}+a_{+}^{2})^{2}}-1\right),\quad
Z_{1},Z_{2}\sim\mathcal{N}(0,1)\text{ iid}.$

The convergence (4.7) reveals that for a fixed input the distribution of the
output of a random with $1-$homogeneous non-linearities at large depth and
width depends in a simple way on the limiting effective network depth $\xi$.
This bolsters the claim that they are all part of the same universality class.
It also means that increasing the network depth $L$ drives it away from the
infinite width Gaussian behavior observed at $\xi=0$ and that the outputs of
deep and wide networks are not well-approximated by a Gaussian at all, unless
$\xi$ is infinitesimal, in which case the log-normal term
$\exp\left[-\mu(\xi,a_{+},a_{-})+\sigma(\xi,a_{+},a_{-})Z_{2}\right]$ is
negligible.

Prior work [20, 21, 23] of the author shows that when $\sigma=\mathrm{ReLU}$
(or any other $1$-homogeneous non-linearity), the distribution at large $n,L$
of not only the network output $z_{i;\alpha}^{{}_{(L+1)}}$ but also is
derivatives with respect to inputs $x_{\alpha}$ and model parameters (e.g.
weights and biases) depends only on the effective depth $L/n.$ We will return
to this theme in our discussion of exploding and vanishing gradients (Theorem
4.5). Finally, we note that it has also been observed that log-normal random
variables describe the structure of gradients in residual networks, even
during/after training [39].

To complete our discussion of the ReLU universality class, we make two final
remarks. First, a direct computation (reviewed briefly in Proposition C.1 of
Appendix C) shows that at criticality for any non-zero inputs
$x_{\alpha_{1}},x_{\alpha_{2}}\in\mathbb{R}^{n_{0}}$ with the same norm we
have

$\lim_{n\rightarrow\infty}\mathrm{Corr}\left(z_{i;\alpha_{1}}^{(\ell)},z_{i;\alpha_{2}}^{(\ell)}\right)=1-\frac{2(a_{+}-a_{-})^{2}}{3\pi(a_{+}^{2}+a_{-}^{2})}\ell^{-2}(1+o(1)).$
(4.8)

The power law exponent $2$ that appears in this estimate is common to all
$1-$homogeneous non-linearities and is another reason to believe they fall
into the same universality class. In contrast, this exponent equals one for
non-linearities in the $K_{*}=0$ universality class presented below. The
estimate (4.8) suggests that in order to define a double scaling limit
$n,L\rightarrow\infty$ and $L/n\rightarrow\xi$ in which the entire field
$x_{\alpha}\mapsto z_{\alpha}^{{}_{(L+1)}}$ is non-degenerate (rather than
just its value at a single input) one must rescale distances in the input
space to prevent the collapse of correlations otherwise guaranteed in (4.8).
We leave this as an interesting direction for future work.

####  4.2.2 The Universality Class of Hyperbolic Tangent

The second class of non-linearities we study is what [53] termed the $K_{*}=0$
universality class, which we take to mean non-linearities $\sigma$ such that

  * •

$\sigma$ is a smooth, odd function satisfying Assumption 2.1.

  * •

$\sigma$ satisfies

$\sigma_{1}\sigma_{3}<0,\qquad\sigma_{j}:=\frac{1}{j!}\frac{d^{j}}{dt^{j}}\bigg{|}_{t=0}\sigma(t).$
(4.9)

  * •

$K_{*}=0$ is the unique fixed point of equation ($*$ ‣ 4.1).

  * •

At criticality, for every non-zero network input
$x_{\alpha}\in\mathbb{R}^{n_{0}}$ and each $\delta\in(0,1)$ we have as
$L\rightarrow\infty$ that

$K_{\alpha\alpha}^{(L)}=\frac{1}{aL}\left(1+O(L^{-1+\delta})\right),$ (4.10)

where the implicit constant depends on $\delta$ and $x_{\alpha}$ and we’ve set

$a:=-6\frac{\sigma_{3}}{\sigma_{1}}.$

This specific value of $a$, which is positive by (4.9), is the only possible
candidate for decay of the form (4.10) that is consistent with the recursion
(2.4).

Some remarks are in order. First, if $K_{*}=0$ is the unique fixed point for
($*$ ‣ 4.1), then a simple computation shows that criticality is achieved if
and only if

$K_{*}=0,\qquad C_{b}=0,\qquad C_{W}=\sigma_{1}^{-2}.$ (4.11)

Next, our definition of the $K_{*}=0$ universality class does not make
apparent whether it is empty. As we will see in Proposition B.1, however, the
$K_{*}=0$ universality class is in fact quite large and contains for example
the hyperbolic tangent and more generally any non-linearity that is
$\tanh$-like in the sense that is smooth with $\sigma_{1}\neq 0$, has the
opposite sign as its second derivative

$\text{for almost every
}z,~{}\mathrm{sgn}\left(\sigma(z)\sigma^{\prime\prime}(z)\right)=-1,$

is sub-linear

$\exists C>0\,\text{ s.t. }\forall
z\in\mathbb{R}\,\left|\sigma(z)\right|\leq\left|\sigma_{1}z\right|,$

and is controlled by its first few non-zero Taylor series coefficients at $0$:

$\exists C\geq 0\text{ s.t. }\forall z\geq
0,\quad\sigma_{1}z+\sigma_{3}z^{3}\leq\sigma(z)\leq\sigma_{1}z+\sigma_{3}z^{3}+Cz^{4}.$

Further, by definition, for the $K_{*}=0$ universality class, the infinite
width variance $K_{\alpha\alpha}^{{}_{(\ell)}}$ of neuron pre-activations
$z_{i;\alpha}^{(\ell)}$ is qualitatively different from that of
$1-$homogeneous non-linearities. Indeed, $K_{\alpha\alpha}^{{}_{(L)}}$ depends
on $L$, decaying polynomially to $0$. Moreover, at large $L$, the value of
$K_{\alpha\alpha}^{{}_{(L)}}$ is independent of the initial condition
$K_{\alpha\alpha}^{{}_{(0)}}$ to leading order in $L$. As a final remark let
us point out that searching for non-linearities $\sigma$ so that $K_{*}=0$ at
criticality is quite natural. Indeed, for any $\sigma$ that is twice
differentiable, we have

$\chi_{||}(K)=\chi_{\perp}(K)+C_{W}\left\langle\sigma(z)\sigma^{\prime\prime}(z)\right\rangle_{K}$

Hence, if $K>0$, then

$\chi_{||}(K)=1,\,\chi_{\perp}(K)=1\qquad\Longrightarrow\qquad\left\langle\sigma(z)\sigma^{\prime\prime}(z)\right\rangle_{K}=0.$

But if $\sigma$ is a sigmoidal function such as $\tanh$, then
$\sigma(z)\sigma^{\prime\prime}(z)<0$ for all $z\neq 0$. Hence,
$\left\langle\sigma(z)\sigma^{\prime\prime}(z)\right\rangle_{K}=0$ can only
occur when $K=0$.

As in the monograph [53], let us now probe the role of network depth by
studying the large $L$ behavior of the cumulants
$\kappa_{2k;\alpha}^{{}_{(L)}},\,k=2,3,4,$ in networks with non-linearities
from the $K_{*}=0$ universality class tuned to criticality. Note that in
(4.10) the limiting behavior of the variance $K_{\alpha\alpha}^{{}_{(L)}}$
depends (mildly) on the non-linearity $\sigma$ in terms of is first few Taylor
coefficients at $0$. As we are about to see, however, the behavior of the
higher cumulants $\kappa_{2k;\alpha}^{{}_{(L)}},\,k=2,3,4,$ when normalized by
the appropriate power of $K_{\alpha\alpha}^{{}_{(L)}}$, is independent of
$\sigma$ at leading order in $n$ and $L$ and depends only on universal
constants and the effective network depth $L/n$.

######  Corollary 4.1.

Suppose $\sigma$ is a non-linearity from the $K_{*}=0$ universality class, and
consider a depth $L$ random neural network with input dimension $n_{0}$,
output dimension $n_{L+1}$, hidden layer widths satisfying

$n_{1},\ldots,n_{L}=n\gg 1,$

and non-linearity $\sigma$ that has been tuned to criticality as in (4.11).
Write $\xi=L/n$ and define the normalized cumulants

$\widehat{\kappa}_{2k;\alpha}^{(L)}:=\frac{\kappa_{2k;\alpha}^{(\ell)}}{(K_{2;\alpha}^{(\ell)})^{k}}.$

We have for $k=2,3,4$ that

$\displaystyle\widehat{\kappa}_{2k;\alpha}^{(L)}$
$\displaystyle=C_{2k}\xi^{k-1}\left(1+O(L^{-1})\right)+O(n^{-k}),$ (4.12)

where

$C_{4}=\frac{2}{3},\qquad C_{6}=\frac{28}{15},\qquad C_{8}=\frac{8756}{315}$

are some positive universal constants. The implicit constants in error terms
$O(L^{-1})$ depend on $\sigma,C_{b},C_{W}$ and constants in
$O(n^{-j}),\,j=2,3,4$ may depend in addition on $L$.

######  Remark 4.2.

Combining estimate (4.12) for $k=2$ with the definition (4.10) of the
$K_{*}=0$ universality class and (3.12) yields that in the setting of
Corollary 3.4 we have to first order in $1/n$ and to leading order in $1/L$
that

$\displaystyle\widehat{\mathrm{Var}}\left[\left(z_{1;\alpha}^{(\ell+1)}\right)^{2}\right]$
$\displaystyle=2(1+\xi),\qquad\mathrm{Corr}\left(\left(z_{1;\alpha}^{(\ell+1)}\right)^{2},\,\left(z_{2;\alpha}^{(\ell+1)}\right)^{2}\right)=\frac{2}{3}\xi,\qquad\xi:=\frac{L}{n},$

where
$\widehat{\mathrm{Var}}[X]=\mathrm{Var}[X]/\mathbb{E}\left[X^{2}\right]$.
Thus, both the fluctuations of a single neuron pre-activation and the
correlation between different neurons is controlled to first order in
$1/n,1/L$ by the effective network depth $\xi$.

We prove Corollary 3.4 in §11.2. The formula (4.12) derived in [61] for $k=2$
at a physics level of rigor. Moreover, Corollary 3.4 represents a partial
analog of the convergence result (4.7), which lends credence to the following

######  Conjecture 4.3 (Existence of Double Scaling Limit for Random Neural
Networks).

Consider a random depth $L$ neural network with input dimension $n_{0}$,
hidden layer widths

$n_{1},\ldots,n_{L}=n\gg 1,$

output dimension $n_{L+1}$ and non-linearity $\sigma$. Suppose further that
this network is tuned to criticality in the sense that (4.4) is satisfied. Fix
a non-zero network input $x_{\alpha}\in\mathbb{R}^{n_{0}}$ and write
$\xi=L/n$. For each $k\geq 1$ there exists $C_{2k}>0$ depending on the
universality class of $\sigma$ so that

$\frac{\kappa_{2k;\alpha}^{(L)}}{\left(K_{2;\alpha}^{(L)}\right)^{k}}=C_{2k}\xi^{k-1}+O\left(\xi^{k}\right).$

Moreover, for each $\xi\in[0,\infty)$ there exists a probability distribution
$\mathbb{P}_{\xi,\sigma}$ on $\mathbb{R}$, depending only on $\xi$ and
$\sigma$, such that in the double scaling limit

$n,L\rightarrow\infty,\qquad\frac{L}{n}\rightarrow\xi,$

the random variable $z_{i;\alpha}^{(\ell)}$ converges in distribution to a
random variable with law $\mathbb{P}_{\xi,\sigma}$.

###  4.3 Gradients in Random Neural Networks

We presented in §3 and §4 a range of results about random fully connected
neural networks. The high-level takeaway was that such networks are succinctly
described when $n,L$ are large by the effective network depth $\xi=L/n$. We
saw that when $\xi=0$ the network outputs are a Gaussian process (Theorem
2.2). When $\xi>0$, however, a different picture emerges. Namely, the higher
order cumulants of the network output grow like powers of $\xi$ (see Corollary
3.4 and Conjecture 4.3). In this section we continue to consider a random
neural network $x_{\alpha}\mapsto z_{\alpha}^{{}_{(L+1)}}$ and study for any
$\ell$ the partial derivatives

$\partial_{x_{j;\alpha}}z_{i;\alpha}^{(\ell)},\qquad j=0,\ldots,n_{0},$

where by definition when $j=0$ we set

$\partial_{x_{0;\alpha}}z_{i;\alpha}^{(\ell)}:=z_{i;\alpha}^{(\ell)}.$

In §10.2 \- §10.4 we obtain recursions with respect to $\ell$ for both the
infinite width covariances

$K_{(ij)}^{(\ell)}:=\lim_{n\rightarrow\infty}\mathrm{Cov}\left(\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{i;\alpha}},\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{j;\alpha}}\right),\quad i,j=1,2.$

and the fourth cumulants

$\displaystyle\kappa_{(ii)(jj)}^{(\ell)}$
$\displaystyle:=(1-\frac{2}{3}\delta_{ij})\kappa\left(\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{i;\alpha}},\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{i;\alpha}},\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{j;\alpha}},\frac{\partial
z_{i;\alpha}^{(\ell+1)}}{\partial x_{j;\alpha}}\right),\qquad i,j=0,1,2.$
$\displaystyle=\kappa\left(\frac{\partial z_{1;\alpha}^{(\ell+1)}}{\partial
x_{i;\alpha}},\frac{\partial z_{1;\alpha}^{(\ell+1)}}{\partial
x_{i;\alpha}},\frac{\partial z_{2;\alpha}^{(\ell+1)}}{\partial
x_{j;\alpha}},\frac{\partial z_{2;\alpha}^{(\ell+1)}}{\partial
x_{j;\alpha}}\right).$

These recursions are valid for any non-linearity and any choice of
$C_{b},C_{W}$. Then, in §11.2 we solve these recursions in the context of a
random neural network with a non-linearity from the $K_{*}=0$ universality
class tuned to criticality. We record several of these solutions in the
following result.

######  Theorem 4.4 (Second and Fourth Joint Cumulants for Values and
Derivatives of Random Neural Networks).

Fix $r\geq 1$ and suppose that $\sigma:\mathbb{R}\rightarrow\mathbb{R}$
satisfies Assumption 2.1 with this value of $r$. Fix also an $x_{\alpha}\neq
0\in\mathbb{R}^{n_{0}}$. Consider a random neural network
$z_{\alpha}^{{}_{(L+1)}}$ with input dimension $n_{0}$, output dimension
$n_{L+1}$, hidden layer widths $n_{1},\ldots,n_{L}=n$, and non-linearity
$\sigma$ belonging to the $K_{*}=0$ universality class that has been tuned to
criticality as in (4.11). We have for any $\delta\in(0,1)$

$\displaystyle K_{(00)}^{(\ell+1)}$
$\displaystyle=\frac{1}{a\ell}+O_{\delta}(\ell^{-2+\delta})$ (4.13)
$\displaystyle K_{(10)}^{(\ell+1)}$
$\displaystyle=\frac{C_{W}e^{-2\gamma}x_{1;\alpha}}{n_{0}\ell^{2}}\left(1+O(\ell^{-1})\right)$
(4.14) $\displaystyle K_{(11)}^{(\ell+1)}$
$\displaystyle=\frac{C_{W}e^{-\gamma}}{n_{0}\ell}\left(1+O(\ell^{-1})\right)+O\left(\frac{(x_{1;\alpha})^{2}}{n_{0}^{2}\ell^{3}}\right)$
(4.15) $\displaystyle K_{(12)}^{(\ell+1)}$
$\displaystyle=O\left(\frac{x_{1;\alpha}x_{2;\alpha}}{n_{0}^{2}\ell^{3}}\right),$
(4.16)

where $\gamma$ is the Euler-Mascheroni constant and the implicit constant in
the $O_{\delta}(\ell^{-2+\delta})$ error depends on $\delta$. Moreover,
denoting by $\xi=L/n$ the effective network depth, we have

$\displaystyle\frac{\kappa_{(11)(00)}^{(\ell+1)}}{K_{(11)}^{(\ell)}K_{(00)}^{(\ell)}}$
$\displaystyle=-\frac{1}{3}\xi(1+O(\ell^{-1}))+O(n^{-2})$ (4.17)
$\displaystyle\frac{\kappa_{(11)(11)}^{(\ell+1)}}{(K_{(11)}^{(\ell)})^{2}}$
$\displaystyle=\frac{8}{3}\xi(1+O(\ell^{-1}))+O(n^{-2})$ (4.18)
$\displaystyle\frac{\kappa_{(11)(22)}^{(\ell+1)}}{K_{(11)}^{(\ell)}K_{(22)}^{(\ell)}}$
$\displaystyle=\frac{2}{3}\xi(1+O(\ell^{-1}))+O(n^{-2}),$ (4.19)

where the implicit constants in the $O(\ell^{-1})$ error terms depend only on
$x_{\alpha},\sigma$ and the implicit constants in the $O(n^{-2})$ error terms
may depend in addition on $\ell$.

Theorem 4.4 reveals several interesting phenomena. First, we find from (4.13)
and (4.14) that at infinite width values and derivatives become uncorrelated
at large $L$:

$\lim_{L\rightarrow\infty}\lim_{n\rightarrow\infty}\mathrm{Corr}\left(z_{1;\alpha}^{(L+1)},\,\partial_{x_{1;\alpha}}z_{1;\alpha}^{(L+1)}\right)=0.$

Moreover, at finite width, the effective network depth $\xi$ controls the
distribution of gradients both in terms of their fluctuations at a single
neuron and their correlations between neurons. For instance, to leading order
in $n,L$:

$\displaystyle\frac{\mathrm{Var}\left[\left(\partial_{x_{1;\alpha}}z_{1;\alpha}^{(L+1)}\right)^{2}\right]}{\mathbb{E}\left[\left(\partial_{x_{1;\alpha}}z_{1;\alpha}^{(L+1)}\right)^{2}\right]}$
$\displaystyle=\frac{8}{3}\xi,\qquad\mathrm{Corr}\left(\left(\frac{\partial
z_{1;\alpha}^{(L+1)}}{\partial x_{1;\alpha}}\right)^{2},\left(\frac{\partial
z_{1;\alpha}^{(L+1)}}{\partial
x_{2;\alpha}}\right)^{2}\right)=\frac{2}{3}\xi.$

In the next section, we will apply such estimates to understand a simple
version of the so-called exploding and vanishing gradient problem.

###  4.4 Exploding and Vanishing Gradient Problem for First Layer Weights

On its face, Theorem 4.4 concerns only derivatives of
$z_{i;\alpha}^{{}_{(L+1)}}$ with respect to the network input $x_{\alpha}$.
However, it is the derivatives of $z_{i;\alpha}^{{}_{(L+1)}}$ with respect to
the model parameters (weights and biases) that are arguably more important. To
explain this more precisely, consider a fully connected neural network
$x_{\alpha}\mapsto z_{\alpha}^{{}_{(L+1)}}$ with input dimension $n_{0}$,
output dimension $n_{L+1}$, and hidden layer widths $n_{1},\ldots,n_{L}$. For
any $m\in\mathbb{N}$ let us set

$[m]:=\left\\{1,\ldots,m\right\\}$

and denote by

$\theta=\left(\theta_{\mu},\,\mu\in[\\#\text{params}]\right)=\left(W_{ij}^{(\ell)},b_{i}^{(\ell)},\,\ell\in[L+1],\,i\in[n_{\ell}],\,j\in[n_{\ell-1}]\right)$

the vector of its trainable parameters, where we’ve set

$\\#\text{params}:=\\#\text{weights and
biases}=\sum_{\ell=1}^{L+1}n_{\ell}\left(1+n_{\ell-1}\right).$

As we briefly discussed in §2.2 the typical use case for a network is to
optimize its parameters $\theta$ by some variant of gradient descent

$\theta_{\mu}(t+1)=\theta_{\mu}(t)-\eta\partial_{\mu}\mathcal{L}(\theta(t)),\qquad\mu\in[\\#\text{params}]$
(4.20)

starting from a random initialization $\theta(0)$ drawn from the distribution
(2.2). The goal of this procedure is to minimize a loss

$\mathcal{L}(\theta)=\mathcal{L}(z_{\mathcal{A}_{\text{train}}}^{(L+1)}(\theta)),\qquad
z_{\mathcal{A}_{\text{train}}}^{(L+1)}:=\left(z_{\alpha}^{(L+1)},\,\alpha\in\mathcal{A}_{\text{train}}\right),$

which often depends only on the network output and measures for each $\theta$
how well the resulting neural net function $x_{\alpha}\mapsto
z_{\alpha}^{(L+1)}(\theta)$ matches given training data
$\left\\{(x_{\alpha},f(x_{\alpha})),\,\alpha\in\mathcal{A}_{\text{train}}\right\\}$
for a function $f:\mathbb{R}^{n_{0}}\rightarrow\mathbb{R}^{n_{L+1}}$. For
example, we might have

$\mathcal{L}(\theta)=\frac{1}{\left|\mathcal{A}_{\text{train}}\right|}\sum_{\alpha\in\mathcal{A}_{\text{train}}}\left|\left|f(x_{\alpha})-z_{\alpha}^{(L+1)}\right|\right|_{2}^{2},$

though many different losses are used in practice. The parameter $\eta$ is
called the learning rate, or step size. An important and well-known [5, 17,
27] numerical stability issue that comes up in the course of using a first
order method such as (4.20) to optimize the parameters of a deep neural
network is called the exploding and vanishing gradient problem (EVGP).
Informally, the EVGP occurs when the update (4.20) is numerically unstable:

$\mathrm{EVGP}\qquad\longleftrightarrow\qquad\text{for many parameters
}\theta_{\mu}\quad\frac{\left|\eta\partial_{\mu}\mathcal{L}(\theta)\right|}{\left|\theta_{\mu}\right|}\approx
0\text{ or }\infty.$

The presence of the EVGP causes optimization to break down since if
$\left|\partial_{\mu}\mathcal{L}(\theta)\right|/\left|\theta_{\mu}\right|\approx
0\text{ or }\infty$, then the relative change in the parameter $\theta_{\mu}$
is either too small to be useful or so large that it amounts to a large random
step in the space of parameters and therefore is unlikely to decrease the
loss.

It has long been known that sufficiently deep fully connected neural networks
are prone to suffer from the EVGP [5, 27, 28]. Indeed, many important
practical innovations such as residual connections [24] were invented in part
to address such numerical issues, though by now they are seen as useful for
many other reasons. To see why deep enough neural nets tend to have
numerically unstable gradients note that since the loss only depends on the
network outputs on a finite training set
$z_{\mathcal{A}_{\text{train}}}^{{}_{(L+1)}}$, we may use the chain rule to
write

$\partial_{\mu}\mathcal{L}(\theta)=\sum_{\alpha\in\mathcal{A}_{\text{train}}}J_{\theta_{\mu}}z_{\alpha}^{(L+1)}J_{z_{\alpha^{(L+1)}}}\mathcal{L}(z_{\mathcal{A}_{\text{train}}}^{(L+1)}),$

where for any function $f$ we denote by $J_{x}f(x)$ its Jacobian with respect
to $x.$ Moreover, suppose that $\theta_{\mu}$ is a weight or bias in layer
$\ell_{0}$. Then, again by the chain rule, the Jacobian of the network output
with respect to $\mu$ equals

$J_{\theta_{\mu}}z_{\alpha}^{(L+1)}=J_{\theta_{\mu}}z_{\alpha}^{(\ell_{0})}J_{z_{\alpha}^{(\ell_{0})}}z_{\alpha}^{(L+1)}.$

Both of these terms may have large fluctuations. Perhaps most importantly, the
second term can be rewritten as a product

$J_{z_{\alpha}^{(\ell_{0})}}z_{\alpha}^{(L+1)}=\prod_{\ell=\ell_{0}}^{L}J_{z_{\alpha}^{(\ell)}}z_{\alpha}^{(\ell+1)},$

of $L-\ell_{0}+1$ layer-to-layer Jacobian matrices (which are random at the
start of training), with each term having the following explicit
representation:

$J_{z_{\alpha}^{(\ell)}}z_{\alpha}^{(\ell+1)}=D^{(\ell+1)}W^{{(\ell+1)}},\qquad
D^{(\ell+1)}:=\mathrm{Diag}\left(\sigma^{\prime}\left(z_{i;\alpha}^{(\ell+1)}\right),\,i=1,\ldots,n_{\ell+1}\right).$

Except in the simple case when $\sigma$ is the identity and we require
$W^{{(\ell+1)}}$ to be orthogonal, the top singular value of each layer-to-
layer Jacobian will tend to deviate from $1$, causing the norms of their
products to either grow or decay exponentially with $L$. This poor
conditioning for large matrix products is a key culprit behind the EVGP.

Often, though not always, the gradients $\partial_{\mu}\mathcal{L}$ have the
largest fluctuations at the start of training. Thus, the EVGP is primarily a
property of random neural networks, and the techniques in this article can be
used to shed light on when it occurs. At the start of training the entries of
the Jacobian $J_{z_{\alpha}^{(\ell_{0})}}z_{\alpha}^{{}_{(L+1)}}$ are
precisely the derivatives of a random neural network with depth $L-\ell_{0}+1$
with respect to its input. One can therefore hope to use Theorems 3.1 and 3.2
to address the following more precise formulation of the EVGP:

$\displaystyle\text{EVGP}\quad\longleftrightarrow\quad$ for a typical random
draw of weights and biases the size of the relative Jacobians of the network
output
$\displaystyle\qquad\qquad\qquad\qquad\qquad|J_{\theta_{\mu}}z_{j;\alpha}^{(L+1)}|/\left|\theta_{\mu}\right|$
$\displaystyle\text{has a large variance over $\mu$ (i.e. over network
parameters)}.$

This characterization of the EVGP says in essence that if we can choose only a
single learning rate $\eta$ for a group of otherwise indistinguishable
parameters, such as the weights in a single layer, then it make sense to set

$\eta=\left(\text{average relative gradient
}|J_{\theta_{\mu}}z_{j;\alpha}^{(L+1)}|/\left|\theta_{\mu}\right|\right)^{-1}.$

The EVGP will then occur when the size of relative gradients have a large
variance across parameters in a typical random draw of weights and biases, so
that for some parameters our setting of $\eta$ is too small whereas for others
it is too large.

We will leave the study of the EVGP in this general form as pertaining to
future work and consider here a special case. Specifically, we compute in
Theorem 4.5 below the average of the empirical variance of the squared
gradients $|J_{\theta_{\mu}}z_{j;\alpha}^{{}_{(L+1)}}|^{2}$ when
$\theta_{\mu}=W_{ij}^{{}_{(1)}}$ varies over all weights in the first layer.
By construction, the typical size of $\left|\theta_{\mu}\right|$ is the same
for all such $\mu$. The EVGP will thus occur for weights in the first layer if
and only if the empirical fluctuations over $\mu$ of the raw Jacobians
$|J_{{W_{ij}^{(1)}}}z_{j;\alpha}^{{}_{(L+1)}}|$ is large. To state the precise
result we fix $x_{\alpha}\neq 0\in\mathbb{R}^{n_{0}}$ and $q\in[n_{L+1}]$. We
then write

$\mathrm{Grad~{}Mean}^{(1)}:=\frac{1}{n_{0}n_{1}}\sum_{\begin{subarray}{c}1\leq
i\leq n_{1}\\\ 1\leq j\leq
n_{0}\end{subarray}}\left(\partial_{W_{ij}^{(1)}}z_{q;\alpha}^{(L+1)}\right)^{2}$

for the empirical mean of the squared gradients of $z_{q;\alpha}^{(L+1)}$ with
respect to the weights $W_{ij}^{{}_{(1)}}$ in the first layer and analogously

$\mathrm{Grad~{}Var}^{(1)}:=\frac{1}{n_{0}n_{1}}\sum_{\begin{subarray}{c}1\leq
i\leq n_{1}\\\ 1\leq j\leq
n_{0}\end{subarray}}\left(\partial_{W_{ij}^{(1)}}z_{q;\alpha}^{(L+1)}\right)^{4}-\left(\mathrm{Grad~{}Mean}^{(1)}\right)^{2}$

for the empirical variance of the squared gradients.

######  Theorem 4.5 (EVGP for First Layer Weights in Fully Connected
Networks).

Let $x_{\alpha}\mapsto z_{\alpha}^{(L+1)}$ be a random neural network with
input dimension $n_{0}$, output dimension $n_{L+1}$, hidden layer widths

$n_{1},\ldots,n_{L}=n\gg 1,$

and non-linearity $\sigma$. Suppose also that

  * •

$\sigma$ satisfies Assumption. 2.1

  * •

$\sigma$ belongs to the $K_{*}=0$ universality class in the sense described in
§4.1 as well as §4.2.2.

  * •

$C_{b},C_{W}$ are tuned to criticality in the sense described as in (4.11).

Then, for any non-zero $x_{\alpha}\in\mathbb{R}^{n_{0}}$ we have

$\frac{\mathbb{E}\left[\mathrm{Grad~{}Var}^{(1)}\right]}{\mathbb{E}\left[\mathrm{Grad~{}Mean}^{(1)}\right]^{2}}=C_{\sigma,\alpha}\bigg{(}1+\frac{8}{3}\xi+O(L^{-1})+O(n^{-1})+O_{L}(n^{-2})\bigg{)},$

where $\xi=\frac{L}{n}$ is the effective network depth and

$C_{\sigma,\alpha}:=3\frac{\left\langle(\sigma^{\prime})^{4}\right\rangle_{K_{\alpha\alpha}^{(1)}}}{\left\langle(\sigma^{\prime})^{2}\right\rangle_{K_{\alpha\alpha}^{(1)}}^{2}}\frac{n_{0}^{-1}\left|\left|x_{\alpha}\right|\right|_{4}^{4}}{(n_{0}^{-1}\left|\left|x_{\alpha}\right|\right|_{2}^{2})^{2}}-1$

is a strictly positive constant depending only on $x_{\alpha},\sigma$, the
implicit constants in the error terms $O(n^{-1}),O(L^{-1})$ depend only on
$x_{\alpha},\sigma$ whereas the implicit constant in the error term
$O_{L}(n^{-2})$ depends on $x_{\alpha},\sigma,L$.

######  Remark 4.6.

Theorem 4.5 shows that, at least to leading order in $1/n$, the exploding and
vanishing gradient problem occurs for first layer weights in a critically
tuned random fully connected with non-linearity from the $K_{*}=0$
universality class if and only if the effective network depth $L/n$ is large.

We prove Theorem 4.5 in §11. To put this result into the context, let us make
several remarks. First, the analog of Theorem 4.5 for $1-$homogeneous non-
linearities was derived in prior work [17, 21]. As we’ve alluded to before,
the behavior of a random fully connected neural network with such non-
linearities can be studied by more specialized methods that reveal not only
the leading order dependence in $1/n$ but actually the full dependence on
$L/n$ plus errors of size $L/n^{2}$ (see Appendix C). The conclusion, just as
for non-linearities from the $K_{*}=0$ universality class, is that the EVGP
occurs if and only if $L/n$ is large.

Second, Theorem 4.5 is the first time the EVGP has been characterized
mathematically in random fully connected networks for a broad class of non-
linearities beyond those that are $1-$homogeneous. However, as throughout, the
author would like to acknowledge the prior work [53] of Roberts, Yaida, and
the author which also studies the EVGP and shows, at a physics level of rigor,
that the variance over random initialization of a single squared gradient
$(\partial_{\mu}z_{q;\alpha}^{{}_{(L+1)}})^{2}$ also scales like $L/n$ at
large $L$. This is conceptually different from showing that in a typical
initialization the empirical variance of
$(\partial_{\mu}z_{q;\alpha}^{{}_{(L+1)}})^{2}$ over $\mu$ in the first layer
is large, which is what Theorem 4.5 establishes. Nonetheless, the conclusions
are the essentially the same and the basic techniques used to derive the
results are similar.

##  5 Overview of Proofs

In this section, we present the essential idea for how we analyze a random
fully connected neural network $x_{\alpha}\mapsto z_{\alpha}^{(L+1)}$ at
finite width. Our approach is based on the following structural properties:

  * •

The sequence of fields $z_{\alpha}^{(\ell)}$ is a Markov Chain with respect to
$\ell$.

  * •

Conditional on the sigma algebra $\mathcal{F}^{(\ell)}$ defined by
$z_{\alpha}^{(\ell)}$ is a Gaussian field with independent components
$z_{i;\alpha}^{(\ell+1)}$. See Lemma 7.1.

  * •

The variance of each component $z_{i;\alpha}^{(\ell+1)}$ depends on
$z_{\alpha}^{(\ell)}$ only through random variables of the form

$\mathcal{O}_{f}^{(\ell)}:=\frac{1}{n_{\ell}}\sum_{j=1}^{n_{\ell}}f(z_{i;\alpha}^{(\ell)}),$

which we refer to as collective observables. See (3.6).

  * •

Collective observables are self-averaging to a similar extent as if the random
variables $f(z_{i;\alpha}^{(\ell)})$ were independent in the sense that for
any $q\geq 0$, we have

$\mathbb{E}\left[\left(\mathcal{O}_{f}^{(\ell)}-\mathbb{E}\left[\mathcal{O}_{f}^{(\ell)}\right]\right)^{q}\right]=O_{q}\left(n^{-\lceil\frac{q}{2}\rceil}\right).$
(5.1)

See Theorem 7.3 and Lemma 7.5.

Let us briefly explain, mostly dispensing with rigor, how these four ideas
come together to obtain a recursive description of the distribution of the
field $z_{\alpha}^{(\ell+1)}$ in terms of that of $z_{\alpha}^{(\ell)}$. To
keep the notation to a minimum, we fix a network input $x_{\alpha}$ and focus
on describing the joint distribution of
$z_{i;\alpha}^{(\ell+1)},\,i=1,\ldots,m$. Extensions to multiple inputs and
derivatives proceed along very similar lines. Denoting by
$\xi=\left(\xi_{1},\ldots,\xi_{m}\right)$ dual variables, consider the
characteristic function

$p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-i\sum_{i=1}^{m}\xi_{i}z_{i;\alpha}^{(\ell+1)}\right]\right].$

Conditioning on $z_{\alpha}^{(\ell)}$ and using (3.6) allows us to write

$p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi\right|\right|^{2}\Sigma_{\alpha\alpha}^{(\ell)}\right]\right],$

where we remind the reader that

$\Sigma_{\alpha\alpha}^{(\ell)}=\mathrm{Var}\left[z_{i;\alpha}^{(\ell+1)}~{}\big{|}~{}\mathcal{F}^{(\ell)}\right]=C_{b}+\frac{C_{W}}{n_{\ell}}\sum_{j=1}^{n_{\ell}}\sigma(z_{j;\alpha}^{(\ell)})^{2}$

is a collective observable at the previous layer. Writing

$\kappa_{\alpha\alpha}^{(\ell)}:=\mathbb{E}\left[\Sigma_{\alpha\alpha}^{(\ell)}\right],\qquad\Delta_{\alpha\alpha}^{(\ell)}:=\Sigma_{\alpha\alpha}^{(\ell)}-\mathbb{E}\left[\Sigma_{\alpha\alpha}^{(\ell)}\right],$

we find

$p^{(\ell+1)}(\xi):=\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi\right|\right|^{2}\Delta_{\alpha\alpha}^{(\ell)}\right]\right]\exp\left[-\frac{1}{2}\left|\left|\xi\right|\right|^{2}\kappa_{\alpha\alpha}^{(\ell)}\right].$

The second term is precisely the characteristic function of a centered
$m$-dimensional Gaussian with iid components of variance
$\kappa_{\alpha\alpha}^{(\ell)}$. Moreover, at least heuristically, the first
term can be written

$\mathbb{E}\left[\exp\left[-\frac{1}{2}\left|\left|\xi\right|\right|^{2}\Delta_{\alpha\alpha}^{(\ell)}\right]\right]=\sum_{q\geq
0}\mathbb{E}\left[\left(\Delta_{\alpha\alpha}^{(\ell)}\right)^{q}\right]\frac{(-1)^{q}}{2^{q}q!}\left|\left|\xi\right|\right|^{2q}.$

The concentration estimates (5.1) ensure that this series converges. Moreover,
since the Fourier transform turns polynomials into derivatives we have

$-\left|\left|\xi\right|\right|^{2}=\text{ Laplacian in the variables
}z_{i;\alpha}^{(\ell+1)}.$

Hence, we obtain for any reasonable test function $f$ that

$\mathbb{E}\left[f(z_{i;\alpha}^{(\ell+1)},\,i=1,\ldots,m)\right]=\sum_{q=0}^{\infty}\frac{1}{2^{q}q!}\mathbb{E}\left[\left(\Delta_{\alpha\alpha}^{(\ell)}\right)^{q}\right]\left\langle\left(\sum_{i=1}^{m}\partial_{z_{i;\alpha}}^{2}\right)^{q}f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}},$

where $(z_{i;\alpha},\,i=1,\ldots,m)$ is a vector of iid centered Gaussians
with variance $\kappa_{\alpha\alpha}^{(\ell)}$. The concentration estimates
(5.1) ensure that this expression is a power series in $1/n$. In particular,

$\displaystyle\mathbb{E}\left[f(z_{i;\alpha}^{(\ell+1)},\,i=1,\ldots,m)\right]$
$\displaystyle=\left\langle
f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}$
(5.2)
$\displaystyle+\frac{\mathbb{E}\left[(\Delta_{\alpha\alpha}^{(\ell)})^{2}\right]}{8}\left\langle\left(\sum_{i=1}^{m}\partial_{z_{i;\alpha}}^{2}\right)^{2}f(z_{i;\alpha},\,i=1,\ldots,m)\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}+O(n^{-2}).$

This is the essence of Theorem 3.2. To derive usable recursions for cumulants
of $z_{i;\alpha}^{(\ell+1)}$, note for instance that, in the notation of
Corollary 3.4,

$\kappa_{4;\alpha}^{(\ell)}:=\frac{1}{3}\kappa\left(z_{i;\alpha}^{(\ell+1)},z_{i;\alpha}^{(\ell+1)},z_{i;\alpha}^{(\ell+1)},z_{i;\alpha}^{(\ell+1)}\right)=\mathbb{E}\left[(\Delta_{\alpha\alpha}^{(\ell)})^{2}\right].$

Writing

$X_{j}:=\sigma(z_{j;\alpha}^{(\ell+1)})^{2}-\mathbb{E}\left[\sigma(z_{j;\alpha}^{(\ell+1)})^{2}\right]$

we thus have

$\kappa_{\alpha\alpha}^{(\ell+1)}=\mathbb{E}\left[\left(\Delta_{\alpha\alpha}^{(\ell)}\right)^{2}\right]=\frac{C_{W}^{2}}{n_{\ell}}\mathbb{E}\left[X_{1}^{2}\right]+C_{W}^{2}\left(1-n_{\ell}^{-1}\right)\mathbb{E}\left[X_{1}X_{2}\right].$

Applying the expansion (5.2) to both these terms and a bit of algebra already
yields

$\displaystyle\kappa_{\alpha\alpha}^{(\ell+1)}$
$\displaystyle=\mathbb{E}\left[\left(\Delta_{\alpha\alpha}^{(\ell)}\right)^{2}\right]$
$\displaystyle=\frac{C_{W}^{2}}{n_{\ell}}\left(\left\langle\sigma^{4}\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}-\left\langle\sigma^{2}\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}^{2}\right)$
$\displaystyle+C_{W}^{2}\left(1-n_{\ell}^{-1}\right)\left(\left(\left\langle\sigma^{2}\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}-\mathbb{E}\left[\sigma(z_{1;\alpha}^{(\ell)})^{2}\right]\right)^{2}+\frac{1}{4}\left\langle\partial^{2}\sigma^{2}\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}^{2}\kappa_{4;\alpha}^{(\ell)}\right)+O(n^{-2}).$

A short argument supplied in §9 shows that
$\left\langle\sigma^{2}\right\rangle_{\kappa_{\alpha\alpha}^{(\ell)}}=\mathbb{E}\left[\sigma(z_{1;\alpha}^{(\ell)})^{2}\right]+O(n^{-1})$
and that we may replace $\kappa_{\alpha\alpha}^{(\ell)}$ by its infinite width
limit $K_{\alpha\alpha}^{(\ell)}$ in all remaining expectations at the cost of
an $O(n^{-1})$ error. This already yields the recursion (3.15) of Corollary
3.4.

##  6 Background

###  6.1 Properties of Cumulants

Recall that, given random variables $X_{1},\ldots,X_{k}$ on the same
probability space, we denote their mixed cumulant by

$\kappa\left(X_{1},\ldots,X_{k}\right):=i^{k}\frac{\partial^{k}}{\partial
t_{1}\cdots\partial
t_{k}}\bigg{|}_{t=0}\log\mathbb{E}\left[\exp\left[-i(t_{1}X_{1}+\cdots+t_{k}X_{k})\right]\right].$

In the following result, we recall the key properties of these mixed cumulants
that we will need.

######  Proposition 6.1 (See Theorem 2.3.1 in [7]).

Mixed cumulants satisfy the following properties.

  1. 1.

Suppose $X=\left(X_{1},\ldots,X_{k}\right)$ is a random vector with finite
moments of all orders. Then

$\displaystyle\kappa\left(X_{1},\ldots,X_{k}\right)=\sum_{\pi=\left(\pi_{1},\ldots\pi_{b}\right)}\kappa\left(\kappa\left(X_{\pi_{1}}~{}|~{}\mathcal{F}\right),\ldots,\kappa\left(X_{\pi_{b}}~{}|~{}\mathcal{F}\right)\right),$
(6.1)

where the sum is over all partitions $\pi$ of $[k]$ and for each
$a=1,\ldots,b$

$X_{\pi_{a}}:=\left(X_{i},\,i\in\pi_{a}\right).$

This is known as the law of total cumulance. See [6].

  2. 2.

Suppose $X=\left(X_{1},\ldots,X_{k}\right)$ is a random vector with finite
moments of all orders. When $X$ can be partitioned into two independent
subsets, the mixed cumulant $\kappa(X)$ vanishes. More precisely, suppose
$I\subseteq[k]$ and that $I,I^{c}\neq\emptyset.$ Then

$X_{I}:=\left(X_{i},\,i\in I\right)\perp X_{I^{c}}=\left(X_{i},\,i\not\in
I\right)\qquad\Longrightarrow\qquad\kappa\left(X_{1},\ldots,X_{k}\right)=0.$
(6.2)

  3. 3.

Mixed cumulants are multi-linear. More precisely if

$\left\\{X_{i,j},\,1\leq j\leq k,\,i\leq T_{j}\right\\}$

are random variables with finite moments defined on the same probability
space, then

$\kappa\left(\sum_{i_{1}=1}^{T_{1}}a_{i_{1},1}X_{i_{1},1},\ldots,\sum_{i_{k}=1}^{T_{k}}a_{i_{k},k}X_{i_{k},k}\right)=\sum_{i_{1}=1}^{T_{1}}\cdots\sum_{i_{k}=1}^{T_{k}}a_{i_{1},1}\cdots
a_{i_{k},k}\kappa\left(X_{i_{1},1},\ldots,X_{i_{k},k}\right)$ (6.3)

for any $a_{i,j}\in\mathbb{R}$.

  4. 4.

Suppose $X=\left(X_{1},\ldots,X_{j}\right)\sim\mathcal{N}(0,\Sigma)$ is a
centered Gaussian with covariance $\Sigma$. Then for any
$i_{1},\ldots,i_{k}\in[j]$

$\kappa\left(X_{i_{1}},\ldots,X_{i_{k}}\right)=\begin{cases}\Sigma_{i_{1}i_{2}},&\quad
k=2\\\ 0,&\quad\text{otherwise}\end{cases}.$ (6.4)

  5. 5.

Moments are polynomials in cumulants. Specifically, suppose
$X=\left(X_{1},\ldots,X_{k}\right)$ is a random vector with finite moments of
all orders. Then,

$\mathbb{E}\left[X_{1}\cdots
X_{k}\right]=\sum_{\pi=\left(\pi_{1},\ldots,\pi_{b}\right)}\prod_{a=1}^{b}\kappa\left(X_{\pi_{a}}\right),$
(6.5)

where the sum is over all partitions of $[k]$ and for each $a\in[b]$ we’ve
written

$X_{\pi_{a}}:=\left(X_{i},\,i\in\pi_{a}\right).$

  6. 6.

Cumulants are polynomials in moments. Specifically,

$\kappa\left(X_{1},\ldots,X_{k}\right)=\sum_{\pi=\left(\pi_{1},\ldots,\pi_{b}\right)}(-1)^{b-1}(b-1)!\prod_{a=1}^{b}\mathbb{E}\left[\prod_{i\in\pi_{a}}X_{i}\right],$
(6.6)

where the sum is over all partitions of $[k]$ and for each $a\in[b]$ we’ve
written

$X_{\pi_{a}}:=\left(X_{i},\,i\in\pi_{a}\right).$

###  6.2 Lemmata

In this section, we collect two simple auxiliary results that we will need.
The first is a Lemma for Solving Certain Recursions.

######  Lemma 6.2.

Fix $C_{1},C_{2},\psi>0$ satisfying

$C_{2}\geq 1,\quad\psi\neq C_{2}+1$

as well as $*\in\left\\{\leq,\,\geq\right\\}$. Suppose also that for each
$\ell\geq 0$ we have

$a_{\ell+1}~{}*~{}\xi_{\ell}+(1-\zeta_{\ell})a_{\ell},\qquad\zeta_{\ell}\in[0,1]$
(6.7)

with $a_{0}\in\mathbb{R}$ given and that there exist
$C_{1}^{\prime},C_{2}^{\prime}>0$ so that

$\left|\xi_{\ell}-C_{1}\ell^{-\psi}\right|\leq
C_{1}^{\prime}\ell^{-1-\psi},\qquad\left|\zeta_{\ell}-C_{2}\ell^{-1}\right|\leq
C_{2}^{\prime}\ell^{-2}.$

Then

$a_{\ell+1}*\frac{\ell^{1-\psi}}{1-\psi+C_{2}}\left(1+O(\ell^{-1})\right)+e^{-C_{2}\gamma}\ell^{-C_{2}}a_{0}\left(1+O(\ell^{-1})\right)$
(6.8)

where $\gamma$ is the Euler-Mascheroni constant and the implied constants
depend only $C_{1},C_{2},C_{1}^{\prime},C_{2}^{\prime}.$

###### Proof.

By unfolding the recursion (6.7) we find

$\displaystyle
a_{\ell+1}~{}*~{}\sum_{\ell^{\prime}=1}^{\ell}\xi_{\ell^{\prime}}\prod_{\ell^{\prime\prime}=\ell^{\prime}+1}^{\ell}\left(1-\zeta_{\ell^{\prime\prime}}\right)+a_{0}\prod_{\ell^{\prime\prime}=0}^{\ell}(1-\zeta_{\ell^{\prime\prime}}).$

We have

$\displaystyle\prod_{\ell^{\prime\prime}=1}^{\ell}(1-\zeta_{\ell^{\prime\prime}})$
$\displaystyle=\exp\left[\sum_{\ell^{\prime\prime}=1}^{\ell}\log\left(1-C_{2}(\ell^{\prime\prime})^{-1}+O(\ell^{-2})\right)\right]$
$\displaystyle=\exp\left[O(\ell^{-1})+\sum_{\ell^{\prime\prime}=1}^{\ell}-C_{2}(\ell^{\prime\prime})^{-1}\right]$
$\displaystyle=\exp\left[O(\ell^{-1})-C_{2}\log(\ell)-C_{2}\gamma\right]$
$\displaystyle=e^{-C_{2}\gamma}\ell^{-C_{2}}\left(1+O(\ell^{-1})\right).$

This gives the second term in (6.8). For the first term, we write

$\displaystyle\sum_{\ell^{\prime}=1}^{\ell}\xi_{\ell^{\prime}}\prod_{\ell^{\prime\prime}=\ell^{\prime}+1}^{\ell}\left(1-\zeta_{\ell}\right)$
$\displaystyle=\sum_{\ell^{\prime}=1}^{\ell}\xi_{\ell^{\prime}}\exp\left[\sum_{\ell^{\prime\prime}=\ell^{\prime}+1}^{\ell}\log\left(1-\zeta_{\ell}\right)\right]$
$\displaystyle=\sum_{\ell^{\prime}=1}^{\ell}\xi_{\ell^{\prime}}\exp\left[\sum_{\ell^{\prime\prime}=\ell^{\prime}+1}^{\ell}-C_{2}(\ell^{\prime\prime})^{-1}+O((\ell^{\prime\prime})^{-2})\right]$
$\displaystyle=\sum_{\ell^{\prime}=1}^{\ell}C_{1}(\ell^{\prime})^{-\psi}(1+O(\ell^{\prime})^{-1})\exp\left[-C_{2}\log\left(\frac{\ell}{\ell^{\prime}}\right)+O((\ell^{\prime})^{-1})\right]$
$\displaystyle=\ell^{-C_{2}}\sum_{\ell^{\prime}=1}^{\ell}C_{1}(\ell^{\prime})^{-\psi+C_{2}}(1+O(\ell^{\prime})^{-1})$
$\displaystyle=\frac{C_{1}}{1+C_{2}-\psi}\ell^{1-\psi}\left(1+O(\ell^{-1})\right).$

This completes the proof of (6.8). ∎

The second result we will need is the following simple Lemma about Gaussian
Integrals.

######  Lemma 6.3.

Fix $r\geq 1$, a $r\times r$ matrix $\Sigma$ and measurable function
$g:\mathbb{R}^{r}\rightarrow\mathbb{R}$ that is polynomially bounded:

$\exists k\geq 1\text{ s.t.
}\sup_{x\in\mathbb{R}^{k}}\left|\left(1+\left|\left|x\right|\right|\right)^{-k}g(x)\right|<\infty.$

If $X$ is a standard Gaussian random vector in $\mathbb{R}^{r}$, then the
function

$\Sigma\mapsto\mathbb{E}\left[g\left(\Sigma^{1/2}X\right)\right]$ (6.9)

is smooth on the open set of strictly positive definite $k\times k$ matrices.
Further, if $g$ is a smooth function and each of its derivatives is
polynomially bounded, then the map (6.9) is extends to a smooth function on
the closed set of positive semi-definite matrices and, in particular,

$\frac{\partial}{\partial\Sigma_{ij}}\mathbb{E}\left[g\left(\Sigma^{1/2}X\right)\right]=\mathbb{E}\left[(\partial_{i}\partial_{j}g)(\Sigma^{1/2}X)\right].$
(6.10)

###### Proof.

On the open set of strictly positive definite matrices, the Gaussian density

$\Sigma\mapsto\exp\left[-\frac{1}{2}x^{T}\Sigma^{-1}x-\frac{1}{2}\log\det(2\pi\Sigma)\right]$

is a smooth function of $\Sigma$ with derivatives that are polynomials in $x$
and the entries of $\Sigma,\Sigma^{-1}$. The assumption that $f$ is
polynomially bounded shows that we may differentiate under the integral sign
and see that that

$\mathbb{E}\left[g(\Sigma^{1/2}X)\right]=\int_{\mathbb{R}^{r}}g(x)\exp\left[-\frac{1}{2}x^{T}\Sigma^{-1}x-\frac{1}{2}\log\det(2\pi\Sigma)\right]dx$

is indeed a smooth function of $\Sigma$. Suppose instead that $g$ is a smooth
function and that it’s derivatives are all polynomially bounded. Suppose first
that $g$ is in fact a Schwartz function. Then, writing $\widehat{g}$ for its
Fourier transform we have

$\mathbb{E}\left[g(\Sigma^{1/2}X)\right]=\int_{\mathbb{R}^{r}}\widehat{g}(\xi)\exp\left[-\frac{1}{2}\xi^{T}\Sigma\xi\right]d\xi.$

Since $\widehat{g}$ is also Schwartz, we may differentiate under the integral
sign to obtain

$\frac{\partial}{\partial\Sigma_{ij}}\mathbb{E}\left[g(\Sigma^{1/2}X)\right]=-\int_{\mathbb{R}^{r}}\xi_{i}\xi_{j}\widehat{g}(\xi)\exp\left[-\frac{1}{2}\xi^{T}\Sigma\xi\right]d\xi=\mathbb{E}\left[\partial_{x_{i}}\partial_{x_{j}}\bigg{|}_{x=\Sigma^{1/2}X}g(x)\right].$
(6.11)

Finally, if $g$ is not Schwartz but is smooth with all derivatives being
polynomially bounded, we consider the convolution

$g_{\epsilon}(x):=(g*\psi_{\epsilon})(x),\qquad\psi_{\epsilon}(y)=\exp\left[-\frac{\left|\left|y\right|\right|^{2}}{2\epsilon}-\frac{1}{2}\log(2\pi\epsilon)\right].$

Then, $g_{\epsilon}$ is Schwartz for all $\epsilon>0$. Moreover, note that
$g_{\epsilon}(\Sigma^{1/2}x)$ is also polynomially bounded for any PSD matrix
$\Sigma.$ Specifically, for any fixed PSD matrix $\Sigma_{0}$ we have for any
$k\geq 1$

$\displaystyle\sup_{\epsilon\in[0,1]}\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\\\ \Sigma\text{
PSD}\end{subarray}}\sup_{x\in\mathbb{R}^{r}}\left|(1+\left|\left|x\right|\right|)^{-k}g_{\epsilon}(\Sigma^{1/2}x)\right|$
$\displaystyle\qquad=\sup_{\epsilon\in[0,1]}\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\\\ \Sigma\text{
PSD}\end{subarray}}\sup_{x\in\mathbb{R}^{r}}\left|(1+\left|\left|x\right|\right|)^{-k}\int_{\mathbb{R}^{r}}g(\Sigma^{1/2}(x-y))\psi_{\epsilon}(y)dy\right|$
$\displaystyle\qquad\leq\sup_{\epsilon\in[0,1]}\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\\\ \Sigma\text{
PSD}\end{subarray}}\sup_{x\in\mathbb{R}^{r}}\left\\{(1+\left|\left|x\right|\right|)^{-k}\int_{\mathbb{R}^{r}}\left(1+\left|\left|\Sigma^{1/2}(x-y)\right|\right|^{k}\right)\psi_{\epsilon}(y)dy\right\\}$
$\displaystyle\qquad<\infty,$ (6.12)

Note that there exists $K>0$ depending only $k,r,\Sigma_{0}$ so that

$\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\end{subarray}}\left|\left|\Sigma^{1/2}(x-y)\right|\right|^{k}\leq
K\left(1+\left|\left|\Sigma_{0}^{1/2}\right|\right|\right)^{k}(\left|\left|x\right|\right|^{k}+\left|\left|y\right|\right|^{k}).$

Hence, since

$\sup_{\epsilon\in[0,1]}\int_{\mathbb{R}^{r}}\left|\left|y\right|\right|^{k}\psi_{\epsilon}(y)dy<\infty$

we find that

$\displaystyle\sup_{\epsilon\in[0,1]}\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\\\ \Sigma\text{
PSD}\end{subarray}}\sup_{x\in\mathbb{R}^{r}}\left|(1+\left|\left|x\right|\right|)^{-k}g_{\epsilon}(\Sigma^{1/2}x)\right|<\infty.$
(6.13)

The estimate above allows us to use dominate convergence to see that for any
PSD $\Sigma$

$\mathbb{E}\left[g(\Sigma^{1/2}X)\right]=\lim_{\epsilon\rightarrow
0}\mathbb{E}\left[g_{\epsilon}(\Sigma^{1/2}X)\right].$ (6.14)

To complete the proof we note that $g_{\epsilon}$ and
$\partial_{i}\partial_{j}g_{\epsilon}$ are both Schwartz for any positive
$\epsilon$. Moreover,
$\partial_{i}\partial_{j}\partial_{k}\partial_{m}g_{\epsilon}$ satisfies
(6.13). Hence, we conclude by applying (6.11) that for any PSD matrix
$\Sigma_{0}$ there exists $C>0$ so that

$\displaystyle\sup_{\begin{subarray}{c}\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1\\\ \Sigma\text{
PSD}\end{subarray}}\sup_{\epsilon\in[0,1]}\frac{\left|\mathbb{E}\left[g_{\epsilon}(\Sigma^{1/2}X)\right]-\mathbb{E}\left[g_{\epsilon}(\Sigma_{0}^{1/2}X)\right]-\sum_{i,j=1}^{r}\mathbb{E}\left[(\partial_{i}\partial_{j}g_{\epsilon})(\Sigma_{0}^{1/2}X)\right]\left(\Sigma-\Sigma_{0}\right)_{ij}\right|}{\left|\left|\Sigma-\Sigma_{0}\right|\right|^{2}}$
$\displaystyle\qquad\leq\sup_{\epsilon\in[0,1]}\sup_{\left|\left|\Sigma-\Sigma_{0}\right|\right|\leq
1}\sum_{i,j,k,m=1,\ldots,r}\left|\mathbb{E}\left[(\partial_{i}\partial_{j}\partial_{k}\partial_{m})g_{\epsilon}(\Sigma^{1/2}X)\right]\right|$
$\displaystyle\qquad\leq C.$

Thus, if
$\Sigma-\Sigma_{0}/\left|\left|\Sigma-\Sigma_{0}\right|\right|\rightarrow\Sigma_{1}$,
we find by applying (6.14) to $\partial_{i}\partial_{j}g$ that

$\displaystyle\lim_{\Sigma\rightarrow\Sigma_{0}}\frac{\mathbb{E}\left[g(\Sigma^{1/2}X)\right]-\mathbb{E}\left[g(\Sigma_{0}^{1/2}X)\right]}{\left|\left|\Sigma-\Sigma_{0}\right|\right|}$
$\displaystyle=\lim_{\Sigma\rightarrow\Sigma_{0}}\lim_{\epsilon\rightarrow
0}\frac{\mathbb{E}\left[g_{\epsilon}(\Sigma^{1/2}X)\right]-\mathbb{E}\left[g_{\epsilon}(\Sigma_{0}^{1/2}X)\right]}{\left|\left|\Sigma-\Sigma_{0}\right|\right|}$
$\displaystyle=\lim_{\Sigma\rightarrow\Sigma_{0}}\lim_{\epsilon\rightarrow
0}\left\\{\sum_{i,j=1}^{r}\mathbb{E}\left[(\partial_{i}\partial_{j}g_{\epsilon}(\Sigma_{0}^{1/2}X))\right]\frac{(\Sigma-\Sigma_{0})_{ij}}{\left|\left|\Sigma-\Sigma_{0}\right|\right|}\right\\}$
$\displaystyle=\sum_{i,j=1}^{r}\mathbb{E}\left[\partial_{i}\partial_{j}g(\Sigma_{0}^{1/2}X)\right](\Sigma_{1})_{ij}.$

This shows that (6.10) holds for any $g$ that is smooth, completing the proof
of Lemma 6.3. ∎

##  7 Proof of Theorem 3.1

Let us first recall the notation. We fix $r\geq 1$ and assume that
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ satisfies assumption 2.1 with this
value of $r$. We also fix a finite collection
$x_{\mathcal{A}}=\left\\{x_{\alpha},\alpha\in\mathcal{A}\right\\}\subseteq\mathbb{R}^{n_{0}}$
of distinct network inputs and $p$ directional derivatives
$d_{1},\ldots,d_{p}$ as in (3.2). We denote by

$N(p,r)=\\#\left\\{J=\left(j_{1},\ldots,j_{p}\right)\in\mathbb{N}^{p}~{}|~{}j_{1}+\cdots+j_{p}\leq
r\right\\},$

which computes the number of possible derivatives of order at most $r$ in the
$p$ directional derivatives $d_{j}$. We also denote by $\mathcal{F}^{(\ell)}$
the sigma algebra generated by the weights and biases in layers up to and
including $\ell$. The starting point for proving Theorem 3.1 is the following
simple but fundamental observation.

######  Lemma 7.1.

For each $\ell\geq 0$, conditional on $\mathcal{F}^{(\ell)}$,

$\left\\{\left(D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)},\,\alpha\in\mathcal{A},\,\left|J\right|\leq
r\right)\right\\}_{i=1}^{n_{\ell+1}}$

is a collection of $n_{\ell+1}$ iid centered Gaussians of dimension $N(p,r)$.

###### Proof.

The defining recursion (2.1) of a fully connected network yields for each
$\alpha,J$

$\displaystyle D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)}$
$\displaystyle=D_{\alpha}^{J}\left\\{b_{i}^{(\ell+1)}+\sum_{j=1}^{n_{\ell}}W_{ij}^{(\ell+1)}\sigma\left(z_{j;\alpha}^{(\ell)}\right)\right\\}=\delta_{\left|J\right|=0}b_{i}^{(\ell+1)}+\sum_{j=1}^{n_{\ell}}W_{ij}^{(\ell+1)}D_{\alpha}^{J}\sigma\left(z_{j;\alpha}^{(\ell)}\right).$
(7.1)

Note that $D_{\alpha}^{J}\sigma(z_{j;\alpha}^{{}_{(\ell)}})$ are measurable
with respect to $\mathcal{F}^{(\ell)}$. The conclusion now follows since the
weights $W_{ij}^{{}_{(\ell+1)}},\,j=1\ldots,n_{\ell}$ and bias
$b_{i}^{{}_{(\ell+1)}}$ are centered Gaussians and are independent for
different $i$. ∎

Thus, the structure of $z_{\alpha}^{{}_{(\ell+1)}}$ and its derivatives is
always that of a Gaussian field after conditioning on
$\mathcal{F}^{{}_{(\ell)}}$. To ease the notation in what comes given
$f:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R}$, let us abbreviate

$f\left(z_{j;\mathcal{A}}^{(\ell)}\right):=f\left(D_{\alpha}^{J}z_{j;\alpha}^{(\ell)},\,\alpha\in\mathcal{A},\,\left|J\right|\leq
r\right),\qquad j\in[n_{\ell}].$

Next, we remind the reader that given
$f:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R}$, which is measurable and polynomially
bounded, the corresponding collective observable $\mathcal{O}_{f}^{(\ell)}$ at
layer $\ell$ is

$\mathcal{O}_{f}^{(\ell)}=\frac{1}{n_{\ell}}\sum_{j=1}^{n_{\ell}}f\left(z_{j;\mathcal{A}}^{(\ell)}\right)$

and that the statement (A.3) in Proposition A.1 ensures

$\sup_{n\geq
1}\mathbb{E}\left[\left|\mathcal{O}_{f}^{(\ell)}\right|\right]<\infty.$ (7.2)

Recall also our notation for the conditional covariance

$\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}:=\mathrm{Cov}\left(z_{i;\alpha_{1}}^{(\ell+1)},z_{i;\alpha_{2}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right)=C_{b}+\frac{C_{W}}{n_{\ell}}\sum_{j=1}^{n_{\ell}}\sigma\left(z_{j;\alpha_{1}}^{(\ell)}\right)\sigma\left(z_{j;\alpha_{2}}^{(\ell)}\right)$

and note that both it and its derivatives

$\displaystyle
D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}$
$\displaystyle=\mathrm{Cov}\left(D_{\alpha_{1}}^{J_{1}}z_{i;\alpha_{1}}^{(\ell+1)},D_{\alpha_{2}}^{J_{2}}z_{i;\alpha_{2}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right)$
$\displaystyle=D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\left(C_{b}+\frac{C_{W}}{n_{\ell}}\sum_{j=1}^{n_{\ell}}\sigma\left(z_{j;\alpha_{1}}^{(\ell)}\right)\sigma\left(z_{j;\alpha_{2}}^{(\ell)}\right)\right)$

are collective observables at layer $\ell$. Our first application of Lemma 7.1
is the following reduction of the study of cumulants of
$D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)}$ to the cumulants of certain collective
observables.

######  Proposition 7.2.

Fix $k,\ell\geq 1$ and $p$-dimensional multi-indices $J_{1},\ldots,J_{k}$ with
$\left|J_{i}\right|\leq r.$ If $k$ is odd, then

$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},\ldots,D_{\alpha_{k}}^{J_{k}}z_{i_{k};\alpha_{k}}^{(\ell+1)}\right)=0$

In contrast, if $k$ is even

$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},\ldots,D_{\alpha_{k}}^{J_{k}}z_{i_{k};\alpha_{k}}^{(\ell+1)}\right)~{}=~{}\text{
finite sums of
}\kappa\left(\mathcal{O}_{f_{1}}^{(\ell)},\ldots,\mathcal{O}_{f_{k/2}}^{(\ell)}\right),$

where $\mathcal{O}_{f_{j}}^{(\ell)}$ are collective observables of the form

$D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)},\qquad\left|J_{1}\right|,\left|J_{2}\right|\leq
r.$ (7.3)

###### Proof.

Using (6.1) and recalling that $\mathcal{F}^{(\ell)}$ is the sigma algebra
generated by weights and biases in layers up to and including $\ell$, we have
that
$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},\ldots,D_{\alpha_{k}}^{J_{k}}z_{i_{k};\alpha_{k}}^{(\ell+1)}\right)$
equals

$\displaystyle\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\end{subarray}}\kappa\left(\kappa\left(\left(D^{J}z^{(\ell+1)}\right)_{\pi_{1}}\big{|}\mathcal{F}^{(\ell)}\right),\ldots,\kappa\left(\left(D^{J}z^{(\ell+1)}\right)_{\pi_{B}}\big{|}\mathcal{F}^{(\ell)}\right)\right),$
(7.4)

where the sum is over partitions $\pi$ of $[k]$ and for $b=1,\ldots,B$ we’ve
abbreviated

$\left(D^{J}z^{(\ell+1)}\right)_{\pi_{b}}:=\left(D_{\alpha_{t}}^{J_{t}}z_{i_{t};\alpha_{t}}^{(\ell+1)},\quad
t\in\pi_{b}\right).$

By Lemma 7.1,
$\\{(D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)},\,\alpha\in\mathcal{A},\,\left|J\right|\leq
d),\,i=1,\ldots,n_{\ell+1}\\}$ are iid centered Gaussians conditional on
$\mathcal{F}^{(\ell)}$. Hence, by the properties (6.2) and (6.3) and (6.4)
from Proposition 6.1, in the sum over partitions above, a term is non-zero
only if

$\forall b\in[B],\qquad\left|\pi_{b}\right|=2\quad\text{and}\quad
i_{\pi_{b}(1)}=i_{\pi_{b}(2)}$

This proves that
$\kappa\left(D^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},\ldots,D^{J_{k}}z_{i_{k};\alpha_{k}}^{(\ell+1)}\right)$
vanishes if $k$ is odd. To treat the case when $k$ is even observe that by
(7.1)

$\displaystyle\kappa\left(D^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},D^{J_{2}}z_{i_{2};\alpha_{2}}^{(\ell+1)}\big{|}\mathcal{F}^{(\ell)}\right)$
$\displaystyle=\delta_{i_{1}i_{2}}D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}.$

Substituting this into (7.4) completes the proof. ∎

When $k=2$, Proposition 7.2 and our assumption (2.1) shows that for each
$\ell\geq 1$, any $i_{1},i_{2}\in[n_{\ell+1}],\,\alpha\in\mathcal{A},$ and
multi-indices $J_{1},J_{2}$ of order at most $d$, there exists a polynomially
bounded function $f:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},d)}\rightarrow\mathbb{R}$ for which

$\kappa\left(D_{\alpha_{1}}^{J_{1}}z_{i_{1};\alpha_{1}}^{(\ell+1)},D_{\alpha_{2}}^{J_{2}}z_{i_{2};\alpha_{2}}^{(\ell+1)}\right)=\mathbb{E}\left[\mathcal{O}_{f}^{(\ell)}\right]$

In light of (7.2) this proves Theorem 3.1 when $k=2$. Further, since the
cumulant of $2$ or more random variables is shift-invariant, we may assume for
$k\geq 3$ that the collective observables
$D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}$
in Proposition 7.2 are replaced by their zero mean versions:

$\Delta_{\alpha_{1}\alpha_{2}}^{J_{1},J_{2},(\ell)}:=D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}-\mathbb{E}\left[D_{\alpha_{1}}^{J_{1}}D_{\alpha_{2}}^{J_{2}}\Sigma_{\alpha_{1}\alpha_{2}}^{(\ell)}\right].$
(7.5)

Hence, Theorem 3.1 is a special case of the following result.

######  Theorem 7.3.

Fix $k,m\geq 1$. Consider any $m$-tuple $F=\left(f_{1},\ldots,f_{m}\right)$
consisting of measurable, functions

$f_{i}:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R},\quad i=1,\ldots,m\qquad$

that are polynomially bounded and satisfy

$\mathbb{E}\left[\mathcal{O}_{f_{i}}^{(\ell)}\right]=\mathbb{E}\left[f_{i}\left(z_{1;\mathcal{A}}^{(\ell)}\right)\right]=0,\qquad
i=1,\ldots,m.$

Define the $m-$tuple of collective observables

$\overrightarrow{\mathcal{O}}_{F}^{(\ell)}:=\left(\mathcal{O}_{f_{i}}^{(\ell)},\,i=1,\ldots,m\right).$

Consider further any measurable polynomially bounded functions

$g_{j}:\mathbb{R}^{m}\rightarrow\mathbb{R},\quad j=1,\ldots,k.$

which are smooth in a neighborhood of $0$. If $f_{i}$ and $\sigma$ are in fact
smooth, then, for every $\ell\geq 1$

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(g_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell)}\right),\ldots,g_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell)}\right)\right)\right|<\infty$
(7.6)

Moreover, (7.6) holds without the assumption that $f_{i},\sigma$ are smooth
provided that for each $\ell$ the vector of iterated directional derivatives
$(D_{\alpha}^{J}z_{i;\alpha}^{(\ell)},\,\left|J\right|\leq
r,\alpha\in\mathcal{A})$ of order at most $r$ is non-degenerate in the sense
of (3.3).

###### Proof.

Our starting point is a reduction of Theorem 7.3 to the case when $g_{j}$ are
polynomials. This is related to a technique called the delta method in some
parts of the mathematical statistics literature [59].

######  Proposition 7.4 (Polynomials are Enough for Theorem 7.3).

Fix $m\geq 1$ and suppose that for each $n\geq 1$ we have an $m-$tuple
$X_{n}=\left(X_{n,1},\ldots,X_{n,m}\right)$ of mean $0$ random variables that
possess bounded moments of all orders:

$\sup_{n\geq 1}\left|\mathbb{E}\left[X_{n,1}^{q_{1}}\cdots
X_{n,m}^{q_{m}}\right]\right|<\infty,\qquad\forall\,q_{1},\ldots,q_{m}\geq 0.$
(7.7)

Suppose for any given polynomials $p_{1},\ldots,p_{k}$ in $m$ variables we
have

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(p_{1}(X_{n}),\ldots,p_{k}(X_{n})\right)\right|<\infty.$
(7.8)

Then, for any measurable, polynomially bounded functions
$g_{j}:\mathbb{R}^{m}\rightarrow\mathbb{R},\,j=1,\ldots,k,$ which are smooth
in some fixed neighborhood of $0$

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(g_{1}\left(X_{n}\right),\ldots,g_{k}\left(X_{n}\right)\right)\right|<\infty.$
(7.9)

###### Proof.

We begin with the following simple Lemma, which allows us to translate between
the cumulants bounds (7.8) and high probability bounds.

######  Lemma 7.5.

For any $q\geq 1$

$\sup_{n\geq 1}\sup_{1\leq i\leq
m}\left|n^{\lceil\frac{q}{2}\rceil}\mathbb{E}\left[X_{i,n}^{q}\right]\right|<\infty.$

###### Proof.

We have by property (6.5) from Proposition 6.1 that

$\displaystyle\mathbb{E}\left[X_{i;n}^{q}\right]=\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\\\
\pi\in
P(m)\end{subarray}}\prod_{b=1}^{B}\kappa\big{(}{\underbrace{X_{i;n},\ldots,X_{i;n}}_{\left|\pi_{b}\right|\text{
times}}}\big{)}.$

Since by assumption $X_{i;n}$ has mean $0$, we have

$\kappa(X_{i;n})=\mathbb{E}\left[X_{i;n}\right]=0.$

Thus, the only partitions $\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\in S(m)$
that give rise to non-zero terms in the expression above must have
$B\leq\lfloor\frac{q}{2}\rfloor$. Moreover, for any such partition, we have

$\left\lceil\frac{q}{2}\right\rceil=q-\left\lfloor\frac{q}{2}\right\rfloor=-\left\lfloor\frac{q}{2}\right\rfloor+\sum_{b=1}^{B}\left|\pi_{b}\right|\leq\sum_{b=1}^{B}\left(\left|\pi_{b}\right|-1\right).$

Hence, we find

$\displaystyle\sup_{n\geq
1}\left|n^{\lceil\frac{q}{2}\rceil}\mathbb{E}\left[X_{i;n}^{q}\right]\right|$
$\displaystyle\leq\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\\\
\pi\in P(m),\,\left|\pi_{b}\right|\geq
2\end{subarray}}\prod_{b=1}^{B}\sup_{n\geq
1}\left|n^{\left|\pi_{b}\right|-1}\kappa\big{(}{\underbrace{X_{i;n},\ldots,X_{i;n}}_{\left|\pi_{b}\right|\text{
times}}}\big{)}\right|<\infty,$

where the final inequality follows from the assumption (7.8). ∎

Applying Markov’s inequality and Lemma 7.5 shows that for any $q\geq 1$ we
have

$\sup_{n\geq 1}n^{q}\mathbb{P}\left(S_{n}^{c}\right)<\infty,\qquad
S_{n}:=\left\\{\left|X_{i;n}\right|\leq n^{-1/4},\qquad
i=1,\ldots,m\right\\}.$ (7.10)

This localization estimate allows us to replace each $g_{i}$ by its Taylor
expansion around $0$. Indeed, note that

$\kappa\left(g_{1}(X_{n}),\ldots,g_{k}(X_{n})\right)=P\left(\mathbb{E}\left[g_{1}(X_{n})^{q_{1}}\cdots
g_{k}(X_{n})^{q_{k}}\right],\quad q_{1}+\cdots+q_{k}\leq k\right)$

for some universal polynomial $P$ evaluated at the mixed moments of $X_{n}$
(the formula for this polynomial is given in (6.6) but is not important).
Moreover, using the growth assumption (7.7) on $X$ and the fact that $g_{i}$
are polynomially bounded we find that

$\sup_{n\geq 1}\mathbb{E}\left[g(X_{n})^{q_{1}}\cdots
g_{k}(X_{n})^{q_{k}}\right]<\infty,\qquad\forall\,q_{1},\ldots,q_{k}\geq 1.$
(7.11)

This, in combination with the localization estimate (7.10) applied with
$q=k-1$ yields

$\kappa\left(g_{1}(X_{n}),\ldots,g_{k}(X_{n})\right)=P\left(\mathbb{E}\left[{\bf
1}_{S_{n}}g_{1}(X_{n})^{q_{1}}\cdots g_{k}(X_{n})^{q_{m}}\right],\quad
q_{1}+\cdots+q_{m}\leq k\right)+O(n^{-k+1}).$

Note that for $n$ sufficiently large, on the event $S_{n}$, the argument
$X_{n}$ is any fixed neighborhood of $0$. Hence, we may write

$g_{j}(X_{n})=p_{j}(X_{n})+O(n^{-k+1}),$

where $p_{j}$ represents the $q-$th order Taylor expansion of $g_{j}$ around
$0$ with $q$ sufficiently large and the constant in the error term is
uniformly bounded. This yields

$\kappa\left(g_{1}(X_{n}),\ldots,g_{k}(X_{n})\right)=P\left(\mathbb{E}\left[{\bf
1}_{S_{n}}p_{1}(X_{n})^{q_{1}}\cdots p_{k}(X_{n})^{q_{m}}\right],\quad
q_{1}+\cdots+q_{m}\leq k\right)+O(n^{-k+1}).$

Finally, using the mixed moment estimates (7.7) and the localization estimate
(7.10), we conclude

$\displaystyle\kappa\left(g_{1}(X_{n}),\ldots,g_{k}(X_{n})\right)$
$\displaystyle=P\left(\mathbb{E}\left[p_{1}(X_{n})^{q_{1}}\cdots
p_{k}(X_{n})^{q_{m}}\right],\quad q_{1}+\cdots+q_{m}\leq k\right)+O(n^{-k+1})$
$\displaystyle=\kappa\left(p_{1}(X_{n}),\ldots,p_{k}(X_{n})\right)+O(n^{-k+1}).$

Recalling (7.8) completes the proof. ∎

Proposition 7.4 shows that, in establishing the conclusion (7.6) of Theorem
7.3, it is sufficient to assume that $g_{j}$ are polynomials. The remainder of
the proof of Theorem 7.3 is by induction on $\ell$, starting with $\ell=1$. In
view of Proposition 7.4, the following result establishes the base case.

######  Proposition 7.6 (Base Case: Theorem 7.3 holds for polynomials at
$\ell=1$).

Fix $k,m\geq 1$ and suppose $f_{i},\,i=1,\ldots,m$ are as in the statement of
Theorem 7.3. Then, if $p_{1},\ldots,p_{k}$ are any polynomials in $m$
variables, we have

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)\right|<\infty.$

###### Proof.

Since cumulants are multi-linear, we may and shall assume that $p_{a}$ are
monomials:

$p_{a}(x)=x^{Q^{(a)}}:=x_{1}^{q_{1}^{(a)}}\cdots x_{m}^{q_{m}^{(a)}},\qquad
x=\left(x_{1},\ldots,x_{m}\right),\quad
Q^{(a)}=\left(q_{1}^{(a)},\ldots,q_{m}^{(a)}\right).$ (7.12)

Recall that

$\mathcal{O}_{f_{i}}^{(1)}:=n_{1}^{-1}\sum_{j=1}^{n_{1}}f_{i}\left(z_{j;\mathcal{A}}^{(1)}\right).$

Therefore, writing $q^{(a)}:=q_{1}^{(a)}+\cdots+q_{m}^{(a)}$ we find

$\displaystyle
p_{a}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)=n_{1}^{-q^{(a)}}\sum_{J^{(a)}}f_{J^{(a)}},\qquad
f_{J^{(a)}}:=\prod_{i=1}^{m}\prod_{q=1}^{q_{i}^{(a)}}f_{i}\big{(}z_{j_{q;i}^{(a)};\mathcal{A}}^{(1)}\big{)},$

where the sum is over tuples of multi-indices

$J^{(a)}=\left(J_{1}^{(a)},\ldots,J_{m}^{(a)}\right),\qquad
J_{i}^{(a)}=\left(j_{q;i}^{(a)}\in[n_{1}],\,i\in[m],\,q\in[q_{i}^{(a)}]\right).$
(7.13)

Hence, using that cumulants are multi-linear (see (6.3)), we obtain

$\displaystyle\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)=n_{1}^{-(q^{(1)}+\cdots+q^{(k)})}\sum_{\begin{subarray}{c}J^{(1)},\ldots,J^{(k)}\end{subarray}}\kappa\left(f_{J^{(1)}},\ldots,f_{J^{(k)}}\right),$

where the sum extends over ordered collections $\left(J^{(a)},\,1\leq a\leq
k\right)$ of multi-indices as in (7.13). The expression on the right hand side
can be interpreted as an average. Namely, we can think of the indices
$j_{q;i}^{(a)}\in[n_{1}]$ are chosen uniformly from $[n_{1}]$ and
independently for all $i,q,a$. Writing $\mathcal{E}$ for the average with
respect to this distribution, we obtain

$\displaystyle\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)=\mathcal{E}\left[\kappa\left(f_{J^{(1)}},\ldots,f_{J^{(k)}}\right)\right].$

Our goal is to show that this average is small. To quantify this, let us
associate to each collection $\left(J^{(a)},\,a\in[k]\right)$ a graph

$\mathcal{G}\left(J^{(a)},\,a\in[k]\right)=\left([k],\,E\left(J^{(a)},\,a\in[k]\right)\right),$
(7.14)

with vertex set $[k]$ and edge set defined by

$(a,a^{\prime})\in\mathcal{E}\left(J^{(a)},\,a\in[k]\right)\quad\Longleftrightarrow\quad\exists
i,i^{\prime}\in[m],\,q\in[q_{i}^{(a)}],\,q^{\prime}\in[q_{i^{\prime}}^{(a^{\prime})}]\text{
s.t. }j_{q;i}^{(a)}=j_{q^{\prime};i^{\prime}}^{(a^{\prime})}.$

The key point is that in light of the vanishing property (6.2) of cumulants
and the fact that neurons at layer $1$ are independent

$\mathcal{G}\left(J^{(a)},\,a\in[k]\right)\text{ disconnected
}\quad\Longrightarrow\quad\kappa\left(f_{J^{(1)}}^{(1)},\ldots,f_{J^{(k)}}^{(1)}\right)=0.$

Hence,

$\displaystyle\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)=\mathcal{E}\left[{\bf
1}_{\left\\{\mathcal{G}\left(J^{(a)},\,a\in[k]\right)\text{
connected}\right\\}}\kappa\left(f_{J^{(1)}}^{(1)},\ldots,f_{J^{(k)}}^{(1)}\right)\right].$

Since $f_{i}$ are assumed to be polynomially bounded and the distribution of
the neuron pre-activations $z_{i;\alpha}^{{}_{(1)}}$ is that of centered
Gaussians with mean $0$ and covariance

$\mathrm{Cov}\left(z_{i_{1};\alpha_{1}}^{{}_{(1)}},z_{i_{2};\alpha_{2}}^{{}_{(1)}}\right)=\delta_{i_{1}i_{2}}\left(C_{b}+\frac{C_{W}}{n_{0}}\sum_{j=1^{n_{0}}}x_{j;\alpha_{1}}x_{j;\alpha_{2}}\right),$

we have for any fixed $k$ that

$\sup_{n\geq
1}\sup_{J^{(1)},\ldots,J^{(a)}}\left|\kappa\left(f_{J^{(1)}}^{(1)},\ldots,f_{J^{(k)}}^{(1)}\right)\right|<\infty.$

Hence,

$\displaystyle\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)=O\left(\mathcal{P}\left(\mathcal{G}\left(J^{(a)},\,a\in[k]\right)\text{
connected}\right)\right),$

where $\mathcal{P}$ is the probability measure associated to our random draw
of $J^{(1)},\ldots,J^{(k)}$. To complete the proof, note that since
$m,q_{i}^{(a)}$ are fixed, by a simple union bound, we obtain

$\displaystyle\mathcal{P}\left(\mathcal{G}\left(J^{(a)},\,a\in[k^{\prime}]\right)\text{
connected}~{}\bigg{|}~{}\mathcal{G}\left(J^{(a)},\,a\in[k^{\prime}-1]\right)\text{
connected}\right)=O(n^{-1}).$

Hence,

$\displaystyle\mathcal{P}\left(\mathcal{G}\left(J^{(a)},\,a\in[k]\right)\text{
connected}\right)$
$\displaystyle\qquad=\prod_{k^{\prime}=2}^{k}\mathcal{P}\left(\mathcal{G}\left(J^{(a)},\,a\in[k^{\prime}]\right)\text{
connected}~{}\bigg{|}~{}\mathcal{G}\left(J^{(a)},\,a\in[k^{\prime}-1]\right)\text{
connected}\right)$ $\displaystyle\qquad=O(n^{-k+1}).$ (7.15)

Thus,

$\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(1)}\right)\right)=O(n^{-k+1}),$

as desired. ∎

Propositions 7.4 and 7.6 together show that the conclusion (7.6) of Theorem
7.3 holds at layer $1$. In conjunction with Proposition 7.4, the following
result establishes that if the conclusion (7.6) of Theorem 7.3 holds at some
layer $\ell\geq 1$ then it also holds at layer $\ell+1$. This will complete
the proof by inductive of Theorem 7.3.

######  Proposition 7.7 (Inductive Step: Reducing polynomial cumulants in
layer $\ell+1$ to smooth cumulants in layer $\ell$).

Fix $\ell\geq 1$.

Case 1: Suppose that $\sigma$ is smooth. Assume that for any collection

$F^{\prime}=\left(f_{i}^{\prime}:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R},\,i=1,\ldots,m\right)$

of smooth and polynomially bounded functions and any $g_{j}$ as in the
statement of Theorem 7.3 the conclusion (7.6) of Theorem 7.3 holds at layer
$\ell$:

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(g_{1}\left(\overrightarrow{\mathcal{O}}_{F^{\prime}}^{(\ell)}\right),\ldots,g_{k}\left(\overrightarrow{\mathcal{O}}_{F^{\prime}}^{(\ell)}\right)\right)\right|<\infty.$

Then, if $p_{1},\ldots,p_{k}$ are any polynomials in $m$ variables, and
$F=\left(f_{i},\,i=1,\ldots,m\right)$ is an arbitrary collection of smooth and
polynomially bounded functions
$f_{i}:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R}$, then

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right)\right)\right|<\infty.$

Case 2: Suppose $\sigma$ is not smooth but satisfies Assumption 2.1 and that
$(D_{\alpha}^{J}z_{i;\alpha}^{(\ell)},\,\alpha\in\mathcal{A},\,\left|J\right|\leq
r)$ is non-degenerate in the infinite width in the sense of (3.3). Assume that
for any collection

$F^{\prime}=\left(f_{i}^{\prime}:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},r)}\rightarrow\mathbb{R},\,i=1,\ldots,m\right)$

of measurable and polynomially bounded functions and any $g_{j}$ as in the
statement of Theorem 7.3 the conclusion (7.6) of Theorem 7.3 holds at layer
$\ell$:

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(g_{1}\left(\overrightarrow{\mathcal{O}}_{F^{\prime}}^{(\ell)}\right),\ldots,g_{k}\left(\overrightarrow{\mathcal{O}}_{F^{\prime}}^{(\ell)}\right)\right)\right|<\infty.$

Then, if $p_{1},\ldots,p_{k}$ are any polynomials in $m$ variables, and
$F=\left(f_{i},\,i=1,\ldots,m\right)$ is an arbitrary collection of measurable
and polynomially bounded functions
$f_{i}:\mathbb{R}^{\left|\mathcal{A}\right|\times
N(n_{0},d)}\rightarrow\mathbb{R}$, then

$\sup_{n\geq
1}\left|n^{k-1}\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right)\right)\right|<\infty.$

###### Proof.

The proof of Proposition 7.7 is similar but somewhat more involved than that
of Proposition 7.6. Moreover, the two cases are proved in essentially the same
way, except that we will employ the different cases in Lemma 6.3. We give the
details in the case when $\sigma$ is smooth and indicate where the proof is
modified slightly to handle the non-smooth case.

To start, as in the proof of Proposition 7.6, note that since cumulants are
multi-linear (see (6.3)), it is enough to assume that $p_{j}$ are monomials.
Thus, borrowing the notation from the proof of Proposition 7.6 (see starting
(7.12)), we find

$\displaystyle\kappa\left(p_{1}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right),\ldots,p_{k}\left(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}\right)\right)=n_{\ell+1}^{-(q^{(1)}+\cdots+q^{(a)})}\sum_{\begin{subarray}{c}J^{(1)},\ldots,J^{(k)}\end{subarray}}\kappa\left(f_{J^{(1)}}^{(\ell+1)},\ldots,f_{J^{(k)}}^{(\ell+1)}\right),$

where

$f_{J^{(a)}}^{(\ell+1)}:=\prod_{i=1}^{m}\prod_{q=1}^{q_{i}^{(a)}}f_{j_{\alpha}^{(a)}}^{(\ell+1)},\qquad
f_{j}^{(\ell+1)}:=f\left(z_{j;\mathcal{A}}^{(\ell+1)}\right).$

Note that, as in Proposition A.1, the polynomially bounded assumption on
$f_{j}$ and the non-linearity $\sigma$ together with the Gaussianity of
weights and biases show that

$\sup_{n\geq
1}\left|\kappa\left(f_{J^{(1)}}^{(\ell+1)},\ldots,f_{J^{(k)}}^{(\ell+1)}\right)\right|<\infty.$
(7.16)

Using the law of total cumulance (6.1), we find that
$\kappa\left(p_{1}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}),\ldots,p_{k}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)})\right)$
equals

$\displaystyle\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\end{subarray}}n_{\ell+1}^{-(q^{(1)}+\cdots+q^{(a)})}\sum_{\begin{subarray}{c}J^{(1)},\ldots,J^{(k)}\end{subarray}}\kappa\left(\kappa\left(f_{J^{(\pi_{1})}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right),\ldots,\kappa\left(f_{J^{(\pi_{B})}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right)\right),$

where $\pi$ is any partition of $[k]$ and

$f_{J^{(\pi_{b})}}:=\left(f_{J^{(a)}},\,a\in\pi_{b}\right).$

Just as in the proof of Proposition 7.6, we may interpret the sum over
$J^{(1)},\ldots,J^{(k)}$ as an average over the distribution in which
$j_{q;i}^{(a)}$ are drawn iid uniformly on $[n_{\ell+1}]$. Writing
$\mathcal{E}$ for averages with respect to this distribution yields

$\kappa\left(p_{1}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}),\ldots,p_{k}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)})\right)=\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{b}\right)\end{subarray}}\mathcal{E}\left[\kappa\left(\kappa\left(f_{J^{(\pi_{1})}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right),\ldots,\kappa\left(f_{J^{(\pi_{b})}}^{(\ell+1)}~{}|~{}\mathcal{F}^{(\ell)}\right)\right)\right]$

As in (7.14), we may associate to each collection $J^{(\pi_{t})}$ the graph
$\mathcal{G}\left(J^{(\pi_{t})}\right)$. Recall that by Lemma 7.1, the neurons
pre-activations $D_{\alpha}^{J}z_{i;\alpha}^{(\ell+1)}$ in layer $\ell+1$ are
independent for different $i$ conditional on $\mathcal{F}^{(\ell)}$. Hence, in
view of the vanishing property (6.2) of cumulants, we obtain that
$\kappa\left(p_{1}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)}),\ldots,p_{k}(\overrightarrow{\mathcal{O}}_{F}^{(\ell+1)})\right)$
equals

$\displaystyle\sum_{\begin{subarray}{c}\pi=\left(\pi_{1},\ldots,\pi_{B}\right)\end{subarray}}\mathcal{E}\left[{\bf
1}_{\left\\{\mathcal{G}\left(J^{(\pi_{b})}\right)\text{ connected }\forall