Title: The asymptotic spectrum of the Hessian of DNN throughout training

Abstract: The dynamics of DNNs during gradient descent is described by the so-called Neural Tangent Kernel (NTK). In this article, we show that the NTK allows one to gain precise insight into the Hessian of the cost of DNNs. When the NTK is fixed during training, we obtain a full characterization of the asymptotics of the spectrum of the Hessian, at initialization and during training. In the so-called mean-field limit, where the NTK is not fixed during training, we describe the first two moments of the Hessian at initialization.
We consider deep fully connected artificial neural networks (DNNs) using the setup and NTK parametrization of (Jacot et al., 2018), taking an arbitrary nonlinearity σ ∈ C 4 b (R) (i.e. σ : R → R that is 4 times continuously differentiable function with all four derivatives bounded). The layers are numbered from 0 (input) to L (output), each containing n neurons for = 0, . . . , L. The P = L-1 =0 (n + 1) n +1 parameters consist of the weight matrices W ( ) ∈ R n +1 ×n and bias vectors b ( ) ∈ R n +1 for = 0, . . . , L -1. We aggregate the parameters into the vector θ ∈ R P .

1. Introduction:

The advent of deep learning has sparked a lot of interest in the loss surface of deep neural networks (DNN), and in particular its Hessian. However to our knowledge, there is still no theoretical description of the spectrum of the Hessian. Nevertheless a number of phenomena have been observed numerically.

The loss surface of neural networks has been compared to the energy landscape of different physical models (Choromanska et al., 2015; Geiger et al., 2018; Mei et al., 2018 ). It appears that the loss surface of DNNs may change significantly depending on the width of the network (the number of neurons in the hidden layer), motivating the distinction between the underand over-parametrized regimes (Baity-Jesi et al., 2018; Geiger et al., 2018; 2019) .

The non-convexity of the loss function implies the existence of a very large number of saddle points, which could slow down training. In particular, in (Pascanu et al., 2014; Dauphin et al., 2014) , a relation between the rank of saddle points (the number of negative eigenvalues of the Hessian) and their loss has been observed.

For overparametrized DNNs, a possibly more important phenomenon is the large number of flat directions (Baity-Jesi et al., 2018) . The existence of these flat minima is conjectured to be related to the generalization of DNNs and may depend on the training procedure (Hochreiter & Schmidhuber, 1997; Chaudhari et al., 2016; Wu et al., 2017) .

In (Jacot et al., 2018 ) it has been shown, using a functional approach, that in the infinitewidth limit, DNNs behave like kernel methods with respect to the so-called Neural Tangent Kernel, which is determined by the architecture of the network. This leads to convergence guarantees for DNNs (Jacot et al., 2018; Du et al., 2019; Allen-Zhu et al., 2018; Huang & Yau, 2019) and strengthens the connections between neural networks and kernel methods (Neal, 1996; Cho & Saul, 2009; Lee et al., 2018) .

Our approach also allows one to probe the so-called mean-field/active limit (studied in (Rotskoff & Vanden-Eijnden, 2018; Chizat & Bach, 2018a; Mei et al., 2018) for shallow networks), where the NTK varies during training.

This raises the question: can we use these new results to gain insight into the behavior of the Hessian of the loss of DNNs, at least in the small region explored by the parameters during training?

1.1. Contributions:

Following ideas introduced in (Jacot et al., 2018) , we consider the training of L + 1-layered DNNs in a functional setting. For a functional cost C, the Hessian of the loss R P θ → C F (L) (θ) is the sum of two P × P matrices I and S. We show the following results for large P and for a fixed number of datapoints N :

• The first matrix I is positive semi-definite and its eigenvalues are given by the (weighted) kernel PCA of the dataset with respect to the NTK. The dominating eigenvalues are the principal components of the data followed by a high number of small eigenvalues. The "flat directions" are spanned by the small eigenvalues and the null-space (of dimension at least P -N when there is a single output).

Because the NTK is asymptotically constant (Jacot et al., 2018) , these results apply at initialization, during training and at convergence.

• The second matrix S can be viewed as residual contribution to H, since it vanishes as the network converges to a global minimum. We compute the limit of the first moment Tr (S) and characterize its evolution during training, of the second moment Tr S 2 which stays constant during training, and show that the higher moments vanish.

• Regarding the sum H = I + S, we show that the matrices I and S are asymptotically orthogonal to each other at initialization and during training. In particular, the moments of the matrices I and S add up: tr(H k ) ≈ tr(I k ) + tr(S k ).

These results give, for any depth and a fairly general non-linearity, a complete description of the spectrum of the Hessian in terms of the NTK at initialization and throughout training.

Our theoretical results are consistent with a number of observations about the Hessian (Hochreiter & Schmidhuber, 1997; Pascanu et al., 2014; Dauphin et al., 2014; Chaudhari et al., 2016; Wu et al., 2017; Pennington & Bahri, 2017; Geiger et al., 2018) , and sheds a new light on them.

1.2. Related works:

The Hessian of the loss has been studied through the decomposition I + S in a number of previous works (Sagun et al., 2017; Pennington & Bahri, 2017; Geiger et al., 2018) .

For least-squares and cross-entropy costs, the first matrix I is equal to the Fisher matrix (Wagenaar, 1998; Pascanu & Bengio, 2013) , whose moments have been described for shallow networks in (Pennington & Worah, 2018) . For deep networks, the first two moments and the operator norm of the Fisher matrix for a least squares loss were computed at initialization in (Karakida et al., 2018) conditionally on a certain independence assumption; our method does not require such assumptions. Note that their approach implicitly uses the NTK.

The second matrix S has been studied in (Pennington & Bahri, 2017; Geiger et al., 2018) for shallow networks, conditionally on a number of assumptions. Note that in the setting of (Pennington & Bahri, 2017) , the matrices I and S are assumed to be freely independent, which allows them to study the spectrum of the Hessian; in our setting, we show that the two matrices I and S are asymptotically orthogonal to each other.

The activations and pre-activations of the layers are defined recursively for an input

x ∈ R n0 , setting α (0) (x; θ) = x : α( +1) (x; θ) = 1 √ n W ( ) α ( ) (x; θ) + βb ( ) , α ( +1) (x; θ) = σ α( +1) (x; θ) .

The parameter β is added to tune the influence of the bias on training1 . All parameters are initialized as iid N (0, 1) Gaussians.

We will in particular study the network function, which maps inputs x to the activation of the output layer (before the last non-linearity):

f θ (x) = α(L) (x; θ).

In this paper, we will study the limit of various objects as n 1 , . . . , n L → ∞ sequentially, i.e. we first take n 1 → ∞, then n 2 → ∞, etc. This greatly simplifies the proofs, but they could in principle be extended to the simultaneous limit, i.e. when n 1 = ... = n L-1 → ∞. All our numerical experiments are done with 'rectangular' networks (with n 1 = ... = n L-1 ) and match closely the predictions for the sequential limit.

In the limit we study in this paper, the NTK is asymptotically fixed, as in (Jacot et al., 2018; Allen-Zhu et al., 2018; Du et al., 2019; Arora et al., 2019; Huang & Yau, 2019) . By rescaling the outputs of DNNs as the width increases, one can reach another limit where the NTK is not fixed (Chizat & Bach, 2018a; b; Rotskoff & Vanden-Eijnden, 2018; Mei et al., 2019) . Some of our results can be extended to this setting, but only at initialization (see Section 3.3). The behavior during training becomes however much more complex.

2.1. Functional viewpoint:

The network function lives in a function space

f θ ∈ F := [R n0 → R n L ]

and we call the function F (L) : R P → F that maps the parameters θ to the network function f θ the realization function. We study the differential behavior of F (L) :

• The derivative DF (L) ∈ R P ⊗ F is a function-valued vector of dimension P . The p-th entry D p F (L) = ∂ θp f θ ∈ F represents how modifying the parameter θ p modifies the function f θ in the space F. • The Hessian HF (L) ∈ R P ⊗ R P ⊗ F is a function-valued P × P matrix.

The network is trained with respect to the cost functional:

C(f ) = 1 N N i=1 c i (f (x i )) ,

for strictly convex c i , summing over a finite dataset x 1 , . . . , x N ∈ R n0 of size N . The parameters are then trained with gradient descent on the composition C • F (L) , which defines the usual loss surface of neural networks.

In this setting, we define the finite realization function Y (L) : R P → R N n L mapping parameters θ to be the restriction of the network function f θ to the training set

y ik = f θ,k (x i ).

The Jacobian DY (L) is hence an N n L × P matrix and its Hessian HY (L) is a L) . For our analysis, we require that the gradient norm DC does not explode during training. The following condition is sufficient: Definition 1. A loss C : R N n L → R has bounded gradients over sublevel sets (BGOSS) if the norm of the gradient is bounded over all sets

P × P × N n L tensor. Defining the restricted cost C(y) = 1 N i c i (y i ), we have C • F (L) = C • Y (

U a = Y ∈ R N n L : C(Y ) ≤ a .

For example, the Mean Square Error (MSE)

C(Y ) = 1 2N Y * -Y 2 for the labels Y * ∈ R N n L has BGOSS because ∇C(Y ) 2 = 1 N Y * -Y 2 = 2C(Y ).

For the binary and softmax cross-entropy the gradient is uniformly bounded, see Proposition 2 in Appendix A.

2.2. Neural Tangent Kernel:

The behavior during training of the network function f θ in the function space F is described by a (multi-dimensional) kernel, the Neural Tangent Kernel (NTK)

Θ (L) k,k (x, x ) = P p=1 ∂ θp f θ,k (x)∂ θp f θ,k (x ).

During training, the function f θ follows the so-called kernel gradient descent with respect to the NTK, which is defined as

∂ t f θ(t) (x) = -∇ Θ (L) C |f θ(t) (x) := - 1 N N i=1 Θ (L) (x, x i )∇c i (f θ(t) (x i )).

In the infinite-width limit (letting n 1 → ∞, . . . , n L-1 → ∞ sequentially) and for losses with BGOSS, the NTK converges to a deterministic limit Θ

(L) → Θ (L)

∞ ⊗ Id n L , which is constant during training, uniformly on finite time intervals [0, T ] (Jacot et al., 2018) . For the MSE loss, the uniform convergence of the NTK was proven for T = ∞ in (Arora et al., 2019) .

The limiting NTK Θ (L):

∞ : R n0 × R n0 → R is constructed as follows:

1. For f, g : R → R and a kernel K :

R n0 × R n0 → R, define the kernel L f,g K : R n0 × R n0 → R by L f,g K (x 0 , x 1 ) = E (a0,a1) [f (a 0 )g(a 1

)] , for (a 0 , a 1 ) a centered Gaussian vector with covariance matrix (K(x i , x j )) i,j=0,1 . For f = g, we denote by L f K the kernel L f,f K . 2. We define the kernels Σ ( ) ∞ for each layer of the network, starting with Σ

(1) ∞ (x 0 , x 1 ) = 1 /n0(x T 0 x 1 ) + β 2 and then recursively by Σ

( +1) ∞ = L σ Σ ( ) ∞

+ β 2 , for = 1, . . . , L -1, where σ is the network non-linearity.

∞ is defined in terms of the kernels Σ ( )

∞ and the kernels Σ( ) ∞ = L σ Σ ( -1) ∞ : Θ (L) ∞ = L =1 Σ ( ) ∞ Σ( +1) ∞ . . . Σ(L) ∞ .

The NTK leads to convergence guarantees for DNNs in the infinite-width limit, and connect their generalization to that of kernel methods (Jacot et al., 2018; Arora et al., 2019) .

2.3. Gram Matrices:

For a finite dataset x 1 , . . . , x N ∈ R n0 and a fixed depth L ≥ 1, we denote by Θ ∈ R N n L ×N n L the Gram matrix of x 1 , . . . , x N with respect to the limiting NTK, defined by Θik,jm = Θ (L) ∞ (x i , x j ) δ km .

It is block diagonal because different outputs k = m are asymptotically uncorrelated.

Similarly, for any (scalar) kernel K (L) (such as the limiting kernels Σ

(L) ∞ , Λ (L) ∞ , Υ (L) ∞ , Φ (L) ∞ , Ξ (L) ∞

introduced later), we denote the Gram matrix of the datapoints by K.

3 Main Theorems

3.1. Hessian as I + S:

Using the above setup, the Hessian H of the loss C • F (L) is the sum of two terms, with the entry H p,p given by

H p,p = HC |f θ (∂ θp F, ∂ θ p F ) + DC |f θ (∂ θp,θ p F ).

For a finite dataset, the Hessian matrix H C • Y (L) is equal to the sum of two matrices I = DY (L) T HCDY (L) and S = ∇C • HY (L) where DY (L) is a N n L × P matrix, HC is a N n L × N n L matrix and HY (L) is a P × P × N n L tensor to which we apply a scalar product (denoted by •) in its last dimension with the N n L vector ∇C to obtain a P × P matrix.

Our main contribution is the following theorem, which describes the limiting moments Tr H k in terms of the moments of I and S:

Theorem 1. For any loss C with BGOSS and σ ∈ C 4 b (R), in the sequential limit

n 1 → ∞, . . . , n L-1 → ∞, we have for all k ≥ 1 Tr H (t) k ≈ Tr I (t) k + Tr S (t) k .

The limits of Tr I (t) k and Tr S (t) k can be expressed in terms of the NTK Θ

(L) ∞ , the kernels Υ (L) ∞ , Ξ (L)

∞ and the non-symmetric kernels

Φ (L) ∞ , Λ (L)

∞ defined in Appendix C:

• The moments Tr I (t)

k converge to the following limits (with the convention that i k+1 = i 1 ):

Tr I (t) k → Tr HC(Y (t)) Θ k = 1 N k N i1,...,i k =1 k m=1 c im (f θ(t) (x im ))Θ (L) ∞ (x im , x im+1 ).

• The first moment Tr (S (t)) converges to the limit:

Tr (S (t)) = (G(t)) T ∇C(Y (t)).

At initialization (G(0), Y (0)) form a Gaussian pair of N n L -vectors, independent for differing output indices k = 1, ..., n L and with covariance

E [G ik (0)G i k (0)] = δ kk Ξ (L) ∞ (x i , x i ) and E [G ik (0)Y i k (0)] = δ kk Φ (L)

∞ (x i , x i ) for the limiting kernel Ξ 

∂ t G(t) = -Λ∇C(Y (t)) ∂ t Y (t) = -Θ∇C(Y (t)).

• The second moment Tr S (t) 2 converges to the following limit defined in terms of the Gram matrix Υ:

Tr S 2 → (∇C(Y (t)))

T Υ∇C(Y (t))

• The higher moments Tr S (t) k for k ≥ 3 vanish.

Proof. The moments of I and S can be studied separately because the moments of their sum is asymptotically equal to the sum of their moments by Proposition 5 below. The limiting moments of I and S are respectively described by Propositions 1 and 4 below.

In the case of a MSE loss

C(Y ) = 1 2N Y -Y * 2

, the first and second derivatives take simple forms ∇C(Y ) = 1 N (Y -Y * ) and HC(Y ) = 1 N Id N n L and the differential equations can be solved to obtain more explicit formulae: 

10 1 10 0 10 1 t 10 2 10 3 Tr[H(t) k ] k = 1 k = 2 k = 3 k = 4 10 1 10 0 10 1 t 10 2 3 × 10 1 4 × 10 1 6 × 10 1 Tr[H(t) k ] k = 1 k = 2 k = 3 k = 4

(R), in the limit n 1 , ..., n L-1 → ∞, we have uniformly over [0, T ] Tr H(t) k → 1 N k Tr Θk + Tr S(t) k where Tr (S(t)) → - 1 N (Y * -Y (0)) T Id N n L + e -t Θ Θ-1 ΛT e -t Θ(Y * -Y (0)) + 1 N G(0) T e -t Θ(Y * -Y (0)) Tr S(t) 2 → 1 N 2 (Y * -Y (0)) T e -t Θ Υe -t Θ(Y * -Y (0)) Tr S(t) k →0 when k > 2.

In expectation we have:

E [Tr (S(t))] → - 1 N T r Id N n L + e -t Θ Θ-1 ΛT e -t Θ Σ + Y * Y * T + 1 N T r e -t Θ ΦT E Tr S(t) 2 → 1 N 2 T r e -t Θ Υe -t Θ Σ + Y * Y * T .

Proof. The moments of I are constant because HC = 1 N Id N n L is constant. For the moments of S, we first solve the differential equation for Y (t):

Y (t) = Y * -e -t Θ(Y * -Y (0)). Noting Y (t) -Y (0) = -Θ t 0 ∇C(s)ds, we have G(t) = G(0) -Λ t 0 ∇C(s)ds = G(0) + Λ Θ-1 (Y (t) -Y (0)) = G(0) + Λ Θ-1 Id N n L + e -t Θ (Y * -Y (0))

The expectation of the first moment of S then follows.

3.2. Mutual Orthogonality of I and S:

A first key ingredient to prove Theorem 1 is the asymptotic mutual orthogonality of the matrices I and S Proposition (Proposition 5 in Appendix D). For any loss C with BGOSS and σ ∈ C 4 b (R), we have uniformly over

[0, T ] lim n L-1 →∞ • • • lim n1→∞ IS F = 0. As a consequence lim n L-1 →∞ • • • lim n1→∞ Tr [I + S] k -Tr I k + Tr S k = 0.

Remark 1. If two matrices A and B are mutualy orthogonal (i.e. AB = 0) the range of A is contained in the nullspace of B and vice versa. The non-zero eigenvalues of the sum A + B are therefore given by the union of the non-zero eigenvalues of A and B. Furthermore the moments of A and B add up:

Tr [A + B] k = Tr A k + Tr B k . Proposition 5 shows that

this is what happens asymptotically for I and S.

Note that both matrices I and S have large nullspaces: indeed assuming a constant width

w = n 1 = ... = n L-1 , we have Rank(I) ≤ N n L and Rank(S) ≤ 2(L -1)wN n L (see Appendix C

), while the number of parameters P scales as w 2 (when L > 2).

Figure 2 illustrates the mutual orthogonality of I and S. All numerical experiments are done for rectangular networks (when the width of the hidden layers are equal) and agree well with our predictions obtained in the sequential limit.

3.3. Mean-field Limit:

For a rectangular network with width w, if the output of the network is divided by √ w and the learning rate is multiplied by w (to keep similar dynamics at initialization), the training dynamics changes and the NTK varies during training when w goes to infinity. The new parametrization of the output changes the scaling of the two matrices:

H C 1 √ w Y (L) = 1 w DY (L) T HCDY (L) + 1 √ w ∇C • HY (L) = 1 w I + 1 √ w S.

The scaling of the learning rate essentially multiplies the whole Hessian by w. In this setting, the matrix I is left unchanged while the matrix S is multiplied by √ w (the k-th moment of S is hence multiplied by w k /2 ). In particular, the two moments of the Hessian are dominated by the moments of S, and the higher moments of S (and the operator norm of S) should not vanish. This suggests that the active regime may be characterised by the fact that S F I F . Under the conjecture that Theorem 1 holds for the infinite-width limit of rectangular networks, the asymptotic of the two first moments of H is given by:

1 / √ wTr (H) → N (0, ∇C T Ξ∇C) 1 /wTr H 2 → ∇C T Υ∇C,

where for the MSE loss we have ∇C = -Y * .

3.4. The matrix S:

The matrix S = ∇C • HY (L) is best understood as a perturbation to I, which vanishes as the network converges because ∇C → 0. To calculate its moments, we note that

Tr ∇C • HY (L) = P p=1 ∂ 2 θ 2 p Y T ∇C = G T ∇C, where the vector G = P k=1 ∂ 2 θ 2 p Y ∈ R N n L is the evaluation of the function g θ (x) = P k=1 ∂ 2 θ 2 p f θ (x) on the training set.

For the second moment we have

Tr ∇C • HY (L) 2 = ∇C T   P p,p =1 ∂ 2 θpθ p Y ∂ 2 θpθ p Y T   ∇C = ∇C T Υ∇C for Υ the Gram matrix of the kernel Υ (L) (x, y) = P p,p =1 ∂ 2 θpθ p f θ (x) ∂ 2 θpθ p f θ (y) T .

The following proposition desribes the limit of the function g θ and the kernel Υ (L) and the vanishing of the higher moments: Proposition (Proposition 4 in Appendix C). For any loss C with BGOSS and σ ∈ C 4 b (R), the first two moments of S take the form

Tr (S(t)) = G(t) T ∇C(t) Tr S(t) 2 = ∇C(t) T Υ(t)∇C(t)

-At initialization, g θ and f θ converge to a (centered) Gaussian pair with covariances

E[g θ,k (x)g θ,k (x )] = δ kk Ξ (L) ∞ (x, x ) E[g θ,k (x)f θ,k (x )] = δ kk Φ (L) ∞ (x, x ) E[f θ,k (x)f θ,k (x )] = δ kk Σ (L)

∞ (x, x ) and during training g θ evolves according to

∂ t g θ,k (x) = N i=1 Λ (L) ∞ (x, x i )∂ ik C(Y (t)) •

-Uniformly over any interval [0, T ], the kernel Υ (L) has a deterministic and fixed limit

lim n L-1 →∞ • • • lim n1→∞ Υ (L) kk (x, x ) = δ kk Υ (L)

∞ (x, x ) with limiting kernel:

Υ (L) ∞ (x, x ) = L-1 =1 Θ ( ) ∞ (x, x ) 2 Σ( ) ∞ (x, x ) + 2Θ ( ) ∞ (x, x ) Σ( ) ∞ (x, x ) Σ( +1) ∞ (x, x ) • • • Σ(L-1) ∞ (x, x ).

-The higher moment k > 2 vanish:

lim n L-1 →∞ • • • lim n1→∞ Tr S k = 0.

This result has a number of consequences for infinitely wide networks:

1. At initialization, the matrix S has a finite Frobenius norm S 2 F = Tr S 2 = ∇C T Υ∇C, because Υ converges to a fixed limit. As the network converges, the derivative of the cost goes to zero ∇C(t) → 0 and so does the Frobenius norm of S. 2. In contrast the operator norm of S vanishes already at initialization (because for all even k, we have S op ≤ k Tr (S k ) → 0). At initialization, the vanishing of S in operator norm but not in Frobenius norm can be explained by the matrix S having a growing number of eigenvalues of shrinking intensity as the width grows. 3. When it comes to the first moment of S, Proposition 4 shows that the spectrum of S is in general not symmetric. For the MSE loss the expectation of the first moment at initialization is

E [Tr(S)] = E (Y -Y * ) T G = E Y T G -(Y * ) T E [G] = Tr Φ -0

which may be positive or negative depending on the choice of nonlinearity: with a smooth ReLU, it is positive, while for the arc-tangent or the normalized smooth ReLU, it can be negative (see Figure 1 ). This is in contrast to the result obtained in (Pennington & Bahri, 2017; Geiger et al., 2018) for the shallow ReLU networks, taking the second derivative of the ReLU to be zero. Under this assumption the spectrum of S is symmetric: if the eigenvalues are ordered from lowest to highest, λ i = -λ P -i and Tr(S) = 0. These observations suggest that S has little influence on the shape of the surface, especially towards the end of training, the matrix I however has an interesting structure.

3.5. The matrix I:

At a global minimizer θ * , the spectrum of I describes how the loss behaves around θ * . Along the eigenvectors of the biggest eigenvalues of I, the loss increases rapidely, while small eigenvalues correspond to flat directions. Numerically, it has been observed that the matrix I features a few dominating eigenvalues and a bulk of small eigenvalues (Sagun et al., 2016; 2017; Gur-Ari et al., 2018; Papyan, 2019) . This leads to a narrow valley structure of the loss around a minimum: the biggest eigenvalues are the 'cliffs' of the valley, i.e. the directions along which the loss grows fastest, while the small eigenvalues form the 'flat directions'or the bottom of the valley.

Note that the rank of I is bounded by N n L and in the overparametrized regime, when N n L < P , the matrix I will have a large nullspace, these are directions along which the value of the function on the training set does not change. Note that in the overparametrized regime, global minima are not isolated: they lie in a manifold of dimension at least P -N n L and the nullspace of I is tangent to this solution manifold.

The matrix I is closely related to the NTK Gram matrix:

Θ = DY (L) DY (L) T and I = DY (L) T HCDY (L) .

As a result, the limiting spectrum of the matrix I can be directly obtained from the NTK2 Proposition 1. For any loss C with BGOSS and σ ∈ C 4 b (R), uniformly over any interval [0, T ], the moments Tr I k converge to the following limit (with the convention that i k+1 = i 1 ):

lim n L-1 →∞ • • • lim n1→∞ Tr I k = Tr HC(Y t ) Θ k = 1 N k N i1,...,i k =1 k m=1 c im (f θ(t) (x im ))Θ (L) ∞ (x im , x im+1 ) Proof. It follows from Tr I k = Tr DY (L) T HCDY (L) k = Tr HC Θ k

and the asymptotic of the NTK (Jacot et al., 2018) .

3.5.1. Mean-Square Error:

When the loss is the MSE, HC is equal to 1 N Id N n L . As a result, Θ and I have the same non-zero eigenvalues up to a scaling of 1 /N. Because the NTK is assymptotically fixed, the spectrum of I is also fixed in the limit.

The eigenvectors of the NTK Gram matrix are the kernel principal components of the data. The biggest principal components are the directions in function space which are most favorised by the NTK. This gives a functional interpretation of the narrow valley structure in DNNs: the cliffs of the valley are the biggest principal components, while the flat directions are the smallest components. Remark 2. As the depth L of the network increases, one can observe two regimes (Poole et al., 2016; Jacot et al., 2019) : Order/Freeze where the NTK converges to a constant and Chaos where the NTK converges to a Kronecker delta. In the Order/Freeze the N n L × N n L Gram matrix approaches a block diagonal matrix with n L constant blocks, and as a result n L eigenvalues of I dominate the other ones, corresponding to constant directions along each outputs (this is in line with the observations of (Papyan, 2019)). This leads to a narrow valley for the loss and slows down training. In contrast, in the Chaos regime, the NTK Gram matrix approaches a scaled identity matrix, and the spectrum of I should hence concentrate around a positive value, hence speeding up training. Figure 3 illustrates this phenomenon: with the smooth ReLU we observe a narrow valley, while with the normalized smooth ReLU (which lies in the Chaos according to (Jacot et al., 2019) ) the narrowness of the loss is reduced. A similar phenomenon may explain why normalization helps smoothing the loss surface and speed up training (Santurkar et al., 2018; Ghorbani et al., 2019) .

3.5.2. Cross-Entropy Loss:

For a binary cross-entropy loss with labels

Y * ∈ {-1, +1} N C(Y ) = 1 N N i=1 log 1 + e -Y * i Yi ,

HC is a diagonal matrix whose entries depend on Y (but not on Y * ):

H ii C(Y ) = 1 N 1 1 + e -Yi + e Yi .

The eigenvectors of I then correspond to the weighted kernel principal component of the data. The positive weights 1 1+e -Y i +e Y i approach 1 /3 as Y i goes to 0, i.e. when it is close to the decision boundary from one class to the other, and as Y i → ±∞ the weight go to zero. The weights evolve in time through Y i , the spectrum of I is therefore not asymptotically fixed as in the MSE case, but the functional interpretation of the spectrum in terms of the kernel principal components remains.

4. Conclusion:

We have given an explicit formula for the limiting moments of the Hessian of DNNs throughout training. We have used the common decomposition of the Hessian in two terms I and S and have shown that the two terms are asymptotically mutually orthogonal, such that they can be studied separately.

The matrix S vanishes in Frobenius norm as the network converges and has vanishing operator norm throughout training. The matrix I is arguably the most important as it describes the narrow valley structure of the loss around a global minimum. The eigendecomposition of I is related to the (weighted) kernel principal components of the data w.r.t. the NTK.

A Proofs:

For the proofs of the theorems and propositions presented in the main text, we reformulate the setup of (Jacot et al., 2018) . For a fixed training set x 1 , ..., x N , we consider a (possibly random) time-varying training direction D(t) ∈ R N n L which describes how each of the outputs must be modified. In the case of gradient descent on a cost C(Y ), the training direction is D(t) = ∇C(Y (t)). The parameters are updated according to the differential equation

∂ t θ(t) = (∂ θ Y (t)) T D(t).

Under the condition that T 0 D(t) 2 dt is stochastically bounded as the width of the network goes to infinity, the NTK Θ (L) converges to its fixed limit uniformly over [0, T ].

The reason we consider a general training direction (and not only a gradient of a loss) is that we can split a network in two at a layer and the training of the smaller network will be according to the training direction D ( )

i (t) given by D ( ) i (t) = diag σ α ( ) (x i ) 1 √ n W ( ) T ...diag σ α (L-1) (x i ) 1 √ n L-1 W (L-1) T D i (t)

because the derivatives σ are bounded and by Lemma 1 of the Appendix of (Jacot et al., 2018) , this training direction satisfies the constraints even though it is not the gradient of a loss. As a consequence, as n 1 → ∞, ..., n -1 → ∞ the NTK of the smaller network Θ ( ) also converges to its limit uniformly over [0, T ]. As we let n → ∞ the pre-activations α( ) i and weights W ( ) ij move at a rate of 1 / √ n . We will use this rate of change to prove that other types of kernels are constant during training.

Published as a conference paper at ICLR 2020 When a network is trained with gradient descent on a loss C with BGOSS, the integral T 0 D(t) 2 dt is stochastically bounded. Because the loss is decreasing during training, the outputs Y (t) lie in the sublevel set U C(Y (0)) for all times t. The norm of the gradient is hence bounded for all times t. Because the distribution of Y (0) converges to a multivariate Gaussian, b(C(Y (0))) is stochastically bounded as the width grows, where b(a) is a bound on the norm of the gradient on U a . We then have the bound

T 0 D(t) 2 dt ≤ T b(C(Y (0))) which is itself stochastically bounded.

For the binary and softmax cross-entropy losses the gradient is uniformly bounded: Proposition 2. For the binary cross-entropy loss C and any

Y ∈ R N , ∇C(Y ) 2 ≤ 1 √ N .

For the softmax cross-entropy loss C on c ∈ N classes and any

Y ∈ R N c , ∇C(Y ) 2 ≤ √ 2c √ N .

Proof. The binary cross-entropy loss with labels

Y * ∈ {0, 1} N is C(Y ) = - 1 N N i=1 log e YiY * i 1 + e Yi = 1 N N i=1 log 1 + e Yi -Y i Y * i

and the gradient at an input i is

∂ i C(Y ) = 1 N e Yi -Y * i (1 + e Yi ) 1 + e Yi which is bounded in absolute value by 1 N for both Y * i = 0, 1 such that ∇C(Y ) 2 ≤ 1 √ N .

The softmax cross-entropy loss over c classes with labels Y * ∈ {1, . . . , c} N is defined by

C(Y ) = - 1 N N i=1 log e Y iY * i c k=1 e Y ik = 1 N N i=1 log c k=1 e Y ik -Y iY * i .

The gradient is at an input i and output class m is

∂ im C(Y ) = 1 N e Yim c k=1 e Y ik -δ Y * i m which is bounded in absolute value by 2 N such that ∇C(Y ) 2 ≤ √ 2c √ N .

B Preliminaries:

To study the moments of the matrix S, we first have to show that two tensors vanish as n 1 , ..., n L-1 → ∞:

Ω (L) k0,k1,k2 (x 0 , x 1 , x 2 ) = (∇f θ,k0 (x 0 )) T Hf θ,k1 (x 1 )∇f θ,k2 (x 2 ) Γ (L) k0,k1,k2,k3 (x 0 , x 1 , x 2 , x 4 ) = (∇f θ,k0 (x 0 )) T Hf θ,k1 (x 1 )Hf θ,k2 (x 2 )∇f θ,k3 (x 3 ).

We study these tensors recursively, for this, we need a recursive definition for the first derivatives ∂ θp f θ,k (x) and second derivatives ∂ 2 θpθ p f θ,k (x). The value of these derivatives depend on the layer the parameters θ p and θ p belong to, and on whether they are connection weights W 

∂ W (L-1) mk f θ,k (x) = 1 √ n L-1 α (L-1) m (x)δ kk ∂ b (L-1) k f θ,k (x) = β 2 δ kk

for parameters θ p which belong to the lower layers the derivatives can be defined recursively by

∂ θp f θ,k (x) = 1 √ n L-1 n L-1 m=1 ∂ θp α(L-1) m (x) σ α(L-1) m (x) W (L-1) mk .

For the second derivatives, we first note that if either of the parameters θ p or θ p are bias of the last layer, or if they are both connection weights of the last layer, then ∂ 2 θpθ p f θ,k (x) = 0. Two cases are left: when one parameter is a connection weight of the last layer and the others belong to the lower layers, and when both belong to the lower layers. Both cases can be defined recursively in terms of the first and second derivatives of α(L-1)

m : ∂ 2 θpW (L) mk f θ,k (x) = 1 √ n L-1 ∂ θp α(L-1) m (x) σ α(L-1) m (x) δ kk ∂ 2 θpθ p f θ,k (x) = 1 √ n L-1 n L-1 m=1 ∂ 2 θpθ p α(L-1) m (x) σ α(L-1) m (x) W (L-1) mk + 1 √ n L-1 n L-1 m=1 ∂ θp α(L-1) m (x)∂ θ p α(L-1) m (x)σ α(L-1) m (x) W (L-1) mk

.

Using these recursive definitions, the tensors Ω (L+1) and Γ (L+1) are given in terms of Θ (L) ,Ω (L) and Γ (L) , in the same manner that the NTK Θ (L+1) is defined recursively in terms of Θ (L) in (Jacot et al., 2018) . Lemma 1. For any loss C with BGOSS and σ ∈ C 4 b (R), we have uniformly over [0, T ]

lim n L-1 →∞ • • • lim n1→∞ Ω (L) k0,k1,k2 (x 0 , x 1 , x 2 ) = 0

Proof. The proof is done by induction. When L = 1 the second derivatives

∂ 2 θpθ p f θ,k (x) = 0 and Ω (L) k0,k1,k2 (x 0 , x 1 , x 2 ) = 0.

For the induction step, we write Ω ( +1) k0,k1,k2 (x 0 , x 1 , x 2 ) recursively as

n -3 /2 m0,m1,m2 Θ ( ) m0,m1 (x 0 , x 1 )Θ ( ) m1,m2 (x 1 , x 2 ) σ(α ( ) m0 (x 0 ))σ(α ( ) m1 (x 1 )) σ(α ( ) m2 (x 2 ))W ( ) m0k0 W ( ) m1k1 W ( ) m2k2 + n -3 /2 m0,m1,m2 Ω ( ) m0,m1,m2 (x 0 , x 1 , x 2 ) σ(α ( ) m0 (x 0 )) σ(α ( ) m1 (x 1 )) σ(α ( ) m2 (x 2 ))W ( ) m0k0 W ( ) m1k1 W ( ) m2k2 + n -3 /2 m0,m1 Θ ( ) m0,m1 (x 0 , x 1 ) σ(α ( ) m0 (x 0 )) σ( α( ) m1 (x 1 ))σ( α( ) m1 (x 2 ))W ( ) m0k0 δ k1k2 + n -3 /2 m1,m2 Θ ( ) m1,m2 (x 1 , x 2 )σ(α ( ) m1 (x 0 )) σ(α ( ) m1 (x 1 )) σ( α( ) m2 (x 2 ))δ k0k1 W ( ) m2k2 .

As n 1 , ..., n -1 → ∞ and for any times t < T , the NTK Θ ( ) converges to its limit while Ω ( ) vanishes. The second summand hence vanishes and the others converge to

n -3 /2 m Θ ( ) ∞ (x 0 , x 1 )Θ ( ) ∞ (x 1 , x 2 ) σ(α ( ) m (x 0 ))σ( α( ) m (x 1 )) σ(α ( ) m (x 2 ))W ( ) mk0 W ( ) mk1 W ( ) mk2 + n -3 /2 m Θ ( ) ∞ (x 0 , x 1 ) σ(α ( ) m (x 0 )) σ( α( ) m (x 1 ))σ(α ( ) m (x 2 ))W ( ) mk0 δ k1k2 + n -3 /2 m Θ ( ) ∞ (x 1 , x 2 )σ(α ( ) m (x 0 )) σ( α( ) m (x 1 )) σ(α ( ) m (x 2 ))δ k0k1 W ( ) mk2 .

At initialization, all terms vanish as n → ∞ because all summands are independent with zero mean and finite variance: in the n 1 → ∞, . . . , n -1 → ∞ limit, the α( ) m (x) are independent for different m, see (Jacot et al., 2018) . During training, the weights W ( ) and preactivations α( ) move at a rate of 1 / √ n (see the proof of convergence of the NTK in (Jacot et al., 2018) ). Since σ is Lipschitz, we obtain that the motion during training of each of the sums is of order n -3 /2+ 1 /2 = n -1 . As a result, uniformly over times t ∈ [0, T ], all the sums vanish.

Similarily, we have:

Lemma 2. For any loss C with BGOSS and σ ∈ C 4 b (R), we have uniformly over [0, T ] lim

n L-1 →∞ • • • lim n1→∞ Γ (L) k0,k1,k2,k3 (x 0 , x 1 , x 2 , x 3 ) = 0

Proof. The proof is done by induction. When L = 1 the hessian HF (1) = 0, such that Γ (L) k0,k1,k2,k3 (x 0 , x 1 , x 2 , x 3 ) = 0.

For the induction step, Γ ( +1) can be defined recursively:

Γ (L+1) k0,k1,k2,k3 (x 0 , x 1 , x 2 , x 3 ) = n -2 L m0,m1,m2,m3 Γ (L) m0,m1,m2,m3 (x 0 , x 1 , x 2 , x 3 ) σ(α (L) m0 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 )) W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2,m3 Θ (L) m0,m1 (x 0 , x 1 )Ω (L) m1,m2,m3 (x 1 , x 2 , x 3 ) σ(α (L) m0 (x 0 ))σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 ))W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2,m3 Ω (L) m0,m1,m2 (x 0 , x 1 , x 2 )Θ (L) m2,m3 (x 2 , x 3 ) σ(α (L) m0 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 ))W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2,m3 Θ (L) m0,m1 (x 0 , x 1 )Θ (L) m1,m2 (x 1 , x 2 )Θ (L) m2,m3 (x 2 , x 3 ) σ(α (L) m0 (x 0 ))σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 ))W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m1,m2,m3 Ω (L) m1,m2,m3 (x 1 , x 2 , x 3 )σ(α (L) m1 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 )) δ k0k1 W (L) m2k2 W (L) m3k3 +n -2 L m1,m2,m3 Θ (L) m1,m2 (x 1 , x 2 )Θ (L) m2,m3 (x 2 , x 3 )σ(α (L) m1 (x 0 )) σ(α (L) m1 (x 1 ))σ(α (L) m2 (x 2 )) σ(α (L) m3 (x 3 )) δ k0k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2 Ω (L) m0,m1,m2 (x 0 , x 1 , x 2 ) σ(α (L) m0 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 ))σ(α (L) m2 (x 3 )) W (L) m0k0 W (L) m1k1 δ k2k3 +n -2 L m0,m1,m2 Θ (L) m0,m1 (x 0 , x 1 )Θ (L) m1,m2 (x 1 , x 2 ) σ(α (L) m0 (x 0 ))σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 ))σ(α (L) m2 (x 3 )) W (L) m0k0 W (L) m1k1 δ k2k3 + n -2 L m1,m2 Θ (L) m1,m2 (x 1 , x 2 )σ(α (L) m1 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m2 (x 2 ))σ(α (L) m2 (x 3 ))δ k0k1 δ k2k3 +n -2 L m0,m1,m3 Θ (L) m0,m1 (x 0 , x 1 )Θ (L) m1,m3 (x 2 , x 3 ) σ(α (L) m0 (x 0 )) σ(α (L) m1 (x 1 )) σ(α (L) m1 (x 2 )) σ(α (L) m3 (x 3 )) W (L) m0k0 δ k1k2 W (L) m3k3

As n 1 , ..., n -1 → ∞ and for any times t < T , the NTK Θ ( ) converges to its limit while Ω ( ) and Γ ( ) vanishes. Γ

(L+1) k0,k1,k2,k3 (x 0 , x 1 , x 2 , x 3 ) therefore converges to: +n -2 L m Θ (L) ∞ (x 0 , x 1 )Θ (L) ∞ (x 1 , x 2 )Θ (L) ∞ (x 2 , x 3 ) σ(α (L) m (x 0 ))σ(α (L) m (x 1 ))σ(α (L) m (x 2 )) σ(α (L) m (x 3 )) W (L) mk0 W (L) mk1 W (L) mk2 W (L) mk3 +n -2 L m Θ (L) ∞ (x 1 , x 2 )Θ (L) ∞ (x 2 , x 3 )σ(α (L) m (x 0 )) σ(α (L) m (x 1 ))σ(α (L) m (x 2 )) σ(α (L) m (x 3 )) δ k0k1 W (L) mk2 W (L) mk3 +n -2 L m Θ (L) ∞ (x 0 , x 1 )Θ (L) ∞ (x 1 , x 2 ) σ(α (L) m (x 0 ))σ(α (L) m (x 1 )) σ(α (L) m (x 2 ))σ(α (L) m (x 3 )) W (L) mk0 W (L) mk1 δ k2k3 +n -2 L m Θ (L) ∞ (x 1 , x 2 )σ(α (L) m (x 0 )) σ(α (L) m (x 1 )) σ(α (L) m (x 2 ))σ(α (L) m (x 3 ))δ k0k1 δ k2k3 +n -2 L m Θ (L) ∞ (x 0 , x 1 )Θ (L) ∞ (x 2 , x 3 ) σ(α (L) m (x 0 )) σ(α (L) m (x 1 )) σ(α (L) m (x 2 )) σ(α (L) m (x 3 )) W (L) mk0 δ k1k2 W (L) mk3

For the convergence during training, we proceed similarily to the proof of Lemma 1. At initialization, all terms vanish as n → ∞ because all summands are independent (after taking the n 1 , . . . , n L-1 → ∞ limit) with zero mean and finite variance. During training, the weights W ( ) and preactivations α( ) move at a rate of 1 / √ n which leads to a change of order n -2+ 1 /2 = n -1.5 , which vanishes for all times t too.

C The Matrix S:

We now have the theoretical tools to describe the moments of the matrix S. We first give a bound for the rank of S:

Proposition 3. Rank(S) ≤ 2(n 1 + ... + n L-1 )N n L

Proof. We first observe that S is given by a sum of N n L matrices:

S pp = N i=1 n L k=1 ∂ ik C∂ 2 θpθp f θ,k (x i ).

It is therefore sufficiant to show that the rank of each matrices Hf

θ,k (x) = ∂ 2 θpθ p f θ,k (x i ) p,p is bounded by 2(n 1 + ... + n L ).

The derivatives ∂ θp f θ,k (x) have different definition depending on whether the parameter θ p is a connection weight W ( )

EQUATION

These formulas only depend on θ through the values α

EQUATION

,i for = 1, ..., L-1 (note that both α

(0) i (x) = x i and ∂ α(L) i (x;θ) f θ,k (x) = δ ik do not depend on θ).

Together there are 2(n 1 + ... + n L-1 ) of them. As a consequence, the map θ → ∂ θp f θ,k (x i ) p can be written as a composition

θ ∈ R P → α ( ) i (x; θ), ∂ α( ) i (x;θ) f θ,k (x) ,i ∈ R 2(n1+...+n L-1 ) → ∂ θp f θ,k (x i ) p ∈ R P

and the matrix Hf θ,k (x) is equal to the Jacobian of this map. By the chain rule, Hf θ,k (x) is the matrix multiplication of the Jacobians of the two submaps, whose rank are bounded by 2(n 1 + ... + n L-1 ), hence bounding the rank of Hf θ,k (x). And because S is a sum of N n L matrices of rank smaller than 2(n 1 + ... + n L-1 ), the rank of S is bounded by 2(n 1 + ... + n L-1 )N n L .

C.1 Moments:

Let us now prove Proposition 4: Proposition 4. For any loss C with BGOSS and σ ∈ C 4 b (R), the first two moments of S take the form

Tr (S(t)) = G(t) T ∇C(t) Tr S(t) 2 = ∇C(t) T Υ(t)∇C(t)

-At initialization, g θ and f θ converge to a (centered) Gaussian pair with covariances

E[g θ,k (x)g θ,k (x )] = δ kk Ξ (L) ∞ (x, x ) E[g θ,k (x)f θ,k (x )] = δ kk Φ (L) ∞ (x, x ) E[f θ,k (x)f θ,k (x )] = δ kk Σ (L)

∞ (x, x ) and during training g θ evolves according to

∂ t g θ,k (x) = N i=1 Λ (L) ∞ (x, x i )∂ ik C(Y (t)) • -Uniformly over any interval [0, T ] where T 0 ∇C(t) 2 dt is stochastically bounded, the kernel Υ (L) has a deterministic and fixed limit lim n L-1 →∞ • • • lim n1→∞ Υ (L) kk (x, x ) = δ kk Υ (L) ∞ (x, x )

with limiting kernel:

Υ (L) ∞ (x, x ) = L-1 =1 Θ ( ) ∞ (x, x ) 2 Σ( ) (x, x ) + 2Θ ( ) ∞ (x, x ) Σ( ) (x, x ) Σ( +1) (x, x ) • • • Σ(L-1) (x, x ).

-The higher moment k > 2 vanish:

lim n L-1 →∞ • • • lim n1→∞ Tr S k = 0.

Proof. The first moment of S takes the form

Tr (S) = p (∇C) T H p,p Y = (∇C) T G

where G is the restriction to the training set of the function g θ (x) = p ∂ 2 θpθp f θ (x). This process is random at initialization and varies during training. Lemma 3 below shows that, in the infinite width limit, it is a Gaussian process at initialization which then evolves according to a simple differential equation, hence describing the evolution of the first moment during training.

The second moment of S takes the form:

Tr(S 2 ) = P p1,p2=1 N i1,i2=1 ∂ 2 θp 1 ,θp 2 f θ,k1 (x 1 )∂ 2 θp 2 ,θp 1 f θ,k2 (x 2 )c i1 (x i1 )c i2 (x i2 ) = (∇C) T Υ∇C where Υ (L) k1,k2 (x 1 , x 2 ) = P p1,p2=1 ∂ 2 θp 1 ,θp 2 f θ,k1 (x 1 )∂ 2 θp 2 ,θp 1 f θ,k2 (x 2

) is a multidimensional kernel and Υ is its Gram matrix. Lemma 4 below shows that in the infinite-width limit, Υ (L) k1,k2 (x 1 , x 2 ) converges to a deterministic and time-independent limit Υ (L)

∞ (x 1 , x 2 )δ k1k2 .

To show that Tr(S k ) → 0 for all k > 2, it suffices to show that S 2

F → 0 as Tr(S k ) < S 2 F S k-2 F

and we know that

S F → (∂ Y C) T Υ∂ Y C is finite. We have that S 2 F = N i0,i1,i2,i3=1 n L k0,k1,k2,k3=1 Ψ (L) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 )∂ f θ,k 0 (xi 0 ) C∂ f θ,k 1 (xi 1 ) C ∂ f θ,k 2 (xi 2 ) C∂ f θ,k 3 (xi 3 ) C = Ψ • (∂ Y C) ⊗4 for Ψ the N n L × N n L × N n L × N n L finite version of Ψ (L) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 ) = P p0,p1,p2,p3=1 ∂ 2 θp 0 ,θp 1 f θ,k0 (x 0 )∂ 2 θp 1 ,θp 2 f θ,k1 (x 1 ) ∂ 2 θp 2 ,θp 3 f θ,k2 (x 2 )∂ 2 θp 3 ,θp 0 f θ,k3 (x 3 ).

which vanishes in the infinite width limit by Lemma 5 below.

Lemma 3. For any loss C with BGOSS and σ ∈ C 4 b (R), at initialization g θ and f θ converge to a (centered) Gaussian pair with covariances

E[g θ,k (x)g θ,k (x )] = δ kk Ξ (L) ∞ (x, x ) E[g θ,k (x)f θ,k (x )] = δ kk Φ (L) ∞ (x, x ) E[f θ,k (x)f θ,k (x )] = δ kk Σ (L)

∞ (x, x ) and during training g θ evolves according to

∂ t g θ (x) = N i=1 Λ (L) ∞ (x, x i )D i (t)

Proof. When L = 1, g θ (x) is 0 for any x and θ.

For the inductive step, the trace g , f θ ) follow a centered Gaussian distribution with (conditioned) covariance

E[g (L+1) θ,k (x)g (L+1) θ,k (x )|g (L) θ,m , α(L) m ] = δ kk n L n L m=1 g (L) θ,m (x) σ α(L) m (x) + Θ (L) ∞ (x, x)σ α(L) m (x) g (L) θ,m (x ) σ α(L) m (x ) + Θ (L) ∞ (x , x )σ α(L) m (x ) E[g (L+1) θ,k (x)f θ,k (x )|g (L) θ,m , α(L) m ] = δ kk n L n L m=1 g (L) θ,m (x) σ α(L) m (x) + Θ (L) ∞ (x, x)σ α(L) m (x) σ α(L) m (x ) E[f θ,k (x)f θ,k (x )|g (L) θ,m , α(L) m ] = δ kk n L n L m=1 σ α(L) m (x) σ α(L) m (x ) + β 2 .

As n L → ∞, by the law of large number, these (random) covariances converge to their expectations which are deterministic, hence the pairs (g

(L+1) k

, f θk ) have asymptotically the same Gaussian distribution independent of g

(L) m , α(L) m : E g (L) θ,k (x)g (L) θ,k (x ) → δ kk Ξ (L) ∞ (x, x ) E g (L) θ,k (x)f (L) θ,k (x ) → δ kk Φ (L) ∞ (x, x) E f (L) θ,k (x)f (L) θ,k (x ) → δ kk Σ (L) ∞ (x, x)

with Ξ

(1)

∞ (x, x ) = Φ (1) ∞ (x, x ) = 0 and Ξ (L+1) ∞ (x, x ) = E [gg σ(α) σ(α )] + Θ (L) ∞ (x , x )E [g σ(α)σ(α )] + Θ (L) ∞ (x, x)E [g σ(α )σ(α)] + Θ (L) ∞ (x, x)Θ (L) ∞ (x , x )E [σ(α )σ(α)] = Ξ (L) ∞ (x, x ) Σ(L) ∞ (x, x ) + Φ (L) ∞ (x, x )Φ (L) ∞ (x , x) + Φ (L) ∞ (x, x)Φ (L) ∞ (x , x ) Σ(L) ∞ (x, x ) + Φ (L) ∞ (x, x )Φ (L) ∞ (x , x )E [ σ(α) ... σ (α )] + Φ (L) ∞ (x, x)Φ (L) ∞ (x , x)E [ ... σ (α) σ(α )] + Θ (L) ∞ (x , x ) Φ (L) ∞ (x, x) Σ(L) ∞ (x, x ) + Φ (L) ∞ (x, x )E [ σ(α) ... σ (α )] + Θ (L) ∞ (x, x) Φ (L) ∞ (x , x ) Σ(L) ∞ (x, x ) + Φ (L) ∞ (x , x)E [ ... σ (α) σ(α )] + Θ (L) ∞ (x, x)Θ (L) ∞ (x , x ) Σ(L) ∞ (x, x ) and Φ (L+1) ∞ (x, x ) = E [g σ(α)σ(α )] + Θ (L) ∞ (x, x)E [σ(α)σ(α )] = Φ (L) ∞ (x, x ) Σ(L+1) (x, x ) + Φ (L) ∞ (x, x) + Θ (L) ∞ (x, x) E [σ(α)σ(α )]

where (g, g , α, α ) is a Gaussian quadruple of covariance

     Ξ (L) ∞ (x, x) Ξ (L) ∞ (x, x ) Φ (L) ∞ (x, x) Φ (L) ∞ (x, x ) Ξ (L) ∞ (x, x ) Ξ (L) ∞ (x , x ) Φ (L) ∞ (x , x) Φ (L) ∞ (x , x ) Φ (L) ∞ (x, x) Φ (L) ∞ (x , x) Σ (L) ∞ (x, x) Σ (L) ∞ (x, x ) Φ (L) ∞ (x, x ) Φ (L) ∞ (x , x ) Σ (L) ∞ (x, x ) Σ (L) ∞ (x , x )      .

During training, the parameters follow the gradient ∂ t θ(t) = (∂ θ Y (t))

T D(t). By the induction hypothesis, the traces g

(L)

θ,m then evolve according to the differential equation

∂ t g (L) θ,m (x) = 1 √ n L N i=1 n L m=1 Λ (L) mm (x, x i ) σ(α (L) m (x)) W (L) m T D i (t)

and in the limit as n 1 , ..., n L-1 → ∞, the kernel Λ (L) mm (x, x i ) converges to a deterministic and fixed limit δ mm Λ (L) ∞ (x, x i ). Note that as n L grows, the g

(L)

θ,m (x) move at a rate of 1 / √ n L just like the pre-activations α(L) m . Even though they move less and less, together they affect the trace g (L+1) θ,k which follows the differential equation

∂ t g (L+1) θ,k (x) = N i=1 n L k =1 Λ (L+1) kk (x, x i )D ik (t) where Λ (L+1) kk (x, x ) = 1 n L m,m Λ (L) mm (x, x ) σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L m,m g (L) θ,m (x)Θ (L) mm (x, x )σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L m g (L) θ,m (x) σ α(L) m (x) σ α(L) m (x ) δ kk + 2 n L m,m Ω (L) m mm (x , x, x)σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L m,m Θ (L) mm (x, x)Θ (L) mm (x, x ) ... σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L m Θ (L) mm (x, x)σ α(L) m (x) σ α(L) m (x ) δ kk .

As n 1 , ..., n L-1 → ∞, the kernels Θ (L) mm (x, x ) and Λ (L) mm (x, x ) converge to their limit and Ω (L) m mm (x , x, x) vanishes:

Λ (L) kk (x, x ) → 1 n L m Λ (L) ∞ (x, x ) σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) mk + 1 n L m g (L) θ,m (x)Θ (L) ∞ (x, x )σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) mk + 1 n L m g (L) θ,m (x) σ α(L) m (x) σ α(L) m (x ) δ kk + 1 n L m Θ (L) ∞ (x, x)Θ (L) ∞ (x, x ) ... σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) mk + 1 n L m Θ (L) ∞ (x, x)σ α(L) m (x) σ α(L) m (x ) δ kk

By the law of large numbers, as n L → ∞, at initialization Λ

(L+1) kk (x, x ) → δ kk Λ (L+1) ∞ (x, x ) where Λ (L+1) ∞ (x, x ) = Λ (L) ∞ (x, x ) Σ(L+1) ∞ (x, x ) + Θ (L) ∞ (x, x )E [gσ (α) σ (α )] + E [g σ (α) σ (α )] + Θ (L) ∞ (x, x)Θ (L) ∞ (x, x )E [ ... σ (α) σ (α )] + Θ (L) ∞ (x, x)E [σ (α) σ (α )] = Λ (L) ∞ (x, x ) Σ(L+1) ∞ (x, x ) + Θ (L) ∞ (x, x ) Φ (L) ∞ (x, x ) Σ(L+1) ∞ (x, x ) + Φ (L) ∞ (x, x)E [ ... σ (α) σ (α )] + Φ (L) ∞ (x, x ) Σ(L+1) ∞ (x, x ) + Φ (L) ∞ (x, x)E [σ (α) σ (α )] + Θ (L) ∞ (x, x)Θ (L) ∞ (x, x )E [ ... σ (α) σ (α )] + Θ (L) ∞ (x, x)E [σ (α) σ (α )] During training Θ (L) ∞ and Λ (L)

∞ are fixed in the limit n 1 , .., n L-1 → ∞, and the values g The next lemma describes the asymptotic limit of the kernel Υ (L) : Lemma 4. For any loss C with BGOSS and σ ∈ C 4 b (R), the second moment of the Hessian of the realization function HF (L) converges uniformly over [0, T ] to a fixed limit as n

1 , ...n L-1 → ∞ Υ (L) kk (x, x ) → δ kk L-1 =1 Θ ( ) ∞ (x, x ) 2 Σ( ) ∞ (x, x ) + 2Θ ( ) ∞ (x, x ) Σ( ) ∞ (x, x ) Σ( +1) ∞ (x, x ) • • • Σ(L-1) ∞ (x, x ).

Proof. The proof is by induction on the depth L. The case L = 1 is trivially true because ∂ 2 θpθ p f θ,k (x) = 0 for all p, p , k, x. For the induction step we observe that Υ

(L) k,k (x, x ) = P p1,p2=1 ∂ 2 θp 1 ,θp 2 f θ,k (x)∂ 2 θp 2 ,θp 1 f θ,k (x ) = 1 n L n L m,m =1 Υ (L) m,m (x, x ) σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L n L m,m =1 Ω (L) m ,m,m (x , x, x ) σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L n L m,m =1 Ω (L) m,m ,m (x, x , x)σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 1 n L n L m,m =1 Θ (L) m,m (x, x )Θ (L) m ,m (x , x)σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) m k + 2 n L n L m=1 Θ (L) m,m (x, x ) σ α(L) m (x) σ α(L) m (x ) δ kk

if we now let the width of the lower layers grow to infinity n 1 , ...n L-1 → ∞, the tensor Ω (L) vanishes and Υ m,m converge to limits which are non-zero only when m = m . As a result, the term above converges to

1 n L n L m=1 Υ (L) ∞ (x, x ) σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) mk + 1 n L n L m=1 Θ (L) ∞ (x, x ) 2 σ α(L) m (x) σ α(L) m (x ) W (L) mk W (L) mk + 2 n L n L m=1 Θ (L) ∞ (x, x ) σ α(L) m (x) σ α(L) m (x ) δ kk

At initialization, we can apply the law of large numbers as n L → ∞ such that it converges to Υ

EQUATION

For the convergence during training, we proceed similarily to the proof of Lemma 1: the activations α(L) m (x) and weights W (L) mk move at a rate of 1 / √ n L and the change to Υ (L+1) kk is therefore of order 1 / √ n L and vanishes as n L → 0.

Finally, the next lemma shows the vanishing of the tensor Ψ (L) k0,k1,k2,k3 to prove that the higher moments of S vanish. Lemma 5. For any loss C with BGOSS and σ

∈ C 4 b (R), uniformly over [0, T ] lim n L-1 →∞ • • • lim n1→∞ Ψ (L) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 ) = 0

Proof. When L = 1 the Hessian is zero and Ψ

(1) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 ) = 0.

For the induction step, we write Ψ (L+1) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 ) recursively, because it contains many terms, we change the notation, writing

x 0 x 1 m 0 m 1 for Θ (L) m0,m1 (x 0 , x 1 ), x 0 x 1 x 2 m 0 m 1 m 2 for Ω (L) m0,m1,m2 (x 0 , x 1 , x 2 ) and x 0 x 1 x 2 x 3 m 0 m 1 m 2 m 3 for Γ (L) m0,m1,m2,m3 (x 0 , x 1 , x 2 , x 3 ). The value Ψ (L+1) k0,k1,k2,k3 (x i0 , x i1 , x i2 , x i3 ) is then equal to n -2 L m0,m1,m2,m3 Ψ (L) m0,m1,m2,m3 (x 0 , x 1 , x 2 , x 3 ) σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m2 (x 2 ) σ α(L) m3 (x 3 ) W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2,m3 x 0 x 1 m 0 m 1 x 1 x 2 m 1 m 2 x 2 x 3 m 2 m 3 x 3 x 0 m 3 m 0 σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m2 (x 2 ) σ α(L) m3 (x 3 ) W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3 +n -2 L m0,m1,m2,m3 x 0 x 1 x 2 m 0 m 1 m 2 x 2 x 3 m 2 m 3 x 3 x 0 m 3 m 0 σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m2 (x 2 ) σ α(L) m3 (x 3 ) W (L) m0k0 W (L) m1k1 W (L) m2k2 W (L) m3k3

+n -2 L m0,m1,m2,m3

x 0 x 1 m 0 m 1

x 1 x 2 x 3 m 1 m 2 m 3

x 3 x 0 m 3 m 0 σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m2 (x 2 ) σ α(L) m3 (x 3 ) W 

:

In our experiments, we take β = 0.1.

This result was already obtained in(Karakida et al., 2018), but without identifying the NTK explicitely and only at initialization.

Figures and Tables:

(x, y) and non-symmetric kernel Φ (L) ∞ (x, y). During training, both vectors follow the differential equations

Figure 1: Comparison of the theoretical prediction of Corollary 1 for the expectation of the first 4 moments (colored lines) to the empirical average over 250 trials (black crosses) for a rectangular network with two hidden layers of finite widths n 1 = n 2 = 5000 (L = 3) with the smooth ReLU (left) and the normalized smooth ReLU (right), for the MSE loss on scaled down 14x14 MNIST with N = 256. Only the first two moments are affected by S at the beginning of training.

Figure 2: Illustration of the mutual orthogonality of I and S. For the 20 first eigenvectors of I (blue) and S (orange), we plot the Rayleigh quotients v T Iv and v T Sv (with L = 3, n 1 = n 2 = 1000 and the normalized ReLU on 14x14 MNIST with N = 256). We see that the directions where I is large are directions where S is small and vice versa.

Figure 3: Plot of the loss surface around a global minimum along the first (along the y coordinate) and fourth (x coordinate) eigenvectors of I. The network has L = 4, width n 1 = n 2 = n 3 = 1000 for the smooth ReLU (left) and the normalized smooth ReLU (right). The data is uniform on the unit disk. Normalizing the non-linearity greatly reduces the narrow valley structure of the loss thus speeding up training.

derivatives with respect to the parameters of the last layer are

∇f θ,m (x) (∇f θ,m (x)) T σ α(L) m (x) W (L) mk First note that Tr ∇f θ,m (x) (∇f θ,m (x)) T = Θ (L)mm (x, x). Now let n 1 , ...n L-1 → ∞, by the induction hypothesis, the pairs (g

mk vary at a rate of 1 / √ n L which induce a change of the same rate to Λ (L) kk (x, x ), which is therefore asymptotically fixed during training as n L → ∞.

m2 (x 2 ) σ α(L) m3 (x 3 ) W ,m2,m3 x 1 x 2 m 1 m 2 x 2 x 3 x 0 x 1 m 2 m 3 m 0 m 1 σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m2 (x 2 ) σ α(L) m3 (x 3 ) W m2 (x 2 ) σ α(L) m3 (x 3 ) W 3 x 0 x 1 m 3 m 0 m σ α(L) m0 (x 0 ) σ α(L) m (x 1 ) σ α(L) m (x 2 ) σ α(L) m3 (x 3 ) W x 3 x 0 x 1 x 2 m m 0 m 1 m σ α(L) m0 (x 0 ) σ α(L) m1 (x 1 ) σ α(L) m (x 2 ) σ α(L) m (x 3 )

Peer review (iclr2020):
This paper uses the Neural Tangent Kernel (NTK) to presents an asymptotic analysis of the evolution of Hessian of the loss (w.r.t. model parameters) throughout training. The authors leverage the Neural Tangent Kernel to analyze the evolution of the Hessian of the loss w.r.t the model parameters. Specifically the authors show that as the width of the neural networks tend to infinity, the Hessian can be decomposed into two components: (1) one that reflects the initialization of the model parameters and reduces as training progresses; and (2) one that captures the principal directions of the data. Technical Soundness: The paper appears to be technically sound, though I did not go through the 14 pages of proofs. Potential Impact: I found the focus of the paper to be quite narrow which I think will negatively affect the impact that this paper could have. The authors take the reader on a notational taxing voyage through the decomposition of the Hessian to finally arrive at a result that is rather intuitive if not entirely obvious. While I learned something about the evolution of the Hessian (in the limit of the infinitely wide NN), I can not say that this paper will have significant impact on how I think about NN training. To address this issue of significance and impact: What novel conclusions about Neural Network learning dynamics can you draw from your analysis? What are the implications for generalization or for future training algorithms? Clarity: Beyond the possibly unduly heavy notation, the paper is rather clear and well written. There is a minor typo in the first equation of Sec. 2.3 (i -> j).