The asymptotic spectrum of the Hessian of DNN throughout training

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 under- and 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 infinite-width 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?

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.

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})\approx 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.

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.

Setup

The parameter $\beta$ is added to tune the influence of the bias on trainingIn our experiments, we take $\beta=0.1$.. All parameters are initialized as iid $\mathcal{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):

In this paper, we will study the limit of various objects as $n_{1},\ldots,n_{L}\to\infty$ sequentially, i.e. we first take $n_{1}\to\infty,$ then $n_{2}\to\infty$, 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}\to\infty$. 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.

The network is trained with respect to the cost functional:

For our analysis, we require that the gradient norm $\left\|\mathcal{D}C\right\|$ does not explode during training. The following condition is sufficient:

2 Neural Tangent Kernel

The behavior during training of the network function $f_{\theta}$ in the function space $\mathcal{F}$ is described by a (multi-dimensional) kernel, the Neural Tangent Kernel (NTK)

During training, the function $f_{\theta}$ follows the so-called kernel gradient descent with respect to the NTK, which is defined as

In the infinite-width limit (letting $n_{1}\to\infty,\ldots,n_{L-1}\to\infty$ sequentially) and for losses with BGOSS, the NTK converges to a deterministic limit $\Theta^{(L)}\to\Theta_{\infty}^{(L)}\otimes Id_{n_{L}}$, which is constant during training, uniformly on finite time intervals $\left[0,T\right]$ (Jacot et al., 2018). For the MSE loss, the uniform convergence of the NTK was proven for $T=\infty$ in (Arora et al., 2019).

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).

3 Gram Matrices

It is block diagonal because different outputs $k\neq m$ are asymptotically uncorrelated.

Main Theorems

𝐼𝑆I+S Using the above setup, the Hessian $H$ of the loss $\mathcal{C}\circ F^{(L)}$ is the sum of two terms, with the entry $H_{p,p^{\prime}}$ given by

For a finite dataset, the Hessian matrix $\mathcal{H}\left(C\circ Y^{(L)}\right)$ is equal to the sum of two matrices

where $\mathcal{D}Y^{(L)}$ is a $Nn_{L}\times P$ matrix, $\mathcal{H}C$ is a $Nn_{L}\times Nn_{L}$ matrix and $\mathcal{H}Y^{(L)}$ is a $P\times P\times Nn_{L}$ tensor to which we apply a scalar product (denoted by $\cdot$) in its last dimension with the $Nn_{L}$ vector $\nabla C$ to obtain a $P\times P$ matrix.

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)=\frac{1}{2N}\left\|Y-Y^{*}\right\|^{2}$, the first and second derivatives take simple forms $\nabla C(Y)=\frac{1}{N}\left(Y-Y^{*}\right)$ and $\mathcal{H}C(Y)=\frac{1}{N}Id_{Nn_{L}}$ and the differential equations can be solved to obtain more explicit formulae:

The moments of $I$ are constant because $\mathcal{H}C=\frac{1}{N}Id_{Nn_{L}}$ is constant. For the moments of $S$, we first solve the differential equation for $Y(t)$:

The expectation of the first moment of $S$ then follows. ∎

2 Mutual Orthogonality of I𝐼I and S𝑆S

A first key ingredient to prove Theorem 1 is the asymptotic mutual orthogonality of the matrices $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)\leq Nn_{L}$ and $Rank(S)\leq 2(L-1)wNn_{L}$ (see Appendix C), while the number of parameters $P$ scales as $w^{2}$ (when $L>2$).

Figure 3 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 Mean-field Limit

For a rectangular network with width $w$, if the output of the network is divided by $\sqrt{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:

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 $\sqrt{w}$ (the $k$-th moment of $S$ is hence multiplied by $w^{\nicefrac{{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}\gg\|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:

where for the MSE loss we have $\nabla C=-Y^{*}$.

4 The matrix S𝑆S

The matrix $S=\nabla C\cdot\mathcal{H}Y^{(L)}$ is best understood as a perturbation to $I$, which vanishes as the network converges because $\nabla C\to 0$. To calculate its moments, we note that

The following proposition desribes the limit of the function $g_{\theta}$ and the kernel $\Upsilon^{(L)}$ and the vanishing of the higher moments:

- At initialization, $g_{\theta}$ and $f_{\theta}$ converge to a (centered) Gaussian pair with covariances

and during training $g_{\theta}$ evolves according to

- Uniformly over any interval $[0,T]$, the kernel $\Upsilon^{(L)}$ has a deterministic and fixed limit $\lim_{n_{L-1}\to\infty}\cdots\lim_{n_{1}\to\infty}\Upsilon_{kk^{\prime}}^{(L)}(x,x^{\prime})=\delta_{kk^{\prime}}\Upsilon_{\infty}^{(L)}(x,x^{\prime})$ with limiting kernel:

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

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

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.

5 The matrix I𝐼I

At a global minimizer $\theta^{*}$, the spectrum of $I$ describes how the loss behaves around $\theta^{*}$. 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 $Nn_{L}$ and in the overparametrized regime, when $Nn_{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-Nn_{L}$ and the nullspace of $I$ is tangent to this solution manifold.

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

As a result, the limiting spectrum of the matrix $I$ can be directly obtained from the NTKThis result was already obtained in (Karakida et al., 2018), but without identifying the NTK explicitely and only at initialization.

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.

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 $Nn_{L}\times Nn_{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).

5.2 Cross-Entropy Loss

For a binary cross-entropy loss with labels $Y^{*}\in\{-1,+1\}^{N}$

$\mathcal{H}C$ is a diagonal matrix whose entries depend on $Y$ (but not on $Y^{*}$):

The eigenvectors of $I$ then correspond to the weighted kernel principal component of the data. The positive weights $\frac{1}{1+e^{-Y_{i}}+e^{Y_{i}}}$ approach $\nicefrac{{1}}{{3}}$ as $Y_{i}$ goes to , i.e. when it is close to the decision boundary from one class to the other, and as $Y_{i}\to\pm\infty$ 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.

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.

Acknowledgements

Clément Hongler acknowledges support from the ERC SG CONSTAMIS grant, the NCCR SwissMAP grant, the NSF DMS-1106588 grant, the Minerva Foundation, the Blavatnik Family Foundation, and the Latsis foundation.

References

Appendix A Proofs

Under the condition that $\int_{0}^{T}\left\|D(t)\right\|_{2}dt$ is stochastically bounded as the width of the network goes to infinity, the NTK $\Theta^{(L)}$ converges to its fixed limit uniformly over $[0,T]$.

When a network is trained with gradient descent on a loss $C$ with BGOSS, the integral $\int_{0}^{T}\left\|D(t)\right\|_{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 $\int_{0}^{T}\left\|D(t)\right\|_{2}dt\leq Tb(C(Y(0)))$ which is itself stochastically bounded.

For the binary and softmax cross-entropy losses the gradient is uniformly bounded:

The binary cross-entropy loss with labels $Y^{*}\in\left\{0,1\right\}^{N}$ is

which is bounded in absolute value by $\frac{1}{N}$ for both $Y_{i}^{*}=0,1$ such that $\left\|\nabla C(Y)\right\|_{2}\leq\frac{1}{\sqrt{N}}$.

The softmax cross-entropy loss over $c$ classes with labels $Y^{*}\in\left\{1,\ldots,c\right\}^{N}$ is defined by

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

which is bounded in absolute value by $\frac{2}{N}$ such that $\left\|\nabla C(Y)\right\|_{2}\leq\frac{\sqrt{2c}}{\sqrt{N}}$. ∎

Appendix 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}\to\infty$:

for parameters $\theta_{p}$ which belong to the lower layers the derivatives can be defined recursively by

Using these recursive definitions, the tensors $\Omega^{(L+1)}$ and $\Gamma^{(L+1)}$ are given in terms of $\Theta^{(L)}$,$\Omega^{(L)}$ and $\Gamma^{(L)}$, in the same manner that the NTK $\Theta^{(L+1)}$ is defined recursively in terms of $\Theta^{(L)}$ in (Jacot et al., 2018).

The proof is done by induction. When $L=1$ the second derivatives $\partial_{\theta_{p}\theta_{p^{\prime}}}^{2}f_{\theta,k}(x)=0$ and $\Omega_{k_{0},k_{1},k_{2}}^{(L)}(x_{0},x_{1},x_{2})=0$.

The proof is done by induction. When $L=1$ the hessian $\mathcal{H}F^{(1)}=0$, such that $\Gamma_{k_{0},k_{1},k_{2},k_{3}}^{(L)}(x_{0},x_{1},x_{2},x_{3})=0$.

Appendix C The Matrix S𝑆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$:

We first observe that $S$ is given by a sum of $Nn_{L}$ matrices:

It is therefore sufficiant to show that the rank of each matrices $\mathcal{H}f_{\theta,k}(x)=\left(\partial_{\theta_{p}\theta_{p^{\prime}}}^{2}f_{\theta,k}(x_{i})\right)_{p,p^{\prime}}$ is bounded by $2(n_{1}+...+n_{L})$.

and the matrix $\mathcal{H}f_{\theta,k}(x)$ is equal to the Jacobian of this map. By the chain rule, $\mathcal{H}f_{\theta,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 $\mathcal{H}f_{\theta,k}(x)$. And because $S$ is a sum of $Nn_{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})Nn_{L}$. ∎

- At initialization, $g_{\theta}$ and $f_{\theta}$ converge to a (centered) Gaussian pair with covariances

and during training $g_{\theta}$ evolves according to

- Uniformly over any interval $[0,T]$ where $\int_{0}^{T}\left\|\nabla C(t)\right\|_{2}dt$ is stochastically bounded, the kernel $\Upsilon^{(L)}$ has a deterministic and fixed limit $\lim_{n_{L-1}\to\infty}\cdots\lim_{n_{1}\to\infty}\Upsilon_{kk^{\prime}}^{(L)}(x,x^{\prime})=\delta_{kk^{\prime}}\Upsilon_{\infty}^{(L)}(x,x^{\prime})$ with limiting kernel:

where $G$ is the restriction to the training set of the function $g_{\theta}(x)=\sum_{p}\partial_{\theta_{p}\theta_{p}}^{2}f_{\theta}(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.

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

and during training $g_{\theta}$ evolves according to

When $L=1$, $g_{\theta}(x)$ is for any $x$ and $\theta$.

For the inductive step, the trace $g_{\theta,k}^{(L+1)}(x)$ is defined recursively as

with $\Xi_{\infty}^{(1)}(x,x^{\prime})=\Phi_{\infty}^{(1)}(x,x^{\prime})=0$ and

where $(g,g^{\prime},\alpha,\alpha^{\prime})$ is a Gaussian quadruple of covariance

During training, the parameters follow the gradient $\partial_{t}\theta(t)=\left(\partial_{\theta}Y(t)\right)^{T}D(t)$. By the induction hypothesis, the traces $g_{\theta,m}^{(L)}$ then evolve according to the differential equation

As $n_{1},...,n_{L-1}\to\infty$, the kernels $\Theta_{mm^{\prime}}^{(L)}(x,x^{\prime})$ and $\Lambda_{mm^{\prime}}^{(L)}(x,x^{\prime})$ converge to their limit and $\Omega_{m^{\prime}mm}^{(L)}(x^{\prime},x,x)$ vanishes:

By the law of large numbers, as $n_{L}\to\infty$, at initialization $\Lambda_{kk^{\prime}}^{(L+1)}(x,x^{\prime})\to\delta_{kk^{\prime}}\Lambda_{\infty}^{(L+1)}(x,x^{\prime})$ where

The next lemma describes the asymptotic limit of the kernel $\Upsilon^{(L)}$:

The proof is by induction on the depth $L$. The case $L=1$ is trivially true because $\partial_{\theta_{p}\theta_{p^{\prime}}}^{2}f_{\theta,k}(x)=0$ for all $p,p^{\prime},k,x$. For the induction step we observe that

if we now let the width of the lower layers grow to infinity $n_{1},...n_{L-1}\to\infty$, the tensor $\Omega^{(L)}$ vanishes and $\Upsilon_{m,m^{\prime}}^{(L)}$ and the NTK $\Theta_{m,m^{\prime}}^{(L)}$ converge to limits which are non-zero only when $m=m^{\prime}$. As a result, the term above converges to

At initialization, we can apply the law of large numbers as $n_{L}\to\infty$ such that it converges to $\Upsilon_{\infty}^{(L+1)}(x,x^{\prime})\delta_{kk^{\prime}}$, for the kernel $\Upsilon_{\infty}^{(L+1)}(x,x^{\prime})$ defined recursively by

and $\Upsilon_{\infty}^{(1)}(x,x^{\prime})=0$.

Finally, the next lemma shows the vanishing of the tensor $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(L)}$ to prove that the higher moments of $S$ vanish.

When $L=1$ the Hessian is zero and $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(1)}(x_{i_{0}},x_{i_{1}},x_{i_{2}},x_{i_{3}})=0$.

For the induction step, we write $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(L+1)}(x_{i_{0}},x_{i_{1}},x_{i_{2}},x_{i_{3}})$ recursively, because it contains many terms, we change the notation, writing \left[\begin{array}[]{cc}x_{0}&x_{1}\\ m_{0}&m_{1}\end{array}\right] for $\Theta_{m_{0},m_{1}}^{(L)}(x_{0},x_{1})$, \left[\begin{array}[]{ccc}x_{0}&x_{1}&x_{2}\\ m_{0}&m_{1}&m_{2}\end{array}\right] for $\Omega_{m_{0},m_{1},m_{2}}^{(L)}(x_{0},x_{1},x_{2})$ and \left[\begin{array}[]{cccc}x_{0}&x_{1}&x_{2}&x_{3}\\ m_{0}&m_{1}&m_{2}&m_{3}\end{array}\right] for $\Gamma_{m_{0},m_{1},m_{2},m_{3}}^{(L)}(x_{0},x_{1},x_{2},x_{3})$. The value $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(L+1)}(x_{i_{0}},x_{i_{1}},x_{i_{2}},x_{i_{3}})$ is then equal to

Even though this is a very large formula one can notice that most terms are “rotation of each other”. Moreover, as $n_{1},...,n_{L-1}\to\infty$, all terms containing either an $\Psi^{(L)}$, an $\Omega^{(L)}$ or a $\Gamma^{(L)}$ vanish. For the remaining terms, we may replace the NTKs $\Theta^{(L)}$ by their limit and as a result $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(L+1)}(x_{i_{0}},x_{i_{1}},x_{i_{2}},x_{i_{3}})$ converges to

And all these sums vanish as $n_{L}\to\infty$ thanks to the prefactor $n_{L}^{-2}$, proving the vanishing of $\Psi_{k_{0},k_{1},k_{2},k_{3}}^{(L+1)}(x_{i_{0}},x_{i_{1}},x_{i_{2}},x_{i_{3}})$ in the infinite width limit.

Appendix D Orthogonality of I𝐼I and S𝑆S

From Lemma 2 and the vanishing of the tensor $\Gamma^{(L)}$ as proven in Lemma 2, we can easily prove the orthogonality of $I$ and $S$ of Proposition 5:

and $\Gamma$ vanishes as $n_{1},...,n_{L-1}\to\infty$ by Lemma 2.

which vanishes in the infinite width limit because $\left\|I\right\|_{F}$ and $\left\|S\right\|_{F}$ are bounded and $\left\|A_{m}A_{m+1}\right\|_{F}=\left\|IS\right\|_{F}$ vanishes. ∎