Year: 2020

Abstract:
Understanding the asymptotic behavior of wide networks is of considerable interest. In this work, we present a general method for analyzing this large width behavior. The method is an adaptation of Feynman diagrams, a standard tool for computing multivariate Gaussian integrals. We apply our method to study training dynamics, improving existing bounds and deriving new results on wide network evolution during stochastic gradient descent. Going beyond the strict large width limit, we present closed-form expressions for higher-order terms governing wide network training, and test these predictions empirically.

heading: 1 INTRODUCTION
text: Neural networks achieve remarkable performance on a wide array of machine learning tasks, yet a complete analytic understanding of deep networks remains elusive. One promising approach is to consider the large width limit, in which the number of neurons in one or several layers is taken to be large. In this limit one can use a mean-field approach to better understand the network’s properties at initialization (Neal (1996); Lee et al. (2018)), as well as its training dynamics (Daniely (2017); Jacot et al. (2018)). Additional related works are cited below.
Suppose that f(x) is the network function evaluated at an input x. Let us denote the vector of model parameters by θ, whose elements are initially chosen to be i.i.d. Gaussian. In this work we consider a class of functions we call correlation functions, obtained by taking the ensemble averages of f , its products, and its derivatives with respect to the parameters θ, evaluated on arbitrary inputs. Here are a few examples of correlation functions.
Eθ [f(x1)f(x2)] , ∑ µ Eθ [ ∂f(x1) ∂θµ ∂f(x2) ∂θµ ] , ∑ µ,ν Eθ [ ∂f(x1) ∂θµ ∂f(x2) ∂θν ∂2f(x3) ∂θµ∂θν f(x4) ] . (1)
Correlation functions often show up in the study of wide networks. For example, the first correlation function in (1) plays a central role in the Gaussian Process picture of wide networks (Lee et al. (2018)), and has been used to diagnose signal propagation in wide networks (Pennington et al. (2017)). The second example in (1) is the ensemble average of the Neural Tangent Kernel (NTK), which controls the evolution of wide networks under gradient flow (Jacot et al. (2018)), and the third example shows up when computing the time derivative of the NTK with MSE loss.
While correlation functions can be computed analytically in some special cases (Cho & Saul (2009)), they are not analytically tractable in general. In this work, we present a method for bounding the asymptotic behavior of such functions at large width. In many cases of interest, the bound we obtain is tight. Derivation of the method relies on Feynman diagrams (Feynman (1949)), a technique for calculating multivariate Gaussian integrals, and specifically on the ’t Hooft expansion (’t Hooft (1974)). However, applying the method is straightforward and does not require any knowledge of Feynman diagrams.
Motivation. The method presented here is a general-purpose tool that can be added to a researcher’s toolbox. It allows a researcher to quickly identify ways in which network behavior and training
∗Authors listed alphabetically.
dynamics simplify at large width, for example by deriving which quantities vanish in this limit. As we demonstrate, using this method one is able to dramatically cut down on the amount of effort required to derive several key results related to wide networks, as well as to go beyond the infinite width limit and investigate ordinary neural networks.

heading: Our contribution.
text: 1. We present a general method for bounding the asymptotic behavior of correlation functions. The method is an adaptation of Feynman diagrams to the case of wide neural networks. The adaptation involves a novel treatment of derivatives of the network function, an element that is not present in the original theoretical physics formulation.
2. We apply the method to the study of wide network evolution under gradient descent. We improve on existing results for gradient flow (Jacot et al. (2018)) by deriving tighter bounds, and extending the analysis to the case of stochastic gradient descent (SGD). Going beyond the infinite-width limit, we present a formalisn for deriving finite-width corrections to network evolution, and present explicit formulas for the first order correction. To our knowledge, this is the first time this correction has been calculated.
3. As additional applications of our method, in Appendix E.2 we show that in the large width limit the SGD updates are linear in the learning rate, and in Appendix E.3 we discuss finite width corrections to the spectrum of the Hessian.
Limitations of our approach. The main result of this paper is a conjecture. We test our predictions extensively using numerical experiments, and prove the conjecture in some cases, but we do not have a proof that applies to all the cases we tested, including for deep networks with general non-linearities. Furthermore, our method can only be used to derive asymptotic bounds at large width; it does not produce the width-independent coefficient, which is often of interest.
Related work. For additional works on wide networks, including relating them to Gaussian processes, see de G. Matthews et al. (2018); Novak et al. (2019); Garriga-Alonso et al. (2019); Yang (2019); Pennington & Bahri (2017); Pennington & Worah (2018); Schoenholz et al. (2017); Xiao et al. (2018); Chen et al. (2018); Daniely et al. (2016); Lee et al. (2019); Du et al. (2018b;a); Allen-Zhu et al. (2018). For additional works discussing the training dynamics of wide networks see Geiger et al. (2019a); Arora et al. (2019). For a previous use of diagrams in this context, see Pennington & Worah (2017). The Neural Tangent Hierarchy presented in Huang & Yau (2019), published during the completion of this version, has significant overlap with the recursive differential equations (10) presented below.
The rest of the paper is organized as follows. In Section 2 we present our main conjecture and supporting evidence. In Section 3 we apply the method to gradient descent evolution of wide networks, and in Section 4 we present details on Feynman diagrams, which is the basic technique used in our proofs. We conclude with a Discussion. Proofs, additional applications, and details can be found in the Appendices.
Note: An earlier version of this work appeared in the ICML 2019 workshop, Theoretical Physics for Deep Learning (Dyer & Gur-Ari (2019)).

heading: 2 CORRELATION FUNCTION ASYMPTOTICS
text: In this section we present our main result: a method for computing asymptotic bounds on correlation functions of wide networks. We present the result as a conjecture, supported by analytic and empirical evidence.

heading: 2.1 NOTATION
text: Let f(x) ∈ R be the network output of a deep network with d hidden layers and input x ∈ RDin , defined by
f(x) = n−1/2V Tσ(n−1/2W d−1 · · ·σ(n−1/2W 1σ(Ux))) . (2)
Here U ∈ Rn×Din , V ∈ Rn, W 1, . . . ,W d−1 are weight matrices of dimension n, and σ : R→ R is the non-linearity.1 We denote the vector of all model parameters by θ. At initialization, the elements of θ are independent Gaussian variables, with each element θµ ∼ N (0, 1). The corresponding distribution of θ is denoted by P0.2
Let us now define correlation functions, the class of functions that is the focus of this work. These functions involve derivative tensors of the network function. We denote the rank-k derivative tensor by Tµ1...µk(x; f) := ∂
kf(x)/∂θµ1 · · · ∂θµk . For k = 0 we define T (x; f) := f(x), and still refer to this as a derivative tensor for consistency. Definition 1. A correlation function is the expectation value of a product of derivative tensors, evaluated at arbitrary inputs, where the tensor indices are summed in pairs over all the model parameters. A general correlation function C takes the form
C(x1, . . . , xm) := ∑
µ1,...,µkm
∆(π)µ1...µkmEθ [ Tµ1...µk1 (x1)Tµk1+1...µk2 (x2) · · ·Tµkm−1+1...µkm (xm) ] .
(3)
Here, 0 ≤ k1 ≤ · · · ≤ km−1 ≤ km are integers,3 m and km are even, π ∈ Skm is a permutation, and ∆(π)µ1...µkm = δµπ(1)µπ(2) · · · δµπ(km−1)µπ(km) . We use δ to denote the Kronecker delta.
If two derivative tensors in a correlation function have matching indices that are summed over, we say that they are contracted. For example, the correlation function ∑ µ Eθ [∂f(x1)/∂θµ · ∂f(x2)/∂θµ] has one pair of contracted tensors. See (1) for additional examples of correlation functions.

heading: 2.2 ASYMPTOTIC BOUNDS ON WIDE NETWORKS
text: We now present our main conjecture, which allows us to place asymptotic bounds on general correlation functions of wide networks. Conjecture 1. Let C(x1, . . . , xm) be a correlation function. The cluster graph GC(V,E) of C is a graph with vertices V = {v1, . . . , vm}, where vi = T (xi) for all i, and edges E = {(vi, vj) | (T (xi), T (xj)) contracted in C}. Suppose that the cluster graph GC has ne connected components with an even size (even number of vertices), and no components of odd size. Then C(x1, . . . , xm) = O(nsC ), where
sC = ne + no 2 − m 2 . (4)
We will refer to the connected components of a cluster graph GC as the clusters of C. Table 1 lists examples of bounds derived using the Conjecture for several correlation functions. The intuition behind Conjecture 1 comes from the following result for deep linear networks. Theorem 1. Conjecture 1 holds for correlation functions of networks with linear activations, including for deep linear networks with biases (affine transformations).
Let us discuss the intuition behind this theorem. Computing correlation functions of deep linear networks amounts to evaluating Gaussian integrals with polynomial integrands in θ. One can evaluate such integrals using Isserlis’ theorem, which tells us how to express moments of multivariate Gaussian variables in terms of their second moments. For example, given centered Gaussian variables z1, ..., z4,
Ez [z1z2z3z4] = Ez [z1z2]Ez [z3z4] + Ez [z1z3]Ez [z2z4] + Ez [z1z4]Ez [z2z3] . (5) Therefore, correlation functions of deep linear networks can be expressed in terms of the covariances Eθ [UiαUjβ ] = δijδαβ , Eθ [ViVj ] = δij , and Eθ [ W (l) ij W (l) kl ] = δikδjl. For example, for a deep linear network with 2 hidden layers, we have
Eθ [f(x1)f(x2)] = 1 n2 Eθ [ V TWUx1V TWUx2 ] = xT1 x2 n2 n∑ i,k δikδik n∑ j,l δjlδjl = x T 1 x2 . (6)
1We take all layers widths to be equal to n for simplicity, but our results hold in the more general case where all widths scale linearly with n.
2We use µ, ν, . . . to denote θ indices, i, j, . . . to denote individual weight matrix and weight vector indices, and α, β, . . . for input dimension indices.
3When ka = ka−1, the tensor T (xa) has no derivatives.
Every correlation function of a deep linear network can be similarly reduced to sums over products of Kronecker delta functions and width-independent functions of the inputs. The asymptotic large width behavior is determined by these sums over delta functions, which are tedious to compute by hand. Feynman diagrams are a graphical tool for computing these sums, allowing us to obtain the asymptotic behavior with minimal effort. This tool, which is described in detail in Section 4, is used to prove Theorem 1.
For networks with non-linear activations we further show the following Theorem 2. Conjecture 1 holds for (1) networks with ReLU activations, where all inputs are set to be equal, and for (2) networks with one hidden layer and smooth activation.
For case (1), the idea behind the proof is to put an asymptotic bound on the ReLU network in terms of a corresponding deep linear network. For case (2), the basic idea is that each network function contains a single sum over the width, and by keeping track of these sums using Feynman diagrams we are able to bound the asymptotic behavior. We refer the reader to Appendix C for details.

heading: 2.3 NUMERICAL EXPERIMENTS
text: Table 1 lists asymptotic bounds on several correlation functions, derived using Conjecture 1. These are compared against the asymptotic behavior computed using numerical experiments. In addition to the results presented here, we performed experiments using the same correlation functions and experimental setup, but with weights sampled uniformly from {±1} instead of from a Gaussian distribution. The results are shown in Appendix A.1.
In all cases that were tested empirically, we found that Conjecture 1 holds. For networks with smooth activations, we found that Conjecture 1 always gives a tight bound. For networks with linear or ReLU activations, we always find that the Conjecture holds as an upper bound, but that sometimes the bound is not tight. One such case is highlighted in Table 1. In such cases, a tight bound can be obtained for deep linear networks using the complete Feynman diagram analysis presented below.4

heading: 3 APPLICATIONS TO TRAINING DYNAMICS
text: In this section we apply Conjecture 1 to study the evolution of wide networks under gradient flow and gradient descent. We begin by briefly reviewing existing results. Let Dtr be a training set of size M ,
4For deep linear networks, the proof of Theorem 1 relies on mapping diagrams to triangulations of 2D surfaces with at least 1 boundary. There are cases, such as the one highlighted in Table 1, in which all diagrams one can write down have at least 2 boundaries. In such cases, the full diagrammatic calculation gives a tight bound, while Conjecture 1 provides a valid upper bound that is not tight.
and let L = ∑
(x,y)∈Dtr `(x, y) be the MSE loss, with single sample loss `(x, y) = 1 2 (f(x) − y)2.
The gradient flow equation is dθdt = −∇θL. The evolution of the network function under gradient flow is given by
df(x)
dt = − ∑ (x′,y′)∈Dtr Θ(x, x′) ∂`(x′, y′) ∂f . (7)
Here, Θ is the Neural Tangent Kernel (NTK), defined by Θ(x1, x2) := ∇θfT (x1)∇θf(x2). The authors of Jacot et al. (2018) showed that the kernel is constant during training up to O(n−1/2) corrections. This leads to a dramatic simplification in training dynamics (Jacot et al. (2018); Lee et al. (2019)). In particular, for MSE loss the network map evaluated on the training data evolves as f(t) = y + e−tΘ (0)
(f (0) − y).5 We will use our technology to derive a tighter bound on finite-width corrections to the kernel during training and present explicit formulas for the leading correction.
The following result is useful in analyzing the behavior of correlation functions under gradient flow.
Lemma 1. Let C(~x) = Eθ [F (~x)] be a correlation function, where ~x = (x1, . . . , xm), and suppose that C = O(nsC ) for sC as defined in Conjecture 1. Then Eθ [ dkF (~x) dtk ] = O(nsC ) for all k.
Here we prove the statement for k = 1. Appendix D.2 contains a proof for the general case.
proof (k = 1). Let C have ne even clusters and no odd clusters. Consider Eθ [ dF (~x) dt ] =
−∑µ∑x′∈Dtr Eθ [∂F (~x)∂θµ ∂f(x′)∂θµ f(x′)]. Denote by n′e (n′o) the number of even (odd) clusters in this correlation function, which has m′ = m + 2 derivative tensors. One can check that either (n′e, n ′ o) = (ne + 1, no) or (n ′ e, n ′ o) = (ne − 1, no + 2), depending on whether the ∂µF derivative is acting on an odd or even cluster in F .6 Therefore, n′e + n′o 2 − m ′ 2 ≤ sC .
With this result, it is easy to understand the constancy of the NTK at large width. The first derivative of the NTK is given by
Eθ [ dΘ(x1, x2)
dt
] = − M∑ x′∈Dtr ∑ µ,ν Eθ [ ∂2f(x1) ∂θµ∂θν ∂f(x2) ∂θµ ∂f(x′) ∂θν f(x′) ] + (x1 ↔ x2) . (8)
Here, (x1 ↔ x2) means “add the same term, but exchange x1 and x2”. This correlation function has ne = 0, no = 2, and m = 4, and is thereforeO(n−1) by Conjecture 1. By Lemma 1, all higher-order time derivatives of the NTK are O(n−1) as well. If we now assume that the time-evolved kernel Θ(t) is analytic in training time t, and that we are free to exchange the Taylor expansion in time with the large width limit, then we find that Eθ [Θ(t)−Θ(0)] = ∑∞ k=1 tk k!Eθ [ dkΘ(0) dtk ] = O(n−1) for any fixed t. This bound, which is tighter than that found in Jacot et al. (2018), was noticed empirically in Lee et al. (2019) as well as in our own experiments, see Figure 1a.
This analysis can be extended to show that Eθ [Θ(t)−Θ(0)] = O(n−1) for SGD as well. The technique is similar, and again relies on Conjecture 1. We refer the reader to Appendix E.1 for details. These results improve on existing ones in several ways. Our method applies to the case of SGD as well as to networks where all layer widths are increased simultaneously — a setup that has proven to be difficult to analyze. In addition, the O(n−1) bound we derive on kernel corrections during SGD is empirically tight, and improves on the existing bound of O(n−1/2) which was derived for gradient flow (Jacot et al. (2018)).
5Here we are using condensed notation: Θ(0) and f (0) are values at initialization, and f , f (0), y are treated as vectors in training set space. The kernel Θ(0) is a square matrix in the same space.
6Here we are extending the use of the term cluster to refer to derivative tensors in the integrand itself.

heading: 3.1 FINITE WIDTH CORRECTIONS
text: Next, we will compute the explicit time dependence of Θ and f at order O(n−1) under gradient flow. This is the leading correction to the infinite width result. We define the functions O1(x) := f(x) and
Os(x1, . . . , xs) := ∑ µ ∂Os−1(x1, . . . , xs−1) ∂θµ ∂f(xs) ∂θµ , s ≥ 2 . (9)
Notice that O2 = Θ is the kernel. It is easy to check that
dOs(x1, . . . , xs)
dt = − ∑ (x′,y′)∈Dtr Os+1(x1, . . . , xs, x ′)(f(x′)− y′) , s ≥ 1 . (10)
We can use equations (9) and (10) to solve for the evolution of the kernel and network map. Notice that each function Os has s derivative tensors and a single cluster. As a result, correlation functions involving operators with larger s are increasingly suppressed in width. In particular, Eθ [dO4/dt] = −∑x∈Dtr Eθ [O5(x)f(x)] = O(n−2) at all times, using Conjecture 1 and Lemma 1. Thus to solve for f and Θ at O(n−1) we can set Os = 0 for all s ≥ 5. Let us denote the time-evolved kernel by Θ(t) = Θ(0) +Θ1(t)+O(n−2), where Θ(0) is the kernel at initialization, and Θ1(t) is the O(n−1) correction we are seeking. Integrating equations (10) starting with s = 4, we find
Θ1(x1, x2; t) =− ∫ t
0 dt′ ∑ x∈Dtr O (0) 3 (x1, x2, x)∆f(x; t ′)
+ ∫ t 0 dt′ ∫ t′ 0 dt′′ ∑
x,x′∈Dtr O
(0) 4 (x1, x2, x, x ′)∆f(x′; t′′)∆f(x; t′) . (11)
Here we have introduced the notation ∆f(x; t) = e−tΘ0(f (0) − y). A detailed derivation can be found in Appendix E.4. There we also evaluate the integrals in (11) in terms of the NTK spectrum.
To obtain the O(n−1) correction to the network map (evaluated for simplicity on the training data), we further integrate (10) for s = 1 and find
f(t) = f0(t)− e−tΘ (0) ∫ t 0 dt′et ′Θ(0)Θ1(t ′)e−t ′Θ(0)(f (0) − y) +O(n−2) . (12)
Here we have denote the infinite width evolution by f0(t) = y + e−tΘ (0) (f (0) − y). Figures 1b and 1c compare these predictions against empirical results.

heading: 4 FEYNMAN DIAGRAMS FOR DEEP LINEAR NETWORKS
text: In this section we present the Feynman diagram technique, and show how it allows us to compute the asymptotic behavior of correlation functions. We end this section with a proof of Theorem 1 for the case of networks with a single hidden layer and linear activations.
Given a correlation function C, we map it to a family of graphs called Feynman diagrams. The graphs are independent of the inputs, and are defined as follows. Definition 2. Let C(x1, . . . , xm) be a correlation function for a network with d hidden layers. The family Γ(C) is the set of all graphs that have the following properties.
1. There are m vertices v1, . . . , vm, each of degree d+ 1.
2. Each edge has a type t ∈ {U,W 1, . . . ,W d−1, V }. Every vertex has one edge of each type. 3. If two derivative tensors Tµ1,...,µ`(xi), Tν1,...,ν`′ (xj) are contracted k times in C, the graph
must have at least k edges (of any type) connecting the vertices vi, vj .
The graphs in Γ(C) are called the Feynman diagrams of C.
Single hidden layer. For the rest of this section we focus on networks with a single hidden layer. We refer the reader to Appendix B for a full treatment of deep linear networks. For networks with one hidden layer and linear activation, the network output is f(x) = n−1/2V TUx. Consider the correlation function C(x1, x2) = Eθ [f(x1)f(x2)]. We have
C(x1, x2) = 1
n n∑ i,j Din∑ α,β EV [ViVj ]EU [Uiα Ujβ ]xα1x β 2 = xT1 x2 n n∑ i,j δijδij = x T 1 x2 . (13)
As the factors of x1 and x2 are independent of n, we see that C(x1, x2) = O(n0). Notice that there are two relevant contributions to this answer: each factor of the network function in the integrand contributes n−1/2, and the summed-over product of Kronecker deltas contributes n. Other details, such as the input dependence, are irrelevant. Feynman diagrams allow us to encode only those details that affect the n scaling, ignoring the rest.
The set Γ(C) for the correlation function (13) consists of a single Feynman diagram, shown in Figure 2a. The asymptotic bound on a correlation function is obtained by the following result, which is due to ’t Hooft (1974). Theorem 3. LetC(x1, . . . , xm) be a correlation function with one hidden layer and linear activation. Then C = O(ns) where s = maxγ∈Γ(C) lγ − m2 , and lγ is the number of loops in γ.7
Let us give some intuition for Theorem 3. Each Feynman diagram γ encodes a subset of the terms contributing to the correlation function. To get the asymptotic bound on the correlation function, we sum over the contributions of individual diagrams. We can compute the asymptotic behavior of a single diagram γ using the following Feynman rules: (1) each vertex contributes a factor of n−1/2, and (2) each loop contributes a factor of n. Therefore, if a diagram has lγ loops, its contribution to
7For networks with a single hidden layer, the number of loops in a diagram is equal to the number of connected components. This is not true for a general deep linear network.
the correlation function scales as nlγ−m/2. Rule (1) is due to the explicit n−1/2 factor in the network definition. Rule (2) follows from applying Isserlis’ theorem, as follows. Each covariance factor (such as the factor E [ViVj ] = δij in eq. (13)) corresponds to an edge in a Feynman diagram. A loop in the diagram corresponds to a sum over a product of Kronecker deltas, yielding a factor of n.
Returning to our example, the graph in Figure 2a has two vertices and one closed loop, so we recover the asymptotic behavior O(n0) from Theorem 3. As another example, consider the correlation function C(x1, x2, x3, x4) = Eθ [f(x1)f(x2)f(x3)f(x4)]. It has two Feynman diagrams, shown in Figures 2b and 2c. Using the Feynman rules, we find that the disconnected graph represents terms that scale as O(n0), while the connected graph represents terms that scale as O(n−1). Therefore, C(x1, x2, x3, x4) = O(n0). This is an example of a more general phenomenon, that connected graphs vanish faster at large n compared with disconnected graphs.
Correlation functions with derivatives. We now extend the Feynman diagram technique to correlation functions that include derivatives of f .8 As a concrete example, consider the correlation function C(x, x′) = Eθ [Θ(x, x′)], where Θ is the kernel defined in Section 3. For a single hidden layer, we have
C(x, x′) = n∑ i=1 Eθ [ ∂f(x) ∂Ui ∂f(x′) ∂Ui + ∂f(x) ∂Vi ∂f(x′) ∂Vi ] . (14)
The two derivative tensors in this correlation function are contracted: their indices are set to be equal and summed over. Therefore, according to Definition 2, Γ(C) includes all diagrams in which the corresponding vertices share at least one edge. The resulting diagrams are shown in Figure 3. The edges forced by the contraction are explicitly marked by dashed lines for clarity, but mathematically they are ordinary edges. We see that in fact there is only one diagram contributing to this correlation function — the same one shown in Figure 2a. Following the Feynman rules, we find that Eθ [Θ(x, x′)] = O(n0). The fact that contracted derivatives should be mapped to forced edges in the Feynman diagrams is proved in Appendix B. The basic reason behind this rule is the relation ∑ k ∂Vi ∂Vk ∂Vj ∂Vk
= δij = E [ViVj ]. Namely, when derivatives act in pairs they yield a Kronecker delta factor (δij), which is equal to the factor obtained from a covariance (E [ViVj ]). While Isserlis’ theorem instructs us to sum over all possible covariance configurations (and therefore over all possible edge configurations), a pair of summed derivative leads to a particular covariance factor. Therefore, we should only consider graphs that include the edge corresponding to this covariance factor.
We are now ready to prove Theorem 1 for the case of single hidden layer with linear activations. A proof for the general case can be found in Appendix B.
Proof (Theorem 1, one hidden layer). Let C(x1, . . . , xm) be a correlation function for a network with a single hidden layer and linear activation. Let GC(V,E) be the cluster graph of C, and let γ(V,E′) ∈ Γ(C) be a Feynman diagram. Notice that E ⊂ E′. Indeed, E′ contains an edge corresponding to each pair of contracted derivative tensors in C, and E = {(vi, vj) | (T (xi), T (xj)) contracted in C}. In addition, notice that γ only contains connected
8We are not aware of a physics application in which such derivatives are included in a Feynman diagram description of correlation functions. Therefore, to our knowledge our treatment of these derivatives is novel both in machine learning and in physics.
components (i.e. loops) with an even number of vertices, because every vertex has exactly one edge of each type. Therefore, nγ ≤ ne+ no2 , where nγ is the number of loops in γ, and ne (no) is the number of even (odd) components in GC . The bound is saturated when every even component of GC belongs to a different component of γ, and pairs of odd components in GC belong to different components of γ. From Theorem 3, we have that C = O(ns) where s = maxγ∈Γ(C) nγ − m2 ≤ ne + no2 − m2 .

heading: 5 DISCUSSION
text: Ensemble averages of the network function and its derivatives are an important class of functions that often show up in the study of wide neural networks. Examples include the ensemble average of the train and test losses, the covariance of the network function, and the Neural Tangent Kernel (Jacot et al. (2018)). In this work we presented Conjecture 1, which allows one to derive the asymptotic behavior of such functions at large width.
For the case of deep linear networks, we presented a complete analytic understanding of the Conjecture based on Feynman diagrams. In addition, we presented empirical and analytic evidence showing that the Conjecture also holds for deep networks with non-linear activations, as well as for networks with non-Gaussian initialization. We found that the Conjecture holds in all cases we tested.
The basic tools presented in this work can be applied to many aspects of wide network research, greatly simplifying theoretical calculations. We presented several applications of our method to the asymptotic behavior of wide networks during stochastic gradient descent, and additional applications are presented in Appendix E. We were able to improve upon known results by tightening existing bounds, and by applying the technique to SGD as well as to gradient flow. In addition, we took a step beyond the infinite width limit, deriving closed-form expressions for the first finite-width correction to the network evolution. These novel results open the door to studying finite-width networks by systematically expanding around the infinite width limit.
A central question in the study of wide networks is whether the infinite width limit is a good model for describing the behavior of realistic deep networks (Chizat et al. (2018); Ghorbani et al. (2019); Geiger et al. (2019b)). In this work we take a step toward answering this question, by working out the next order in a perturbative expansion around the infinite width limit, potentially bringing us closer to an analytic description of finite-width networks. We hope that the techniques presented here provide a basis to systematically answering these and other questions about the behavior of wide networks.

heading: ACKNOWLEDGEMENTS
text: The authors would like to thank Alex Alemi, Yasaman Bahri, Boris Hanin, Jared Kaplan, Jaehoon Lee, Behanm Neyshabur, Sam Schoenholz, Sylvia Smullin, and Jascha Sohl-Dickstein for useful discussion. The authors would especially like to thank Ying Xiao for extensive comments on early versions of this manuscript.

heading: A EXPERIMENTAL DETAILS AND ADDITIONAL RESULTS

heading: A.1 NON-GAUSSIAN INITIALIZATION
text: In Table 1 we listed asymptotic bounds on several correlation functions, where the model parameters were initialized from a Gaussian distribution. Table 2 shows additional results using the same experimental setup, but with weights sampled uniformly from {±1}. We again find good agreement with the predictions of Conjecture 1.

heading: A.2 EXPERIMENTAL DETAILS
text: The experiments in Figure 1 were performed on two-class MNIST, computing a single randomlychosen component of the kernel Θ. Sub-figure (a) uses networks trained for 1024 steps with learning rate 1.0 and 1000 samples per class, averaged over 100 initializations. Each curve in figure (b) represents a single instance of the network map evaluated on a random image over the corse of training. The models were trained with 10 samples per class and learning rate 0.1. The input to the network is normalized by the square root of the input dimension as in (Jacot et al. (2018))
f(x) = n−1/2V Tσ(n−1/2W d−1 · · ·σ(n−1/2W 1σ(D−1/2in Ux))) . (15)

heading: B FEYNMAN DIAGRAMS FOR DEEP LINEAR NETWORKS
text: Feynman diagrams can be used to derive asymptotic upper bounds on deep linear networks in the large width limit. In this section we describe the method in detail, and use it to prove Theorem 1.

heading: B.1 FEYNMAN DIAGRAMS AND DOUBLE-LINE DIAGRAMS
text: In this section we build on the results of Section 4 and consider correlation functions of deep linear networks with d hidden layers. The network function was defined in (2), and here we set the activation σ to be the identity. Definition 2 describes how to map a correlation function C to Γ(C), a family of graphs called Feynman diagrams. The Feynman diagram method relies on Isserlis’ theorem, which allows us to express arbitrary moments of multivariate Gaussian variables in terms of their covariance.
Theorem 4 (Isserlis). Let z = (z1, . . . , zl) be a centered multivariate Gaussian variable. For any positive integer k,
Ez [zi1 · · · zi2k ] = 1
2kk! ∑ π∈S2k E [ ziπ(1)ziπ(2) ] E [ ziπ(3)ziπ(4) ] · · ·E [ ziπ(2k−1)ziπ(2k) ] , (16)
Ez [ zi1 · · · zi2k−1 ] = 0 . (17)

heading: In particular, if the covariance matrix of z is the identity then
text: Ez [zi1 · · · zi2k ] = 1
2kk! ∑ π∈S2k δiπ(1)iπ(2)δiπ(3)iπ(4) · · · δiπ(2k−1)iπ(2k) . (18)
Using this theorem, a correlation function C can be expressed as a sum over permutations as in (16). Each term in this sum maps to a Feynman diagram in Γ(C).
As an example, consider the correlation function C(x, x′) = Eθ [f(x)f(x′)] for a network with 2 hidden layers. An explicit calculation gives
C(x, x′) = 1
n2 n∑ i,j,k,l Din∑ α,β Eθ [ ViW 1 ijUjαVkW 1 klUlβ ] xαx ′ β (19)
= 1
n2 n∑ i,j,k,l Din∑ α,β EV [ViVk]EW [ W 1ijW 1 kl ] EU [Ujα Ulβ ]xαx′β (20)
= xTx′
n2 n∑ i,k=1 δikδik n∑ j,l=1 δjlδjl = x Tx′ = O(n0) . (21)
To get from (19) to (20), we applied Isserlis’ theorem for every choice of indices i, j, k, l, α, β. We find that there is at most one permutation π of the network parameters such that the covariances do not vanish. This is because the covariance of parameters across different layers vanishes identically (for example Eθ [ ViW 1 jk ] = 0). Correspondingly, this correlation function has a single Feynman diagram, shown in Figure 4a.
For networks with one hidden layer, we saw in Section 4 that the asymptotic behavior is determined by the number of loops in a graph. This observation does not immediately generalize to networks with general depth, because it is not obvious how to count loops in diagrams such as the one in Figure 4a. The problem can be traced back to the fact that weight matrices have covariances of the form E [WijWkl] = δikδjl involving two Kronecker deltas, but the procedure described so far assumes that each covariance (and each edge in the graph) corresponds to a single Kronecker delta.
A similar question appeared in the context of theoretical physics, and the correct generalization is due to ’t Hooft (1974). The idea is to treat each Feynman diagram as the triangulation of a Riemann surface, and to define the number of loops in a graph to be the number of faces of the triangulation. In practice, this involves mapping each Feynman diagram to a new double-line diagram: A graph in which each edge corresponds to a single Kronecker delta factor, and loops correspond to triangulation faces of the original diagram. Definition 3. Let γ ∈ Γ(C) be a Feynman diagram for a correlation functionC involving k derivative tensors for a network of depth d. Its double-line graph, DL(γ) is a graph with kd vertices of degree 2, defined by the following blow-up procedure.
• Each vertex v(i) in γ is mapped to d vertices v(i)1 , . . . , v (i) d in DL(γ).
• Each edge (v(i), v(j)) in γ of type U is mapped to a single edge (v(i)1 , v (j) 1 ).
• Each edge (v(i), v(j)) in γ of type W l is mapped to two edges (v(i)l , v (j) l ), (v (i) l+1, v (j) l+1).
• Each edge (v(i), v(j)) in γ of type V is mapped to a single edge (v(i)d , v (j) d ).
The number of faces in γ is given by the number of loops in the double-line graph DL(γ).
Figure 4 shows the Feynman diagram and corresponding double-line diagram for Eθ [f(x)f(x′)] with 2 hidden layers. We can interpret this Feynman diagram as a triangulation of the disc: a 2-dimensional
surface with a single boundary. The triangulation has 2 vertices, 3 edges corresponding to the edges of the Feynman diagram, and 2 faces correponding to the loops of the double-line diagram. Figure 5 shows additional examples of double-line diagrams, and Figure 6 shows the double-line diagrams of a correlation function with derivatives. As explained in Section 4, contracted derivative tensors in a correlation function C map to forced edges in Γ(C), and these are marked with dashed lines on the diagrams.
Generally, every Feynman diagram of a deep linear network can be interpreted as the triangulation of a 2-dimensional manifold with at least one boundary. Intuitively, the presence of the boundary is due to the fact that the U, V weights at both ends of network have only one dimension that scales linearly with the width. As a result, in the double-line diagrams the U, V edges become single lines rather than double lines. These ‘missing’ lines translate into fewer faces in the triangulation, and the ‘missing’ faces can be geometrically interpreted as boundaries in the corresponding surface.9
The discussion above is summarized by the following result, due to ’t Hooft (1974), that describes how the asymptotic behavior of a general correlation function can be computed using the Feynman rules for deep linear networks.
Theorem 5. Let C(x1, . . . , xm) be a correlation function of a deep linear network with d hidden layers, and let γ ∈ Γ(C) be a Feynman diagram. The diagram represents a subset of terms that contribute to C, and its asymptotic behavior is determined by the Feynman rules: the subset isO(nsγ ) where sγ = lγ − dm2 , and lγ is the number of loops in the double-line diagram DL(γ). Furthremore, the correlation function is C = O(ns), where s = maxγ∈Γ(C) sγ .
The intuition for the formula sγ = lγ − dm2 is similar to the single hidden layer case, Theorem 3. The term lγ counts the number of factors of the form ∑ i1,...,ik
δi1i2δi2i3 · · · δiki1 = n that appear in the Correlation function after applying Isserlis’ theorem. The term ( −dm2 ) is due to the explicit n−d/2 normalization of the network function.
9One can consider the cyclic model f̃(x) = n−(d+1)/2 Tr ( Ṽ W d−1 · · ·W 1Ũ ) x with 1D input x, in which
all the weight tensors Ũ , Ṽ ,W 1, . . . ,W d−1 are n× n matrices. The Feynman diagram construction for this model is similar to the deep linear case, except that all edges in the Feynman diagram map to double edges in the double-line diagrams. Such diagrams can be interpreted as triangulations of surfaces with no boundaries. Unlike the deep network diagrams, they have no ‘missing’ loops.

heading: B.2 ASYMPTOTICS OF DEEP LINEAR NETWORKS
text: We now prove Theorem 1. The theorem follows from the following lemma, again due to ’t Hooft (1974), that relates the asymptotic behavior to the number of connected components in a Feynman diagram.
Lemma 2. Let C(x1, . . . , xm) be a correlation function for a deep linear network. Let cγ be the number of connected components of a graph γ ∈ Γ(C). Then C = O(ns), where
s = max γ∈Γ(C)
cγ − m
2 . (22)
Proof. We prove the result for 1D inputs (Din = 1), and it is easy to generalize to arbitrary input dimension. Let γ ∈ Γ(C) be a Feynman diagram and let γ′ be a connected component of γ with vγ′ vertices. Notice that we can apply the Feynman rules of Theorem 5 separately to each component γ′ and find a bound O(nsγ′ ). Then, γ = O(nsγ ) where sγ = ∑ γ′ sγ′ , and the sum runs over the connected components of sγ . We will show below that sγ′ ≤ 1− vγ′2 , and therefore sγ ≤ cγ − m2 , which is sufficient to prove (22).
Let us prove the remaining statement about sγ′ . The Euler character of the graph γ′ with v = vγ′ vertices, e edges and f faces is χ = v − e+ f . The degree of each vertex in the graph is d+ 1, and therefore e = (d+1)v2 . Using Theorem 5 the graph is O(nsγ′ ) where sγ′ = f − dv2 . We therefore find that χ = v2 + sγ′ . The graph γ
′ is a triangulation of some connected surface with at least one boundary. The Euler character for such a surface is bounded by χ ≤ 1, and therefore sγ′ ≤ 1− v2 .
Let us now prove Theorem 1.
Proof. Let C(x1, . . . , xm) be a correlation function for a deep linear network. Suppose that the cluster graph GC has ne even size components and no odd size components. Let γ ∈ Γ(C) be a Feynman diagram with cγ connected components. We will show that cγ ≤ ne + no2 . It then follows immediately from Lemma 2 that C = O(ns) where s = ne + no2 − m2 , concluding the proof. Let us derive the bound on cγ . First, all vertices that belong to a given cluster (a component of GC) will also belong to the same connected component in γ. This is because every edge in GC is also an edge in γ (note that GC and γ have the same set of vertices). Therefore, cγ ≤ ne + no. Second, note that every connected component of the graph γ has an even number of vertices. Indeed, each edge has a type t, and each vertex has exactly one edge of each type. Therefore, a connected component with v vertices has v2 edges of each type, and so v must be even. It follows that the vertices of even clusters can form their own connected components in a Feynman diagrams, while odd clusters must be connected in sets of 2 or more to form connected components. The bound on cγ then follows.
For deep linear networks with bias the proof is the same, except we use Lemma 3 below instead of Lemma 2.

heading: B.3 DEEP LINEAR NETWORKS WITH BIAS
text: In this section we argue that Conjecture 1 holds for deep linear networks with bias. Let us add bias to the definition (2) of a deep linear network (setting σ to be the identity). We choose x = 1 without loss of generality. The activations zl at layer l and the network function f are given by
z1 = U , (23)
zl+1 = n−1/2W lzl + bl , l = 1, . . . , d− 1 , (24) f = n−1/2V T zd . (25)
Here b1, . . . , bd−1 ∈ Rn are biases whose elements are chosen independently from N (0, 1) at initialization. The network function can be written as
f = d+1∑ k=2 n−(k−1)/2V TW d−1 · · ·W d+2−kbd+1−k , (26)
where we defined b0 = U . In the rest of this Section we extend our definitions of Feynman diagrams and double-line diagrams to deep linear networks with bias, and then extend several of our results to this case. Definition 4. Let C(x1, . . . , xm) be a correlation function for a network with d hidden layers and bias. The family Γ(C) is the set of all graphs that have the following properties.
1. There are m vertices v1, . . . , vm, where each vi has degree ki ∈ {2, . . . , d+ 1}. 2. Each edge has a type t. Every vertex of degree k has one edge of each type t ∈ {bd+1−k,W d+2−k, . . . ,W d−1, V }.
3. If two derivative tensors Tµ1,...,µ`(xi), Tν1,...,ν`′ (xj) are contracted s times in C, the graph must have at least s edges (of any type) connecting the vertices vi, vj .
The intuition behind this definition is that term k in the sum in (26) is equal to the network function of a deep linear network with depth k − 1. Correlation functions of such a function can be represented diagramatically by including vertices of degree k. To represent the entire sum we must consider all diagrams with vertices of varying degree.
We now extend Definition 3 of double-line graphs to accommodate vertices with varying degree. Definition 5. Let γ ∈ Γ(C) be a Feynman diagram for a correlation function C involving m derivative tensors for a network of depth d with biases. Its double-line graph, DL(γ) is a graph with vertices of degree 2, defined by the following blow-up procedure.
• Each vertex v(i) in γ of degree ki is mapped to ki − 1 vertices v(i)1 , . . . , v (i) ki−1 in DL(γ).
• Each edge (v(i), v(j)) in γ of type W l is mapped to two edges, (v(i)l+ki−d−1, v (j) l+kj−d−1) and
(v (i) l+ki−d, v (j) l+kj−d).
• Each edge (v(i), v(j)) in γ of type bd+1−ki is mapped to a single edge (v(i)1 , v (j) 1 ) (notice
that ki = kj).
• Each edge (v(i), v(j)) in γ of type V is mapped to a single edge (v(i)ki−1, v (j) kj−1).
As in Definition 3, the number of faces in γ is given by the number of loops in the double-line graph DL(γ).
Figure 7 is an example of the new single line and double line diagrams that appear in the correlation function Eθ [f(x1)f(x2)f(x3)f(x4)] for a two-hidden-layer network with biases.
Let us now extend Theorem 5 (the Feynman rules) to accomodate bias. The only difference is in the explicit normalization. Previously, each vertex had degree d+ 1 and contributed a factor of n−d/2 due to the explicit normalization of the function map. With biases, as a result of (26), a vertex with degree k contributes a factor of n(k−1)/2. Theorem 6. Let C(x1, . . . , xm) be a correlation function of a deep linear network with d hidden layers and with biases. Let γ ∈ Γ(C) be a Feynman diagram. The diagram represents a subset of
terms that contribute to C, and its asymptotic behavior is determined by the Feynman rules: the subset isO(nsγ ) where sγ = lγ− ∑m i=1 ki−1 2 , with ki the degree of the i-th vertex, and lγ the number of loops in the double-line diagram DL(γ). Furthremore, the correlation function is C = O(ns), where s = maxγ∈Γ(C) sγ .
We now introduce the following result, which generalizes Lemma 2 to the case of networks with bias. This Lemma is used above in the proof of Theorem 1. Lemma 3. Let C(x1, . . . , xm) be a correlation function for a deep linear network with biases. Let cγ be the number of connected components of a graph γ ∈ Γ(C). Then C = O(ns), where
s = max γ∈Γ(C)
cγ − m
2 . (27)
Proof. It is enough to show that each connected component γ′ in γ is bounded as O(nsγ′ ) where sγ′ ≤ 1− vγ′2 . The graph γ′ is a triangulation of a Riemann surface with f faces e edges, and vγ′ vertices of degrees k1, . . . , kvγ′ . The Feynman rules give
sγ′ = f − vγ′∑ i=1 ki − 1 2 . (28)
Using the relation e = ∑vγ′ i=1 ki 2 and the definition of the Euler character, χ = vγ′ − e+ f we have
sγ′ = χ− vγ′
2 . (29)
The diagram γ′ is a triangulation of a Riemann surface with at least one boundary, thus χ ≤ 1 and sγ′ ≤ 1− vγ′2 .

heading: C NON-LINEARITIES
text: In previous sections we presented the Feynman diagram method for computing the large width asymptotics of correlation functions. In this section we show that the method applies as-is for deep networks with ReLU non-linearities and all-equal inputs, as well as to networks with a single hidden layer, a broader class of non-linearities, and arbitrary inputs. Theorem 2 follows immediately from the results presented in this section.

heading: C.1 DEEP NETWORKS WITH RELU NON-LINEARITIES
text: The following result guarantees that the presence of ReLU non-linearities does not change the asymptotic upper bound compared with linear activations, when all inputs are the same. Theorem 7. Let fnl be a network function of the form (2) with ReLU activation. Let f be a network with the same architecture but with linear activation. Let C(x1, . . . , xm; fnl) be a correlation function of the deep network fnl with width n, and let f be a deep linear network with the same width and depth. Suppose that all inputs are the same, x1 = x2 = · · · = xm. Then
C(x1, . . . , xm; fnl) = O(C(x1, . . . , xm; f)) . (30)
Intuitively, we will rely on the fact that for ReLU networks we can, in some cases, treat the binary neuron activations as being statistically independent of the weights. This result is due to (Hanin & Nica (2018)). Given this result, we can bound the contribution of the binary activations, and the remaining Gaussian integral is equivalent to that found in a deep linear network. The proof does not work for correlation functions with non-equal inputs, because in that case the independence result of Hanin & Nica (2018) no longer holds.
Proof. We can write the network function as
fnl(x) = n −d/2V TDd(x)W d−1Dd−1(x)W d−2 · · ·D2(x)W 1D1(x)Ux , (31) Dj(x) := H(W j−1σ(W j−2 · · ·σ(Ux))) , j = 1, . . . , d . (32)
Here, H is the Heaviside step function acting elementwise on its vector argument. We now introduce the construction from Hanin & Nica (2018). Let ξj , ηj , j = 1, . . . , d be diagonal matrices of dimension n, whose diagonal elements are ±1-Bernoulli(p) variables with p = 12 . We define the new variables
Û := ξ1U , V̂ := ηdV , Ŵ j := ξj+1W jηj , j = 1, . . . , d . (33)
Let f̂nl be a network function with the same architecture as fnl but using the re-defined weights. We define
f̂nl(x) = n −d/2V̂ T D̂d(x)Ŵ d−1 · · · D̂2(x)Ŵ 1D̂1(x)Ûx (34)
= n−d/2V T ρdD̂d(x)W d−1ρd−1 · · · ρ2D̂2(x)W 1ρ1D̂1(x)Ux , (35) D̂j(x) := H(Ŵ j−1σ(Ŵ j−2 · · ·σ(Ûx))) , j = 1, . . . , d , (36)
ρj := ηjξj , j = 1, . . . , d . (37) The following was shown in Hanin & Nica (2018).
1. fnl and f̂nl are equal in distribution.
2. {D̂j(x), ρj , j = 1, . . . , d} are independent of {U, V,W 1, . . . ,W d−1}.
3. {D̂j(x), j = 1, . . . , d} are independent of each other for fixed x. The diagonal entries of each diagonal matrix D̂j(x) are independent, and take the values {0, 1} with probability 12 .
Now, the correlation function can be written as C(x1, . . . , xm; fnl) = Eθ [F (x1, . . . , xm; fnl)] = Eθ,η,ξ [ F (x1, . . . , xm; f̂nl) ] . (38)
The second equality follows from the fact that fnl d = f̂nl. Let us assume for now that the correlation function has no contracted derivatives. The we have
C(x1, . . . , xm; fnl) = Eθ,η,ξ [ f̂nl(x1) · · · f̂nl(xm) ] (39)
= 1
ndm/2 n∑ ~i Din∑ ~α E [ m∏ l=1 Vil,d ]( d−1∏ l′=1 E [ m∏ l′′=1 W l ′ il′′,l′+1,il′′,l′ ]) E [ m∏ l=1 Uil,1,αl ] E~i,~α . (40)
Here,~i = {i1,1, . . . , im,d} and ~α = {α1, . . . , αm}, and
E~i,~α := E [ d∏ l=1 m∏ l′=1 ρlil′,lD̂ l il′,l (xl′) ] m∏ l=1 xlαl . (41)
In writing the equation (40) we used the facts that D̂, ρ are independent of the parameters, and that parameters in different layers are independent. We now use Theorem 4 (Isserlis), which says that each of the expectation values over products of V , W j , and U elements is equal to a sum over permutations, where each term is a product over Kronecker delta functions — the covariance matrices of the parameters.
C(x1, . . . , xm; fnl) = 1
ndm/2 ∑ σ1,...,σd+1∈Sm n∑ ~i Din∑ ~α ∆̃~σ~i,~αE~i,~α . (42)
Here, ∆̃~σ~i,~α is a product of Kronecker delta functions. The precise form of this object will not be important for our purpose, though it can be easily derived using Theorem 4.
We can now bound the correlation function as follows.
|C(x1, . . . , xm; fnl)| ≤ 1
ndm/2 ∑ σ1,...,σd+1∈Sm n∑ ~i Din∑ ~α ∆̃~σ~i,~α ∣∣∣E~i,~α∣∣∣ (43) ≤ Emax ndm/2 ∑ σ1,...,σd+1∈Sm n∑ ~i Din∑ ~α ∆̃~σ~i,~α (44) = Emax|C(v, . . . , v; f)| . (45)
Here, Emax = max~i,~α ∣∣∣E~i,~α∣∣∣, and v ∈ RDin is the all-ones vector. The diagonal elements of
ρl, D̂l(x) are identical variables. Therefore, E~i,~α only takes an O(1) number of values in the large width limit. Note that each of these values are independent of n, and therefore Emax = O(1). We managed to bound the correlation function for fnl in terms of the correlation function for the corresponding linear network f with fixed inputs. The asymptotics of the linear-network correlation function do not depend on the inputs, and therefore this concludes the proof.

heading: C.2 SINGLE HIDDEN LAYER NETWORKS
text: For networks with a single hidden layer, defined by
fnl(x) = 1√ n n∑ i=1 Viσ(U T i x) , (46)
we can extend our asymptotic analysis to smooth non-linearities. We will show in Theorem 8 that for any correlation function C, we have
C(x1, . . . , xm; fnl) = O(C(x1, . . . , xm; f)) , (47) where f is a deep linear network of equal width and sufficient depth. Therefore, computing the asymptotics using Feynman diagrams for deep linear networks yields a bound on networks with a single hidden layer and smooth non-linearities.
Before delving into the proof of this claim, we consider a few simple examples. Let us begin with the correlation function Eθ [fnl(x)fnl(x′)].
Eθ [fnl(x)fnl(x′)] = 1
n n∑ i,j Eθ [ ViVjσ(U T i x)σ(U T j x ′) ]
(48)
= 1
n n∑ i Eθ [ σ(UTi x)σ(U T i x ′) ] = O(1) . (49)
Here we used two facts. First, the weights Vi are unaffected by the non-linearity, so we can carry out the V integral. Second, the summand in the last equation, E [ σ(UTi x)σ(U T i x ′) ] , is independent of both i and n because Ui are i.i.d. variables.
Next, consider the following correlation function. For simplicity, here we set the input dimension to be Din = 1, and we set all inputs to 1; this does not change the asymptotics.
Eθ [ dΘ
dt
] = n∑ j,k Eθ [ fnl ∂2fnl ∂Uj∂Vk ∂fnl ∂Uj ∂fnl ∂Vk ]
= 1
n2 n∑ i1,i2 Eθ [Vi1Vi2σ(Ui1)σ′(Ui2)σ′(Ui2)σ(Ui2)]
= 1
n2 n∑ i1 Eθ [ σ(Ui1) 2σ′(Ui1) 2 ]
= O(n−1) . (50) In the last line, we again used the fact that the summands are equal and independent of n.
In general, a correlation function C(x1, . . . , xm; fnl) can be reduced to the form
C = 1
nm/2 K∑ α=1 n∑ i1,...,irα S(α)i1...irα , (51)
S(α)i1...irα := EU [( σ(` α 1 )(UTi1xσ̄α(1)) · · ·σ (`αk1 )(UTi1xσ̄α(kα1 )) ) × · · ·×(
σ (`αm+1−kαrα ) (UTirαxσ̄α(m+1−kαrα )) · · ·σ (`αm)(UTirαxσ̄α(m)) ) F (α)(x1, . . . , xm) ] .
(52)
This form is obtained by carrying out the Vi integrals, as well as all the sums that can be trivially carried out due to the presence of Kronecker deltas. Here, the α sum represents a sum over all K terms that appear from performing the Vi integrals using Isserlis’ theorem; σ(`) denotes the `-th derivative of the non-linearity; kαs is the number of σ
(·)(UTisx) factors sharing the is index; σ̄α ∈ Sm is a permutation; and F (α) is a function whose exact form is not important for us. We will denote by rmax the maximum number of index sums appearing in (51), namely rmax := maxα rα.
Our approach to establishing the asymptotic scaling will be to first bound the maximum number of index sums appearing in any term in our correlation function, written in the form (51), and then to argue that the summands are bounded by an n-independent constant.
Let us introduce a family of diagrams, Γ′(C), which are different in general than the Feynman diagrams. A given diagram g ∈ Γ′(C) is constructed as follows. Definition 6. Let C(x1, . . . , xm; fnl) be a correlation function, where fnl is the output of a network with one hidden layer, defined in (46). The family Γ′(C) is the set of all graphs that have the following properties.
• Each derivative tensor in C is mapped to a vertex in the graph. • Each edge has a type that corresponds to one of the weight vectors U or V . • Each vertex has exactly one edge of V type. • If two derivative tensors are contracted in C, the graph must have at least one edge (of any
type) connecting the corresponding vertices for each contraction.
The reason for introducing this graphical structure is that it allows us to make two important statements. If we define c̃g the number of connected components in a graph g ∈ Γ′(C) and c̃max := maxg∈Γ′(C)c̃g then the following holds.
1. The maximal number of sums appearing in a correlation function of the form (51) is c̃max.
2. c̃max = cmax, where cmax is the maximal number of connected components in the family of Feynman diagrams corresponding to a correlation function C(x1, . . . , xm; fd), where fd is a deep linear network with d hidden layers, and d is large enough that none of the derivative tensors vanish.10
These two results, combined with a bound on the summands occurring in C(x1, . . . , xm; fnl) will establish the bound (47).
Lemma 4. A correlation function C has rmax = c̃max.
Proof. Each vertex corresponds to a derivative tensor, Tµ1...µk , which contains a single sum over paired Ui and Vi indices. If two vertices are connected by one or more edges, there is a Kronecker delta factor that sets the corresponding indices to be equal (due either to a derivative contraction, or to a V covariance), reducing the number of sums by one. The result is that any connected component of a graph g ∈ Γ′(C) corresponds to a single sum, and the total numer of index sums corresponding to a particular graph is c̃g , the number of connected components. As a result, the maximum number of sums corresponding to the collection of graphs Γ′(C) is c̃max.
Lemma 5. The maximal number of connected components, c̃max, over the collection of graphs Γ′(C) corresponding to C(x1, . . . , xm; fnl) is bounded by the maximal number of components, cmax, over the collection Γ(C) corresponding to C(x1, . . . , xm; fd) of Feynman diagrams, where fd is a deep linear network of sufficient depth d.
Proof. Let ne(no) be the number of even(odd) clusters in the cluster graphGC of C(x1, . . . , xm; fd). The cluster graph, GC is a subgraph of any graph g ∈ Γ′(C). We can thus think about the embedding of even and odd clusters into g. In any graph g ∈ Γ′(C), an even cluster may belong to its own connected component, while for odd clusters there must be an even number of them in any connected
10Note that any derivative tensor of f that has rank greater than d vanishes.
component. This is because an even (odd) cluster contains an even (odd) number of factors of V , which must be paired up in any connected component. We find that
c̃max = ne + no 2 ≤ cmax . (53)
The last inequality was used below Lemma 2 in the proof of Theorem 1.
Theorem 8. Let C(x1, . . . , xm; fnl) be a correlation function for a one hidden layer network. Let cmax be the maximal number of connected components over the collection of graphs Γ(C) corresponding to C(x1, . . . , xm; fd), where fd(x) is a deep linear network map, with with depth d greater than or equal to the maximum number of derivatives appearing in any single derivative tensor in C. Then C = O(ns) where s = cmax − m2 . Furthermore,
C(x1, . . . , xm; fnl) = O(C(x1, . . . , xm; f)) . (54)
Proof. The correlation function in (51) can be bound as
|C(x1, . . . , xm; fnl)| ≤ 1
nm/2 K∑ α=1 n∑ i1,...,ir ∣∣∣S(α)i1,...,irα ∣∣∣ . (55) For fixed inputs, S(α)i1,...,irα can take at most O(1) different values as a function of its indices, and the values are independent of n. This is because the variables Ui are identical. We define smax as the maximum value of |S(α)i1,...,irα | as a function of α and the indices. Combining this with the above lemmas we can then write.
|C(x1, . . . , xm; fnl)| ≤ Ksmaxncmax−m/2 . (56) The result of the theorem follows from Lemma 2.

heading: D CORRELATION FUNCTION ASYMPTOTICS
text: In this section we prove several general results about correlation function asymptotics in the large width limit. Throughout this section, we assume that Conjecture 1 holds.
D.1 VARIANCE ASYMPTOTICS
Conjecture 1 can be used to bound the variance of the integrands that appear inside correlation functions.
Lemma 6. Let C̃(x1, . . . , xm) = Eθ [Fθ(x1, . . . , xm)] be a correlation function, and let GC be the cluster graph of C with ne even components and no odd components. Assume Conjecture 1 holds such that C = O(nsC ), where sC = ne + no2 − m2 . Then Varθ [Fθ(x1, . . . , xm)] = O(n2sC ).
Proof. To bound the variance, it is enough to bound the correlation function
C̃(x1, . . . , x2m) := Eθ [Fθ(x1, . . . , xm)Fθ(xm+1, . . . , x2m)] , (57)
because Varθ [Fθ(x1, . . . , xm)] ≤ C̃(x1, . . . , xm, x1, . . . , xm). The correlation function C̃ has 2ne even clusters and 2no odd clusters. It follows from Conjecture 1 that C̃ = O(n2sC ).
As a corrolary, notice that if C = O(ns) according to Conjecture 1, then typical realizations of the integrand will also be O(ns). In other words, the asymptotic bound of Conjecture 1 holds for typical initializations, not just in expectation.
Table 3 shows empirical results for the variance of several functions. In all cases we find that the bound of Lemma 6 holds. For ∑ µ Varθ [∂µf(x1)∂µf(x2)] we prove a tight bound below in Appendix E.1 for deep linear networks.

heading: D.2 GRADIENT FLOW
text: The following results are used in the gradient flow calculations of Section 3.
Lemma 7. Let G′ be a graph with m′ vertices, n′e even components, and n′o odd components. Let G be a subgraph of G′ with m vertices, ne even components, and no odd components. Then s(ne, no,m) ≥ s(n′e, n′o,m′) where s(a, b, c) := a+ b−c2 .
Proof. It is enough to show that s(ne, no,m) does not increase if we (1) add a vertex to G, or (2) add an edge to G, because G′ can be obtained from G by performing such operations finitely many times. Adding a vertex to G changes ne 7→ ne, no 7→ no + 1, and m 7→ m+ 1, leaving s(ne, no,m) unchanged. Next, if we add an edge to G then m does not change, and there are 4 possibilities for how ne and no change.
1. The edge connects two even components. Then ne 7→ ne − 1, no 7→ no, and s(ne, no,m) decreases by 1.
2. The edge connects two odd components. Then ne 7→ ne+1, no 7→ no−2, and s(ne, no,m) does not change.
3. The edge connects an even component and an odd component. Then ne 7→ ne−1, no 7→ no, and s(ne, no,m) decreases by 1.
4. The edge connects two vertices that belong to the same component. In this case ne, no, and s(ne, no,m) do not change.
We now prove Lemma 1 giving the scaling of time derivatives of correlation functions at initialization. We prove the result for polynomial loss functions. Extension to more general loss functions requires interchanging the large width limit and the Taylor expansion of the loss, which we do not discuss.
proof (Lemma 1). Notice that
Eθ [ dkF (x1, . . . , xm)
dtk
] = Eθ  ∑
µ ∑ (x′,y′)∈Dtr ∂`(x′, y′) ∂f ∂f(x′) ∂θµ ∂ ∂θµ
k F (x1, . . . , xm)  .
(58) On the right-hand side we have a linear combination ∑ A αACA of correlation functions CA, where the coefficients depend on the training set labels. Here we used the polynomial loss assumption. For each correlation function CA, its integrand can be obtained from F by finitely many operations of the form (1) multiply the integrand by f(x) for some input x, and (2) act with a pair of contracted
derivatives on two of the derivative tensors. In the cluster graph representation, these two operations correspond to (1) adding a vertex, and (2) adding an edge. Therefore, the cluster graph GC of C = Eθ [F ] is a subgraph of the cluster graph GCA of each one of the correlation functions CA. By Lemma 7 we have that CA = O(nsC ) for all A, and therefore the bound applies to Eθ [ dkF/dtk ] as well.

heading: D.3 STOCHASTIC GRADIENT DESCENT
text: In this section we show that the asymptotics of a correlation function do not change under stochastic gradient descent (SGD) updates. Let C(x1, . . . , xm) = Eθ∼P0 [Fθ(x1, . . . , xm)] be a correlation function, where Fθ is the integrand which explicitly depends on the parameters θ. Under a single SGD step, the parameters are updated as θt+1 = θt − η∇Lt(θt), where L is the mini-batch loss at time t. Let Pt denote the distribution of parameters at step t, where P0 is the initial distribution. We define the evolved correlation function at step t to be
Ct(x1, . . . , xm) := Eθ∼Pt [Fθ(x1, . . . , xm)] . (59)
We have the following
Theorem 9. Let C be a correlation function for a network with linear activations, and assume that Conjecture 1 holds, namely C = O(nsC ) where sC is defined in the Conjecture. If Ct is the evolved correlation function after t SGD steps, then Ct = O(nsC ).
Proof. Notice that
Ct+1 = Eθ∼Pt+1 [Fθ] = Eθ∼Pt [ Fθ−η∇L(θ) ] . (60)
The integrand Fθ can be written as a product of derivative tensors of the form Tµ1...µk(x; θ), with contracted derivatives. Suppose that the network has d hidden layers. Then under an SGD step, we have
Tµ1...µl(x; θ − η∇Lt(θ)) = d+1∑ k=0 (−η)k k! ∑ ν1...νk ∂Lt ∂θν1 · · · ∂Lt ∂θνk Tµ1...µlν1...νk(x; θ) (61)
= d+1∑ k=0 (−η)k k! ∑ x′1,...,x ′ k∈DB [ ∂Lt ∂f(x′1) · · · ∂Lt ∂f(x′k)
× ∑
ν1,...,νk
Tν1(x ′ 1) · · ·Tνk(x′k)Tµ1...µlν1...νk(x)
] . (62)
Here DB is the mini-batch, and x′a are mini-batch samples. The k sum is truncated because higherorder derivatives of f vanish.
We can now see how taking a gradient descent step affects the cluster graph. After taking a step, each derivative tensor in Ct is replaced by a sum (over k, x′1, . . . , x ′ k). Each term in the combination of these sums is a correlation function, whose cluster graph is a subgraph of C. Therefore, by Lemma 7, Ct = O(nsC ).
We note that for general activation functions the sum in (62) may be infinite. In this case, to complete the proof we would need to show that the infinite sum obeys the same bound as each individual term in the sum. We leave this for future work.

heading: E APPLICATIONS
text: Here we present several applications of our Feynman diagram method for computing large width asymptotics. We assume throughout this section that Conjecture 1 holds.

heading: E.1 NEURAL TANGENT KERNEL
text: In this section we prove two results regarding the NTK at large width. We show that the kernel converges in probability, and that during gradient descent it is constant up to O(n−1) corrections. Theorem 10. The Neural Tangent Kernel Θ of a deep linear network converges in probability in the large width limit, and its variance is O(n−1).
Conjecture 1 is not sufficient for proving this theorem, as we need to use a more detailed Feynman diagram argument. Therefore, we only prove the theorem for the case of deep linear networks.
Proof. First, we will show that Var[Θ] = O(n−1). Given a model function f , the variance is given by
Var[Θ(x, x′)] = A(x, x′)−B(x, x′) , (63)
A(x, x′) = ∑ µ,ν E [ ∂f (1)(x) ∂θµ ∂f (2)(x′) ∂θµ ∂f (3)(x) ∂θν ∂f (4)(x′) ∂θν ] , (64)
B(x, x′) = ∑ µ,ν E [ ∂f (1)(x) ∂θµ ∂f (2)(x′) ∂θµ ] E [ ∂f (3)(x) ∂θν ∂f (4)(x′) ∂θν ] . (65)
Here f (1) = · · · = f (4) = f ; we introduced the superscripts so we can easily refer to individual factors. The crux of the argument is that the set of Feynman diagrams representing the expression (63) includes only connected graphs. Assuming this is the case, according to the bound (22) these diagrams will all scale as O(n(2−4)/2) = O(n−1) and so will the variance, completing the proof. To see why only connected graphs contribute, notice that
• A(x, x′) includes all diagrams in which the vertices corresponding to f (1), f (2) share an edge, and also the vertices f (3), f (4) share an edge (due to the explicit derivatives);
• B(x, x′) includes all diagrams in which all edges are either between f (1), f (2) or between f (3), f (4).
Therefore, in the full expression (63), the only remaining diagrams (i.e. the diagrams that do not cancel between the two terms) are those that include
• an edge connecting f (1), f (2),
• an edge connecting f (3), f (4), and
• an edge connecting one of f (1), f (2) to one of f (3), f (4).
These diagrams are all connected graphs, and this completes the proof that Var[Θ] = O(n−1). Here and above we have established that Var[Θ] = O(n−1) and E [Θ] = O(n0). More generally, let us consider a random variable O equal to the product of f and its derivatives, where the derivatives indices are fully summed in pairs. As we will now show, if E [O] = O(n0) and Var[O] = O(n−1) then (1) the limit limn→∞ E [O] exists, and (2) O converges in probability to this limit. In particular, the NTK is an example of such a random variable, and so this will conclude the proof that the NTK converges in probability.
First, consider the mean. There is a finite number of diagrams contributing to E [O], and each has a well-defined n scaling. Therefore, we can write
E [O] = kmax∑ k=0 ck nk
(66)
for some values of kmax ≥ 0 (integer) and c0, c1, . . . ∈ R. We see that the mean has a well-defined large n limit,
lim n→∞
E [O] = c0 . (67)
Next, let us show that O − E [O] converges in probability using Chebyshev’s inequality. Indeed, by the variance assumption there exists c̃ > 0 such that
P (|O − E [O] | > ) ≤ Var[O] 2 ≤ c̃ n 2 → 0 . (68)
Combining the facts that O − E [O]→p 0 and E [O]→ c0, we find that O converges in probability to c0.
Next, we show that the large width NTK is constant during training, and compute the asymptotics of the higher-order terms. The following argument is phrased for deep linear networks. More generally, the same argument holds for deep networks with smooth non-linear activations under the additional assumption that the large width limit and Taylor series can be exchanged (note that for such networks, the network function is analytic in the weights).
Theorem 11. Let f(x) be the network output of a deep linear network with MSE loss L. Let Θt(x, x′) be the Neural Tangent Kernel at SGD step t, for some inputs x, x′. Then in the large width limit, the kernel is constant in t in expectation, and
Eθ∼P0 [Θt(x, x′)−Θ0(x, x′)] = O(n−1) , (69) Varθ∼P0 [Θt(x, x ′)−Θ0(x, x′)] = O(n−2) . (70)
Proof. It is enough to show that Eθ∼P0 [Θ1(x, x′)−Θ0(x, x′)] = O(n−1) . It then follows from Theorem 9 that Eθ∼P0 [Θt+1(x, x′)−Θt(x, x′)] = O(n−1) for all t, and therefore
Eθ∼P0 [Θt(x, x′)−Θ0(x, x′)] = t−1∑ t′=0 Eθ∼P0 [Θt′+1(x, x′)−Θt′(x, x′)] = O(n−1) , (71)
concluding the proof. We have Eθ∼P0 [Θ1(x, x′)−Θ0(x, x′)] = d+1∑ k,l=0 k+l≥1 (−η)k+l k!l! ∑ µ Eθ∼P0 [ ∑ µ1,...,µk ∂L ∂θµ1 · · · ∂L ∂θµk Tµµ1...µk(x)
∑ ν1,...,νl ∂L ∂θν1 · · · ∂L ∂θνl Tµν1...νl(x ′)
] .
(72)
The fact that the k, l sums are truncated at d+ 1 follows from using linear activations, as higher-order derivatives of the network function vanish in this case. All terms in the sum over k, l include a tensor product of the form ∑ µ,~µ,~ν Tµµ1...µk(x)Tµν1...νl(x
′)(· · · ) with either k ≥ 1 or l ≥ 1, where (· · · ) stands for additional derivative tensor factors. Therefore, all terms in the k, l sum are correlation functions that have a cluster of size at least 3, including T (x), T (x′), and at least one other tensor contracted through the µ1 or ν1 index. It follows from the Conjecture that each term in the sum is O(n−1). Lemma 6 then implies that the variance of these updates is O(n−2).

heading: E.2 WIDE NETWORK EVOLUTION IS LINEAR IN THE LEARNING RATE
text: In this section we prove that, at large width, the NTK determines the evolution of the network function not just for continuous-time gradient descent but also for discrete-time gradient descent. A similar result holds for stochastic gradient descent, using a stochastic kernel.11 Again we prove the deep linear case explicitly, but the result holds for deep networks with smooth non-linear activations under the additional assumption that the large width limit and Taylor series can be exchanged.
11The perspective presented here helps understand the results of (Lee et al. (2019)) where it was observed empirically that linearized evolution is a good description of wide networks even for relatively large learning rates.
Theorem 12. Let f(x) be the network output of a deep linear network, and let ft(x) be the evolved function after t gradient descent steps, defined by θt+1 = θt − η∇L(θt). In the large width limit, each gradient descent step update of ft is linear in the learning rate η. Furthermore,
Eθ [ft+1(x)− ft(x)] = −η ∑ (x′,y′)∈Dtr Eθ [ Θ0(x, x ′) ∂`t(x ′, y′) ∂f ] +O(n−1) . (73)
Here, Θ0 is the NTK at initialization and `t is the single-sample loss at time t.
Proof. For a deep linear network, under a single gradient descent step we have
Eθ [ft+1(x)− ft(x)] = d+1∑ k=1 (−η)k k! ∑ µ1,...,µk Eθ [ ∂kft(x) ∂θµ1 · · · ∂θµk ∂Lt ∂θµ1 · · · ∂Lt ∂θµk ] (74)
= −η ∑ (x′,y′)∈Dtr Eθ [ Θt(x, x ′) ∂`t(x ′, y′) ∂f ]
+ d+1∑ k=2 (−η)k k! ∑ µ1,...,µk Eθ [ ∂kft(x) ∂θµ1 · · · ∂θµk ∂Lt ∂θµ1 · · · ∂Lt ∂θµk ] . (75)
First, consider the sum over k in (73). Each term in the sum is a correlation function for which the cluster graph contains a connected subgraph of size at least 3, and is therefore O(n−1) by Lemma 7 and Theorem 1. In the remaining O(η) term, by the same argument as Theorem 11 we can replace Θt = Θ0 +O(n−1) in the correlation function. The result is equation (73).
As mentioned above, we note that the proof goes through when using stochastic gradient descent updates, with the difference that in (73) we should sum over mini-batch samples instead of over the entire training set.

heading: E.3 SPECTRAL PROPERTIES OF THE NTK AND THE HESSIAN
text: With an eye towards understanding the structure of the loss landscape at large width and as another example use case of our approach, we investigate the relation between the spectra of the Hessian, and the NTK.
Among other observations, we present an argument that for a network with d hidden layers and MSE loss, the top Din ordered eigenvalues of the Hessian, λ (H) i , are related to those of the kernel λ (Θ) i by
Eθ [ λ (H) i − λ (Θ) i ] = { O(n−1/2) , d = 1 O(n−1) , d > 1 . (76)
The remaining Hessian eigenvalues vanish at large width as O(n−1/2) for one hidden layer networks and O(n−1) for deeper networks.12 This is experimentally corroborated in Figure 9. The Hessian of a general loss takes the form,
Hµν = ∑
(x,y)∈Dtr
[ ∂2`(x, y) ∂f2 ∂f(x) ∂θµ ∂f(x)
∂θν︸ ︷︷ ︸ Aµν
+ ∂`(x, y)
∂f
∂2f(x)
∂θµ∂θν︸ ︷︷ ︸ Bµν
] . (77)
The moments of H are a useful way to understand its spectrum. As in other examples that we have seen, traces of powers of B lead to connected diagrams with higher and higher numbers of vertices, and so scale as increasingly negative powers of n. More explicitly, traces involving powers of B contain multiple contracted derivative tensors. As a result, taking the average of these traces over the
12Empirically we actually find the even stronger bound Eθ [ λ (H) i − λ (Θ) i ] = O(n−1) for the top Din eigenvalue differences and O(n−1/2) for the remaining eigenvalues in the one hidden layer case. We can gain insight into this improved scaling through the perspective of degenerate eigenvalue perturbation theory (Trefethen & Bau (1997)), but this is outside the scope of the current presentation.
Figure 8: A prototypical diagram corresponding to Eθ
[ Tr ( B2 )]
weight distribution leads to correlation functions where the associated cluster graph, GC , contain subgraphs with at least three connected vertices. By Lemma 7 and Conjecture 1 these correlation functions vanish at infinite width. For example Eθ [Tr (AB)] = O(1/n). There is a single exception to this, as can be easily seen in figure 8, Eθ [ Tr(B2) ] is O(1). Thus most moments of the Hessian are equal to moments of A in expectation.13
Eθ [Tr (Hm)] = { Eθ [Tr (Am)] ,m 6= 2 Eθ [ Tr ( A2 )] + Eθ [ Tr ( B2 )] ,m = 2 +O(1/n) . (78)
What’s more, we can relate moments of A to those of the kernel, as both are built out of two logit derivatives. Explicitly,
Tr(Am) = Tr((MΘ)m) , (79) and so the moments of the Hessian are also related to those of the kernel.14
Eθ [Tr (Hm)] = { Eθ [Tr ((MΘ)m)] ,m 6= 2 Eθ [ Tr ( (MΘ)2 )] + Eθ [ Tr ( B2 )] ,m = 2 +O(1/n) . (80)
Here, we have defined,
M(xa, xb) = δab ∂2`(xa, ya)
∂f2 : (xa, ya) ∈ Dtr . (81)
For the case of MSE loss, we can go even further. In that case, M(xa, xb) = δab and B decays to zero during training. We thus have,
Initially: Eθ [Tr (Hm)] = { Eθ [Tr (Θm)] ,m 6= 2 Eθ [ Tr ( Θ2 )] + Eθ [ Tr ( B2 )] ,m = 2 +O(1/n) At late times: Eθ [Tr (Hm)]→ Eθ [Tr (Θm)] ∀m. (82)
These results indicate that the only difference between the spectra of the Hessian and the NTK come from B and that B must have eigenvalues which scale as 1/ √ dim(B). As dim(B) = O(n) for one hidden layer networks and O(n2) for deep networks we are left with the relation (76) between the eigenvalues of Θ and H . As the network trains, the difference between these eigenvalues gets even smaller. These results are confirmed experimentally in Figure 9 and Figure 10.
Multi-class. This story can be generalized to multi-class neural network maps. In this case the logit is a map fA : RDin → RNclass , and the NTK has two class indices.
ΘAB(x, x′) = ∂fA(x) ∂θµ ∂fB(x′) ∂θµ . (83)
Expression (80) relating the moments of the Hessian to the NTK still holds in this context provided we take Tr ((MΘ)m)→ Tr ( (MABΘ AB)m ) , with the more general matrix,
MAB(xa, xb) = δab ∂`(xa, ya)
∂fA∂fB : (xa, ya) ∈ Dtr . (84)
13For the first moment in linear or ReLU networks, TrB = 0. 14Here we have argued for relations relating the mean of moments. It is not too difficult to see that these
relations will also hold for typical realizations. This follows from Lemma 6.
At large width, the NTK approaches the identity matrix in class space, ΘAB → N−1classδAB Tr (Θ) (Jacot et al. (2018)). This implies that the NTK spectrum consists of Nclass repeated copies. This has consequences for the spectrum of the Hessian at large width. For instance, for the case of MSE error it also implies Nclass repeated copies of the Hessian spectrum (See Figure 9b). It is intriguing to think this could serve as a path towards understanding the emergence of the Nclass (or Nclass − 1) large eigenvalues and corresponding subspace observed in the Hessian spectrum (Sagun et al. (2016; 2017); Gur-Ari et al. (2018); Ghorbani et al. (2019); Papyan (2019)).

heading: E.4 HIGHER-ORDER NETWORK EVOLUTION
text: In this section we will explain how to compute higher-order corrections to training dynamics. In principle, this prescription allows one to compute model dynamics as an expansion in 1/n to arbitrary order. We apply this explicitly to compute the O(1/n) training dynamics. In Figure 11, we present experimental confirmation of our predictions for the evolution of the NTK. In the main text, we presented these results for the special case of gradient flow with MSE loss.
So far we have mostly been using diagrammatic techniques to understanding the leading order scaling of correlation functions. It is interesting to try and understand finite width networks by asking how training dynamics are modified beyond the leading order asymptotics. There are three clear sources of corrections to the leading behavior.
1. The initial kernel, Θ0, receives finite width corrections.
2. The linear (in learning rate) update for the network function, equation (73), is modified away from infinite width.
3. The kernel is not constant at finite width.
The first source of corrections is automatically taken into account in typical empirical settings, as the finite-width Θ0 can be computed explicitly. The other sources of corrections are non-trivial, and will be the focus of this section. We begin by explaining how to take into account the non-constancy of the kernel order by order in 1/n, while maintaining the continuous time approximation. Next, we back away from the continuous time limit and write down the full discrete evolution. We introduce a method to compute arbitrary order corrections, and explicitly work out the corrections at order 1/n for MSE loss.

heading: E.4.1 CONTINUOUS TIME.
text: Continuous time evolution of the model function in neural networks is governed by the differential equation
df(x; t)
dt = − ∑ (x′,y′)∈Dtr Θ(x, x′; t) ∂`(x′, y′) ∂f . (85)
As we have discussed at length, in the large width limit, this equation simplifies and Θ(t) is asymptotically constant (Jacot et al. (2018)). Our goal is to move beyond this leading behavior at large width and solve (85) order-by-order in a 1/n expansion. In Section 3 we described how to compute the O(n−1) corrections for the case of MSE loss. Here we explain how to handle corrections more generally as well as giving a more detailed discussion of the MSE case.
The network map f and kernel, Θ are members of the following family of operators.
Os(x1, . . . , xs; t) := ∑ µ ∂Os−1(x1, . . . , xs−1; t) ∂θµ ∂f(xs; t) ∂θµ , s ≥ 2 , (86)
with O1 := f . Here, as above, O2 = Θ is the kernel. With a general loss function, these operators satisfy.
dOs(x1, . . . , xs; t)
dt = − ∑ (x′,y′)∈Dtr Os+1(x1, . . . , xs, x ′; t) ∂`(x′, y′; t) ∂f , s ≥ 1 . (87)
Equations (86) and (87). Give an infinite tower of first order ODEs, the solution of which gives the time evolution of the network map and the kernel.
df(x1; t)
dt = − ∑ (x,y)∈Dtr Θ(x1, x; t) ∂`(x, y) ∂f
dΘ(x1, x2; t)
dt = − ∑ (x,y)∈Dtr O3(x1, x2, x; t) ∂`(x, y) ∂f
dO(1)(x, x′, x′′; t) dt
= − ∑
(x,y)∈Dtr O4(x1, x2, x3, x; t)
∂`(x, y)
∂f (88)
...
Solving this infinite tower is not feasible. If we wish to work to a given order in an expansion in 1/n, however, there is a dramatic simplification, which makes a solution possible. Firstly, we can truncate these equations. To see this note that the operators Os contain s contracted derivative tensors. As a result, by Lemma 7 and Conjecture 1, correlation functions involving the operators, Os satisfy
Eθ [Os(x1, . . . , xs; t)F (~x; t)] =  O ( n 2−s 2 ) , s even O ( n 1−s 2 ) , s odd , (89)
where F (~x; t) is arbitrary additional contribution to the integrand.
Thus, if we wish to work to solve for the time evolution up to corrections which scale as O(n−r) we can truncate the tower at s = 2r and set O2r(t) to be equal to its initial value. Note that the leading order solution, (85), is the result of this procedure with r = 1 and the results presented in the main text for the O(n−1) evolution correspond to r = 2. The truncation provides a dramatic simplification, however it is not immediately clear how to solve even the truncated differential equations in (87). We now describe how to organize the perturbative expansion of the operators Os(t) (including f and Θ) in such a way that the differential equations become tractable.
The central idea is to write each operator, Os(t), as an expansion.
Os(x1, . . . , xs; t) = ∞∑ r=b s−12 c O(r)s (x1, . . . , xs; t) (90)
where each order O(r)s captures the O(n−r) evolution of Os. For example, f(x1; t) = f (0)(x1; t) + f (1)(x1; t) + · · ·
Θ(x1, x2; t) = Θ (0)(x1, x2; t) + Θ (1)(x1, x2; t) + · · · O3(x1, x2, x3; t) = O (1) 3 (x1, x2, x3; t) +O (2) 3 (x1, x2, x3; t) + · · · (91)
...
The notation O(r)s means both that any correlation function containing O (r) s is O(n−r) and, by Lemma 6, that typical realizations of the operators scale as O(n−r).15
Once the operators Os are organized in this way, solving the differential equations (85) is tractable. As the differential equations describing the evolution of O(r)s for r > 0 only depend on the time dependent solutions of O(r−1)s , we can iteratively solve for the O (r) s order by order. For example, in Section 3, we used the leading order solution for f and Θ to solve for Θ(1).
In principle this procedure can be extended to arbitrary order in 1/n. Before going onto explain the finite step corrections to this procedure, we reproduce the results of section 3 in more detail.
15Note that in Section 3 we used the alternate notation Θ1 for Θ(1).
MSE continuous time example. The MSE loss is,
LMSE = 1
2 ∑ (x,y)∈Dtr (f(x)− ya)2 (92)
As such the update equation simplifies to,
df(x; t)
dt = − ∑ (x′,y′)∈Dtr Θ(x, x′; t) (f(x′; t)− y′) , (93)
The leading order solution to equation (93) is exponential, kernel evolution,
f (0)(t) = y + e−tΘ0(f0 − y) (94) Θ(0)(t) = Θ0 . (95)
Here we are using a condensed notation where f0 := f (0)(0), Θ0 := Θ(0)(0). Equation (94) is a vector equation with f0, f (0)(t), and y are vectors over the training set and the leading order kernel Θ0 is a square matrix over the same space.
As discussed above to study the O(1/n) evolution, we can truncate the set of equations in (87) at s = 4, and set O4(t) = O (1) 4 (t) = O4(0) up to corrections which scale as O(n−2).
Moving up a rung in (87), we have
dO3(x1, x2, x3; t)
dt = − ∑ (x,y)∈Dtr O4(x1, x2, x3, x; t) (f(x; t)− y) . (96)
To solve for O(1)3 (t), i.e. neglecting O(n−2) corrections, we can plug in the leading order approximations, to O4(t) and f(t).
dO (1) 3 (x1, x2, x3; t)
dt = − ∑ (x,y)∈Dtr O4(x1, x2, x3, x; 0) ( f (0)(x; t)− y ) . (97)
This gives,
O (1) 3 (x1, x2, x3; t) = O3(x1, x2, x3; 0)− ∫ t 0 dt′ ∑
(x,y)∈Dtr O4(x1, x2, x3, x; 0)
( f (0)(x; t)− y ) .
(98)
Using the explicit form of f (0)(t) we can write,
O (1) 3 (~x; t) = O3(~x; 0)− ∑ a O4(~x; 0) · ( Θ−10 (1− e−tΘ0)(f0 − y) ) . (99)
Here we are again using a condensed notation. ~x = (x1, x2, x3). f0 and y are vectors over the training set while Θ0 is a square matrix, and O4(~x; 0) is also a vector over the training set, with the value O4(~x, x′; 0) on the point x′ ∈ Dtr. Moving one more step up the ladder gives
Θ(x1, x2; t) = Θx1,x2;0 + Θ (1)(x1, x2; t) +O(n−2) , Θ(1)(x1, x2; t) = − ∫ t
0
dt′ ∑
(x,y)∈Dtr O
(1) 3 (x1, x2, x; t
′) ( f (0)(x; t′)− y ) . (100)
Plugging in O(1)3 , and using the eigen-decomposition of Θ0 to perform the integrals gives, Θ(1)(~x; t) =− ∑ i 1 λi (O3(~x; 0) · êi)(∆f0 · êi) [ 1− e−tλi ] + ∑ ij 1 λj (êTi O4(~x; 0)êj)(∆f0 · êi)(∆f0 · êj) [ 1− e−tλi λi − 1− e −t(λi+λj) λi + λj ] .
(101)
Here, we have introduced the eigenvectors, êi, which are vectors over the training dataset and eigenvalues λi of the initial kernel, Θ0. The vector ~x = (x1, x2) as Θ depends on two inputs. O3(~x; 0) is a vector over the dataset while O4(~x; 0) is a square matrix and ∆f0 := f0 − y is a vector over the training data.
Finally, we can plug in the expression (100) into (93) and give the sub-leading behavior for f(t)
df(x; t)
dt = − ∑ (x,y)∈Dtr ( Θ(x, x′; 0) + Θ(1)(x, x′; t) ) (f(x′; t)− y) +O(n−2) (102)
f(t) = y + e−Θ0t ( 1− ∫ t
0
dt′et ′Θ0Θ(1)(t′)e−t ′Θ0 ) (f0 − y) . (103)
This completes the full O(n−1) time dependance.This prediction for the O(1/n) time dependence is confirmed experimentally in Figure 11.
Note, one consequence of (101) is that the O(1/n) corrections to Θ go to a constant at late times.
Θ(1)∞ := lim t→∞ Θ(1)(t) = − ∑ i 1 λi (O3(~x; 0) · êi)(∆f0 · êi) + ∑ ij
1
λi(λi + λj) (êTi O4(~x; 0)êj)(∆f0 · êi)(∆f0 · êj) .
(104)
And the asymptotic evolution of ∆f(t) is again constant kernel evolution, but with a corrected kernel.
f(t)→ y + e−t(Θ0+Θ(1)∞ ) (f0 − y) . (105)
It is worth noting that these predictions are quite detailed, giving the full time dependence, rather than just the overall scaling with width. In deriving these results and in the discrete time analysis below, we implicitly make an assumption that we can control the discrete time corrections or in the continuous time case that f(t), Θ(t) and all the Os(t) are differentiable. This relies on being able to differentiate through the activation function. For non-smooth activations such as ReLU, this assumption is suspect, and indeed the evolution is more subtle in the ReLU case and a full analysis is left to future work. For this reason we present experimental results for networks with smooth activation functions.

heading: E.4.2 DISCRETE TIME.
text: With the 1/n corrections to continuous time evolution under our belts, it is possible to also keep track of corrections coming from the discrete update step. The essential point, that it is possible to solve iteratively, order by order in 1/n is not changed. To see this we can look at how the update equations are modified. Beginning with the network output update equation, we have,
ft+1(x)− ft(x) = −η ∑
(x′,y′)∈Dtr Θt(x, x
′) ∂`t(x ′, y′) ∂f
+ η2
2 ∑ µν ∂2ft(x) ∂θµ∂θν gµt g ν t − η3 3! ∑ µνρ ∂3ft(x) ∂θµ∂θν∂θρ gµt g ν t g ρ t + · · · . (106)
In Section E.2 we argued that the O(η2) terms vanish at large width. In more detail, each factor of the gradient, gµ, contains a contracted derivative tensor. Thus, the cluster graph for correlation functions involving the η2 term will always contain three contracted vertices and thus the correlation functions will scale as O(n−1). The η3 term will give four vertices, and thus also scale as O(1/n), the η4 contribution will be O(1/n2) and so on. Just as in the continuous time case this equation is still solvable order by order in 1/n. What we mean by solvable is the following. All terms on the RHS appearing with ft are already higher order as a result of their derivative structure. This means, just as in the continuous case, we can use the lower order solution of ft in the gradient terms on the RHS to solve for ft on the LHS.
Similarly, the update equation for Θ receives discrete time modifications. Θt+1(x1, x2)−Θt(x1, x2) = −η ∑ a O3;t(x1, x2, x) ∂`t(x, y) ∂f + η2 2 ∑ µν ∂2Θt(x1, x2) ∂θµ∂θν gµt g ν t + · · · .
(107)
Here we have adopted a notation where, ~x = (x1, . . . , xs). Just as for f , increased powers in η are increasingly suppressed in n and we can iteratively solve this equation using solutions at leading orders in n to solve for sub-leading behavior.
More generally the tower of differential equations defined recursively in (87) becomes Os;t+1(~x)−Os;t(~x) =− ∑
(x′,y′)∈Dtr Os+1;t(~x, x
′) ∂`(x′, y′; t)
∂f
+ ∞∑ k=2 ηk k! ∑ µ1,µ2,...,µk
∂kOs;t+1(~x)
∂θµ1∂θµ2 . . . ∂θµk gµ1gµ2 · · · gµk , s ≥ 1 . (108)
As in the continuous time case, at a given order in n, we can truncate this tower and use lower order solutions to solve for the time evolution up to the desired order. To ground this discussion in a concrete use case, we again walk through the 1/n corrections for MSE loss, now in discrete time.

heading: MSE discrete time example.
text: For discrete time evolution, the leading order behavior of the model map is
f (0) t = (1− ηΘ0)t (f0 − y) (109)
Θ (0) t = Θ0 . (110)
At next order we can proceed by solving the truncated set of equations in (108). The first equation, for O4;t still gives a constant solution,
O (1) 4;t = O4;0 . (111)
Here, both the discrete and continuous time terms that would appear on the RHS are O(1/n2). Next, we must solve for O3;t. The discrete time update is,
O3;t+1(~x)−O3;t(~x) = −η ∑
(x′,y′)∈Dtr O4;t(~x, x
′) (ft(x ′)− y′) + η
2
2 ∑ µν ∂2O3;t(~x) ∂θµ∂θν gµt g ν t + · · · .
(112)
To order 1/n we can drop the order η2 and higher terms, as they contain expressions with 5 or more contracted derivative tensors. Thus at this order, we are left with the discrete analogue of (96).
O3;t+1(~x)−O3;t(~x) = − ∑
(x′,y′)∈Dtr O4;t(~x, x
′) (ft(x ′)− y′) . (113)
which we can sum to get the discrete version of (98),
O (1) 3;t (~x) = O3;0(~x)− η t∑ t′=1 ∑ (x′,y′)∈Dtr O4;0(~x, x ′) (f0;t(x ′)− y′) . (114)
Explicitly, this takes the form
O (1) 3;t (~x) = O3;0(~x)−O4;0(~x) · ( Θ−10 (1− ηΘ0)(1− (1− ηΘ0)t) (f0 − y) ) . (115)
Here we have adopted notation similar to above, where O4;0(~x) is a vector over the training dataset.
So far, this procedure has been a discrete analogue of what we have done in the continuous time case, however as we move on to compute Θ and f we will have to keep track of the novel corrections, which vanish in continuous time. Explicitly, at order 1/n the discrete update for Θt is given by,
Θ (1) t+1(~x)−Θ (1) t (~x) = −η ∑ a O (1) 3;t (~x, x ′) ( f (0) t (x ′)− y′ )
+ η2
2 ∑ x′,x′′,y′,y′′∈Dtr ∑ µν ∂2Θ0(~x) ∂θµ∂θν ∂f0(x ′) ∂θµ ∂f0(x ′′) ∂θν ( f (0) t (x ′)− y′ )( f (0) t (x ′′)− y′′ ) .
(116)
This can be summed to give Θ(1)t . To move up to the neural network map itself there is one additional complication. The O(η2/n) term,∑ µν ∂2ft(x) ∂θµ∂θν ∂ft(x ′) ∂θµ ∂ft(x ′′) ∂θν , in (106) has non trivial time dependence. We can deal with this just as we have been doing, by taking an extra time derivative and integrating. Defining,
O1,1;t(x1, x2, x2) := ∑ µν ∂2ft(x1) ∂θµ∂θν ∂ft(x2) ∂θµ ∂ft(x3) ∂θν . (117)
We have
O (1) 1,1;t+1(~x)−O (1) 1,1;t(~x) = −η ∑ (x′,y′)∈Dtr ∑ µ ∂O1,1;0(~x) ∂θµ ∂f (0) t (x ′) ∂θµ ( f (0) t (x ′)− y′ ) , (118)
with the solution,
O (1) 1,1;t(~x) = O1,1;0(~x)− ∑ µ ∂O1,1;0(~x) ∂θµ ∂f0 ∂θµ · ( Θ−10 (1− ηΘ0)(1− (1− ηΘ0)t)∆f0 ) a . (119)
This can in turn be plugged into the update equation for ft. ft+1(x)− ft(x) = −η ∑ a ( Θ0(x, xa) + Θ (1) t (x, xa) ) (ft(xa)− ya)
+ η2
2 ∑ a,b O (1) 1,1;t(x, xa, xb) ( f (0) t (xa)− ya )( f (0) t (xb)− yb ) − η 3
3! ∑ a,b,c ∂3f0(x) ∂θµ∂θν∂θρ ∂f0(xa) ∂θµ ∂f0(xb) ∂θν ∂f0(xc) ∂θρ (ft(xa)− ya) (ft(xb)− yb) (ft(xc)− yc) .
(120)
Here (xa, ya) are elements of the training set Dtr and summed over. This equation can be solved to give f (1)t .
f (1) t = (1− ηΘ0)t t∑ t′=1 (1− ηΘ0)−(t ′+1) [ −ηΘNLO, t′(1− ηΘ0)t ′ ∆f0 + Disct′ ] . (121)
Where Disct contains the discrete derivative updates at O(1/n).
Disct := η2
2 ∑ a,b O (1) 1,1; t(xa, xb) (ft(xa)− ya) (ft(xb)− yb)
− η 3
3! ∑ a,b,c ∂3f0 ∂θµ∂θν∂θρ ∂f0(xa) ∂θµ ∂f0(xb) ∂θν ∂f0(xc) ∂θρ (ft(xa)− ya) (ft(xb)− yb) (ft(xc)− yc) .
(122)
These expressions may look fairly intimidating. The key point is that all terms in the summand in (116) and thus (121) are known functions of time and initial data, just as in the continuous time setting.

Review summaries:
The paper investigates the asymptotic behavior of correlation functions of wide networks. The main part of the work revolves around Conjecture 1 which assesses how the correlation function scales depending on the number of connected components of even/odd size of the constructed cluster graph. The authors then study training dynamics of wide networks under gradient flow and (stochastic) gradient descent and present empirical evidence to support their claims.