Wide Neural Networks of Any Depth Evolve as Linear Models Under Gradient Descent

TL;DR

This paper demonstrates that wide neural networks behave like linear models during training, with their dynamics governed by a first-order Taylor expansion, and their predictions linked to Gaussian processes, providing a simplified theoretical understanding.

Contribution

The work shows that in the infinite width limit, neural networks' training dynamics are linear and can be accurately approximated by a first-order Taylor expansion, bridging neural networks and Gaussian processes.

Findings

Linearized models accurately predict training dynamics.

Wide neural networks' predictions match Gaussian process priors.

Empirical results show strong agreement even for finite networks.

Abstract

A longstanding goal in deep learning research has been to precisely characterize training and generalization. However, the often complex loss landscapes of neural networks have made a theory of learning dynamics elusive. In this work, we show that for wide neural networks the learning dynamics simplify considerably and that, in the infinite width limit, they are governed by a linear model obtained from the first-order Taylor expansion of the network around its initial parameters. Furthermore, mirroring the correspondence between wide Bayesian neural networks and Gaussian processes, gradient-based training of wide neural networks with a squared loss produces test set predictions drawn from a Gaussian process with a particular compositional kernel. While these theoretical results are only exact in the infinite width limit, we nevertheless find excellent empirical agreement between the predictions of the original network and those of the linearized version even for finite practically-sized networks. This agreement is robust across different architectures, optimization methods, and loss functions.

Tables

Table 1. Table S1 : Wide Residual Network architecture from  Zagoruyko and Komodakis [ 14 ] . In the residual block, we follow Batch Normalization-ReLU-Conv ordering.

group name	output size	block type
conv1	32 $\times$ 32	[3$\times$3, channel size]
conv2	32 $\times$ 32	$\times$ N
conv3	16 $\times$ 16	$\times$ N
conv4	8 $\times$ 8	$\times$ N
avg-pool	1 $\times$ 1	[8 $\times$ 8]

Equations

$$\{\begin{array}{ll}{h}^{l+1}& =x^l{W}^{l+1}+b^{l+1}\\ x^{l+1}& =\varphi \left(h^{l+1}\right)\end{array}\text{and}\{\begin{array}{ll}{W}_{i,j}^l& =\frac{{\sigma }_{\omega }}{\sqrt{n_l}}{\omega }_{ij}^l\\ b_j^l& ={\sigma }_b{\beta }_j^l\end{array},$$

$${\dot{\theta }}_t=-\eta {{\nabla }_{\theta }{f}_t(\mathcal{X})}^T{\nabla }_{f_t(\mathcal{X})}\mathcal{L}$$

$${\dot{f}}_t(\mathcal{X})={\nabla }_{\theta }{f}_t(\mathcal{X}){\dot{\theta }}_t=-\eta {\hat{\Theta }}_t(\mathcal{X},\mathcal{X}){\nabla }_{f_t(\mathcal{X})}\mathcal{L}$$

$${\hat{\Theta }}_t$$

$$f_t^{\text{lin}}(x)\equiv f_0(x)+{{\nabla }_{\theta }{f}_0(x)|}_{\theta ={\theta }_0}{\omega }_t,$$

$${\dot{\omega }}_t=-\eta {{\nabla }_{\theta }{f}_0(\mathcal{X})}^T{\nabla }_{f_t^{\text{lin}}(\mathcal{X})}\mathcal{L}$$

$${\dot{f}}_t^{\text{lin}}(x)=-\eta {\hat{\Theta }}_0(x,\mathcal{X}){\nabla }_{f_t^{\text{lin}}(\mathcal{X})}\mathcal{L}.$$

$${\omega }_t=-{{\nabla }_{\theta }{f}_0(\mathcal{X})}^T{\hat{\Theta }}_0^{-1}\left(I-e^{-\eta {\hat{\Theta }}_{0}t}\right)\left(f_0(\mathcal{X})-\mathcal{Y}\right),$$

$$f_t^{\text{lin}}(\mathcal{X})=(I-e^{-\eta {\hat{\Theta }}_{0}t})\mathcal{Y}+e^{-\eta {\hat{\Theta }}_{0}t}{f}_0(\mathcal{X}).$$

$${\mu }_t(x)={\hat{\Theta }}_0(x,\mathcal{X}){\hat{\Theta }}_0^{-1}\left(I-e^{-\eta {\hat{\Theta }}_{0}t}\right)\mathcal{Y}$$

$${\gamma }_t(x)=f_0(x)-{\hat{\Theta }}_0\left(x,\mathcal{X}\right){\hat{\Theta }}_0^{-1}\left(I-e^{-\eta {\hat{\Theta }}_{0}t}\right)f_0(\mathcal{X}).$$

$${\mathcal{K}}^{i,j}(x,x^{\prime })=\underset{\min \left(n_1,\dots ,n_L\right)\to \infty }{\lim }\mathbb{E}\left[f_0^i(x)\cdot f_0^j(x^{\prime })\right],$$

$$\mu (x)=\mathcal{K}(x,\mathcal{X}){\mathcal{K}}^{-1}\mathcal{Y},\Sigma (x)=\mathcal{K}(x,x)-\mathcal{K}(x,\mathcal{X}){\mathcal{K}}^{-1}\mathcal{K}(x,\mathcal{X}{)}^T,$$

$$\mu ({\mathcal{X}}_T)=\Theta \left({\mathcal{X}}_T,\mathcal{X}\right){\Theta }^{-1}\left(I-e^{-\eta \Theta t}\right)\mathcal{Y},$$

$$\Sigma ({\mathcal{X}}_T,{\mathcal{X}}_T)=\mathcal{K}\left({\mathcal{X}}_T,{\mathcal{X}}_T\right)+\Theta ({\mathcal{X}}_T,\mathcal{X}){\Theta }^{-1}\left(I-e^{-\eta \Theta t}\right)\mathcal{K}\left(I-e^{-\eta \Theta t}\right){\Theta }^{-1}\Theta \left(\mathcal{X},{\mathcal{X}}_T\right)$$

$$-\left(\Theta ({\mathcal{X}}_T,\mathcal{X}){\Theta }^{-1}\left(I-e^{-\eta \Theta t}\right)\mathcal{K}\left(\mathcal{X},{\mathcal{X}}_T\right)+h.c.\right).$$

$$\displaystyle\mathcal{N}\big{(}$$

$$\displaystyle\mathcal{K}\left(\mathcal{X}_{T},\mathcal{X}_{T}\right)+\Theta(\mathcal{X}_{T},\mathcal{X})\Theta^{-1}\mathcal{K}\Theta^{-1}\Theta\left(\mathcal{X},\mathcal{X}_{T}\right)-\left(\Theta(\mathcal{X}_{T},\mathcal{X})\Theta^{-1}\mathcal{K}\left(\mathcal{X},\mathcal{X}_{T}\right)+h.c.\right)\big{)}.$$

$$\underset{t\ge 0}{\sup }{\Vert f_t(x)-f_t^{\text{lin}}(x)\Vert }_2,\underset{t\ge 0}{\sup }\frac{{\Vert {\theta }_t-{\theta }_0\Vert }_2}{\sqrt{n}},\underset{t\ge 0}{\sup }{\Vert {\hat{\Theta }}_t-{\hat{\Theta }}_0\Vert }_F=\mathcal{O}(n^{-\frac{1}{2}}),\mathrm{as}n\to \infty .$$

$${\theta }_{i+1}={\theta }_i+\beta ({\theta }_i-{\theta }_{i-1})-\eta {\nabla }_{\theta }\mathcal{L}{|}_{\theta ={\theta }_i}.$$

$$f_{i+1}^{\text{lin}}(x)=f_i^{\text{lin}}(x)-\eta {\hat{\Theta }}_0(x,\mathcal{X}){\nabla }_{f_i^{\text{lin}}(\mathcal{X})}\mathcal{L}+\beta (f_i^{\text{lin}}(x)-f_{i-1}^{\text{lin}}(x))$$

$${\ddot{\omega }}_t$$

$${\ddot{f}}_t{}^{\text{lin}}(x)$$

$$\ell (f,y)=-\sum _i{y}^i\log \sigma (f^i),\sigma (f^i)\equiv \frac{\exp (f^i)}{{\sum }_j\exp (f^j)}.$$

$${\dot{f}}_t^i(x)={\nabla }_{\theta }{f}_t^i(x)\frac{d\theta }{dt}$$

$$=-\eta \sum _{(z,y)\in \mathcal{D}}\sum _j{\nabla }_{\theta }{f}_t^i(x){\nabla }_{\theta }{f}_t^j(z{)}^T\left(\sigma (f_t^j(z))-y^j\right)$$

$$\dot{f_t}(\mathcal{X})$$

$$\dot{f_t}(x)$$

$${\dot{f_t}}^{\text{lin}}(\mathcal{X})$$

$${\dot{f_t}}^{\text{lin}}(x)$$

$$k_n(x,y)=\int {\varphi }^n(x\cdot w){\varphi }^n(y\cdot w)e^{-\Vert w{\Vert }^2/2}dw\cdot (2\pi {)}^{-d/2}=\frac{1}{2\pi }\Vert x{\Vert }^n\Vert y{\Vert }^n{J}_n(\theta )$$

$$\varphi (x)$$

$$J_0(\theta )$$

$$\mathcal{T}(\Sigma )$$

$$\dot{\mathcal{T}}(\Sigma )$$

$$\mathcal{T}(\Sigma )$$

$$\dot{\mathcal{T}}(\Sigma )$$

Full text

1 Introduction

Machine learning models based on deep neural networks have achieved unprecedented performance across a wide range of tasks [1, 2, 3]. Typically, these models are regarded as complex systems for which many types of theoretical analyses are intractable. Moreover, characterizing the gradient-based training dynamics of these models is challenging owing to the typically high-dimensional non-convex loss surfaces governing the optimization. As is common in the physical sciences, investigating the extreme limits of such systems can often shed light on these hard problems. For neural networks, one such limit is that of infinite width, which refers either to the number of hidden units in a fully-connected layer or to the number of channels in a convolutional layer. Under this limit, the output of the network at initialization is a draw from a Gaussian process (GP); moreover, the network output remains governed by a GP after exact Bayesian training using squared loss [4, 5, 6, 7, 8]. Aside from its theoretical simplicity, the infinite-width limit is also of practical interest as wider networks have been found to generalize better [5, 7, 9, 10, 11].

In this work, we explore the learning dynamics of wide neural networks under gradient descent and find that the weight-space description of the dynamics becomes surprisingly simple: as the width becomes large, the neural network can be effectively replaced by its first-order Taylor expansion with respect to its parameters at initialization. For this linear model, the dynamics of gradient descent become analytically tractable. While the linearization is only exact in the infinite width limit, we nevertheless find excellent agreement between the predictions of the original network and those of the linearized version even for finite width configurations. The agreement persists across different architectures, optimization methods, and loss functions.

For squared loss, the exact learning dynamics admit a closed-form solution that allows us to characterize the evolution of the predictive distribution in terms of a GP. This result can be thought of as an extension of “sample-then-optimize" posterior sampling [12] to the training of deep neural networks. Our empirical simulations confirm that the result accurately models the variation in predictions across an ensemble of finite-width models with different random initializations.

Here we summarize our contributions:

• •

Parameter space dynamics: We show that wide network training dynamics in parameter space are equivalent to the training dynamics of a model which is affine in the collection of all network parameters, the weights and biases. This result holds regardless of the choice of loss function. For squared loss, the dynamics admit a closed-form solution as a function of time.

• •

Sufficient conditions for linearization: We formally prove that there exists a threshold learning rate $\eta_{{\rm critical}}$ (see Theorem 2.1), such that gradient descent training trajectories with learning rate smaller than $\eta_{{\rm critical}}$ stay in an $\mathcal{O}\left(n^{-1/2}\right)$-neighborhood of the trajectory of the linearized network when $n$, the width of the hidden layers, is sufficiently large.

• •

Output distribution dynamics: We formally show that the predictions of a neural network throughout gradient descent training are described by a GP as the width goes to infinity (see Theorem 2.2), extending results from Jacot et al. [13]. We further derive explicit time-dependent expressions for the evolution of this GP during training. Finally, we provide a novel interpretation of the result. In particular, it offers a quantitative understanding of the mechanism by which gradient descent differs from Bayesian posterior sampling of the parameters: while both methods generate draws from a GP, gradient descent does not generate samples from the posterior of any probabilistic model.

• •

Large scale experimental support: We empirically investigate the applicability of the theory in the finite-width setting and find that it gives an accurate characterization of both learning dynamics and posterior function distributions across a variety of conditions, including some practical network architectures such as the wide residual network [14].

• •

Parameterization independence: We note that linearization result holds both in standard and NTK parameterization (defined in §2.1), while previous work assumed the latter, emphasizing that the effect is due to increase in width rather than the particular parameterization.

• •

Analytic $\operatorname{ReLU}$ and $\operatorname{erf}$ neural tangent kernels: We compute the analytic neural tangent kernel corresponding to fully-connected networks with $\operatorname{ReLU}$ or $\operatorname{erf}$ nonlinearities.

• •

Source code: Example code investigating both function space and parameter space linearized learning dynamics described in this work is released as open source code within [15].111Note that the open source library has been expanded since initial submission of this work. We also provide accompanying interactive Colab notebooks for both parameter space222colab.sandbox.google.com/github/google/neural-tangents/blob/master/notebooks/weight_space_linearization.ipynb** and function space333colab.sandbox.google.com/github/google/neural-tangents/blob/master/notebooks/function_space_linearization.ipynb linearization.**

1.1 Related work

We build on recent work by Jacot et al. [13] that characterize the exact dynamics of network outputs throughout gradient descent training in the infinite width limit. Their results establish that full batch gradient descent in parameter space corresponds to kernel gradient descent in function space with respect to a new kernel, the Neural Tangent Kernel (NTK). We examine what this implies about dynamics in parameter space, where training updates are actually made.

Daniely et al. [16] study the relationship between neural networks and kernels at initialization. They bound the difference between the infinite width kernel and the empirical kernel at finite width $n$, which diminishes as $\mathcal{O}(1/\sqrt{n})$. Daniely [17] uses the same kernel perspective to study stochastic gradient descent (SGD) training of neural networks.

Saxe et al. [18] study the training dynamics of deep linear networks, in which the nonlinearities are treated as identity functions. Deep linear networks are linear in their inputs, but not in their parameters. In contrast, we show that the outputs of sufficiently wide neural networks are linear in the updates to their parameters during gradient descent, but not usually their inputs.

Du et al. [19], Allen-Zhu et al. [20, 21], Zou et al. [22] study the convergence of gradient descent to global minima. They proved that for i.i.d. Gaussian initialization, the parameters of sufficiently wide networks move little from their initial values during SGD. This small motion of the parameters is crucial to the effect we present, where wide neural networks behave linearly in terms of their parameters throughout training.

Mei et al. [23], Chizat and Bach [24], Rotskoff and Vanden-Eijnden [25], Sirignano and Spiliopoulos [26] analyze the mean field SGD dynamics of training neural networks in the large-width limit. Their mean field analysis describes distributional dynamics of network parameters via a PDE. However, their analysis is restricted to one hidden layer networks with a scaling limit $\left(1/n\right)$ different from ours $\left(1/\sqrt{n}\right)$, which is commonly used in modern networks [2, 27].

Chizat et al. [28]444We note that this is a concurrent work and an expanded version of this note is presented in parallel at NeurIPS 2019. argued that infinite width networks are in ‘lazy training’ regime and maybe too simple to be applicable to realistic neural networks. Nonetheless, we empirically investigate the applicability of the theory in the finite-width setting and find that it gives an accurate characterization of both the learning dynamics and posterior function distributions across a variety of conditions, including some practical network architectures such as the wide residual network [14].

2 Theoretical results

2.1 Notation and setup for architecture and training dynamics

Let $\mathcal{D}\subseteq\mathbb{R}^{n_{0}}\times\mathbb{R}^{k}$ denote the training set and $\mathcal{X}=\left\{x:(x,y)\in\mathcal{D}\right\}$ and $\mathcal{Y}=\left\{y:(x,y)\in\mathcal{D}\right\}$ denote the inputs and labels, respectively. Consider a fully-connected feed-forward network with $L$ hidden layers with widths $n_{l}$, for $l=1,...,L$ and a readout layer with $n_{L+1}=k$. For each $x\in\mathbb{R}^{n_{0}}$, we use $h^{l}(x),x^{l}(x)\in\mathbb{R}^{n_{l}}$ to represent the pre- and post-activation functions at layer $l$ with input $x$. The recurrence relation for a feed-forward network is defined as

[TABLE]

where $\phi$ is a point-wise activation function, $W^{l+1}\in\mathbb{R}^{n_{l}\times n_{l+1}}$ and $b^{l+1}\in\mathbb{R}^{n_{l+1}}$ are the weights and biases, $\omega_{ij}^{l}$ and $b_{j}^{l}$ are the trainable variables, drawn i.i.d. from a standard Gaussian $\omega_{ij}^{l},\beta_{j}^{l}\sim\mathcal{N}(0,1)$ at initialization, and $\sigma_{\omega}^{2}$ and $\sigma_{b}^{2}$ are weight and bias variances. Note that this parametrization is non-standard, and we will refer to it as the NTK parameterization. It has already been adopted in several recent works [29, 30, 13, 19, 31]. Unlike the standard parameterization that only normalizes the forward dynamics of the network, the NTK-parameterization also normalizes its backward dynamics. We note that the predictions and training dynamics of NTK-parameterized networks are identical to those of standard networks, up to a width-dependent scaling factor in the learning rate for each parameter tensor. As we derive, and support experimentally, in Supplementary Material (SM) §F and §G, our results (linearity in weights, GP predictions) also hold for networks with a standard parameterization.

We define $\theta^{l}\equiv\operatorname{vec}\left(\{W^{l},b^{l}\}\right)$, the $\left((n_{l-1}+1)n_{l}\right)\times 1$ vector of all parameters for layer $l$. $\theta=\operatorname{vec}\left(\cup_{l=1}^{L+1}{\theta^{l}}\right)$ then indicates the vector of all network parameters, with similar definitions for $\theta^{\leq l}$ and $\theta^{>l}$. Denote by $\theta_{t}$ the time-dependence of the parameters and by $\theta_{0}$ their initial values. We use $f_{t}(x)\equiv h^{L+1}(x)\in\mathbb{R}^{k}$ to denote the output (or logits) of the neural network at time $t$. Let $\ell(\hat{y},y):\mathbb{R}^{k}\times\mathbb{R}^{k}\to\mathbb{R}$ denote the loss function where the first argument is the prediction and the second argument the true label. In supervised learning, one is interested in learning a $\theta$ that minimizes the empirical loss555To simplify the notation for later equations, we use the total loss here instead of the average loss, but for all plots in §3, we show the average loss., $\mathcal{L}=\sum_{(x,y)\in\mathcal{D}}\ell(f_{t}(x,\theta),y).$

Let $\eta$ be the learning rate666Note that compared to the conventional parameterization, $\eta$ is larger by factor of width [31]. The NTK parameterization allows usage of a universal learning rate scale irrespective of network width.. Via continuous time gradient descent, the evolution of the parameters $\theta$ and the logits $f$ can be written as

[TABLE]

where $f_{t}(\mathcal{X})=\operatorname{vec}\left(\left[f_{t}\left(x\right)\right]_{x\in\mathcal{X}}\right)$, the $k|\mathcal{D}|\times 1$ vector of concatenated logits for all examples, and $\nabla_{f_{t}(\mathcal{X})}\mathcal{L}$ is the gradient of the loss with respect to the model’s output, $f_{t}(\mathcal{X})$. $\hat{\Theta}_{t}\equiv\hat{\Theta}_{t}(\mathcal{X},\mathcal{X})$ is the tangent kernel at time $t$, which is a $k|\mathcal{D}|\times k|\mathcal{D}|$ matrix

[TABLE]

One can define the tangent kernel for general arguments, e.g. $\hat{\Theta}_{t}(x,\mathcal{X})$ where $x$ is test input. At finite-width, $\hat{\Theta}$ will depend on the specific random draw of the parameters and in this context we refer to it as the empirical tangent kernel.

The dynamics of discrete gradient descent can be obtained by replacing $\dot{\theta}_{t}$ and $\dot{f}_{t}(\mathcal{X})$ with $(\theta_{i+1}-\theta_{i})$ and $(f_{i+1}(\mathcal{X})-f_{i}(\mathcal{X}))$ above, and replacing $e^{-\eta\hat{\Theta}_{0}t}$ with $(1-(1-\eta\hat{\Theta}_{0})^{i})$ below.

2.2 Linearized networks have closed form training dynamics for parameters and outputs

In this section, we consider the training dynamics of the linearized network. Specifically, we replace the outputs of the neural network by their first order Taylor expansion,

[TABLE]

where $\omega_{t}\equiv\theta_{t}-\theta_{0}$ is the change in the parameters from their initial values. Note that $f^{\textrm{lin}}_{t}$ is the sum of two terms: the first term is the initial output of the network, which remains unchanged during training, and the second term captures the change to the initial value during training. The dynamics of gradient flow using this linearized function are governed by,

[TABLE]

As ${\nabla_{\theta}f_{0}(x)}$ remains constant throughout training, these dynamics are often quite simple. In the case of an MSE loss, i.e., $\ell(\hat{y},y)=\frac{1}{2}\|\hat{y}-y\|_{2}^{2}$, the ODEs have closed form solutions

[TABLE]

For an arbitrary point $x$, $f^{\textrm{lin}}_{t}(x)=\mu_{t}(x)+\gamma_{t}(x)$, where

[TABLE]

Therefore, we can obtain the time evolution of the linearized neural network without running gradient descent. We only need to compute the tangent kernel $\hat{\Theta}_{0}$ and the outputs $f_{0}$ at initialization and use Equations 8, 10, and 11 to compute the dynamics of the weights and the outputs.

2.3 Infinite width limit yields a Gaussian process

As the width of the hidden layers approaches infinity, the Central Limit Theorem (CLT) implies that the outputs at initialization $\left\{f_{0}(x)\right\}_{x\in\mathcal{X}}$ converge to a multivariate Gaussian in distribution. Informally, this occurs because the pre-activations at each layer are a sum of Gaussian random variables (the weights and bias), and thus become a Gaussian random variable themselves. See [32, 33, 5, 34, 35] for more details, and [36, 7] for a formal treatment.

Therefore, randomly initialized neural networks are in correspondence with a certain class of GPs (hereinafter referred to as NNGPs), which facilitates a fully Bayesian treatment of neural networks [5, 6]. More precisely, let $f_{t}^{i}$ denote the $i$-th output dimension and $\mathcal{K}$ denote the sample-to-sample kernel function (of the pre-activation) of the outputs in the infinite width setting,

[TABLE]

then $f_{0}(\mathcal{X})\sim\mathcal{N}(0,\mathcal{K}(\mathcal{X},\mathcal{X}))$, where $\mathcal{K}^{i,j}(x,x^{\prime})$ denotes the covariance between the $i$-th output of $x$ and $j$-th output of $x^{\prime}$, which can be computed recursively (see Lee et al. [5, §2.3] and SM §E). For a test input $x\in\mathcal{X}_{T}$, the joint output distribution $f\left([x,\mathcal{X}]\right)$ is also multivariate Gaussian. Conditioning on the training samples777 This imposes that $h^{L+1}$ directly corresponds to the network predictions. In the case of softmax readout, variational or sampling methods are required to marginalize over $h^{L+1}$. , $f(\mathcal{X})=\mathcal{Y}$, the distribution of $\left.f(x)\right|\mathcal{X},\mathcal{Y}$ is also a Gaussian $\mathcal{N}\left(\mu(x),\Sigma(x)\right)$,

[TABLE]

and where $\mathcal{K}=\mathcal{K}(\mathcal{X},\mathcal{X})$. This is the posterior predictive distribution resulting from exact Bayesian inference in an infinitely wide neural network.

2.3.1 Gaussian processes from gradient descent training

If we freeze the variables $\theta^{\leq L}$ after initialization and only optimize $\theta^{L+1}$, the original network and its linearization are identical. Letting the width approach infinity, this particular tangent kernel $\hat{\Theta}_{0}$ will converge to $\mathcal{K}$ in probability and Equation 10 will converge to the posterior Equation 13 as $t\to\infty$ (for further details see SM §D). This is a realization of the “sample-then-optimize" approach for evaluating the posterior of a Gaussian process proposed in Matthews et al. [12].

If none of the variables are frozen, in the infinite width setting, $\hat{\Theta}_{0}$ also converges in probability to a deterministic kernel $\Theta$ [13, 37], which we sometimes refer to as the analytic kernel, and which can also be computed recursively (see SM §E). For $\operatorname{ReLU}$ and $\operatorname{erf}$ nonlinearity, $\Theta$ can be exactly computed (SM §C) which we use in §3. Letting the width go to infinity, for any $t$, the output $f^{\textrm{lin}}_{t}(x)$ of the linearized network is also Gaussian distributed because Equations 10 and 11 describe an affine transform of the Gaussian $[f_{0}(x),f_{0}(\mathcal{X})]$. Therefore

Corollary 1**.**

For every test points in $x\in\mathcal{X}_{T}$, and $t\geq 0$, $f^{\textrm{lin}}_{t}(x)$ converges in distribution as width goes to infinity to a Gaussian with mean and covariance given by888Here “+h.c.” is an abbreviation for “plus the Hermitian conjugate”.

[TABLE]

Therefore, over random initialization, $\lim_{t\to\infty}\lim_{n\to\infty}f^{\textrm{lin}}_{t}(x)$ has distribution

[TABLE]

Unlike the case when only $\theta^{L+1}$ is optimized, Equations 14 and 15 do not admit an interpretation corresponding to the posterior sampling of a probabilistic model.999One possible exception is when the NNGP kernel and NTK are the same up to a scalar multiplication. This is the case when the activation function is the identity function and there is no bias term. We contrast the predictive distributions from the NNGP, NTK-GP (i.e. Equations 14 and 15) and ensembles of NNs in Figure 2.

Infinitely-wide neural networks open up ways to study deep neural networks both under fully Bayesian training through the Gaussian process correspondence, and under GD training through the linearization perspective. The resulting distributions over functions are inconsistent (the distribution resulting from GD training does not generally correspond to a Bayesian posterior). We believe understanding the biases over learned functions induced by different training schemes and architectures is a fascinating avenue for future work.

2.4 Infinite width networks are linearized networks

Equation 2 and 3 of the original network are intractable in general, since $\hat{\Theta}_{t}$ evolves with time. However, for the mean squared loss, we are able to prove formally that, as long as the learning rate $\eta<\eta_{{\rm critical}}:=2({\lambda_{\rm{min}}(\Theta)+\lambda_{\rm{max}}(\Theta)})^{-1}$, where ${\lambda_{\textrm{min/max}}}(\Theta)$ is the min/max eigenvalue of $\Theta$, the gradient descent dynamics of the original neural network falls into its linearized dynamics regime.

Theorem 2.1** (Informal).**

Let $n_{1}=\dots=\ n_{L}=n$ and assume $\lambda_{\rm{min}}(\Theta)>0$. Applying gradient descent with learning rate $\eta<\eta_{{\rm critical}}$ (or gradient flow), for every $x\in\mathbb{R}^{n_{0}}$ with $\|x\|_{2}\leq 1$, with probability arbitrarily close to 1 over random initialization,

[TABLE]

Therefore, as $n\to\infty$, the distributions of $f_{t}(x)$ and $f^{\textrm{lin}}_{t}(x)$ become the same. Coupling with Corollary 1, we have

Theorem 2.2**.**

If $\eta<\eta_{\rm critical}$, then for every $x\in\mathbb{R}^{n_{0}}$ with $\|x\|_{2}\leq 1$, as $n\to\infty$, $f_{t}(x)$ converges in distribution to the Gaussian with mean and variance given by Equation 14 and Equation 15.

We refer the readers to Figure  2 for empirical verification of this theorem. The proof of Theorem 2.1 consists of two steps. The first step is to prove the global convergence of overparameterized neural networks [19, 20, 21, 22] and stability of the NTK under gradient descent (and gradient flow); see SM §G. This stability was first observed and proved in [13] in the gradient flow and sequential limit (i.e. letting $n_{1}\to\infty$, …, $n_{L}\to\infty$ sequentially) setting under certain assumptions about global convergence of gradient flow. In §G, we show how to use the NTK to provide a self-contained (and cleaner) proof of such global convergence and the stability of NTK simultaneously. The second step is to couple the stability of NTK with Grönwall’s type arguments [38] to upper bound the discrepancy between $f_{t}$ and $f^{\textrm{lin}}_{t}$, i.e. the first norm in Equation 17. Intuitively, the ODE of the original network (Equation 3) can be considered as a $\|\hat{\Theta}_{t}-\hat{\Theta}_{0}\|_{F}$-fluctuation from the linearized ODE (Equation 7). One expects the difference between the solutions of these two ODEs to be upper bounded by some functional of $\|\hat{\Theta}_{t}-\hat{\Theta}_{0}\|_{F}$; see SM §H. Therefore, for a large width network, the training dynamics can be well approximated by linearized dynamics.

Note that the updates for individual weights in Equation 6 vanish in the infinite width limit, which for instance can be seen from the explicit width dependence of the gradients in the NTK parameterization. Individual weights move by a vanishingly small amount for wide networks in this regime of dynamics, as do hidden layer activations, but they collectively conspire to provide a finite change in the final output of the network, as is necessary for training. An additional insight gained from linearization of the network is that the individual instance dynamics derived in [13] can be viewed as a random features method,101010We thank Alex Alemi for pointing out a subtlety on correspondence to a random features method. where the features are the gradients of the model with respect to its weights.

2.5 Extensions to other optimizers, architectures, and losses

Our theoretical analysis thus far has focused on fully-connected single-output architectures trained by full batch gradient descent. In SM §B we derive corresponding results for: networks with multi-dimensional outputs, training against a cross entropy loss, and gradient descent with momentum.

In addition to these generalizations, there is good reason to suspect the results to extend to much broader class of models and optimization procedures. In particular, a wealth of recent literature suggests that the mean field theory governing the wide network limit of fully-connected models [32, 33] extends naturally to residual networks [35], CNNs [34], RNNs [39], batch normalization [40], and to broad architectures [37]. We postpone the development of these additional theoretical extensions in favor of an empirical investigation of linearization for a variety of architectures.

3 Experiments

In this section, we provide empirical support showing that the training dynamics of wide neural networks are well captured by linearized models. We consider fully-connected, convolutional, and wide ResNet architectures trained with full- and mini- batch gradient descent using learning rates sufficiently small so that the continuous time approximation holds well. We consider two-class classification on CIFAR-10 (horses and planes) as well as ten-class classification on MNIST and CIFAR-10. When using MSE loss, we treat the binary classification task as regression with one class regressing to $+1$ and the other to $-1$.

Experiments in Figures 1, 4, S2, S3, S4, S5 and S6, were done in JAX [41]. The remaining experiments used TensorFlow [42]. An open source implementation of this work providing tools to investigate linearized learning dynamics is available at www.github.com/google/neural-tangents  [15].

Predictive output distribution: In the case of an MSE loss, the output distribution remains Gaussian throughout training. In Figure 2, the predictive output distribution for input points interpolated between two training points is shown for an ensemble of neural networks and their corresponding GPs. The interpolation is given by $x(\alpha)=\alpha x^{(1)}+(1-\alpha)x^{(2)}$ where $x^{(1,2)}$ are two training inputs with different classes. We observe that the mean and variance dynamics of neural network outputs during gradient descent training follow the analytic dynamics from linearization well (Equations 14, 15). Moreover the NNGP predictive distribution which corresponds to exact Bayesian inference, while similar, is noticeably different from the predictive distribution at the end of gradient descent training. For dynamics for individual function draws see SM Figure S1.

Comparison of training dynamics of linearized network to original network: For a particular realization of a finite width network, one can analytically predict the dynamics of the weights and outputs over the course of training using the empirical tangent kernel at initialization. In Figures  3, 4 (see also S2, S3), we compare these linearized dynamics (Equations 8, 9) with the result of training the actual network. In all cases we see remarkably good agreement. We also observe that for finite networks, dynamics predicted using the empirical kernel $\hat{\Theta}$ better match the data than those obtained using the infinite-width, analytic, kernel $\Theta$. To understand this we note that $\|\hat{\Theta}^{(n)}_{T}-\hat{\Theta}^{(n)}_{0}\|_{F}=\mathcal{O}(1/n)\leq\mathcal{O}(1/{\sqrt{n}})=\|\hat{\Theta}^{(n)}_{0}-\Theta\|_{F}$, where $\hat{\Theta}^{(n)}_{0}$ denotes the empirical tangent kernel of width $n$ network, as plotted in Figure 1.

One can directly optimize parameters of $f^{\textrm{lin}}$ instead of solving the ODE induced by the tangent kernel $\hat{\Theta}$. Standard neural network optimization techniques such as mini-batching, weight decay, and data augmentation can be directly applied. In Figure 4 (S2, S3), we compared the training dynamics of the linearized and original network while directly training both networks.

With direct optimization of linearized model, we tested full ($|\mathcal{D}|=50,000$) MNIST digit classification with cross-entropy loss, and trained with a momentum optimizer (Figure S3). For cross-entropy loss with softmax output, some logits at late times grow indefinitely, in contrast to MSE loss where logits converge to target value. The error between original and linearized model for cross entropy loss becomes much worse at late times if the two models deviate significantly before the logits enter their late-time steady-growth regime (See Figure S4).

Linearized dynamics successfully describes the training of networks beyond vanilla fully-connected models. To demonstrate the generality of this procedure we show we can predict the learning dynamics of subclass of Wide Residual Networks (WRNs) [14]. WRNs are a class of model that are popular in computer vision and leverage convolutions, batch normalization, skip connections, and average pooling. In Figure 4, we show a comparison between the linearized dynamics and the true dynamics for a wide residual network trained with MSE loss and SGD with momentum, trained on the full CIFAR-10 dataset. We slightly modified the block structure described in Table S1 so that each layer has a constant number of channels (1024 in this case), and otherwise followed the original implementation. As elsewhere, we see strong agreement between the predicted dynamics and the result of training.

Effects of dataset size: The training dynamics of a neural network match those of its linearization when the width is infinite and the dataset is finite. In previous experiments, we chose sufficiently wide networks to achieve small error between neural networks and their linearization for smaller datasets. Overall, we observe that as the width grows the error decreases (Figure S5). Additionally, we see that the error grows in the size of the dataset. Thus, although error grows with dataset this can be counterbalanced by a corresponding increase in the model size.

4 Discussion

We showed theoretically that the learning dynamics in parameter space of deep nonlinear neural networks are exactly described by a linearized model in the infinite width limit. Empirical investigation revealed that this agrees well with actual training dynamics and predictive distributions across fully-connected, convolutional, and even wide residual network architectures, as well as with different optimizers (gradient descent, momentum, mini-batching) and loss functions (MSE, cross-entropy). Our results suggest that a surprising number of realistic neural networks may be operating in the regime we studied. This is further consistent with recent experimental work showing that neural networks are often robust to re-initialization but not re-randomization of layers (Zhang et al. [43]).

In the regime we study, since the learning dynamics are fully captured by the kernel $\hat{\Theta}$ and the target signal, studying the properties of $\hat{\Theta}$ to determine trainability and generalization are interesting future directions. Furthermore, the infinite width limit gives us a simple characterization of both gradient descent and Bayesian inference. By studying properties of the NNGP kernel $\mathcal{K}$ and the tangent kernel $\Theta$, we may shed light on the inductive bias of gradient descent.

Some layers of modern neural networks may be operating far from the linearized regime. Preliminary observations in Lee et al. [5] showed that wide neural networks trained with SGD perform similarly to the corresponding GPs as width increase, while GPs still outperform trained neural networks for both small and large dataset size. Furthermore, in Novak et al. [7], it is shown that the comparison of performance between finite- and infinite-width networks is highly architecture-dependent. In particular, it was found that infinite-width networks perform as well as or better than their finite-width counterparts for many fully-connected or locally-connected architectures. However, the opposite was found in the case of convolutional networks without pooling. It is still an open research question to determine the main factors that determine these performance gaps. We believe that examining the behavior of infinitely wide networks provides a strong basis from which to build up a systematic understanding of finite-width networks (and/or networks trained with large learning rates).

Acknowledgements

We thank Greg Yang and Alex Alemi for useful discussions and feedback. We are grateful to Daniel Freeman, Alex Irpan and anonymous reviewers for providing valuable feedbacks on the draft. We thank the JAX team for developing a language which makes model linearization and NTK computation straightforward. We would like to especially thank Matthew Johnson for support and debugging help.

Appendix A Additional figures

Appendix B Extensions

B.1 Momentum

One direction is to go beyond vanilla gradient descent dynamics. We consider momentum updates111111Combining the usual two stage update into a single equation.

[TABLE]

The discrete update to the function output becomes

[TABLE]

where $f^{\textrm{lin}}_{t}(x)$ is the output of the linearized network after $t$ steps. One can take the continuous time limit as in Qian [44], Su et al. [45] and obtain

[TABLE]

where continuous time relates to steps $t=i\sqrt{\eta}$ and $\tilde{\beta}=(\beta-1)/\sqrt{\eta}$. These equations are also amenable to analytic treatment for MSE loss. See Figure S2, S3 and 4 for experimental agreement.

B.2 Multi-dimensional output and cross-entropy loss

One can extend the loss function to general functions with multiple output dimensions. Unlike for squared error, we do not have a closed form solution to the dynamics equation. However, the equations for the dynamics can be solved using an ODE solver as an initial value problem.

[TABLE]

Recall that $\frac{\partial\ell}{\partial\hat{y}^{i}}=\sigma(\hat{y}^{i})-y^{i}$. For general input point $x$ and for an arbitrary parameterized function $f^{i}(x)$ parameterized by $\theta$, gradient flow dynamics is given by

[TABLE]

Let $\hat{\Theta}^{ij}(x,\mathcal{X})=\nabla_{\theta}f^{i}(x)\nabla_{\theta}f^{j}(\mathcal{X})^{T}$. The above is

[TABLE]

The linearization is

[TABLE]

For general loss, e.g. cross-entropy with softmax output, we need to rely on solving the ODE Equations S10 and S11. We use the dopri5 method for ODE integration, which is the default integrator in TensorFlow (tf.contrib.integrate.odeint).

Appendix C Neural Tangent kernel for $\operatorname{ReLU}$ and $\operatorname{erf}$

For $\operatorname{ReLU}$ and $\operatorname{erf}$ activation functions, the tangent kernel can be computed analytically. We begin with the case $\phi=$ $\operatorname{ReLU}$; using the formula from Cho and Saul [46], we can compute $\mathcal{T}$ and $\dot{\mathcal{T}}$ in closed form. Let $\Sigma$ be a $2\times 2$ PSD matrix. We will use

[TABLE]

where

[TABLE]

Let $d=2$ and $u=(x\cdot w,y\cdot w)^{T}$. Then $u$ is a mean zero Gaussian with $\Sigma=[[x\cdot x,x\cdot y];[x\cdot y,y\cdot y]]$. Then

[TABLE]

For $\phi=\operatorname{erf}$, let $\Sigma$ be the same as above. Following Williams [47], we get

[TABLE]

Appendix D Gradient flow dynamics for training only the readout-layer

The connection between Gaussian processes and Bayesian wide neural networks can be extended to the setting when only the readout layer parameters are being optimized. More precisely, we show that when training only the readout layer, the outputs of the network form a Gaussian process (over an ensemble of draws from the parameter prior) throughout training, where that output is an interpolation between the GP prior and GP posterior.

Note that for any $x,x^{\prime}\in\mathbb{R}^{n_{0}}$, in the infinite width limit $\bar{x}(x)\cdot\bar{x}(x^{\prime})=\hat{\mathcal{K}}(x,x^{\prime})\to\mathcal{K}(x,x^{\prime})$ in probability, where for notational simplicity we assign $\bar{x}(x)=\left[\frac{\sigma_{w}x^{L}(x)}{\sqrt{n_{L}}},{\sigma_{b}}\right]$. The regression problem is specified with mean-squared loss

[TABLE]

and applying gradient flow to optimize the readout layer (and freezing all other parameters),

[TABLE]

where $\eta$ is the learning rate. The solution to this ODE gives the evolution of the output of an arbitrary $x^{*}$. So long as the empirical kernel $\bar{x}(\mathcal{X})\bar{x}(\mathcal{X})^{T}$ is invertible, it is

[TABLE]

For any $x,x^{\prime}\in\mathbb{R}^{n_{0}}$, letting $n_{l}\to\infty$ for $l=1,\dots,L$, one has the convergence in distribution in probability and distribution respectively

[TABLE]

Moreover $\bar{x}(\mathcal{X})\theta_{0}^{L+1}$ and the term containing $f_{0}(\mathcal{X})$ are the only stochastic term over the ensemble of network initializations, therefore for any $t$ the output $f(x^{*})$ throughout training converges to a Gaussian distribution in the infinite width limit, with

[TABLE]

Thus the output of the neural network is also a GP and the asymptotic solution (i.e. $t\to\infty$) is identical to the posterior of the NNGP (Equation 13). Therefore, in the infinite width case, the optimized neural network is performing posterior sampling if only the readout layer is being trained. This result is a realization of sample-then-optimize equivalence identified in Matthews et al. [12].

Appendix E Computing NTK and NNGP Kernel

For completeness, we reproduce, informally, the recursive formula of the NNGP kernel and the tangent kernel from [5] and [13], respectively. Let the activation function $\phi:\mathbb{R}\to\mathbb{R}$ be absolutely continuous. Let $\mathcal{T}$ and $\dot{\mathcal{T}}$ be functions from $2\times 2$ positive semi-definite matrices $\Sigma$ to $\mathbb{R}$ given by

[TABLE]

In the infinite width limit, the NNGP and tangent kernel can be computed recursively. Let $x,x^{\prime}$ be two inputs in $\mathbb{R}^{n_{0}}$. Then $h^{l}(x)$ and $h^{l}(x^{\prime})$ converge in distribution to a joint Gaussian as $\min\{n_{1},\dots,n_{l-1}\}$. The mean is zero and the variance $\mathcal{K}^{l}(x,x^{\prime})$ is

[TABLE]

[TABLE]

with base case

[TABLE]

Using this one can also derive the tangent kernel for gradient descent training. We will use induction to show that

[TABLE]

where

[TABLE]

with $\tilde{\Theta}^{1}=\tilde{\mathcal{K}}^{1}$. Let

[TABLE]

Then

[TABLE]

Letting $n_{1},\dots,n_{l-1}\to\infty$ sequentially, the first term converges to the NNGP kernel $\mathcal{K}^{l}(x,x^{\prime})$. By applying the chain rule and the induction step (letting $n_{1},\dots,n_{l-2}\to\infty$ sequentially), the second term is

[TABLE]

Appendix F Results in function space for NTK parameterization transfer to standard parameterization

In this Section we present a sketch for why the function space linearization results, derived in [13] for NTK parameterized networks, also apply to networks with a standard parameterization. We follow this up with a formal proof in §G of the convergence of standard parameterization networks to their linearization in the limit of infinite width. A network with standard parameterization is described as:

[TABLE]

The NTK parameterization in Equation 1 is not commonly used for training neural networks. While the function that the network represents is the same for both NTK and standard parameterization, training dynamics under gradient descent are generally different for the two parameterizations. However, for a particular choice of layer-dependent learning rate training dynamics also become identical. Let $\eta^{l}_{\text{NTK},w}$ and $\eta^{l}_{\text{NTK},b}$ be layer-dependent learning rate for $W^{l}$ and $b^{l}$ in the NTK parameterization, and $\eta_{\text{std}}=\frac{1}{n_{\text{max}}}\eta_{0}$ be the learning rate for all parameters in the standard parameterization, where $n_{\text{max}}=\max_{l}n_{l}$. Recall that gradient descent training in standard neural networks requires a learning rate that scales with width like $\frac{1}{n_{\text{max}}}$, so $\eta_{0}$ defines a width-invariant learning rate [31]. If we choose

[TABLE]

then learning dynamics are identical for networks with NTK and standard parameterizations. With only extremely minor modifications, consisting of incorporating the multiplicative factors in Equation S37 into the per-layer contributions to the Jacobian, the arguments in §2.4 go through for an NTK network with learning rates defined in Equation S37. Since an NTK network with these learning rates exhibits identical training dynamics to a standard network with learning rate $\eta_{\text{std}}$, the result in §2.4 that sufficiently wide NTK networks are linear in their parameters throughout training also applies to standard networks.

We can verify this property of networks with the standard parameterization experimentally. In Figure S7, we see that for different choices of dataset, activation function and loss function, final performance of two different parameterization leads to similar quality model for similar value of normalized learning rate $\eta_{\textrm{std}}=\eta_{\textrm{NTK}}/n$. Also, in Figure S8, we observe that our results is not due to the parameterization choice and holds for wide networks using the standard parameterization.

Appendix G Convergence of neural network to its linearization, and stability of NTK under gradient descent

In this section, we show that how to use the NTK to provide a simple proof of the global convergence of a neural network under (full-batch) gradient descent and the stability of NTK under gradient descent. We present the proof for standard parameterization. With some minor changes, the proof can also apply to the NTK parameterization. To lighten the notation, we only consider the asymptotic bound here. The neural networks are parameterized as in Equation S36. We make the following assumptions:

Assumptions [1-4]:

The widths of the hidden layers are identical, i.e. $n_{1}=\dots=n_{L}=n$ (our proof extends naturally to the setting $\frac{n_{l}}{n_{l^{\prime}}}\to\alpha_{l,l^{\prime}}\in(0,\infty)$ as $\min\{n_{1},\dots,n_{L}\}\to\infty$.) 2. 2.

The analytic NTK $\Theta$ (defined in Equation S42) is full-rank, i.e. $0<\lambda_{\rm{min}}:=\lambda_{\rm{min}}(\Theta)\leq\lambda_{\rm{max}}:=\lambda_{\rm{max}}(\Theta)<\infty.$ We set $\eta_{{\rm critical}}=2(\lambda_{\rm{min}}+\lambda_{\rm{max}})^{-1}$ . 3. 3.

The training set $(\mathcal{X},\mathcal{Y})$ is contained in some compact set and $x\neq\tilde{x}$ for all $x,\tilde{x}\in\mathcal{X}$. 4. 4.

The activation function $\phi$ satisfies

[TABLE]

Assumption 2 indeed holds when $\mathcal{X}\subseteq\{x\in\mathbb{R}^{n_{0}}\}:\|x\|_{2}=1\}$ and $\phi(x)$ grows non-polynomially for large $x$ [13]. Throughout this section, we use $C>0$ to denote some constant whose value may depend on $L$, $|\mathcal{X}|$ and $(\sigma_{w}^{2},\sigma_{b}^{2})$ and may change from line to line, but is always independent of $n$.

Let $\theta_{t}$ denote the parameters at time step $t$. We use the following short-hand

[TABLE]

where $|\mathcal{X}|$ is the cardinality of the training set and $k$ is the output dimension of the network. The empirical and analytic NTK of the standard parameterization is defined as

[TABLE]

Note that the convergence of the empirical NTK in probability is proved rigorously in [37]. We consider the MSE loss

[TABLE]

Since $f(\theta_{t})$ converges in distribution to a mean zero Guassian with covariance $\mathcal{K}$, one can show that for arbitrarily small $\delta_{0}>0$, there are constants $R_{0}>0$ and $n_{0}$ (both may depend on $\delta_{0}$, $|\mathcal{X}|$ and $\mathcal{K}$) such that for every $n\geq n_{0}$, with probability at least $(1-\delta_{0})$ over random initialization,

[TABLE]

The gradient descent update with learning rate $\eta$ is

[TABLE]

and the gradient flow equation is

[TABLE]

We prove convergence of neural network training and the stability of NTK for both discrete gradient descent and gradient flow. Both proofs rely on the local lipschitzness of the Jacobian $J(\theta)$.

Lemma 1** **(Local Lipschitzness of the Jacobian).

There is a $K>0$ such that for every $C>0$, with high probability over random initialization (w.h.p.o.r.i.) the following holds

[TABLE]

where

[TABLE]

The following are the main results of this section.

Theorem G.1** **(Gradient descent).

Assume Assumptions [1-4]. For $\delta_{0}>0$ and $\eta_{0}<\eta_{{\rm critical}}$, there exist $R_{0}>0$, $N\in\mathbb{N}$ and $K>1$, such that for every $n\geq N$, the following holds with probability at least $(1-\delta_{0})$ over random initialization when applying gradient descent with learning rate $\eta=\frac{\eta_{0}}{n}$,

[TABLE]

and

[TABLE]

Theorem G.2** **(Gradient Flow).

Assume Assumptions[1-4]. For $\delta_{0}>0$, there exist $R_{0}>0$, $N\in\mathbb{N}$ and $K>1$, such that for every $n\geq N$, the following holds with probability at least $(1-\delta_{0})$ over random initialization when applying gradient flow with “learning rate" $\eta=\frac{\eta_{0}}{n}$

[TABLE]

and

[TABLE]

See the following two subsections for the proof.

Remark 1**.**

One can extend the results in Theorem G.1 and Theorem G.2 to other architectures or functions as long as

The empirical NTK converges in probability and the limit is positive definite. 2. 2.

Lemma 1 holds, i.e. the Jacobian is locally Lipschitz.

G.1 Proof of Theorem G.1

As discussed above, there exist $R_{0}$ and $n_{0}$ such that for every $n\geq n_{0}$, with probability at least $(1-\delta_{0}/10)$ over random initialization,

[TABLE]

Let $C=\frac{3KR_{0}}{\lambda_{\rm{min}}}$ in Lemma 1. We first prove Equation S49 by induction. Choose $n_{1}>n_{0}$ such that for every $n\geq n_{1}$ Equation S47 and Equation S53 hold with probability at least $(1-\delta_{0}/5)$ over random initialization. The $t=0$ case is obvious and we assume Equation S49 holds for $t=t$. Then by induction and the second estimate of Equation S47

[TABLE]

which gives the first estimate of Equation S49 for $t+1$ and which also implies $\|\theta_{j}-\theta_{0}\|_{2}\leq\frac{3KR_{0}}{\lambda_{\rm{min}}}n^{-\frac{1}{2}}$ for $j=0,\dots,t+1$. To prove the second one, we apply the mean value theorem and the formula for gradient decent update at step $t+1$

[TABLE]

where $\tilde{\theta}_{t}$ is some linear interpolation between $\theta_{t}$ and $\theta_{t+1}$. It remains to show with probability at least $(1-\delta_{0}/2)$,

[TABLE]

This can be verified by Lemma 1. Because $\hat{\Theta}_{0}\to\Theta$ [37] in probability, one can find $n_{2}$ such that the event

[TABLE]

has probability at least $(1-\delta_{0}/5)$ for every $n\geq n_{2}$. The assumption $\eta_{0}<\frac{2}{\lambda_{\rm{min}}+\lambda_{\rm{max}}}$ implies

[TABLE]

Thus

[TABLE]

with probability as least $(1-\delta_{0}/2)$ if

[TABLE]

Therefore, we only need to set

[TABLE]

To verify Equation S50, notice that

[TABLE]

where we have applied the second estimate of Equation S49 and Equation S47.

G.2 Proof of Theorem G.2

The first step is the same. There exist $R_{0}$ and $n_{0}$ such that for every $n\geq n_{0}$, with probability at least $(1-\delta_{0}/10)$ over random initialization,

[TABLE]

Let $C=\frac{3KR_{0}}{\lambda_{\rm{min}}}$ in Lemma 1. Using the same arguments as in Section G.1, one can show that there exists $n_{1}$ such that for all $n\geq n_{1}$, with probability at least $(1-\delta_{0}/10)$

[TABLE]

Let

[TABLE]

We claim $t_{1}=\infty$. If not, then for all $t\leq t_{1}$, $\theta_{t}\in B(\theta_{0},Cn^{-\frac{1}{2}})$ and

[TABLE]

Thus

[TABLE]

and

[TABLE]

Note that

[TABLE]

which implies, for all $t\leq t_{1}$

[TABLE]

This contradicts to the definition of $t_{1}$ and thus $t_{1}=\infty$. Note that Equation S78 is the same as the first equation of Equation S51.

G.3 Proof of Lemma 1

The proof relies on upper bounds of operator norms of random Gaussian matrices.

Theorem G.3** (Corollary 5.35 [48]).**

Let $A=A_{N,n}$ be an $N\times n$ random matrix whose entries are independent standard normal random variables. Then for every $t\geq 0$, with probability at least $1-2\exp(-t^{2}/2)$ one has

[TABLE]

For $l\geq 1$, let

[TABLE]

Let $\theta=\{W^{l},b^{l}\}$ and $\tilde{\theta}=\{\tilde{W}^{l},\tilde{b}^{l}\}$ be any two points in $B(\theta_{0},\frac{C}{\sqrt{n}})$. By the above theorem and the triangle inequality, w.h.p. over random initialization,

[TABLE]

Using this and the assumption on $\phi$ Equation S38, it is not difficult to show that there is a constant $K_{1}$, depending on $\sigma_{\omega}^{2},\sigma_{b}^{2},|\mathcal{X}|$ and $L$ such that with high probability over random initialization121212These two estimates can be obtained via induction. To prove bounds relating to $x^{l}$ and $\delta^{l}$, one starts with $l=1$ and $l=L$, respectively.

[TABLE]

Lemma 1 follows from these two estimates. Indeed, with high probability over random initialization

[TABLE]

and similarly

[TABLE]

G.4 Remarks on NTK parameterization

For completeness, we also include analogues of Theorem G.1 and Lemma 1 with NTK parameterization.

Theorem G.4** (NTK parameterization).**

Assume Assumptions [1-4]. For $\delta_{0}>0$ and $\eta_{0}<\eta_{{\rm critical}}$, there exist $R_{0}>0$, $N\in\mathbb{N}$ and $K>1$, such that for every $n\geq N$, the following holds with probability at least $(1-\delta_{0})$ over random initialization when applying gradient descent with learning rate $\eta={\eta_{0}}$,

[TABLE]

and

[TABLE]

Lemma 2** (NTK parameterization: Local Lipschitzness of the Jacobian).**

There is a $K>0$ such that for every $C>0$, with high probability over random initialization the following holds

[TABLE]

Appendix H Bounding the discrepancy between the original and the linearized network: MSE loss

We provide the proof for the gradient flow case. The proof for gradient descent can be obtained similarly. To simplify the notation, let $g^{\textrm{lin}}(t)\equiv f^{\textrm{lin}}_{t}(\mathcal{X})-\mathcal{Y}$ and $g(t)\equiv f_{t}(\mathcal{X})-\mathcal{Y}$. The theorem and proof apply to both standard and NTK parameterization. We use the notation $\lesssim$ to hide the dependence on uninteresting constants.

Theorem H.1**.**

Same as in Theorem G.2. For every $x\in\mathbb{R}^{n_{0}}$ with $\|x\|_{2}\leq 1$, for $\delta_{0}>0$ arbitrarily small, there exist $R_{0}>0$ and $N\in\mathbb{N}$ such that for every $n\geq N$, with probability at least $(1-\delta_{0})$ over random initialization,

[TABLE]

Proof.

[TABLE]

Integrating both sides and using the fact $g^{\textrm{lin}}(0)=g(0)$,

[TABLE]

Let $\lambda_{0}>0$ be the smallest eigenvalue of $\hat{\Theta}_{0}$ (with high probability $\lambda_{0}>\frac{1}{3}\lambda_{\rm{min}}$). Taking the norm gives

[TABLE]

Let

[TABLE]

The above can be written as

[TABLE]

Note that $\alpha(t)$ is non-decreasing. Applying an integral form of the Grönwall’s inequality (see Theorem 1 in [38]) gives

[TABLE]

Note that

[TABLE]

Then

[TABLE]

Let $\sigma_{t}=\sup_{0\leq s\leq t}\|\hat{\Theta}_{s}-\hat{\Theta}_{0}\|_{op}$. Then

[TABLE]

As it is proved in Theorem G.1, for every $\delta_{0}>0$, with probability at least $(1-\delta_{0})$ over random initialization,

[TABLE]

when $n_{1}=\dots=n_{L}=n\to\infty$. Thus for large $n$ and any polynomial $P(t)$ (we use $P(t)=t$ here)

[TABLE]

Therefore

[TABLE]

as $n\to\infty$.

Now we control the discrepancy on a test point $x$. Let $y$ be its true label. Similarly,

[TABLE]

Integrating over $[0,t]$ and taking the norm imply

[TABLE]

Similarly, Lemma 1 implies

[TABLE]

This gives

[TABLE]

Using Equation S118 and Equation S119,

[TABLE]

∎

Appendix I Convergence of empirical kernel

As in Novak et al. [7], we can use Monte Carlo estimates of the tangent kernel (Equation 4) to probe convergence to the infinite width kernel (analytically computed using Equations S26, S29). For simplicity, we consider random inputs drawn from ${\mathcal{N}}(0,1)$ with $n_{0}=1024$. In Figure S9, we observe convergence as both width $n$ increases and the number of Monte Carlo samples $M$ increases. For both NNGP and tangent kernels we observe $\|\hat{\Theta}^{(n)}-\Theta\|_{F}=\mathcal{O}\left(1/\sqrt{n}\right)$ and $\|\hat{\mathcal{K}}^{(n)}-\mathcal{K}\|_{F}=\mathcal{O}\left({1}/\sqrt{n}\right)$, as predicted by a CLT in Daniely et al. [16].

Appendix J Details on Wide Residual Network