Wide Neural Networks of Any Depth Evolve as Linear Models Under Gradient Descent
Jaehoon Lee, Lechao Xiao, Samuel S. Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, Jeffrey Pennington
WWide Neural Networks of Any Depth Evolve asLinear Models Under Gradient Descent
Jaehoon Lee ∗ , Lechao Xiao ∗ , Samuel S. Schoenholz, Yasaman BahriRoman Novak, Jascha Sohl-Dickstein, Jeffrey Pennington Google Brain { jaehlee, xlc, schsam, yasamanb, romann, jaschasd, jpennin}@google.com Abstract
A longstanding goal in deep learning research has been to precisely characterizetraining and generalization. However, the often complex loss landscapes of neuralnetworks have made a theory of learning dynamics elusive. In this work, we showthat for wide neural networks the learning dynamics simplify considerably andthat, in the infinite width limit, they are governed by a linear model obtained fromthe first-order Taylor expansion of the network around its initial parameters. Fur-thermore, mirroring the correspondence between wide Bayesian neural networksand Gaussian processes, gradient-based training of wide neural networks with asquared loss produces test set predictions drawn from a Gaussian process with aparticular compositional kernel. While these theoretical results are only exact in theinfinite width limit, we nevertheless find excellent empirical agreement betweenthe predictions of the original network and those of the linearized version evenfor finite practically-sized networks. This agreement is robust across differentarchitectures, optimization methods, and loss functions.
Machine learning models based on deep neural networks have achieved unprecedented performanceacross a wide range of tasks [1, 2, 3]. Typically, these models are regarded as complex systems forwhich many types of theoretical analyses are intractable. Moreover, characterizing the gradient-basedtraining dynamics of these models is challenging owing to the typically high-dimensional non-convexloss surfaces governing the optimization. As is common in the physical sciences, investigating theextreme 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 ofthe network at initialization is a draw from a Gaussian process (GP); moreover, the network outputremains governed by a GP after exact Bayesian training using squared loss [4, 5, 6, 7, 8]. Aside fromits theoretical simplicity, the infinite-width limit is also of practical interest as wider networks havebeen found to generalize better [5, 7, 9, 10, 11].In this work, we explore the learning dynamics of wide neural networks under gradient descent andfind that the weight-space description of the dynamics becomes surprisingly simple: as the widthbecomes large, the neural network can be effectively replaced by its first-order Taylor expansion withrespect to its parameters at initialization. For this linear model, the dynamics of gradient descentbecome analytically tractable . While the linearization is only exact in the infinite width limit, wenevertheless find excellent agreement between the predictions of the original network and those of ∗ Both authors contributed equally to this work. Work done as a member of the Google AI Residency program(https://g.co/airesidency).33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada. a r X i v : . [ s t a t . M L ] D ec he linearized version even for finite width configurations. The agreement persists across differentarchitectures, optimization methods, and loss functions.For squared loss, the exact learning dynamics admit a closed-form solution that allows us to charac-terize the evolution of the predictive distribution in terms of a GP. This result can be thought of as anextension 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 acrossan 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 spaceare equivalent to the training dynamics of a model which is affine in the collection of all networkparameters, the weights and biases. This result holds regardless of the choice of loss function. Forsquared 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 learningrate η critical (see Theorem 2.1), such that gradient descent training trajectories with learning ratesmaller than η critical stay in an O (cid:0) n − / (cid:1) -neighborhood of the trajectory of the linearized networkwhen n , the width of the hidden layers, is sufficiently large. • Output distribution dynamics : We formally show that the predictions of a neural networkthroughout gradient descent training are described by a GP as the width goes to infinity (seeTheorem 2.2), extending results from Jacot et al. [13]. We further derive explicit time-dependentexpressions for the evolution of this GP during training. Finally, we provide a novel interpretationof the result. In particular, it offers a quantitative understanding of the mechanism by whichgradient descent differs from Bayesian posterior sampling of the parameters: while both methodsgenerate draws from a GP, gradient descent does not generate samples from the posterior of anyprobabilistic model. • Large scale experimental support : We empirically investigate the applicability of the theory inthe finite-width setting and find that it gives an accurate characterization of both learning dynamicsand posterior function distributions across a variety of conditions, including some practical networkarchitectures such as the wide residual network [14]. • Parameterization independence : We note that linearization result holds both in standard andNTK parameterization (defined in §2.1), while previous work assumed the latter, emphasizing thatthe effect is due to increase in width rather than the particular parameterization. • Analytic
ReLU and erf neural tangent kernels : We compute the analytic neural tangent kernelcorresponding to fully-connected networks with
ReLU or erf nonlinearities. • Source code : Example code investigating both function space and parameter space linearizedlearning dynamics described in this work is released as open source code within [15]. We alsoprovide accompanying interactive Colab notebooks for both parameter space and functionspace linearization. We build on recent work by Jacot et al. [13] that characterize the exact dynamics of network outputsthroughout gradient descent training in the infinite width limit. Their results establish that full batchgradient descent in parameter space corresponds to kernel gradient descent in function space withrespect to a new kernel, the Neural Tangent Kernel (NTK). We examine what this implies aboutdynamics in parameter space, where training updates are actually made.Daniely et al. [16] study the relationship between neural networks and kernels at initialization. Theybound the difference between the infinite width kernel and the empirical kernel at finite width n ,which diminishes as O (1 / √ n ) . Daniely [17] uses the same kernel perspective to study stochasticgradient descent (SGD) training of neural networks.Saxe et al. [18] study the training dynamics of deep linear networks, in which the nonlinearitiesare treated as identity functions. Deep linear networks are linear in their inputs, but not in their Note that the open source library has been expanded since initial submission of this work. colab.sandbox.google.com/github/google/neural-tangents/blob/master/notebooks/weight_space_linearization.ipynb colab.sandbox.google.com/github/google/neural-tangents/blob/master/notebooks/function_space_linearization.ipynb (1 /n ) different from ours (1 / √ n ) , which is commonly used in modern networks [2, 27].Chizat et al. [28] argued that infinite width networks are in ‘lazy training’ regime and maybe toosimple to be applicable to realistic neural networks. Nonetheless, we empirically investigate theapplicability of the theory in the finite-width setting and find that it gives an accurate characterizationof 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]. Let
D ⊆ R n × R k denote the training set and X = { x : ( x, y ) ∈ D} and Y = { y : ( x, y ) ∈ D} 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 ∈ R n ,we use h l ( x ) , x l ( x ) ∈ 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 (cid:26) h l +1 = x l W l +1 + b l +1 x l +1 = φ (cid:0) h l +1 (cid:1) and (cid:40) W li,j = σ ω √ n l ω lij b lj = σ b β lj , (1)where φ is a point-wise activation function, W l +1 ∈ R n l × n l +1 and b l +1 ∈ R n l +1 are the weights andbiases, ω lij and b lj are the trainable variables, drawn i.i.d. from a standard Gaussian ω lij , β lj ∼ N (0 , at initialization, and σ ω and σ b 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 severalrecent works [29, 30, 13, 19, 31]. Unlike the standard parameterization that only normalizes theforward 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 identicalto those of standard networks, up to a width-dependent scaling factor in the learning rate for eachparameter tensor. As we derive, and support experimentally, in Supplementary Material (SM) §Fand §G, our results (linearity in weights, GP predictions) also hold for networks with a standardparameterization.We define θ l ≡ vec (cid:0) { W l , b l } (cid:1) , the (( n l − + 1) n l ) × vector of all parameters for layer l . θ =vec (cid:0) ∪ L +1 l =1 θ l (cid:1) then indicates the vector of all network parameters, with similar definitions for θ ≤ l and θ >l . Denote by θ t the time-dependence of the parameters and by θ their initial values. Weuse f t ( x ) ≡ h L +1 ( x ) ∈ R k to denote the output (or logits) of the neural network at time t . Let (cid:96) (ˆ y, y ) : R k × R k → R denote the loss function where the first argument is the prediction and thesecond argument the true label. In supervised learning, one is interested in learning a θ that minimizesthe empirical loss , L = (cid:80) ( x,y ) ∈D (cid:96) ( f t ( x, θ ) , y ) . We note that this is a concurrent work and an expanded version of this note is presented in parallel atNeurIPS 2019. To simplify the notation for later equations, we use the total loss here instead of the average loss, but for allplots in §3, we show the average loss. η be the learning rate . Via continuous time gradient descent, the evolution of the parameters θ and the logits f can be written as ˙ θ t = − η ∇ θ f t ( X ) T ∇ f t ( X ) L (2) ˙ f t ( X ) = ∇ θ f t ( X ) ˙ θ t = − η ˆΘ t ( X , X ) ∇ f t ( X ) L (3)where f t ( X ) = vec (cid:0) [ f t ( x )] x ∈X (cid:1) , the k |D| × vector of concatenated logits for all examples, and ∇ f t ( X ) L is the gradient of the loss with respect to the model’s output, f t ( X ) . ˆΘ t ≡ ˆΘ t ( X , X ) is thetangent kernel at time t , which is a k |D| × k |D| matrix ˆΘ t = ∇ θ f t ( X ) ∇ θ f t ( X ) T = L +1 (cid:88) l =1 ∇ θ l f t ( X ) ∇ θ l f t ( X ) T . (4)One can define the tangent kernel for general arguments, e.g. ˆΘ t ( x, X ) where x is test input. Atfinite-width, ˆΘ will depend on the specific random draw of the parameters and in this context werefer to it as the empirical tangent kernel.The dynamics of discrete gradient descent can be obtained by replacing ˙ θ t and ˙ f t ( X ) with ( θ i +1 − θ i ) and ( f i +1 ( X ) − f i ( X )) above, and replacing e − η ˆΘ t with (1 − (1 − η ˆΘ ) i ) below. In this section, we consider the training dynamics of the linearized network. Specifically, we replacethe outputs of the neural network by their first order Taylor expansion, f lin t ( x ) ≡ f ( x ) + ∇ θ f ( x ) | θ = θ ω t , (5)where ω t ≡ θ t − θ is the change in the parameters from their initial values. Note that f lin t is thesum of two terms: the first term is the initial output of the network, which remains unchanged duringtraining, and the second term captures the change to the initial value during training. The dynamicsof gradient flow using this linearized function are governed by, ˙ ω t = − η ∇ θ f ( X ) T ∇ f lin t ( X ) L (6) ˙ f lin t ( x ) = − η ˆΘ ( x, X ) ∇ f lin t ( X ) L . (7)As ∇ θ f ( x ) remains constant throughout training, these dynamics are often quite simple. In the caseof an MSE loss, i.e., (cid:96) (ˆ y, y ) = (cid:107) ˆ y − y (cid:107) , the ODEs have closed form solutions ω t = −∇ θ f ( X ) T ˆΘ − (cid:16) I − e − η ˆΘ t (cid:17) ( f ( X ) − Y ) , (8) f lin t ( X ) = ( I − e − η ˆΘ t ) Y + e − η ˆΘ t f ( X ) . (9)For an arbitrary point x , f lin t ( x ) = µ t ( x ) + γ t ( x ) , where µ t ( x ) = ˆΘ ( x, X ) ˆΘ − (cid:16) I − e − η ˆΘ t (cid:17) Y (10) γ t ( x ) = f ( x ) − ˆΘ ( x, X ) ˆΘ − (cid:16) I − e − η ˆΘ t (cid:17) f ( X ) . (11)Therefore, we can obtain the time evolution of the linearized neural network without running gradientdescent. We only need to compute the tangent kernel ˆΘ and the outputs f at initialization and useEquations 8, 10, and 11 to compute the dynamics of the weights and the outputs. As the width of the hidden layers approaches infinity, the Central Limit Theorem (CLT) impliesthat the outputs at initialization { f ( x ) } x ∈X converge to a multivariate Gaussian in distribution. Note that compared to the conventional parameterization, η is larger by factor of width [31]. The NTKparameterization allows usage of a universal learning rate scale irrespective of network width. f it denote the i -th output dimension and K denote the sample-to-samplekernel function (of the pre-activation) of the outputs in the infinite width setting, K i,j ( x, x (cid:48) ) = lim min( n ,...,n L ) →∞ E (cid:104) f i ( x ) · f j ( x (cid:48) ) (cid:105) , (12)then f ( X ) ∼ N (0 , K ( X , X )) , where K i,j ( x, x (cid:48) ) denotes the covariance between the i -th outputof x and j -th output of x (cid:48) , which can be computed recursively (see Lee et al. [5, §2.3] and SM§E). For a test input x ∈ X T , the joint output distribution f ([ x, X ]) is also multivariate Gaussian.Conditioning on the training samples , f ( X ) = Y , the distribution of f ( x ) | X , Y is also a Gaussian N ( µ ( x ) , Σ( x )) , µ ( x ) = K ( x, X ) K − Y , Σ( x ) = K ( x, x ) − K ( x, X ) K − K ( x, X ) T , (13)and where K = K ( X , X ) . This is the posterior predictive distribution resulting from exact Bayesianinference in an infinitely wide neural network. If we freeze the variables θ ≤ L after initialization and only optimize θ L +1 , the original network and itslinearization are identical. Letting the width approach infinity, this particular tangent kernel ˆΘ willconverge to K in probability and Equation 10 will converge to the posterior Equation 13 as t → ∞ (for further details see SM §D). This is a realization of the “sample-then-optimize" approach forevaluating the posterior of a Gaussian process proposed in Matthews et al. [12].If none of the variables are frozen, in the infinite width setting, ˆΘ also converges in probability toa deterministic kernel Θ [13, 37], which we sometimes refer to as the analytic kernel, and whichcan also be computed recursively (see SM §E). For ReLU and erf nonlinearity, Θ can be exactlycomputed (SM §C) which we use in §3. Letting the width go to infinity, for any t , the output f lin t ( x ) of the linearized network is also Gaussian distributed because Equations 10 and 11 describe an affinetransform of the Gaussian [ f ( x ) , f ( X )] . Therefore Corollary 1.
For every test points in x ∈ X T , and t ≥ , f lin t ( x ) converges in distribution as widthgoes to infinity to a Gaussian with mean and covariance given by µ ( X T ) = Θ ( X T , X ) Θ − (cid:16) I − e − η Θ t (cid:17) Y , (14) Σ( X T , X T ) = K ( X T , X T ) + Θ( X T , X )Θ − (cid:16) I − e − η Θ t (cid:17) K (cid:16) I − e − η Θ t (cid:17) Θ − Θ ( X , X T ) − (cid:16) Θ( X T , X )Θ − (cid:16) I − e − η Θ t (cid:17) K ( X , X T ) + h.c. (cid:17) . (15) Therefore, over random initialization, lim t →∞ lim n →∞ f lin t ( x ) has distribution N (cid:0) Θ ( X T , X ) Θ − Y , K ( X T , X T ) + Θ( X T , X )Θ − K Θ − Θ ( X , X T ) − (cid:0) Θ( X T , X )Θ − K ( X , X T ) + h.c. (cid:1) (cid:1) . (16)Unlike the case when only θ L +1 is optimized, Equations 14 and 15 do not admit an interpretationcorresponding to the posterior sampling of a probabilistic model. We contrast the predictivedistributions 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 Bayesiantraining through the Gaussian process correspondence, and under GD training through the lineariza-tion perspective. The resulting distributions over functions are inconsistent (the distribution resulting 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 . Here “ +h.c.” is an abbreviation for “plus the Hermitian conjugate”. One possible exception is when the NNGP kernel and NTK are the same up to a scalar multiplication. Thisis the case when the activation function is the identity function and there is no bias term. t -12 -10 -8 -6 -4 -2 k ˆΘ ( n ) t − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F n -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k θ lT − θ l k F / k θ l k F / p n /n layer 1 (input)layer 2layer 3layer 4 (output) n -7 -6 -5 -4 -3 -2 -1 k ˆΘ ( n ) T − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F / p n /n NN n -7 -6 -5 -4 -3 -2 -1 k ˆ K ( n ) T − ˆ K ( n )0 k F / k ˆ K ( n )0 k F / p n /n NN Figure 1:
Relative Frobenius norm change during training.
Three hidden layer
ReLU net-works trained with η = 1 . on a subset of MNIST ( |D| = 128 ). We measure changes of (in-put/output/intermediary) weights, empirical ˆΘ , and empirical ˆ K after T = 2 steps of gradientdescent updates for varying width. We see that the relative change in input/output weights scales as / √ n while intermediate weights scales as /n , this is because the dimension of the input/outputdoes not grow with n . The change in ˆΘ and ˆ K is upper bounded by O (1 / √ n ) but is closer to O (1 /n ) . See Figure S6 for the same experiment with 3-layer tanh and 1-layer ReLU networks. SeeFigures S9 and S10 for additional comparisons of finite width empirical and analytic kernels.from GD training does not generally correspond to a Bayesian posterior). We believe understand-ing the biases over learned functions induced by different training schemes and architectures is afascinating avenue for future work.
Equation 2 and 3 of the original network are intractable in general, since ˆΘ t evolves with time.However, for the mean squared loss, we are able to prove formally that, as long as the learning rate η < η critical := 2( λ min (Θ) + λ max (Θ)) − , where λ min/max (Θ) is the min/max eigenvalue of Θ , thegradient descent dynamics of the original neural network falls into its linearized dynamics regime. Theorem 2.1 (Informal) . Let n = · · · = n L = n and assume λ min (Θ) > . Applying gradientdescent with learning rate η < η critical (or gradient flow), for every x ∈ R n with (cid:107) x (cid:107) ≤ , withprobability arbitrarily close to 1 over random initialization, sup t ≥ (cid:13)(cid:13) f t ( x ) − f lin t ( x ) (cid:13)(cid:13) , sup t ≥ (cid:107) θ t − θ (cid:107) √ n , sup t ≥ (cid:13)(cid:13)(cid:13) ˆΘ t − ˆΘ (cid:13)(cid:13)(cid:13) F = O ( n − ) , as n → ∞ . (17)Therefore, as n → ∞ , the distributions of f t ( x ) and f lin t ( x ) become the same. Coupling withCorollary 1, we have Theorem 2.2. If η < η critical , then for every x ∈ R n with (cid:107) x (cid:107) ≤ , as n → ∞ , f t ( x ) convergesin 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.1consists of two steps. The first step is to prove the global convergence of overparameterized neuralnetworks [19, 20, 21, 22] and stability of the NTK under gradient descent (and gradient flow); seeSM §G. This stability was first observed and proved in [13] in the gradient flow and sequential limit(i.e. letting n → ∞ , . . . , n L → ∞ sequentially) setting under certain assumptions about globalconvergence of gradient flow. In §G, we show how to use the NTK to provide a self-contained (andcleaner) proof of such global convergence and the stability of NTK simultaneously. The second stepis to couple the stability of NTK with Grönwall’s type arguments [38] to upper bound the discrepancybetween f t and f lin t , i.e. the first norm in Equation 17. Intuitively, the ODE of the original network(Equation 3) can be considered as a (cid:107) ˆΘ t − ˆΘ (cid:107) F -fluctuation from the linearized ODE (Equation 7).One expects the difference between the solutions of these two ODEs to be upper bounded by somefunctional of (cid:107) ˆΘ t − ˆΘ (cid:107) F ; see SM §H. Therefore, for a large width network, the training dynamicscan be well approximated by linearized dynamics.Note that the updates for individual weights in Equation 6 vanish in the infinite width limit, which forinstance 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 finaloutput of the network, as is necessary for training. An additional insight gained from linearization6f the network is that the individual instance dynamics derived in [13] can be viewed as a randomfeatures method, where the features are the gradients of the model with respect to its weights. Our theoretical analysis thus far has focused on fully-connected single-output architectures trainedby 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 muchbroader class of models and optimization procedures. In particular, a wealth of recent literaturesuggests 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], andto broad architectures [37]. We postpone the development of these additional theoretical extensionsin favor of an empirical investigation of linearization for a variety of architectures. α O u t p u t V a l u e t=0 α t=16 α t=256 α t=65536 NNsNNGPNTK
Figure 2:
Dynamics of mean and variance of trained neural network outputs follow analyticdynamics from linearization . Black lines indicate the time evolution of the predictive outputdistribution from an ensemble of 100 trained neural networks (NNs). The blue region indicates theanalytic prediction of the output distribution throughout training (Equations 14, 15). Finally, the redregion indicates the prediction that would result from training only the top layer, corresponding to anNNGP (Equations S22, S23). The trained network has 3 hidden layers of width 8192, tanh activationfunctions, σ w = 1 . , no bias, and η = 0 . . The output is computed for inputs interpolated betweentwo training points (denoted with black dots) x ( α ) = αx (1) + (1 − α ) x (2) . The shaded region anddotted lines denote 2 standard deviations ( ∼ quantile) from the mean denoted in solid lines.Training was performed with full-batch gradient descent with dataset size |D| = 128 . For dynamicsfor individual function initializations, see SM Figure S1. In this section, we provide empirical support showing that the training dynamics of wide neuralnetworks are well captured by linearized models. We consider fully-connected, convolutional, andwide ResNet architectures trained with full- and mini- batch gradient descent using learning ratessufficiently small so that the continuous time approximation holds well. We consider two-classclassification on CIFAR-10 (horses and planes) as well as ten-class classification on MNIST andCIFAR-10. When using MSE loss, we treat the binary classification task as regression with one classregressing to +1 and the other to − .Experiments in Figures 1, 4, S2, S3, S4, S5 and S6, were done in JAX [41]. The remaining experi-ments used TensorFlow [42]. An open source implementation of this work providing tools to inves-tigate linearized learning dynamics is available at [15]. Predictive output distribution : In the case of an MSE loss, the output distribution remains Gaussianthroughout training. In Figure 2, the predictive output distribution for input points interpolatedbetween two training points is shown for an ensemble of neural networks and their correspondingGPs. The interpolation is given by x ( α ) = αx (1) + (1 − α ) x (2) where x (1 , are two training inputs We thank Alex Alemi for pointing out a subtlety on correspondence to a random features method. -2 -1 t Test output value
Neural NetworkLinearized Model 10 -2 -1 t -0.1-0.08-0.06-0.04-0.020.00.020.040.06 Weight change ( ω ) -2 -1 t Loss
TrainTrain f lin TestTest f lin -2 -1 t -3 -2 -1 q › ( f ( x ) − f lin ( x )) fi ΘˆΘ
Figure 3:
Full batch gradient descent on a model behaves similarly to analytic dynamics onits linearization, both for network outputs, and also for individual weights.
A binary CIFARclassification task with MSE loss and a
ReLU fully-connected network with 5 hidden layers of width n = 2048 , η = 0 . , |D| = 256 , k = 1 , σ w = 2 . , and σ b = 0 . . Left two panes show dynamics fora randomly selected subset of datapoints or parameters. Third pane shows that the dynamics of lossfor training and test points agree well between the original and linearized model. The last pane showsthe dynamics of RMSE between the two models on test points. We observe that the empirical kernel ˆΘ gives more accurate dynamics for finite width networks.with different classes. We observe that the mean and variance dynamics of neural network outputsduring gradient descent training follow the analytic dynamics from linearization well (Equations14, 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 descenttraining. For dynamics for individual function draws see SM Figure S1. Comparison of training dynamics of linearized network to original network : For a particularrealization of a finite width network, one can analytically predict the dynamics of the weights andoutputs over the course of training using the empirical tangent kernel at initialization. In Figures3, 4 (see also S2, S3), we compare these linearized dynamics (Equations 8, 9) with the result oftraining the actual network. In all cases we see remarkably good agreement. We also observethat for finite networks, dynamics predicted using the empirical kernel ˆΘ better match the datathan those obtained using the infinite-width, analytic, kernel Θ . To understand this we note that (cid:107) ˆΘ ( n ) T − ˆΘ ( n )0 (cid:107) F = O (1 /n ) ≤ O (1 / √ n ) = (cid:107) ˆΘ ( n )0 − Θ (cid:107) F , where ˆΘ ( n )0 denotes the empirical tangentkernel of width n network, as plotted in Figure 1.One can directly optimize parameters of f lin instead of solving the ODE induced by the tangentkernel ˆΘ . Standard neural network optimization techniques such as mini-batching, weight decay, anddata augmentation can be directly applied. In Figure 4 (S2, S3), we compared the training dynamicsof the linearized and original network while directly training both networks.With direct optimization of linearized model, we tested full ( |D| = 50 , ) MNIST digit classifica-tion with cross-entropy loss, and trained with a momentum optimizer (Figure S3). For cross-entropyloss with softmax output, some logits at late times grow indefinitely, in contrast to MSE loss wherelogits converge to target value. The error between original and linearized model for cross entropyloss becomes much worse at late times if the two models deviate significantly before the logits entertheir late-time steady-growth regime (See Figure S4).Linearized dynamics successfully describes the training of networks beyond vanilla fully-connectedmodels. To demonstrate the generality of this procedure we show we can predict the learningdynamics of subclass of Wide Residual Networks (WRNs) [14]. WRNs are a class of model that arepopular in computer vision and leverage convolutions, batch normalization, skip connections, andaverage pooling. In Figure 4, we show a comparison between the linearized dynamics and the truedynamics for a wide residual network trained with MSE loss and SGD with momentum, trained onthe full CIFAR-10 dataset . We slightly modified the block structure described in Table S1 so thateach layer has a constant number of channels (1024 in this case), and otherwise followed the originalimplementation. As elsewhere, we see strong agreement between the predicted dynamics and theresult of training. Effects of dataset size : The training dynamics of a neural network match those of its linearizationwhen the width is infinite and the dataset is finite. In previous experiments, we chose sufficientlywide networks to achieve small error between neural networks and their linearization for smaller8 t Neural NetworkLinearized Model t t TrainTrain f lin TestTest f lin t Figure 4:
A wide residual network and its linearization behave similarly when both are trainedby SGD with momentum on MSE loss on CIFAR-10.
We adopt the network architecturefrom Zagoruyko and Komodakis [14]. We use N = 1 , channel size , η = 1 . , β = 0 . , k = 10 , σ w = 1 . , and σ b = 0 . . See Table S1 for details of the architecture. Both the linearizedand original model are trained directly on full CIFAR-10 ( |D| = 50 , ), using SGD with batch size8. Output dynamics for a randomly selected subset of train and test points are shown in the first twopanes. Last two panes show training and accuracy curves for the original and linearized networks.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 canbe counterbalanced by a corresponding increase in the model size. We showed theoretically that the learning dynamics in parameter space of deep nonlinear neuralnetworks are exactly described by a linearized model in the infinite width limit. Empirical investiga-tion revealed that this agrees well with actual training dynamics and predictive distributions acrossfully-connected, convolutional, and even wide residual network architectures, as well as with differentoptimizers (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 theregime we studied. This is further consistent with recent experimental work showing that neuralnetworks 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 ˆΘ and the targetsignal, studying the properties of ˆΘ to determine trainability and generalization are interesting futuredirections. Furthermore, the infinite width limit gives us a simple characterization of both gradientdescent and Bayesian inference. By studying properties of the NNGP kernel K and the tangent kernel Θ , 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. Preliminaryobservations in Lee et al. [5] showed that wide neural networks trained with SGD perform similarlyto the corresponding GPs as width increase, while GPs still outperform trained neural networks forboth small and large dataset size. Furthermore, in Novak et al. [7], it is shown that the comparisonof performance between finite- and infinite-width networks is highly architecture-dependent. Inparticular, it was found that infinite-width networks perform as well as or better than their finite-widthcounterparts for many fully-connected or locally-connected architectures. However, the opposite wasfound in the case of convolutional networks without pooling. It is still an open research question todetermine the main factors that determine these performance gaps. We believe that examining thebehavior of infinitely wide networks provides a strong basis from which to build up a systematicunderstanding 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 toDaniel Freeman, Alex Irpan and anonymous reviewers for providing valuable feedbacks on thedraft. We thank the JAX team for developing a language which makes model linearization and NTKcomputation straightforward. We would like to especially thank Matthew Johnson for support anddebugging help. 9 eferences [1] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deepconvolutional neural networks. In
Advances in Neural Information Processing Systems . 2012.[2] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for imagerecognition. In
Conference on Computer Vision and Pattern Recognition , pages 770–778, 2016.[3] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. Bert: Pre-training ofdeep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 ,2018.[4] Radford M. Neal. Priors for infinite networks (tech. rep. no. crg-tr-94-1).
University of Toronto ,1994.[5] Jaehoon Lee, Yasaman Bahri, Roman Novak, Sam Schoenholz, Jeffrey Pennington, and JaschaSohl-dickstein. Deep neural networks as gaussian processes. In
International Conference onLearning Representations , 2018.[6] Alexander G. de G. Matthews, Jiri Hron, Mark Rowland, Richard E. Turner, and ZoubinGhahramani. Gaussian process behaviour in wide deep neural networks. In
InternationalConference on Learning Representations , 4 2018. URL https://openreview.net/forum?id=H1-nGgWC- .[7] Roman Novak, Lechao Xiao, Jaehoon Lee, Yasaman Bahri, Greg Yang, Jiri Hron, Daniel A.Abolafia, Jeffrey Pennington, and Jascha Sohl-Dickstein. Bayesian deep convolutional net-works with many channels are gaussian processes. In
International Conference on LearningRepresentations , 2019.[8] Adrià Garriga-Alonso, Laurence Aitchison, and Carl Edward Rasmussen. Deep convolutionalnetworks as shallow gaussian processes. In
International Conference on Learning Representa-tions , 2019.[9] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro. In search of the real inductive bias:On the role of implicit regularization in deep learning. In
International Conference on LearningRepresentations workshop track , 2015.[10] Roman Novak, Yasaman Bahri, Daniel A. Abolafia, Jeffrey Pennington, and Jascha Sohl-Dickstein. Sensitivity and generalization in neural networks: an empirical study. In
InternationalConference on Learning Representations , 2018.[11] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro. Therole of over-parametrization in generalization of neural networks. In
International Conferenceon Learning Representations , 2019.[12] Alexander G. de G. Matthews, Jiri Hron, Richard E. Turner, and Zoubin Ghahramani. Sample-then-optimize posterior sampling for bayesian linear models. In
NeurIPS Workshop on Advancesin Approximate Bayesian Inference , 2017. URL http://approximateinference.org/2017/accepted/MatthewsEtAl2017.pdf .[13] Arthur Jacot, Franck Gabriel, and Clement Hongler. Neural tangent kernel: Convergence andgeneralization in neural networks. In
Advances in Neural Information Processing Systems ,2018.[14] Sergey Zagoruyko and Nikos Komodakis. Wide residual networks. In
British Machine VisionConference , 2016.[15] Roman Novak, Lechao Xiao, Jiri Hron, Jaehoon Lee, Alexander A. Alemi, Jascha Sohl-Dickstein, and Samuel S. Schoenholz. Neural tangents: Fast and easy infinite neural networks inpython. https://github.com/google/neural-tangents , https://arxiv.org/abs/1912.02803 , 2019.[16] Amit Daniely, Roy Frostig, and Yoram Singer. Toward deeper understanding of neural networks:The power of initialization and a dual view on expressivity. In Advances In Neural InformationProcessing Systems , 2016. 1017] Amit Daniely. SGD learns the conjugate kernel class of the network. In
Advances in NeuralInformation Processing Systems , 2017.[18] Andrew M Saxe, James L McClelland, and Surya Ganguli. Exact solutions to the nonlineardynamics of learning in deep linear neural networks. In
International Conference on LearningRepresentations , 2014.[19] Simon S Du, Jason D Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent findsglobal minima of deep neural networks. In
International Conference on Machine Learning ,2019.[20] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning viaover-parameterization. In
International Conference on Machine Learning , 2019.[21] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. On the convergence rate of training recurrentneural networks. arXiv preprint arXiv:1810.12065 , 2018.[22] Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Stochastic gradient descent optimizesover-parameterized deep relu networks.
Machine Learning , 2019.[23] Song Mei, Andrea Montanari, and Phan-Minh Nguyen. A mean field view of the landscapeof two-layer neural networks.
Proceedings of the National Academy of Sciences , 115(33):E7665–E7671, 2018.[24] Lenaic Chizat and Francis Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. In
Advances in neural information processingsystems , 2018.[25] Grant M Rotskoff and Eric Vanden-Eijnden. Parameters as interacting particles: long timeconvergence and asymptotic error scaling of neural networks. In
Advances in neural informationprocessing systems , 2018.[26] Justin Sirignano and Konstantinos Spiliopoulos. Mean field analysis of neural networks. arXivpreprint arXiv:1805.01053 , 2018.[27] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforwardneural networks. In
International Conference on Artificial Intelligence and Statistics , pages249–256, 2010.[28] Lenaic Chizat, Edouard Oyallon, and Francis Bach. On lazy training in differentiable program-ming. arXiv preprint arXiv:1812.07956 , 2018.[29] Twan van Laarhoven. L2 regularization versus batch and weight normalization. arXiv preprintarXiv:1706.05350 , 2017.[30] Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen. Progressive growing of GANs forimproved quality, stability, and variation. In
International Conference on Learning Representa-tions , 2018.[31] Daniel S. Park, Jascha Sohl-Dickstein, Quoc V. Le, and Samuel L. Smith. The effect of networkwidth on stochastic gradient descent and generalization: an empirical study. In
InternationalConference on Machine Learning , 2019.[32] Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli.Exponential expressivity in deep neural networks through transient chaos. In
Advances InNeural Information Processing Systems , pages 3360–3368, 2016.[33] Samuel S Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein. Deep informa-tion propagation.
International Conference on Learning Representations , 2017.[34] Lechao Xiao, Yasaman Bahri, Jascha Sohl-Dickstein, Samuel Schoenholz, and Jeffrey Penning-ton. Dynamical isometry and a mean field theory of CNNs: How to train 10,000-layer vanillaconvolutional neural networks. In
International Conference on Machine Learning , 2018.1135] Ge Yang and Samuel Schoenholz. Mean field residual networks: On the edge of chaos. In
Advances in Neural Information Processing Systems . 2017.[36] Alexander G de G Matthews, Mark Rowland, Jiri Hron, Richard E Turner, and ZoubinGhahramani. Gaussian process behaviour in wide deep neural networks. arXiv preprintarXiv:1804.11271 , 9 2018.[37] Greg Yang. Scaling limits of wide neural networks with weight sharing: Gaussian pro-cess behavior, gradient independence, and neural tangent kernel derivation. arXiv preprintarXiv:1902.04760 , 2019.[38] Sever Silvestru Dragomir.
Some Gronwall type inequalities and applications . Nova SciencePublishers New York, 2003.[39] Minmin Chen, Jeffrey Pennington, and Samuel Schoenholz. Dynamical isometry and a meanfield theory of RNNs: Gating enables signal propagation in recurrent neural networks. In
International Conference on Machine Learning , 2018.[40] Greg Yang, Jeffrey Pennington, Vinay Rao, Jascha Sohl-Dickstein, and Samuel S. Schoen-holz. A mean field theory of batch normalization. In
International Conference on LearningRepresentations , 2019.[41] Roy Frostig, Peter Hawkins, Matthew Johnson, Chris Leary, and Dougal Maclaurin. JAX:Autograd and XLA. , 2018.[42] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, MatthieuDevin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. Tensorflow: A system forlarge-scale machine learning. In , 2016.[43] Chiyuan Zhang, Samy Bengio, and Yoram Singer. Are all layers created equal? arXiv preprintarXiv:1902.01996 , 2019.[44] Ning Qian. On the momentum term in gradient descent learning algorithms.
Neural networks ,12(1):145–151, 1999.[45] Weijie Su, Stephen Boyd, and Emmanuel Candes. A differential equation for modeling nes-terov’s accelerated gradient method: Theory and insights. In
Advances in Neural InformationProcessing Systems , pages 2510–2518, 2014.[46] Youngmin Cho and Lawrence K Saul. Kernel methods for deep learning. In
Advances in neuralinformation processing systems , 2009.[47] Christopher KI Williams. Computing with infinite networks. In
Advances in neural informationprocessing systems , pages 295–301, 1997.[48] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXivpreprint arXiv:1011.3027 , 2010. 12 upplementary MaterialA Additional figures α O u t p u t V a l u e t=0 α t=16 α t=256 α t=65536 NNsNNGPNTK
Figure S1:
Sample of neural network outputs.
The lines correspond to the functions learned for100 different initializations. The configuration is the same as in Figure 2. O u t p u t V a l u e Train output
Neural NetworkLinearized Model 1.00.50.00.51.01.52.0
Test output -0.008-0.006-0.004-0.0020.00.0020.004
Weight change( ω ) t Loss
TrainTrain f lin TestTest f lin t Accuracy t -8 -7 -6 -5 -4 -3 -2 -1 q › ( f ( x ) − f lin ( x )) fi Figure S2:
A convolutional network and its linearization behave similarly when trained usingfull batch gradient descent with a momentum optimizer.
Binary CIFAR classification task withMSE loss, tanh convolutional network with 3 hidden layers of channel size n = 512 , × sizefilters, average pooling after last convolutional layer, η = 0 . , β = 0 . , |D| = 128 , σ w = 2 . and σ b = 0 . . The linearized model is trained directly by full batch gradient descent with momentum,rather than by integrating its continuous time analytic dynamics. Panes are the same as in Figure 3.1 O u t p u t V a l u e Train output
Neural NetworkLinearized Model 105051015
Test output -0.06-0.04-0.020.00.020.040.06
Weight change( ω ) t Loss
TrainTrain f lin TestTest f lin t Accuracy t -5 -4 -3 -2 -1 q › ( f ( x ) − f lin ( x )) fi Figure S3:
A neural network and its linearization behave similarly when both are trainedvia SGD with momentum on cross entropy loss on MNIST.
Experiment is for 10 class MNISTclassification using a
ReLU fully connected network with 2 hidden layers of width n = 2048 , η = 1 . , β = 0 . , |D| = 50 , , k = 10 , σ w = 2 . , and σ b = 0 . . Both models are trained usingstochastic minibatching with batch size 64. Panes are the same as in Figure 3, except that the top rowshows all ten logits for a single randomly selected datapoint. -1 t n = 32 Neural NetworkLinearized Model -1 t n = 64 -1 t n = 128 -1 t n = 256 -1 t n = 512 -1 t n = 1024 -1 t n = 2048 -1 t n = 4096 Figure S4:
Logit deviation for cross entropy loss.
Logits for models trained with cross entropy lossdiverge at late times. If the deviation between the logits of the linearized model and original modelare large early in training, as shown for the narrower networks (first row), logit deviation at latetimes can be significantly large. As the network becomes wider (second row), the logit deviates at alater point in training. Fully connected tanh network L = 4 trained on binary CIFAR classificationproblem. 2 Width2 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 FC L =1 . L =2 . L =4 . L =8 . L =16 . Width2 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 FC, Dataset Size | D | =32 . | D | =64 . | D | =128 . | D | =256 . | D | =512 . | D | =1024 . | D | =2048 . | D | =4096 . Channels2 -3 -2 -1 CNN L =4 . L =8 . L =16 . L =32 . Channels2 -2 -1 Wide ResNet L =10 . L =16 . Figure S5:
Error dependence on depth, width, and dataset size.
Final value of the RMSE forfully-connected, convolutional, wide residual network as networks become wider for varying depthand dataset size. Error in fully connected networks as the depth is varied from 1 to 16 (first) and thedataset size is varied from 32 to 4096 (last). Error in convolutional networks as the depth is variedbetween 1 and 32 (second), and WRN for depths 10 and 16 corresponding to N=1,2 described inTable S1 (third). Networks are critically initialized σ w = 2 . , σ b = 0 . , trained with gradient descenton MSE loss. Experiments in the first three panes used |D| = 128 . t -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k ˆΘ ( n ) t − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F n -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k θ lT − θ l k F / k θ l k F / p n / p n layer 1 (input) layer 2 (output) n -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k ˆΘ ( n ) T − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F / p n /n NN n -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k ˆ K ( n ) T − ˆ K ( n )0 k F / k ˆ K ( n )0 k F / p n /n NN t -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 k ˆΘ ( n ) t − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F n -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 k θ lT − θ l k F / k θ l k F / p n /n layer 1 (input)layer 2layer 3layer 4 (output) n -8 -7 -6 -5 -4 -3 -2 -1 k ˆΘ ( n ) T − ˆΘ ( n )0 k F / k ˆΘ ( n )0 k F / p n /n NN n -7 -6 -5 -4 -3 -2 -1 k ˆ K ( n ) T − ˆ K ( n )0 k F / k ˆ K ( n )0 k F / p n /n NN Figure S6:
Relative Frobenius norm change during training. (top)
One hidden layer,
ReLU networks trained with η = 1 . , on a 2-class CIFAR10 subset of size |D| = 128 . We measure changesof (read-out/non read-out) weights, empirical ˆΘ and empirical ˆ K after T = 2 steps of gradientdescent updates for varying width. (bottom) Networks with three layer tanh nonlinearity and otherdetails are identical to Figure 1.
B Extensions
B.1 Momentum
One direction is to go beyond vanilla gradient descent dynamics. We consider momentum updates θ i +1 = θ i + β ( θ i − θ i − ) − η ∇ θ L| θ = θ i . (S1)The discrete update to the function output becomes f lin i +1 ( x ) = f lin i ( x ) − η ˆΘ ( x, X ) ∇ f lin i ( X ) L + β ( f lin i ( x ) − f lin i − ( x )) (S2) Combining the usual two stage update into a single equation. f lin t ( x ) is the output of the linearized network after t steps. One can take the continuous timelimit as in Qian [44], Su et al. [45] and obtain ¨ ω t = ˜ β ˙ ω t − ∇ θ f lin ( X ) T ∇ f lin t ( X ) L (S3) ¨ f t lin ( x ) = ˜ β ˙ f lin t ( x ) − ˆΘ ( x, X ) ∇ f lin t ( X ) L (S4)where continuous time relates to steps t = i √ η and ˜ β = ( β − / √ η . These equations are alsoamenable 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. Unlikefor squared error, we do not have a closed form solution to the dynamics equation. However, theequations for the dynamics can be solved using an ODE solver as an initial value problem. (cid:96) ( f, y ) = − (cid:88) i y i log σ ( f i ) , σ ( f i ) ≡ exp( f i ) (cid:80) j exp( f j ) . (S5)Recall that ∂(cid:96)∂ ˆ y i = σ (ˆ y i ) − y i . For general input point x and for an arbitrary parameterized function f i ( x ) parameterized by θ , gradient flow dynamics is given by ˙ f it ( x ) = ∇ θ f it ( x ) dθdt = − η ∇ θ f it ( x ) (cid:88) j (cid:88) ( z,y ) ∈D (cid:20) ∇ θ f jt ( z ) T ∂(cid:96) ( f t , y ) ∂ ˆ y j (cid:21) (S6) = − η (cid:88) ( z,y ) ∈D (cid:88) j ∇ θ f it ( x ) ∇ θ f jt ( z ) T (cid:16) σ ( f jt ( z )) − y j (cid:17) (S7)Let ˆΘ ij ( x, X ) = ∇ θ f i ( x ) ∇ θ f j ( X ) T . The above is ˙ f t ( X ) = − η ˆΘ t ( X , X ) ( σ ( f t ( X )) − Y ) (S8) ˙ f t ( x ) = − η ˆΘ t ( x, X ) ( σ ( f t ( X )) − Y ) . (S9)The linearization is ˙ f t lin ( X ) = − η ˆΘ ( X , X ) (cid:0) σ ( f lin t ( X )) − Y (cid:1) (S10) ˙ f t lin ( x ) = − η ˆΘ ( x, X ) (cid:0) σ ( f lin t ( X )) − Y (cid:1) . (S11)For general loss, e.g. cross-entropy with softmax output, we need to rely on solving the ODEEquations S10 and S11. We use the dopri5 method for ODE integration, which is the defaultintegrator in TensorFlow ( tf.contrib.integrate.odeint ). C Neural Tangent kernel for
ReLU and erf
For
ReLU and erf activation functions, the tangent kernel can be computed analytically. We beginwith the case φ = ReLU ; using the formula from Cho and Saul [46], we can compute T and ˙ T inclosed form. Let Σ be a × PSD matrix. We will use k n ( x, y ) = (cid:90) φ n ( x · w ) φ n ( y · w ) e −(cid:107) w (cid:107) / dw · (2 π ) − d/ = 12 π (cid:107) x (cid:107) n (cid:107) y (cid:107) n J n ( θ ) (S12)where φ ( x ) = max( x, , θ ( x, y ) = arccos (cid:18) x · y (cid:107) x (cid:107)(cid:107) y (cid:107) (cid:19) ,J ( θ ) = π − θ , J ( θ ) = sin θ + ( π − θ ) cos θ = (cid:115) − (cid:18) x · y (cid:107) x (cid:107)(cid:107) y (cid:107) (cid:19) + ( π − θ ) (cid:18) x · y (cid:107) x (cid:107)(cid:107) y (cid:107) (cid:19) . (S13)4et d = 2 and u = ( x · w, y · w ) T . Then u is a mean zero Gaussian with Σ = [[ x · x, x · y ]; [ x · y, y · y ]] .Then T (Σ) = k ( x, y ) = 12 π (cid:107) x (cid:107)(cid:107) y (cid:107) J ( θ ) (S14) ˙ T (Σ) = k ( x, y ) = 12 π J ( θ ) (S15)For φ = erf , let Σ be the same as above. Following Williams [47], we get T (Σ) = 2 π sin − (cid:32) x · y (cid:112) (1 + 2 x · x )(1 + 2 y · y ) (cid:33) (S16) ˙ T (Σ) = 4 π det( I + 2Σ) − / (S17) D Gradient flow dynamics for training only the readout-layer
The connection between Gaussian processes and Bayesian wide neural networks can be extended tothe setting when only the readout layer parameters are being optimized. More precisely, we show thatwhen training only the readout layer, the outputs of the network form a Gaussian process (over anensemble of draws from the parameter prior) throughout training, where that output is an interpolationbetween the GP prior and GP posterior.Note that for any x, x (cid:48) ∈ R n , in the infinite width limit ¯ x ( x ) · ¯ x ( x (cid:48) ) = ˆ K ( x, x (cid:48) ) → K ( x, x (cid:48) ) inprobability, where for notational simplicity we assign ¯ x ( x ) = (cid:104) σ w x L ( x ) √ n L , σ b (cid:105) . The regression problemis specified with mean-squared loss L = 12 (cid:107) f ( X ) − Y(cid:107) = 12 (cid:107) ¯ x ( X ) θ L +1 − Y(cid:107) , (S18)and applying gradient flow to optimize the readout layer (and freezing all other parameters), ˙ θ L +1 = − η ¯ x ( X ) T (cid:0) ¯ x ( X ) θ L +1 − Y (cid:1) , (S19)where η 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 ¯ x ( X )¯ x ( X ) T is invertible, it is f t ( x ∗ ) = f ( x ∗ ) + ˆ K ( x, X ) ˆ K ( X , X ) − (cid:16) exp (cid:16) − ηt ˆ K ( X , X ) (cid:17) − I (cid:17) ( f ( X ) − Y ) (S20)For any x, x (cid:48) ∈ R n , letting n l → ∞ for l = 1 , . . . , L , one has the convergence in distribution inprobability and distribution respectively ¯ x ( x )¯ x ( x (cid:48) ) → K ( x, x (cid:48) ) and ¯ x ( X ) θ L +10 → N (0 , K ( X , X )) . (S21)Moreover ¯ x ( X ) θ L +10 and the term containing f ( X ) are the only stochastic term over the ensembleof network initializations, therefore for any t the output f ( x ∗ ) throughout training converges to aGaussian distribution in the infinite width limit, with E [ f t ( x ∗ )] = K ( x ∗ , X ) K − ( I − e − η K t ) Y , (S22) Var[ f t ( x ∗ )] = K ( x ∗ , x ∗ ) − K ( x ∗ , X ) K − ( I − e − η K t ) K ( x ∗ , X ) T . (S23)Thus the output of the neural network is also a GP and the asymptotic solution (i.e. t → ∞ ) isidentical to the posterior of the NNGP (Equation 13). Therefore, in the infinite width case, theoptimized 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]. E Computing NTK and NNGP Kernel
For completeness, we reproduce, informally, the recursive formula of the NNGP kernel and thetangent kernel from [5] and [13], respectively. Let the activation function φ : R → R be absolutely5ontinuous. Let T and ˙ T be functions from × positive semi-definite matrices Σ to R given by (cid:26) T (Σ) = E [ φ ( u ) φ ( v )]˙ T (Σ) = E [ φ (cid:48) ( u ) φ (cid:48) ( v )] ( u, v ) ∼ N (0 , Σ) . (S24)In the infinite width limit, the NNGP and tangent kernel can be computed recursively. Let x, x (cid:48) be two inputs in R n . Then h l ( x ) and h l ( x (cid:48) ) converge in distribution to a joint Gaussian as min { n , . . . , n l − } . The mean is zero and the variance K l ( x, x (cid:48) ) is K l ( x, x (cid:48) ) = ˜ K l ( x, x (cid:48) ) ⊗ Id n l (S25) ˜ K l ( x, x (cid:48) ) = σ ω T (cid:18)(cid:20) ˜ K l − ( x, x ) ˜ K l − ( x, x (cid:48) )˜ K l − ( x, x (cid:48) ) ˜ K l − ( x (cid:48) , x (cid:48) ) (cid:21)(cid:19) + σ b (S26)with base case K ( x, x (cid:48) ) = σ ω · n x T x (cid:48) + σ b . (S27)Using this one can also derive the tangent kernel for gradient descent training. We will use inductionto show that Θ l ( x, x (cid:48) ) = ˜Θ l ( x, x (cid:48) ) ⊗ Id n l (S28)where ˜Θ l ( x, x (cid:48) ) = ˜ K l ( x, x (cid:48) ) + σ ω ˜Θ l − ( x, x (cid:48) ) ˙ T (cid:18)(cid:20) ˜ K l − ( x, x ) ˜ K l − ( x, x (cid:48) )˜ K l − ( x, x (cid:48) ) ˜ K l − ( x (cid:48) , x (cid:48) ) (cid:21)(cid:19) (S29)with ˜Θ = ˜ K . Let J l ( x ) = ∇ θ ≤ l h l ( x ) = [ ∇ θ l h l ( x ) , ∇ θ In this Section we present a sketch for why the function space linearization results, derived in [13] forNTK parameterized networks, also apply to networks with a standard parameterization. We followthis up with a formal proof in §G of the convergence of standard parameterization networks to theirlinearization in the limit of infinite width. A network with standard parameterization is described as: (cid:26) h l +1 = x l W l +1 + b l +1 x l +1 = φ (cid:0) h l +1 (cid:1) and (cid:40) W li,j = ω lij ∼ N (cid:16) , σ ω n l (cid:17) b lj = β lj ∼ N (cid:0) , σ b (cid:1) . (S36)6 -9 -7 -5 -3 -1 η A cc u r a c y relu-xent-ntkrelu-mse-ntktanh-xent-ntktanh-mse-ntkrelu-xent-standardrelu-mse-standardtanh-xent-standardtanh-mse-standard (a) MNIST -13 -11 -9 -7 -5 -3 -1 η A cc u r a c y relu-xent-ntkrelu-mse-ntktanh-xent-ntktanh-mse-ntkrelu-xent-standardrelu-mse-standardtanh-xent-standardtanh-mse-standard (b) CIFAR Figure S7: NTK vs Standard parameterization. Across different choices of dataset, activationfunction and loss function, models obtained from (S)GD training for both parameterization (circleand triangle denotes NTK and standard parameterization respectively) get similar performance.The NTK parameterization in Equation 1 is not commonly used for training neural networks. Whilethe 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. How-ever, for a particular choice of layer-dependent learning rate training dynamics also become identical.Let η l NTK ,w and η l NTK ,b be layer-dependent learning rate for W l and b l in the NTK parameterization,and η std = n max η be the learning rate for all parameters in the standard parameterization, where n max = max l n l . Recall that gradient descent training in standard neural networks requires a learningrate that scales with width like n max , so η defines a width-invariant learning rate [31]. If we choose η l NTK, w = n l n max σ ω η , and η l NTK, b = 1 n max σ b η , (S37)then learning dynamics are identical for networks with NTK and standard parameterizations. Withonly extremely minor modifications, consisting of incorporating the multiplicative factors in EquationS37 into the per-layer contributions to the Jacobian, the arguments in §2.4 go through for an NTKnetwork with learning rates defined in Equation S37. Since an NTK network with these learning ratesexhibits identical training dynamics to a standard network with learning rate η std , the result in §2.4that sufficiently wide NTK networks are linear in their parameters throughout training also applies tostandard networks.We can verify this property of networks with the standard parameterization experimentally. InFigure S7, we see that for different choices of dataset, activation function and loss function, finalperformance of two different parameterization leads to similar quality model for similar value ofnormalized learning rate η std = η NTK /n . Also, in Figure S8, we observe that our results is not due tothe parameterization choice and holds for wide networks using the standard parameterization. G Convergence of neural network to its linearization, and stability of NTKunder gradient descent In this section, we show that how to use the NTK to provide a simple proof of the global convergenceof 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 alsoapply to the NTK parameterization. To lighten the notation, we only consider the asymptotic boundhere. The neural networks are parameterized as in Equation S36. We make the following assumptions: Assumptions [1-4]: 1. The widths of the hidden layers are identical, i.e. n = · · · = n L = n (our proof extendsnaturally to the setting n l n l (cid:48) → α l,l (cid:48) ∈ (0 , ∞ ) as min { n , . . . , n L } → ∞ .)7 O u t p u t V a l u e Train output Neural NetworkLinearized Model 210123 Test output -2e-05-1.5e-05-1e-05-5e-060.05e-061e-051.5e-052e-052.5e-05 Weight change( ω ) -2 -1 t Loss TrainTrain f lin TestTest f lin -2 -1 t Accuracy -2 -1 t -3 -2 -1 q › ( f ( x ) − f lin ( x )) fi Figure S8: Exact and experimental dynamics are nearly identical for network outputs, andare similar for individual weights (Standard parameterization). Experiment is for an MSE loss, ReLU network with 5 hidden layers of width n = 2048 , η = 0 . / |D| = 256 , k = 1 , σ w = 2 . , and σ b = 0 . . All three panes in the first row show dynamics for a randomly selectedsubset of datapoints or parameters. First two panes in the second row show dynamics of loss andaccuracy for training and test points agree well between original and linearized model. Bottom rightpane shows the dynamics of RMSE between the two models on test points using empirical kernel.2. The analytic NTK Θ (defined in Equation S42) is full-rank, i.e. < λ min := λ min (Θ) ≤ λ max := λ max (Θ) < ∞ . We set η critical = 2( λ min + λ max ) − .3. The training set ( X , Y ) is contained in some compact set and x (cid:54) = ˜ x for all x, ˜ x ∈ X .4. The activation function φ satisfies | φ (0) | , (cid:107) φ (cid:48) (cid:107) ∞ , sup x (cid:54) =˜ x | φ (cid:48) ( x ) − φ (cid:48) (˜ x ) | / | x − ˜ x | < ∞ . (S38)Assumption 2 indeed holds when X ⊆ { x ∈ R n } : (cid:107) x (cid:107) = 1 } and φ ( x ) grows non-polynomiallyfor large x [13]. Throughout this section, we use C > to denote some constant whose value maydepend on L , |X | and ( σ w , σ b ) and may change from line to line, but is always independent of n .Let θ t denote the parameters at time step t . We use the following short-hand f ( θ t ) = f ( X , θ t ) ∈ R |X |× k (S39) g ( θ t ) = f ( X , θ t ) − Y ∈ R |X |× k (S40) J ( θ t ) = ∇ θ f ( θ t ) ∈ R ( |X | k ) ×| θ | (S41)where |X | is the cardinality of the training set and k is the output dimension of the network. Theempirical and analytic NTK of the standard parameterization is defined as (cid:40) ˆΘ t := ˆΘ t ( X , X ) = n J ( θ t ) J ( θ t ) T Θ := lim n →∞ ˆΘ in probability . (S42)Note that the convergence of the empirical NTK in probability is proved rigorously in [37]. Weconsider the MSE loss L ( t ) = 12 (cid:107) g ( θ t ) (cid:107) . (S43)8ince f ( θ t ) converges in distribution to a mean zero Guassian with covariance K , one can show thatfor arbitrarily small δ > , there are constants R > and n (both may depend on δ , |X | and K )such that for every n ≥ n , with probability at least (1 − δ ) over random initialization, (cid:107) g ( θ ) (cid:107) < R . (S44)The gradient descent update with learning rate η is θ t +1 = θ t − ηJ ( θ t ) T g ( θ t ) (S45)and the gradient flow equation is ˙ θ t = − J ( θ t ) T g ( θ t ) . (S46)We prove convergence of neural network training and the stability of NTK for both discrete gradientdescent and gradient flow. Both proofs rely on the local lipschitzness of the Jacobian J ( θ ) . Lemma 1 ( Local Lipschitzness of the Jacobian). There is a K > such that for every C > ,with high probability over random initialization (w.h.p.o.r.i.) the following holds √ n (cid:107) J ( θ ) − J (˜ θ ) (cid:107) F ≤ K (cid:107) θ − ˜ θ (cid:107) √ n (cid:107) J ( θ ) (cid:107) F ≤ K , ∀ θ, ˜ θ ∈ B ( θ , Cn − ) (S47) where B ( θ , R ) := { θ : (cid:107) θ − θ (cid:107) < R } . (S48)The following are the main results of this section. Theorem G.1 ( Gradient descent). Assume Assumptions [1-4] . For δ > and η < η critical , thereexist R > , N ∈ N and K > , such that for every n ≥ N , the following holds with probability atleast (1 − δ ) over random initialization when applying gradient descent with learning rate η = η n , (cid:107) g ( θ t ) (cid:107) ≤ (cid:16) − η λ min (cid:17) t R (cid:80) tj =1 (cid:107) θ j − θ j − (cid:107) ≤ η KR √ n (cid:80) tj =1 (1 − η λ min ) j − ≤ KR λ min n − (S49) and sup t (cid:107) ˆΘ − ˆΘ t (cid:107) F ≤ K R λ min n − . (S50) Theorem G.2 ( Gradient Flow). Assume Assumptions[1-4] . For δ > , there exist R > , N ∈ N and K > , such that for every n ≥ N , the following holds with probability at least (1 − δ ) over random initialization when applying gradient flow with “learning rate" η = η n (cid:107) g ( θ t ) (cid:107) ≤ e − η λ min3 t R (cid:107) θ t − θ (cid:107) ≤ KR λ min (1 − e − η λ min t ) n − (S51) and sup t (cid:107) ˆΘ − ˆΘ t (cid:107) F ≤ K R λ min n − . (S52)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 orfunctions as long as1. The empirical NTK converges in probability and the limit is positive definite.2. Lemma 1 holds, i.e. the Jacobian is locally Lipschitz. .1 Proof of Theorem G.1 As discussed above, there exist R and n such that for every n ≥ n , with probability at least (1 − δ / over random initialization, (cid:107) g ( θ ) (cid:107) < R . (S53)Let C = KR λ min in Lemma 1. We first prove Equation S49 by induction. Choose n > n such thatfor every n ≥ n Equation S47 and Equation S53 hold with probability at least (1 − δ / overrandom initialization. The t = 0 case is obvious and we assume Equation S49 holds for t = t . Thenby induction and the second estimate of Equation S47 (cid:107) θ t +1 − θ t (cid:107) ≤ η (cid:107) J ( θ t ) (cid:107) op (cid:107) g ( θ t ) (cid:107) ≤ Kη √ n (cid:18) − η λ min (cid:19) t R , (S54)which gives the first estimate of Equation S49 for t +1 and which also implies (cid:107) θ j − θ (cid:107) ≤ KR λ min n − for j = 0 , . . . , t + 1 . To prove the second one, we apply the mean value theorem and the formula forgradient decent update at step t + 1 (cid:107) g ( θ t +1 ) (cid:107) = (cid:107) g ( θ t +1 ) − g ( θ t ) + g ( θ t ) (cid:107) (S55) = (cid:107) J (˜ θ t )( θ t +1 − θ t ) + g ( θ t ) (cid:107) (S56) = (cid:107) − ηJ (˜ θ t ) J ( θ t ) T g ( θ t ) + g ( θ t ) (cid:107) (S57) ≤ (cid:107) − ηJ (˜ θ t ) J ( θ t ) T (cid:107) op (cid:107) g ( θ t ) (cid:107) (S58) ≤ (cid:107) − ηJ (˜ θ t ) J ( θ t ) T (cid:107) op (cid:18) − η λ min (cid:19) t R , (S59)where ˜ θ t is some linear interpolation between θ t and θ t +1 . It remains to show with probability atleast (1 − δ / , (cid:107) − ηJ (˜ θ t ) J ( θ t ) T (cid:107) op ≤ − η λ min . (S60)This can be verified by Lemma 1. Because ˆΘ → Θ [37] in probability, one can find n such that theevent (cid:107) Θ − ˆΘ (cid:107) F ≤ η λ min (S61)has probability at least (1 − δ / for every n ≥ n . The assumption η < λ min + λ max implies (cid:107) − η Θ (cid:107) op ≤ − η λ min . (S62)Thus (cid:107) − ηJ (˜ θ t ) J ( θ t ) T (cid:107) op (S63) ≤(cid:107) − η Θ (cid:107) op + η (cid:107) Θ − ˆΘ (cid:107) op + η (cid:107) J ( θ ) J ( θ ) T − J (˜ θ t ) J ( θ t ) T (cid:107) op (S64) ≤ − η λ min + η λ min η K ( (cid:107) θ t − θ (cid:107) + (cid:107) ˜ θ t − θ (cid:107) ) (S65) ≤ − η λ min + η λ min η K KR λ min √ n ≤ − η λ min (S66)with probability as least (1 − δ / if n ≥ (cid:18) K R λ (cid:19) . (S67)Therefore, we only need to set N = max (cid:40) n , n , n , (cid:18) K R λ (cid:19) (cid:41) . (S68)10o verify Equation S50, notice that (cid:107) ˆΘ − ˆΘ t (cid:107) F = 1 n (cid:107) J ( θ ) J ( θ ) T − J ( θ t ) J ( θ t ) T (cid:107) F (S69) ≤ n (cid:0) (cid:107) J ( θ ) (cid:107) op (cid:107) J ( θ ) T − J ( θ t ) T (cid:107) F + (cid:107) J ( θ t ) − J ( θ ) (cid:107) op (cid:107) J ( θ t ) T (cid:107) F (cid:1) (S70) ≤ K (cid:107) θ − θ t (cid:107) (S71) ≤ K R λ min √ n , (S72)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 and n such that for every n ≥ n , with probability at least (1 − δ / over random initialization, (cid:107) g ( θ ) (cid:107) < R . (S73)Let C = KR λ min in Lemma 1. Using the same arguments as in Section G.1, one can show that thereexists n such that for all n ≥ n , with probability at least (1 − δ / n J ( θ ) J ( θ ) T (cid:31) λ min Id ∀ θ ∈ B ( θ , Cn − ) (S74)Let t = inf (cid:26) t : (cid:107) θ t − θ (cid:107) ≥ KR λ min n − (cid:27) (S75)We claim t = ∞ . If not, then for all t ≤ t , θ t ∈ B ( θ , Cn − ) and ˆΘ t (cid:31) λ min Id . (S76)Thus ddt (cid:0) (cid:107) g ( t ) (cid:107) (cid:1) = − η g ( t ) T ˆΘ t g ( t ) ≤ − η λ min (cid:107) g ( t ) (cid:107) (S77)and (cid:107) g ( t ) (cid:107) ≤ e − η λ min t (cid:107) g (0) (cid:107) ≤ e − η λ min t R . (S78)Note that ddt (cid:107) θ t − θ (cid:107) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ddt θ t (cid:13)(cid:13)(cid:13)(cid:13) = η n (cid:107) J ( θ t ) g ( t ) (cid:107) ≤ η KR e − η λ min t n − / (S79)which implies, for all t ≤ t (cid:107) θ t − θ (cid:107) ≤ KR λ min (1 − e − η λ min t ) n − ≤ KR λ min (1 − e − η λ min t ) n − < KR λ min n − . (S80)This contradicts to the definition of t and thus t = ∞ . Note that Equation S78 is the same as thefirst 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 × n random matrix whose entries areindependent standard normal random variables. Then for every t ≥ , with probability at least − − t / one has √ N − √ n − t ≤ λ min ( A ) ≤ λ max ( A ) ≤ √ N + √ n + t. (S81)11or l ≥ , let δ l ( θ, x ) := ∇ h l ( θ,x ) f L +1 ( θ, x ) ∈ R kn (S82) δ l ( θ, X ) := ∇ h l ( θ, X ) f L +1 ( θ, X ) ∈ R ( k ×|X | ) × ( n ×X ) (S83)Let θ = { W l , b l } and ˜ θ = { ˜ W l , ˜ b l } be any two points in B ( θ , C √ n ) . By the above theorem and thetriangle inequality, w.h.p. over random initialization, (cid:107) W (cid:107) op , (cid:107) ˜ W (cid:107) op ≤ σ ω √ n √ n , (cid:107) W l (cid:107) op , (cid:107) ˜ W l (cid:107) op ≤ σ ω for 2 ≤ l ≤ L + 1 (S84)Using this and the assumption on φ Equation S38, it is not difficult to show that there is a constant K , depending on σ ω , σ b , |X | and L such that with high probability over random initialization n − (cid:107) x l ( θ, X ) (cid:107) , (cid:107) δ l ( θ, X ) (cid:107) ≤ K , (S85) n − (cid:107) x l ( θ, X ) − x l (˜ θ, X ) (cid:107) , (cid:107) δ l ( θ, X ) − δ l (˜ θ, X ) (cid:107) ≤ K (cid:107) ˜ θ − θ (cid:107) (S86)Lemma 1 follows from these two estimates. Indeed, with high probability over random initialization (cid:107) J ( θ ) (cid:107) F = (cid:88) l (cid:107) J ( W l ) (cid:107) F + (cid:107) J ( b l ) (cid:107) F (S87) = (cid:88) l (cid:88) x ∈X (cid:107) x l − ( θ, x ) δ l ( θ, x ) T (cid:107) F + (cid:107) δ l ( θ, x ) T (cid:107) F (S88) ≤ (cid:88) l (cid:88) x ∈X (1 + (cid:107) x l − ( θ, x ) (cid:107) F ) (cid:107) δ l ( θ, x ) T (cid:107) F (S89) ≤ (cid:88) l (1 + K n ) (cid:88) x (cid:107) δ l ( θ, x ) T (cid:107) F (S90) ≤ (cid:88) l K (1 + K n ) (S91) ≤ L + 1) K n, (S92)and similarly (cid:107) J ( θ ) − J (˜ θ ) (cid:107) F (S93) = (cid:88) l (cid:88) x ∈X (cid:107) x l − ( θ, x ) δ l ( θ, x ) T − x l − (˜ θ, x ) δ l (˜ θ, x ) T (cid:107) F + (cid:107) δ l ( θ, x ) T − δ l (˜ θ, x ) T (cid:107) F (S94) ≤ (cid:32)(cid:88) l (cid:0) K n + K n (cid:1) + K (cid:33) (cid:107) θ − ˜ θ (cid:107) (S95) ≤ L + 1) K n (cid:107) θ − ˜ θ (cid:107) . (S96) G.4 Remarks on NTK parameterization For completeness, we also include analogues of Theorem G.1 and Lemma 1 with NTK parameteriza-tion. Theorem G.4 (NTK parameterization) . Assume Assumptions [1-4] . For δ > and η < η critical ,there exist R > , N ∈ N and K > , such that for every n ≥ N , the following holds withprobability at least (1 − δ ) over random initialization when applying gradient descent with learningrate η = η , (cid:107) g ( θ t ) (cid:107) ≤ (cid:16) − η λ min (cid:17) t R (cid:80) tj =1 (cid:107) θ j − θ j − (cid:107) ≤ Kη (cid:80) tj =1 (1 − η λ min ) j − R ≤ KR λ min (S97) These two estimates can be obtained via induction. To prove bounds relating to x l and δ l , one starts with l = 1 and l = L , respectively. nd sup t (cid:107) ˆΘ − ˆΘ t (cid:107) F ≤ K R λ min n − . (S98) Lemma 2 (NTK parameterization: Local Lipschitzness of the Jacobian) . There is a K > such thatfor every C > , with high probability over random initialization the following holds (cid:107) J ( θ ) − J (˜ θ ) (cid:107) F ≤ K (cid:107) θ − ˜ θ (cid:107) (cid:107) J ( θ ) (cid:107) F ≤ K , ∀ θ, ˜ θ ∈ B ( θ , C ) (S99) H Bounding the discrepancy between the original and the linearizednetwork: MSE loss We provide the proof for the gradient flow case. The proof for gradient descent can be obtainedsimilarly. To simplify the notation, let g lin ( t ) ≡ f lin t ( X ) − Y and g ( t ) ≡ f t ( X ) − Y . The theoremand proof apply to both standard and NTK parameterization. We use the notation (cid:46) to hide thedependence on uninteresting constants. Theorem H.1. Same as in Theorem G.2. For every x ∈ R n with (cid:107) x (cid:107) ≤ , for δ > arbitrarilysmall, there exist R > and N ∈ N such that for every n ≥ N , with probability at least (1 − δ ) over random initialization, sup t (cid:13)(cid:13) g lin ( t ) − g ( t ) (cid:13)(cid:13) , sup t (cid:13)(cid:13) g lin ( t, x ) − g ( t, x ) (cid:13)(cid:13) (cid:46) n − R . (S100) Proof. ddt (cid:16) exp( η ˆΘ t )( g lin ( t ) − g ( t )) (cid:17) (S101) = η (cid:16) ˆΘ exp( η ˆΘ t )( g lin ( t ) − g ( t )) + exp( η ˆΘ t )( − ˆΘ g lin ( t ) + ˆΘ t g ( t )) (cid:17) (S102) = η (cid:16) exp( η ˆΘ t )( ˆΘ t − ˆΘ ) g ( t ) (cid:17) (S103)Integrating both sides and using the fact g lin (0) = g (0) , ( g lin ( t ) − g ( t )) = − (cid:90) t η (cid:16) exp( η ˆΘ ( s − t ))( ˆΘ s − ˆΘ )( g lin ( s ) − g ( s )) (cid:17) ds (S104) + (cid:90) t η (cid:16) exp( η ˆΘ ( s − t ))( ˆΘ s − ˆΘ ) g lin ( s ) (cid:17) ds (S105)Let λ > be the smallest eigenvalue of ˆΘ (with high probability λ > λ min ). Taking the normgives (cid:107) g lin ( t ) − g ( t ) (cid:107) ≤ η (cid:16) (cid:90) t (cid:107) exp( ˆΘ η ( s − t )) (cid:107) op (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op (cid:107) g lin ( s ) − g ( s ) (cid:107) ds (S106) + (cid:90) t (cid:107) exp( ˆΘ η ( s − t )) (cid:107) op (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op (cid:107) g lin ( s ) (cid:107) ds (cid:17) (S107) ≤ η (cid:16) (cid:90) t e η λ ( s − t ) (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op (cid:107) g lin ( s ) − g ( s ) (cid:107) ds (S108) + (cid:90) t e η λ ( s − t ) (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op (cid:107) g lin ( s ) (cid:107) ds (cid:17) (S109)Let u ( t ) ≡ e λ η t (cid:107) g lin ( t ) − g ( t ) (cid:107) (S110) α ( t ) ≡ η (cid:90) t e λ η s (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op (cid:107) g lin ( s ) (cid:107) ds (S111) β ( t ) ≡ η (cid:107) ( ˆΘ t − ˆΘ ) (cid:107) op (S112)13he above can be written as u ( t ) ≤ α ( t ) + (cid:90) t β ( s ) u ( s ) ds (S113)Note that α ( t ) is non-decreasing. Applying an integral form of the Grönwall’s inequality (seeTheorem 1 in [38]) gives u ( t ) ≤ α ( t ) exp (cid:18)(cid:90) t β ( s ) ds (cid:19) (S114)Note that (cid:107) g lin ( t ) (cid:107) = (cid:107) exp (cid:16) − η ˆΘ t (cid:17) g lin (0) (cid:107) ≤ (cid:107) exp (cid:16) − η ˆΘ t (cid:17) (cid:107) op (cid:107) g lin (0) (cid:107) = e − λ η t (cid:107) g lin (0) (cid:107) . (S115)Then (cid:107) g lin ( t ) − g ( t ) (cid:107) ≤ η e − λ η t (cid:90) t e λ η s (cid:107) ˆΘ s − ˆΘ (cid:107) op (cid:107) g lin ( s ) (cid:107) ds exp (cid:18)(cid:90) t η (cid:107) ˆΘ s − ˆΘ (cid:107) op ds (cid:19) (S116) ≤ η e − λ η t (cid:107) g lin (0) (cid:107) (cid:90) t (cid:107) ( ˆΘ s − ˆΘ ) (cid:107) op ds exp (cid:18)(cid:90) t η (cid:107) ˆΘ s − ˆΘ (cid:107) op ds (cid:19) (S117)Let σ t = sup ≤ s ≤ t (cid:107) ˆΘ s − ˆΘ (cid:107) op . Then (cid:107) g lin ( t ) − g ( t ) (cid:107) (cid:46) (cid:0) η tσ t e − λ η t + σ t η t (cid:1) (cid:107) g lin (0) (cid:107) (S118)As it is proved in Theorem G.1, for every δ > , with probability at least (1 − δ ) over randominitialization, sup t σ t ≤ sup t (cid:107) ˆΘ − ˆΘ t (cid:107) F (cid:46) n − / R → (S119)when n = · · · = n L = n → ∞ . Thus for large n and any polynomial P ( t ) (we use P ( t ) = t here) sup t e − λ η t + σ t η t η P ( t ) = O (1) (S120)Therefore sup t (cid:107) g lin ( t ) − g ( t ) (cid:107) (cid:46) sup t σ t R (cid:46) n − / R → , (S121)as n → ∞ .Now we control the discrepancy on a test point x . Let y be its true label. Similarly, ddt (cid:0) g lin ( t, x ) − g ( t, x ) (cid:1) = − η (cid:16) ˆΘ ( x, X ) − ˆΘ t ( x, X ) (cid:17) g lin ( t ) + η ˆΘ t ( x, X )( g ( t ) − g lin ( t )) . (S122)Integrating over [0 , t ] and taking the norm imply (cid:13)(cid:13) g lin ( t, x ) − g ( t, x ) (cid:13)(cid:13) (S123) ≤ η (cid:90) t (cid:13)(cid:13)(cid:13) ˆΘ ( x, X ) − ˆΘ s ( x, X ) (cid:13)(cid:13)(cid:13) (cid:107) g lin ( s ) (cid:107) ds + η (cid:90) t (cid:107) ˆΘ s ( x, X ) (cid:107) (cid:107) g ( s ) − g lin ( s ) (cid:107) ds (S124) ≤ η (cid:107) g lin (0) (cid:107) (cid:90) t (cid:13)(cid:13)(cid:13) ˆΘ ( x, X ) − ˆΘ s ( x, X ) (cid:13)(cid:13)(cid:13) e − η λ s ds (S125) + η (cid:90) t ( (cid:107) ˆΘ ( x, X ) (cid:107) + (cid:107) ˆΘ s ( x, X ) − ˆΘ ( x, X ) (cid:107) ) (cid:107) g ( s ) − g lin ( s ) (cid:107) ds (S126)14 Width ( n ) -9 -8 -7 -6 -5 -4 -3 -2 -1 L = 1 , || ˆ K ( M,n ) − K|| F / ||K|| F M=1M=2M=4 M=8M=16M=32 M=64M=128M=256 Width ( n ) L = 3 , || ˆ K ( M,n ) − K|| F / ||K|| F Width ( n ) L = 1 , || ˆΘ ( M,n ) − Θ || F / || Θ || F Width ( n ) L = 3 , || ˆΘ ( M,n ) − Θ || F / || Θ || F Figure S9: Kernel convergence. Kernels computed from randomly initialized ReLU networks withone and three hidden layers converge to the corresponding analytic kernel as width n and number ofMonte Carlo samples M increases. Colors indicate averages over different numbers of Monte Carlosamples.Similarly, Lemma 1 implies sup t (cid:13)(cid:13)(cid:13) ˆΘ ( x, X ) − ˆΘ t ( x, X ) (cid:13)(cid:13)(cid:13) (cid:46) n − R (S127)This gives ( S (cid:46) n − R . (S128)Using Equation S118 and Equation S119,(S126) (cid:46) (cid:107) ˆΘ ( x, X ) (cid:107) (cid:90) t (cid:0) η sσ s e − λ η s + σ s η s (cid:1) (cid:107) g lin (0) (cid:107) dt (cid:46) n − . (S129) I Convergence of empirical kernel As in Novak et al. [7], we can use Monte Carlo estimates of the tangent kernel (Equation 4) to probeconvergence to the infinite width kernel (analytically computed using Equations S26, S29). Forsimplicity, we consider random inputs drawn from N (0 , with n = 1024 . In Figure S9, we observeconvergence as both width n increases and the number of Monte Carlo samples M increases. For bothNNGP and tangent kernels we observe (cid:107) ˆΘ ( n ) − Θ (cid:107) F = O (1 / √ n ) and (cid:107) ˆ K ( n ) − K(cid:107) F = O (1 / √ n ) ,as predicted by a CLT in Daniely et al. [16]. J Details on Wide Residual Network Table S1: Wide Residual Network architecture from Zagoruyko and Komodakis [14] . In theresidual block, we follow Batch Normalization-ReLU-Conv ordering.group name output size block typeconv1 32 × 32 [3 × 3, channel size]conv2 32 × (cid:20) × , channel size × , channel size (cid:21) × Nconv3 16 × (cid:20) × , channel size × , channel size (cid:21) × Nconv4 8 × (cid:20) × , channel size × , channel size (cid:21) × Navg-pool 1 × × -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 | K ∞ − K M , n | F | K ∞ | F l og ( n ) Number of Samples ( M )2 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 | Θ ∞ − Θ M , n | F | Θ ∞ | F n=1n=2n=4n=8n=16 n=32n=64n=128n=256 n=512n=1024n=2048n=4096 log ( M ) l og ( n ) Figure S10: Kernel convergence. Kernels from single hidden layer randomly initialized