Year: 2019

Abstract:
Are neural networks biased toward simple functions? Does depth always help learn more complex features? Is training the last layer of a network as good as training all layers? How to set the range for learning rate tuning? These questions seem unrelated at face value, but in this work we give all of them a common treatment from the spectral perspective. We will study the spectra of the Conjugate Kernel, CK, (also called the Neural Network-Gaussian Process Kernel), and the Neural Tangent Kernel, NTK. Roughly, the CK and the NTK tell us respectively “what a network looks like at initialization” and “what a network looks like during and after training.” Their spectra then encode valuable information about the initial distribution and the training and generalization properties of neural networks. By analyzing the eigenvalues, we lend novel insights into the questions put forth at the beginning, and we verify these insights by extensive experiments of neural networks. We believe the computational tools we develop here for analyzing the spectra of CK and NTK serve as a solid foundation for future studies of deep neural networks. We have open-sourced the code for it and for generating the plots in this paper at github.com/jxVmnLgedVwv6mNcGCBy/NNspectra.

heading: 1 INTRODUCTION
text: Understanding the behavior of neural networks and why they generalize has been a central pursuit of the theoretical deep learning community. Recently, Valle-Pérez et al. (2018) observed that neural networks have a certain “simplicity bias” and proposed this as a solution to the generalization question. One of the ways with which they argued that this bias exists is the following experiment: they drew a large sample of boolean functions by randomly initializing neural networks and thresholding the output. They observed that there is a bias toward some "simple" functions which get sampled disproportionately more often. However, their experiments were only done for relu networks. Can one expect this “simplicity bias” to hold universally, for any architecture?
A priori, this seems difficult, as the nonlinear nature seems to present an obstacle in reasoning about the distribution of random networks. However, this question turns out to be more easily treated if we allow the width to go to infinity. A long line of works starting with Neal (1995) and extended recently by Lee et al. (2018); Novak et al. (2018); Yang (2019) have shown that randomly initialized, infinite-width networks are distributed as Gaussian processes. These Gaussian processes also describe finite width random networks well (Valle-Pérez et al., 2018). We will refer to the corresponding kernels as the Conjugate Kernels (CK), following the terminology of Daniely et al. (2016). Given the CK K, the simplicity bias of a wide neural network can be read off quickly from the spectrum of K: If the largest eigenvalue of K accounts for most of trK, then a typical random network looks like a function from the top eigenspace of K.
In this paper, we will use this spectral perspective to probe not only the simplicity bias, but more generally, questions regarding how hyperparameters affect the generalization of neural networks.
Via the usual connection between Gaussian processes and linear models with features, the CK can be thought of as the kernel matrix associated to training only the last layer of a wide randomly initialized network. It is a remarkable recent advance (Jacot et al., 2018; Allen-Zhu et al., 2018a;c; Du et al., 2018) that, under a certain regime, a wide neural network of any depth evolves like a linear model even when training all parameters. The associated kernel is call the Neural Tangent Kernel, which is typically different from CK. While its theory was initially derived in the infinite width setting, Lee et al. (2019) confirmed with extensive experiment that this limit is predictive of finite width neural networks as well. Thus, just as the CK reveals information about what a network looks like at
initialization, NTK reveals information about what a network looks like after training. As such, if we can understand how hyperparameters change the NTK, we can also hope to understand how they affect the performance of the corresponding finite-width network.
Our Contributions In this paper, in addition to showing that the simplicity bias is not universal, we will attempt a first step at understanding the effects of the hyperparameters on generalization from a spectral perspective.
At the foundation is a spectral theory of the CK and the NTK on the boolean cube. In Section 3, we show that these kernels, as integral operators on functions over the boolean cube, are diagonalized by the natural Fourier basis, echoing similar results for over the sphere (Smola et al., 2001). We also partially diagonalize the kernels over standard Gaussian, and show that, as expected, the kernels over the different distributions (boolean cube, sphere, standard Gaussian) behave very similarly in high dimensions. However, the spectrum is much easier to compute over the boolean cube: while the sphere and Gaussian eigenvalues would require integration against a kind of polynomials known as the Gegenbauer polynomials, the boolean ones only require calculating a linear combination of a small number of terms. For this reason, in the rest of the paper we focus on analyzing the eigenvalues over the boolean cube.
Just as the usual Fourier basis over R has a notion of frequency that can be interpreted as a measure of complexity, so does the boolean Fourier basis (this is just the degree; see Section 3.1). While not perfect, we adopt this natural notion of complexity in this work; a “simple” function is then one that is well approximated by “low frequencies.”
This spectral perspective immediately yields that the simplicity bias is not universal (Section 4). In particular, while it seems to hold more or less for relu networks, for sigmoidal networks, the simplicity bias can be made arbitrarily weak by changing the weight variance and the depth. In the extreme case, the random function obtained from sampling a deep erf network with large weights is distributed like a “white noise.” However, there is a very weak sense in which the simplicity bias does hold: the eigenvalues of more “complex” eigenspaces cannot be bigger than those of less “complex” eigenspaces (Thm 4.1).
Next, we examine how hyperparameters affect the performance of neural networks through the lens of NTK and its spectrum. To do so, we first need to understand the simpler question of how a kernel affects the accuracy of the function learned by kernel regression. A coarse-grained theory, concerned with big-O asymptotics, exists from classical kernel literature (Yao et al., 2007; Raskutti et al., 2013; Wei et al.; Lin and Rosasco; Schölkopf and Smola, 2002). However, the fine-grained details, required for discerning the effect of hyperparameters, have been much less studied. We make a first attempt at a heuristic, fractional variance (i.e. what fraction of the trace of the kernel does an eigenspace contribute), for understanding how a minute change in kernel effects a change in performance. Intuitively, if an eigenspace has very large fractional variance, so that it accounts for most of the trace, then a ground truth function from this eigenspace should be very easy to learn.
Using this heuristic, we make two predictions about neural networks, motivated by observations in the spectra of NTK and CK, and verify them with extensive experiments.
• Deeper networks learn more complex features, but excess depth can be detrimental as well. Spectrally, depth can increase fractional variance of an eigenspace, but past an optimal depth, it will also decrease it. (Section 5) Thus, deeper is not always better.
• Training all layers is better than training just the last layer when it comes to more complex features, but the opposite is true for simpler features. Spectrally, fractional variances of more “complex” eigenspaces for the NTK are larger than the correponding quantities of the CK. (Section 6)
Finally, we use our spectral theory to predict the maximal nondiverging learning rate (“max learning rate”) of SGD (Section 7).
In general, we will not only verify our theory with experiments on the theoretically interesting distributions, i.e. uniform measures over the boolean cube and the sphere, or the standard Gaussian, but also confirm these findings on real data like MNIST and CIFAR10 1.
1The code for computing the eigenvalues and for reproducing the plots of this paper is available at github. com/jxVmnLgedVwv6mNcGCBy/NNspectra, which will be open sourced upon publication.
For space concerns, we review relevant literature along the flow of the main text, and relegate a more complete discussion of the related research landscape in Appendix A.

heading: 2 KERNELS ASSOCIATED TO NEURAL NETWORKS
text: As mentioned in the introduction, we now know several kernels associated to infinite width, randomly initialized neural networks. The most prominent of these are the neural tangent kernel (NTK) (Jacot et al., 2018) and the conjugate kernel (CK) (Daniely et al., 2016), which is also called the NNGP kernel (Lee et al., 2018). We briefly review them below. First we introduce the following notation that we will repeatedly use. Definition 2.1. For φ : R→ R, write Vφ for the function that takes a PSD (positive semidefinite) kernel function to a PSD kernel of the same domain by the formula
Vφ(K)(x, x ′) = E f∼N (0,K) φ(f(x))φ(f(x′)).
Conjugate Kernel Neural networks are commonly thought of as learning a high-quality embedding of inputs to the latent space represented by the network’s last hidden layer, and then using its final linear layer to read out a classification given the embedding. The conjugate kernel is just the kernel associated to the embedding induced by a random initialization of the neural network. Consider an MLP with widths {nl}l, weight matrices {W l ∈ Rn l×nl−1}l, and biases {bl ∈ Rn l}l, l = 1, . . . , L. For simplicity of exposition, in this paper, we will only consider scalar output nL = 1. Suppose it is parametrized by the NTK parametrization, i.e. its computation is given recursively as
h1(x) = σw√ n0 W 1x+ σbb 1 and hl(x) = σw√ nl−1 W lφ(hl−1(x)) + σbb l (MLP)
with some hyperparameters σw, σb that are fixed throughout training2. At initialization time, suppose W lαβ , b l α ∼ N (0, 1) for each α ∈ [nl], β ∈ [nl−1]. It can be shown that, for each α ∈ [nl], hlα is a Gaussian process with zero mean and kernel function Σl in the limit as all hidden layers become infinitely wide (nl →∞, l = 1, . . . , L− 1), where Σl is defined inductively on l as
Σ1(x, x′) def = σ2w(n 0)−1〈x, x′〉+ σ2b , Σl def = σ2wVφ(Σ l−1) + σ2b (CK)
The kernel ΣL corresponding the the last layer L is the network’s conjugate kernel, and the associated Gaussian process limit is the reason for its alternative name Neural Network-Gaussian process kernel. In short, if we were to train a linear model with features given by the embedding x 7→ hL−1(x) when the network parameters are randomly sampled as above, then the CK is the kernel of this linear model. See Daniely et al. (2016); Lee et al. (2018) and Appendix F for more details.
Neural Tangent Kernel On the other hand, the NTK corresponds to training the entire model instead of just the last layer. Intuitively, if we let θ be the entire set of parameters {W l}l ∪ {bl}l of Eq. (MLP), then for θ close to its initialized value θ0, we expect
hL(x; θ)− hL(x; θ0) ≈ 〈∇θhL(x; θ0), θ − θ0〉
via a naive first-order Taylor expansion. In other words, hL(x; θ)− hL(x; θ0) behaves like a linear model with feature of x given by the gradient taken w.r.t. the initial network, ∇θhL(x; θ0), and the weights of this linear model are the deviation θ− θ0 of θ from its initial value. It turns out that, in the limit as all hidden layer widths tend to infinity, this intuition is correct (Jacot et al., 2018; Lee et al., 2018; Yang, 2019), and the following inductive formula computes the corresponding infinite-width kernel of this linear model:
Θ1 def = Σ1, Θl(x, x′) def = Σl(x, x′) + σ2wΘ l−1(x, x′)Vφ′(Σ l−1)(x, x′). (NTK)
Computing CK and NTK While in general, computing Vφ and Vφ′ requires evaluating a multivariate Gaussian expectation, in specific cases, such as when φ = relu or erf , there exists explicit, efficient formulas that only require pointwise evaluation of some simple functions (see Facts F.1 and F.2). This allows us to evaluate CK and NTK on a set X of inputs in only time O(|X |2L).
2SGD with learning rate α in this parametrization is roughly equivalent to SGD with learning rate α/width in the standard parametrization with Glorot initialization; see Lee et al. (2018)
What Do the Spectra of CK and NTK Tell Us? In summary, the CK governs the distribution of a randomly initialized neural network and also the properties of training only the last layer of a network, while the NTK governs the dynamics of training (all parameters of) a neural network. A study of their spectra thus informs us of the “implicit prior” of a randomly initialized neural network as well as the “implicit bias” of GD in the context of training neural networks.
In regards to the implicit prior at initialization, we know from Lee et al. (2018) that a randomly initialized network as in Eq. (MLP) is distributed as a Gaussian process N (0,K), where K is the corresponding CK, in the infinite-width limit. If we have the eigendecomposition
K = ∑ i≥1 λiui ⊗ ui (1)
with eigenvalues λi in decreasing order and corresponding eigenfunctions ui, then each sample from this GP can be obtained as ∑
i≥1
√ λiωiui, ωi ∼ N (0, 1).
If, for example, λ1 ∑ i≥2 λi, then a typical sample function is just a very small perturbation of u1. We will see that for relu, this is indeed the case (Section 4), and this explains the “simplicity bias” in relu networks found by Valle-Pérez et al. (2018).
Training the last layer of a randomly initialized network via full batch gradient descent for an infinite amount of time corresponds to Gaussian process inference with kernel K (Lee et al., 2018; 2019). A similar intuition holds for NTK: training all parameters of the network (Eq. (MLP)) for an infinite amount of time yields the mean prediction of the GPN (0,NTK) in expectation; see Lee et al. (2019) and Appendix F.4 for more discussion.
Thus, the more the GP prior (governed by the CK or the NTK) is consistent with the ground truth function f∗, the more we expect the Gaussian process inference and GD training to generalize well. We can measure this consistency in the “alignment” between the eigenvalues λi and the squared coefficients a2i of f
∗’s expansion in the {ui}i basis. The former can be interpreted as the expected magnitude (squared) of the ui-component of a sample f ∼ N (0,K), and the latter can be interpreted as the actual magnitude squared of such component of f∗. In this paper, we will investigate an even cleaner setting where f∗ = ui is an eigenfunction. Thus we would hope to use a kernel whose ith eigenvalue λi is as large as possible.
Neural Kernels From the forms of the equation Eqs. (CK) and (NTK) and the fact that Vφ(K)(x, x
′) only depends on K(x, x),K(x, x′), and K(x′, x′), we see that CK or NTK of MLPs takes the form
K(x, y) = Φ ( 〈x, y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) (2)
for some function Φ : R3 → R. We will refer to this kind of kernel as Neural Kernel in this paper.
Kernels as Integral Operators We will consider input spaces of various forms X ⊆ Rd equipped with some probability measure. Then a kernel function K acts as an integral operator on functions f ∈ L2(X ) by
Kf(x) = (Kf)(x) = E y∼X K(x, y)f(y).
We will use the “juxtaposition syntax” Kf to denote this application of the integral operator. 3 Under certain assumptions, it then makes sense to speak of the eigenvalues and eigenfunctions of the integral operator K. While we will appeal to an intuitive understanding of eigenvalues and eigenfunctions in the main text below, we include a more formal discussion of Hilbert-Schmidt operators and their spectral theory in Appendix G for completeness. In the next section, we investigate the eigendecomposition of neural kernels as integral operators over different distributions.
3In cases when X is finite, K can be also thought of as a big matrix and f as a vector — but do not confuse Kf with their multiplication! If we use · to denote matrix multiplication, then the operator application Kf is the same as the matrix multiplication K ·D · f where D is the diagonal matrix encoding the probability values of each point in X .

heading: 3 THE SPECTRA OF NEURAL KERNELS

heading: 3.1 BOOLEAN CUBE
text: We first consider a neural kernelK on the boolean cubeX = ddef= {±1}d, equipped with the uniform measure. In this case, since each x ∈ X has the same norm, K(x, y) = Φ ( 〈x,y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) effectively only depends on 〈x, y〉, so we will treat Φ as a single variate function in this section, Φ(c) = Φ(c, 1, 1).
Brief review of basic Fourier analysis on the boolean cube d (O’Donnell (2014)). The space of real functions on d forms a 2d-dimensional space. Any such function has a unique expansion into a multilinear polynomial (polynomials whose monomials do not contain xpi , p ≥ 2, of any variable xi). For example, the majority function over 3 bits has the following unique multilinear expansion
maj3 : 3 → 1, maj3(x1, x2, x3) =
1 2 (x1 + x2 + x3 − x1x2x3).
In the language of Fourier analysis, the 2d multilinear monomial functions
χS(x) def = xS def = ∏ i∈S xi, for each S ⊆ [d] (3)
form a Fourier basis of the function space L2( d) = {f : d → R}, in the sense that their inner products satisfy
E x∼ d χS(x)χT (x) = I(S = T ).
Thus, any function f : d → R can be always written as f(x) = ∑ S⊆[d] f̂(S)χX(x)
for a unique set of coefficients {f̂(S)}S⊆[d].
It turns out that K is always diagonalized by this Fourier basis {χS}S⊆[d]. Theorem 3.1. On the d-dimensional boolean cube d, for every S ⊆ [d], χS is an eigenfunction of K with eigenvalue
µ|S| def = E x∈ d xSK(x,1) = E x∈ d xSΦ (∑ i xi/d ) , (4)
where 1 = (1, . . . , 1) ∈ d. This definition of µ|S| does not depend on the choice S, only on the cardinality of S. These are all of the eigenfunctions of K by dimensionality considerations.4
Define T∆ to be the shift operator on functions over [−1, 1] that sends Φ(·) to Φ(· −∆). Then we can re-express the eigenvalue as follows. Lemma 3.2. With µk as in Thm 3.1,
µk = 2 −d(I − T∆)k(I + T∆)d−kΦ(1) (5)
= 2−d d∑ r=0 Cd−k,kr Φ (( d 2 − r ) ∆ ) (6)
where
Cd−k,kr def = ∑ j=0 (−1)r+j ( d− k j )( k r − j ) . (7)
Eq. (5) will be important for computational purposes, and we will come back to discuss this more in Section 3.5. It also turns out µk affords a pretty expression via the Fourier series coefficients of Φ. As this is not essential to the main text, we relegate its exposition to Appendix H.1.
4Readers familiar with boolean Fourier analysis may be reminded of the noise operator Tρ, ρ ≤ 1 (O’Donnell, 2014, Defn 2.46). In the language of this work, Tρ is a neural kernel with eigenvalues µk = ρk.

heading: 3.2 SPHERE
text: Now let’s consider the case when X = √ dSd−1 is the radius√ d sphere in Rd equipped with the uniform measure. Again, because x ∈ X all have the same norm, we will treat Φ as a univariate function with K(x, y) = Φ(〈x, y〉/‖x‖‖y‖) = Φ(〈x, y〉/d). As is long known (Schoenberg, 1942; Gneiting, 2013; Xu and Cheney, 1992; Smola et al., 2001), K is diagonalized by spherical harmonics, and the eigenvalues are given by the coefficients of Φ against a system of orthogonal polynomials called Gegenbuaer polynomials. We relegate a complete review of this topic to Appendix H.2.

heading: 3.3 ISOTROPIC GAUSSIAN
text: Now let’s consider X = Rd equipped with standard isotropic Gaussian N (0, I), so that K behaves like
Kf(x) = E y∼N (0,I) K(x, y)f(y) = E y∼N (0,I) Φ ( 〈x, y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) f(y)
for any f ∈ L2(N (0, I)). In contrast to the previous two sections, K will essentially depend on the effect of the norms ‖x‖ and ‖y‖ on Φ. Nevertheless, because an isotropic Gaussian vector can be obtained by sampling its direction uniformly from the sphere and its magnitude from a chi distribution, K can still be partially diagonalized into a sum of products between spherical harmonics and kernels on R equipped with a chi distribution (Thm H.14). In certain cases, we can obtain complete eigendecompositions, for example when Φ is positive homogeneous. See Appendix H.3 for more details.
3.4 KERNEL IS SAME OVER BOOLEAN CUBE, SPHERE, OR GAUSSIAN WHEN d 1
The reason we have curtailed a detailed discussion of neural kernels on the sphere and on the standard Gaussian is because, in high dimension, the kernel behaves the same under these distributions as under uniform distribution over the boolean cube. Indeed, by intuition along the lines of the central limit theorem, we expect that uniform distribution over a high dimension boolean cube should approximate high dimensional standard Gaussian. Similarly, by concentration of measure, most of the mass of a Gaussian is concentrated around a thin shell of radius √ d. Thus, morally, we expect the same kernel function K induces approximately the same integral operator on these three distributions in high dimension, and as such, their eigenvalues should also approximately coincide. We verify empirically and theoretically this is indeed the case in Appendix H.4.

heading: 3.5 COMPUTING THE EIGENVALUES
text: As the eigenvalues of K over the different distributions are very close, we will focus in the rest of this paper on eigenvalues over the boolean cube. This has the additional benefit of being much easier to compute.
Each eigenvalue over the sphere and the standard Gaussian requires an integration of Φ against a Gegenbauer polynomial. In high dimension d, these Gegenbauer polynomials varies wildly in a sinusoidal fashion, and blows up toward the boundary (see Fig. 15 in the Appendix). As such, it is difficult to obtain a numerically stable estimate of this integral in an efficient manner when d is large.
In contrast, we have multiple ways of computing boolean cube eigenvalues, via Eqs. (5) and (6). In either case, we just take some linear combination of the values of Φ at a grid of points on [−1, 1], spaced apart by ∆ = 2/d. While the coefficients Cd−k,kr (defined in Eq. (7)) are relatively efficient to compute, the change in the sign of Cd−k,kr makes this procedure numerically unstable for large d. Instead, we use Eq. (5) to isolate the alternating part to evaluate in a numerically stable way: Since µk = ( I+T∆
2 )d−k ( I−T∆ 2 )k Φ(1), we can evaluate Φ̃ def= ( I−T∆ 2 )k Φ via k finite differences, and then
compute
( I + T∆
2
)d−k Φ̃(1) = 1
2d−k d−k∑ r=0 ( d− k r ) Φ̃(1− r∆). (8)
When Φ arises from the CK or the NTK of an MLP, all derivatives of Φ at 0 are nonnegative (Thm I.3). Thus intuitively, the finite difference Φ̃ should be also all nonnegative, and this sum can be evaluated without worry about floating point errors from cancellation of large terms.
A slightly more clever way to improve the numerical stability when 2k ≤ d is to note that (I + T∆)d−k (I − T∆)k Φ(1) = (I + T∆)d−2k ( I − T 2∆ )k Φ(1) = (I + T∆)d−2k (I − T2∆)k Φ(1).
So an improved algorithm is to first compute the kth finite difference (I − T2∆)k with the larger step size 2∆, then compute the sum (I + T∆)d−2k as in Eq. (8).

heading: 4 CLARIFYING THE “SIMPLICITY BIAS” OF RANDOM NEURAL NETWORKS
text: As mentioned in the introduction, Valle-Pérez et al. (2018) claims that neural networks are biased toward simple functions. We show that this phenomenon depends crucially on the nonlinearity, the sampling variances, and the depth of the network. In Fig. 1(a), we have repeated their experiment for 104 random functions obtained by sampling relu neural networks with 2 hidden layers, 40 neurons each, following Valle-Pérez et al. (2018)’s architectural choices5. We also do the same for erf networks of the same depth and width, varying as well the sampling variances of the weights and biases, as shown in the legend. As discussed in Valle-Pérez et al. (2018), for relu, there is indeed this bias, where a single function gets sampled more than 10% of the time. However, for erf, as we increase σ2w, we see this bias disappear, and every function in the sample gets sampled only once.
This phenomenon can be explained by looking at the eigendecomposition of the CK, which is the Gaussian process kernel of the distribution of the random networks as their hidden widths tend to infinity. In Fig. 1(b), we plot the normalized eigenvalues {µk/ ∑7 i=0 ( 7 i ) µi}7k=0 for the CKs corresponding to the networks sampled in Fig. 1(a). Immediately, we see that for relu and σ2w = σ 2 b = 2, the degree 0 eigenspace, corresponding to constant functions, accounts for more than 80% of the variance. This means that a typical infinite-width relu network of 2 layers is expected to be almost constant, and this should be even more true after we threshold the network to be a boolean function. On the other hand, for erf and σb = 0, the even degree µks all vanish, and most of the variance comes from degree 1 components (i.e. linear functions). This concentration in degree 1 also lessens as σ2w increases. But because this variance is spread across a dimension 7 eigenspace, we don’t see duplicate function samples nearly as much as in the relu case. As σw increases, we also see the eigenvalues become more equally distributed, which corresponds to the flattening of
5Valle-Pérez et al. (2018) actually performed their experiments over the {0, 1}7 cube, not the {±1}7 cube we are using here. This does not affect our conclusion. See Appendix J for more discussion
the probability-vs-rank curve in Fig. 1(a). Finally, we observe that a 32-layer erf network with σ2w = 4 has all its nonzero eigenvalues (associated to odd degrees) all equal (see points marked by ∗ in Fig. 1(b)). This means that its distribution is a "white noise" on the space of odd functions, and the distribution of boolean functions obtained by thresholding the Gaussian process samples is the uniform distribution on odd functions. This is the complete lack of simplicity bias modulo the oddness constraint.
However, from the spectral perspective, there is a weak sense in which a simplicity bias holds for all neural network-induced CKs and NTKs. Theorem 4.1 (Weak Spectral Simplicity Bias). Let K be the CK or NTK of an MLP on a boolean cube d. Then the eigenvalues µk, k = 0, . . . , d, satisfy
µ0 ≥ µ2 ≥ · · · ≥ µ2k ≥ · · · , µ1 ≥ µ3 ≥ · · · ≥ µ2k+1 ≥ · · · . (9)
Even though it’s not true that the fraction of variance contributed by the degree k eigenspace is decreasing with k, the eigenvalue themselves will be in a nonincreasing pattern across even and odd degrees. In fact, if we fix k and let d→∞, then we can show that (Thm I.6)
µk = Θ(d −k).
Of course, as we have seen, this is a very weak sense of simplicity bias, as it doesn’t prevent “white noise” behavior as in the case of erf CK with large σ2w and large depth.

heading: 5 DEEPER NETWORKS LEARN MORE COMPLEX FEATURES
text: In the rest of this work, we compute the eigenvalues µk over the 128-dimensional boolean cube ( d, with d = 128) for a large number of different hyperparameters, and analyze how the latter affect the former. We vary the degree k ∈ [0, 8], the nonlinearity between relu and erf, the depth (number of hidden layers) from 1 to 128, and σ2b ∈ [0, 4]. We fix σ2w = 2 for relu kernels, but additionally vary σ2w ∈ [1, 5] for erf kernels. Comprehensive contour plots of how these hyperparameters affect the kernels are included in Appendix D, but in the main text we summarize several trends we see.
We will primarily measure the change in the spectrum by the degree k fractional variance, which is just
degree k fractional variance def=
( d k ) µk∑d
i=0 ( d i ) µi .
This terminology comes from the fact that, if we were to sample a function f from a Gaussian process with kernel K, then we expect that r% of the total variance of f comes from degree k components of f , where r% is the degree k fractional variance. If we were to try to learn a homogeneous degree-k
polynomial using a kernel K, intuitively we should try to choose K such that its µk is maximized, relative to other eigenvalues. Fig. 3(a) shows that this is indeed the case even with neural networks: over a large number of different hyperparameter settings, degree k fractional variance is inversely related to the validation loss incurred when learning a degree k polynomial. However, this plot also shows that there does not seem like a precise, clean relationship between fractional variance and validation loss. Obtaining a better measure for predicting generalization is left for future work.
Before we continue, we remark that the fractional variance of a fixed degree k converges to a fixed value as the input dimension d→∞: Theorem 5.1 (Asymptotic Fractional Variance). Let K be the CK or NTK of an MLP on a boolean cube d. ThenK can be expressed asK(x, y) = Φ(〈x, y〉/d) for some analytic function Φ : R→ R. If we fix k and let the input dimension d→∞, then the fractional variance of degree k converges to
(k!)−1Φ(k)(0)/Φ(1) = (k!)−1Φ(k)(0)∑ j≥0(j!) −1Φ(j)(0)
where Φ(k) denotes the kth derivative of Φ.
For the fractional variances we compute in this paper, their values at d = 128 are already very close to their d→∞ limit, so we focus on the d = 128 case experimentally. If K were to be the CK or NTK of a relu or erf MLP, then we find that for higher k, the depth of the network helps increase the degree k fractional variance. In Fig. 2(a) and (b), we plot, for each degree k, the depth that (with some combination of other hyperparameters like σ2b ) achieves this maximum, for respectively relu and erf kernels. Clearly, the maximizing depths are increasing with k for relu, and also for erf when considering either odd k or even k only. The slightly differing behavior between even and odd k is expected, as seen in the form of Thm 4.1. Note the different scales of y-axes for relu and erf — the depth effect is much stronger for erf than relu.
For relu NTK and CK, σ2b = 0 maximizes fractional variance in general, and the same holds for erf NTK and CK in the odd degrees (see Appendix D). In Fig. 2(c) and Fig. 2(d) we give a more fine-grained look at the σ2b = 0 slice, via heatmaps of fractional variance against degree and depth.
Brighter color indicates higher variance, and we see the optimal depth for each degree k clearly increases with k for relu NTK, and likewise for odd degrees of erf NTK. However, note that as k increases, the difference between the maximal fractional variance and those slightly suboptimal becomes smaller and smaller, reflected by suppressed range of color moving to the right. The heatmaps for relu and erf CKs look similar and are omitted.
We verify this increase of optimal depth with degree in Fig. 3(b). There we have trained relu networks of varying depth against a ground truth multilinear polynomial of varying degree. We see clearly that the optimal depth is increasing with degree. We also verify this phenomenon when the input distribution changes to the standard Gaussian or the uniform distribution over the sphere √ dSd−1; see Fig. 6.
Note that implicit in our results here is a highly nontrivial observation: Past some point (the optimal depth), high depth can be detrimental to the performance of the network, beyond just the difficulty to train, and this detriment can already be seen in the corresponding NTK or CK. In particular, it’s not true that the optimal depth is infinite. We confirm the existence of such an optimal depth even in real distributions like MNIST and CIFAR10; see Fig. 7. This adds significant nuance to the folk wisdom that “depth increases expressivity and allows neural networks to learn more complex features.”

heading: 6 NTK FAVORS MORE COMPLEX FEATURES THAN CK
text: We generally find the degree k fractional variance of NTK to be higher than that of CK when k is large, and vice versa when k is small, as shown in Fig. 4. This means that, if we train only the last layer of a neural network (i.e. CK dynamics), we intuitively should expect to learn simpler features better, while, if we train all parameters of the network (i.e. NTK dynamics), we should expect to learn more complex features better. Similarly, if we were to sample a function from a Gaussian process with the CK as kernel (recall this is just the distribution of randomly initialized infinite width MLPs (Lee et al., 2018)), this function is more likely to be accurately approximated by low degree polynomials than the same with the NTK.
We verify this intuition by training a large number of neural networks against ground truth functions of various homogeneous polynomials of different degrees, and show a scatterplot of how training the last layer only measures against training all layers (Fig. 3(c)). This phenomenon remains true over the standard Gaussian or the uniform distribution on the sphere (Fig. 8). Consistent with our theory, the only place training the last layer works meaningfully better than training all layers is when the ground truth is a constant function. However, we reiterate that fractional variance is an imperfect indicator of performance. Even though for erf neural networks and k ≥ 1, degree k fractional variance of NTK is not always greater than that of the CK, we do not see any instance where training the last layer of an erf network is better than training all layers. We leave an investigation of this discrepancy to future work.

heading: 7 PREDICTING THE MAXIMUM LEARNING RATE
text: In any setup that tries to push deep learning benchmarks, learning rate tuning is a painful but indispensable part. In this section, we show that our spectral theory can accurately predict the maximal nondiverging learning rate over real datasets as well as toy input distributions, which would help set the correct upper limit for a learning rate search.
By Jacot et al. (2018), in the limit of large width and infinite data, the function g : X → R represented by our neural network evolves like
gt+1 = gt − 2αK(gt − g∗), t = 0, 1, 2, . . . , (10) when trained under full batch GD (with the entire population) with L2 loss L(f, g) = Ex∼X (f(x)− g(x))2, ground truth g∗, and learning rate α, starting from randomly initialization. If we train only the last layer, then K is the CK; if we train all layers, then K is the NTK. Given an eigendecomposition of K as in Eq. (1), if g0 − g∗ = ∑ i aiui is the decomposition of g
0 in the eigenbasis {ui}i, then one can easily deduce that gt − g∗ = ∑ i ai(1− 2αλi)tui.
Consequently, we must have α < (maxi λi) −1 in order for Eq. (10) to converge 6
When the input distribution is the uniform distribution over d, the maximum learning rate is max(µ0, µ1) by Thm 4.1. By Thm 5.1, as long as the Φ function corresonding to K has Φ(0) 6= 0, when d is large, we expect µ0 ≈ Φ(0) but µ1 ∼ d−1Φ′(0) µ0. Therefore, we should predict 1Φ(0) for the maximal learning rate when training on the boolean cube. However, as Fig. 5 shows, this prediction is accurate not only for the boolean cube, but also over the sphere, the standard Gaussian, and even MNIST and CIFAR10!

heading: 8 CONCLUSION
text: In this work, we have taken a first step at studying how hyperparameters change the initial distribution and the generalization properties of neural networks through the lens of neural kernels and their spectra. We obtained interesting insights by computing kernel eigenvalues over the boolean cube and relating them to generalization through the fractional variance heuristic. While it inspired valid predictions that are backed up by experiments, fractional variance is clearly just a rough indicator. We hope future work can refine on this idea to produce a much more precise prediction of test loss. Nevertheless, we believe the spectral perspective is the right line of research that will not only shed light on mysteries in deep learning but also inform design choices in practice.

heading: A RELATED WORKS
text: The Gaussian process behavior of neural networks was found by Neal (1995) for shallow networks and then extended over the years to different settings and architectures (Williams, 1997; Le Roux and Bengio, 2007; Hazan and Jaakkola, 2015; Daniely et al., 2016; Lee et al., 2018; Matthews et al., 2018; Novak et al., 2018). This connection was exploited implicitly or explicitly to build new models (Cho and Saul, 2009; Lawrence and Moore, 2007; Damianou and Lawrence, 2013; Wilson et al., 2016a;b; Bradshaw et al., 2017; van der Wilk et al., 2017; Kumar et al., 2018; Blomqvist et al., 2018; Borovykh, 2018; Garriga-Alonso et al., 2018; Novak et al., 2018; Lee et al., 2018). The Neural Tangent Kernel is a much more recent discovery by Jacot et al. (2018) and later Allen-Zhu et al. (2018a;c;b); Du et al. (2018); Arora et al. (2019b); Zou et al. (2018) came upon the same reasoning independently. Like CK, NTK has also been applied toward building new models or algorithms (Arora et al., 2019a; Achiam et al., 2019).
Closely related to the discussion of CK and NTK is the signal propagation literature, which tries to understand how to prevent pathological behaviors in randomly initialized neural networks when they are deep (Poole et al., 2016; Schoenholz et al., 2017; Yang and Schoenholz, 2017; 2018; Hanin, 2018; Hanin and Rolnick, 2018; Chen et al., 2018; Yang et al., 2019; Pennington et al., 2017a; Hayou et al., 2018; Philipp and Carbonell, 2018). This line of work can trace its original at least to the advent of the Glorot and He initialization schemes for deep networks (Glorot and Bengio, 2010; He et al., 2015). The investigation of forward signal propagation, or how random neural networks change with depth, corresponds to studying the infinite-depth limit of CK, and the investigation of backward signal propagation, or how gradients of random networks change with depth, corresponds to studying the infinite-depth limit of NTK. Some of the quite remarkable results from this literature includes how to train a 10,000 layer CNN (Xiao et al., 2017) and that, counterintuitively, batch normalization causes gradient explosion (Yang et al., 2019).
This signal propagation perspective can be refined via random matrix theory (Pennington et al., 2017a; 2018). In these works, free probability is leveraged to compute the singular value distribution of the input-output map given by the random neural network, as the input dimension and width tend to infinity together. Other works also investigate various questions of neural network training and generalization from the random matrix perspective (Pennington and Worah, 2017; Pennington and Bahri, 2017; Pennington and Worah, 2018).
Yang (2019) presents a common framework, known as Tensor Programs, unifying the GP, NTK, signal propagation, and random matrix perspectives, as well as extending them to new scenarios, like recurrent neural networks. It proves the existence of and allows the computation of a large number of infinite-width limits (including ones relevant to the above perspectives) by expressing the quantity of interest as the output of a computation graph and then manipulating the graph mechanically.
Several other works also adopt a spectral perspective on neural networks (Candès, 1999; Sonoda and Murata, 2017; Eldan and Shamir, 2016; Barron, 1993; Xu et al., 2018; Zhang et al., 2019; Xu et al., 2019; Xu, 2018); here we highlight a few most relevant to us. Rahaman et al. (2018) studies the real Fourier frequencies of relu networks and perform experiments on real data as well as synthetic ones. They convincingly show that relu networks learn low frequencies components first. They also investigate the subtleties when the data manifold is low dimensional and embedded in various ways in the input space. In contrast, our work focuses on the spectra of the CK and NTK (which indirectly informs the Fourier frequencies of a typical network). Nevertheless, our results are complementary to theirs, as they readily explain the low frequency bias in relu that they found. Karakida et al. (2018) studies the spectrum of the Fisher information matrix, which share the nonzero eigenvalues with the NTK. They compute the mean, variance, and maximum of the eigenvalues Fisher eigenvalues (taking the width to infinity first, and then considering finite amount of data sampled iid from a Gaussian). In comparison, our spectral results yield all eigenvalues of the NTK (and thus also all nonzero eigenvalues of the Fisher) as well as eigenfunctions.
Finally, we note that several recent works (Xie et al., 2016; Bietti and Mairal, 2019; Basri et al., 2019; Ghorbani et al., 2019) studied one-hidden layer neural networks over the sphere, building on Smola et al. (2001)’s observation that spherical harmonics diagonalize dot product kernels, with the latter two concurrent to us. This is in contrast to the focus on boolean cube here, which allows us to study the fine-grained effect of hyperparameters on the spectra, leading to a variety of insights into neural networks’ generalization properties.

heading: B UNIVERSALITY OF OUR BOOLEAN CUBE OBSERVATIONS IN OTHER INPUT DISTRIBUTIONS
text: Using the spectral theory we developed in this paper, we made three observations, that can be roughly summarized as follows: 1) the simplicity bias noted by Valle-Pérez et al. (2018) is not universal; 2) for each function of fixed “complexity” there is an optimal depth such that networks shallower or deeper will not learn it as well; 3) training last layer only is better than training all layers when learning “simpler” features, and the opposite is true for learning “complex” features.
In this section, we discuss the applicability of these observations to distributions that are not uniform over the boolean cube: in particular, the uniform distribution over the sphere √ dSd−1, the standard Gaussian N (0, Id), as well as realistic data distributions such as MNIST and CIFAR10.
Simplicity bias The simplicity bias noted by Valle-Pérez et al. (2018), in particular Fig. 1, depends on the finiteness of the boolean cube as a domain, so we cannot effectively test this on the distributions above, which all have uncountable support.
Optimal depth With regard to the second observation, we can test whether an optimal depth exists for learning functions over the distributions above. Since polynomial degrees remain the natural indicator of complexity for the sphere and the Gaussian (see Appendices H.2 and H.3 for the relevant spectral theory), we replicated the experiment in Fig. 3(b) for these distributions, using the same ground truth functions of polynomials of different degrees. The results are shown in Fig. 6. We see the same phenomenon as in the boolean cube case, with an optimal depth for each degree, and with the optimal depth increasing with degree.
For MNIST and CIFAR10, the notion of “feature complexity” is less clear, so we will not test the hypothesis that “optimal depth increases with degree” for these distributions but only test for the existence of the optimal depth for the ground truth marked by the labels of the datasets. We do so by training a large number of MLPs of varying depth on these datasets until convergence, and plot the results in Fig. 7. This figure clearly shows that such an optimal depth exists, such that shallower or deeper networks do monotonically worse as the depth diverge away from this optimal depth.
Again, the existence of the optimal depth is not obvious at all, as conventional deep learning wisdom would have one believe that adding depth should always help.
Training last layer only vs training all layers Finally, we repeat the experiment in Fig. 3(c) for the sphere and the standard Gaussian, with polynomials of different degrees as ground truth functions. The results are shown in Fig. 8. We see the same phenomenon as in the boolean cube case: for degree 0 polynomials, training last layer works better in general, but for higher degree polynomials, training all layers fares better.
Note that, unlike the sphere and the Gaussian, whose spectral theory tells us that (harmonic) polynomial degree is a natural notion of complexity, for MNIST and CIFAR10 we have much less clear idea of what a “complex” or a “simple” feature is. Therefore, we did not attempt a similar experiment on these datasets.

heading: C THEORETICAL VS EMPIRICAL MAX LEARNING RATES UNDER DIFFERENT PREPROCESSING FOR MNIST AND CIFAR10
text: In the main text Fig. 5, on the MNIST and CIFAR10 datasets, we preprocessed the data by centering and normalizing to the sphere (see Appendix E.2 for a precise description). With this preprocessing, our theory accurately predicts the max learning rate in practice.
In general, if we go by another preprocessing, such as PCA or ZCA, or no preprocessing, our theoretical max learning rate 1/Φ(0) is less accurate but still correlated in general. The only exception seems to be relu networks on PCA- or ZCA- preprocessed CIFAR10. See Fig. 9.
Theoretical vs empirical max learning rate under different preprocessing
D VISUALIZING THE SPECTRAL EFFECTS OF σ2w, σ 2 b , AND DEPTH
While in the main text, we summarized several trends of interest kn several plots, they do not give the entire picture of how eigenvalues and fractional variances vary with σ2w, σ 2 b , and depth. Here we try to present this relationship more completely in a series of contour plots. Fig. 10 shows how varying depth and σ2b changes the fractional variances of each degree, for relu CK and NTK. We are fixing σ2w = 2 in the CK plots, as the fractional variances only depend on the ratio σ 2 b/σ 2 w; even though this is not true for relu NTK, we fix σ2w = 2 as well for consistency. For erf, however, the fractional variance will crucially depend on both σ2w and σ 2 b , so we present 3D contour plots of how σ 2 w, σ 2 b , and depth changes fractional variance in Fig. 13. Complementarily, we also show in Figs. 11 and 12 a few slices of these 3D contour plots for different fixed values of σ2b , for erf CK and NTK.

heading: E EXPERIMENTAL DETAILS
text: E.1 FIG. 3
Fig. 3(a), (b) and (c) differ in the set of hyperparameters they involve (to be specified below), but in all of them, we train relu networks against a randomly generated ground truth multilinear polynomial, with input space 128 and L2 loss L(f) = Ex∈ d(f(x)− f∗(x))2.
Training We perform SGD with batch size 1000. In each iteration, we freshly sample a new batch, and we train for a total of 100,000 iterations, so the network potentially sees 108 different examples. At every 1000 iterations, we validate the current network on a freshly drawn batch of 10,000 examples. We thus record a total of 100 validation losses, and we take the lowest to be the “best validation loss.”
Generating the Ground Truth Function The ground truth function f∗(x) is generated by first sampling 10 monomials m1, . . . ,m10 of degree k, then randomly sampling 10 coefficients a1, . . . , a10 for them. The final function is obtained by normalizing {ai} such that the sum of their squares is 1:
f∗(x) def = 10∑ i=1 aimi/ 10∑ j=1 a2j . (11)

heading: Hyperparameters for Fig. 3(a)
text: • The learning rate is half the theoretical maximum learning rate7 12 max(µ0, µ1) −1
• Ground truth degree k ∈ {0, 1, 2, 3} • Depth ∈ {0, . . . , 10} • activation = relu • σ2w = 2 • σ2b = 0 • width = 1000 • 10 random seeds per hyperparameter combination • training last layer (marked “ck”), or all layers (marked “ntk”). In the latter case, we use the
NTK parametrization of the MLP (Eq. (MLP)).

heading: Hyperparameters for Fig. 3(b)
text: • The learning rate is half the theoretical maximum learning rate 12 max(µ0, µ1) −1
• Ground truth degree k ∈ {0, 1, 2, 3} • Depth ∈ {0, . . . , 10} • activation = relu • σ2w = 2 • σ2b = 0 • width = 1000 • 100 random seeds per hyperparameter combination • training last layer weight and bias only
7Note that, because the L2 loss here is L(f) = Ex∈ d(f(x) − f∗(x))2, the maximum learning rate is λ−1max = max(µ0, µ1) −1 (see Thm 4.1). If we instead adopt the convention L(f) = Ex∈ d 12 (f(x)− f ∗(x))2, then the maximum learning rate would be 2λ−1max = 2max(µ0, µ1)−1
Algorithm 1 Binary Search for Empirical Max Learning Rate upper ← 16× theoretical max lr lower ← 0 tol← 0.01× theoretical max lr while |upper − lower| > tol do
α← (upper + lower)/2 Run SGD with learning rate α for 1000 iterations if loss diverges then
upper ← α else
lower ← α end if
end while Output: upper

heading: Hyperparameters for Fig. 3(c)
text: • The learning rate ∈ {0.05, 0.1, 0.5} • Ground truth degree k ∈ {0, 1, . . . , 6} • Depth ∈ {1, . . . , 5} • activation ∈ {relu, erf} • σ2w = 2 for relu, but σ2w ∈ {1, 2, . . . , 5} for erf • σ2b ∈ {0, 1, . . . , 4} • width = 1000 • 1 random seed per hyperparameter combination • Training all layers, using the NTK parametrization of the MLP (Eq. (MLP))

heading: E.2 MAX LEARNING RATE EXPERIMENTS
text: Here we describe the experimental details for the experiments underlying Figs. 5 and 9.
Theoretical max learning rate For a fixed setup, we compute Φ according to Eq. (CK) (if only last layer is trained) or Eq. (NTK) (if all layers are trained). For ground truth problems where the output is n-dimensional, the theoretical max learning rate is nΦ(0)−1; in particular, the max learning rates for MNIST and CIFAR10 are 10 times those for boolean cube, sphere, and Gaussian. This is because the kernel for an multi-output problem effectively becomes
1 n K⊕n = 1 n K 0 0 0
. . . 0 0 0 K  where the 1n factor is due to the 1 n factor in the scaled square loss L(f, f ∗) = Ex∼X 1n ∑n i=1(f(x)i− f∗(x)i) 2. The top eigenvalue for 1nK ⊕n is just 1n times the top eigenvalue for K.
Empirical max learning rate For a fixed setup, we perform binary search for the empirical max learning rate as in Algorithm 1.
Preprocessing In Fig. 5, for MNIST and CIFAR10, we center and project each image onto the sphere √ dSd−1, where d = 28× 28 = 784 for MNIST and d = 3× 32× 32 = 3072 for CIFAR10. More precisely, we compute the average image x̄ over the entire dataset, and we preprocess each image x as √ d x−x̄‖x−x̄‖ . In Fig. 9, there are three different preprocessing schemes. For “no preprocessing,” we load the MNIST and CIFAR10 data as is. In “PCA128,” we take the top 128 eigencomponents of the data, so that the data has only 128 dimensions. In “ZCA128,” we take the top 128 eigencomponents but rotate it back to the original space, so that the data still has dimension d, where d = 28×28 = 784 for MNIST and d = 3× 32× 32 = 3072 for CIFAR10.

heading: Hyperparameters
text: • Target function: For boolean cube, sphere, and standard Gaussian, we randomly sample a degree 1 polynomial as in Eq. (11). For MNIST and CIFAR10, we just use the label in the dataset, encoded as a one-hot vector for square-loss regression. • Depth ∈ {1, 2, 4, 8, 16} • activation ∈ {relu, erf} • σ2w = 2 for relu, but σ2w ∈ {1, 2, . . . , 5} for erf • σ2b ∈ {1, . . . , 4} • width = 1000 • 1 random seed per hyperparameter combination • Training last layer (CK) or all layers (NTK). In the latter case, we use the NTK parametriza-
tion of the MLP (Eq. (MLP)).

heading: F REVIEW OF THE THEORY OF NEURAL TANGENT KERNELS

heading: F.1 CONVERGENCE OF INFINITE-WIDTH KERNELS AT INITIALIZATION
text: Conjugate Kernel Via a central-limit-like intuition, each unit hl(x)α of Eq. (MLP) should behave like a Gaussian as width nl−1 →∞, as it is a sum of a large number of roughly independent random variables (Poole et al., 2016; Schoenholz et al., 2017; Yang and Schoenholz, 2017). The devil, of course, is in what “roughly independent” means and how to apply the central limit theorem (CLT) to this setting. It can be done, however, (Lee et al., 2018; Matthews et al., 2018; Novak et al., 2018), and in the most general case, using a “Gaussian conditioning” technique, this result can be rigorously generalized to almost any architecture Yang (2019). In any case, the consequence is that, for any finite set S ⊆ X ,
{hlα(x)}x∈S converges in distribution to N (0,Σl(S, S)),
as min{n1, . . . , nl−1} → ∞, where Σl is the CK as given in Eq. (CK).
Neural Tangent Kernel By a slightly more involved version of the “Gaussian conditioning” technique, Yang (2019) also showed that, for any x, y ∈ X ,
〈∇θhL(x),∇θhL(y)〉 converges almost surely to ΘL(x, y)
as the widths tend to infinity, where Θl is the NTK as given in Eq. (NTK).

heading: F.2 FAST EVALUATIONS OF CK AND NTK
text: For certain φ like relu or erf, Vφ and V′φ can be evaluated very quickly, so that both the CK and NTK can be computed in O(|X |2L) time, where X is the set of points we want to compute the kernel function over, and L is the number of layers. Fact F.1 (Cho and Saul (2009)). For any kernel K
Vrelu(K)(x, x ′) =
1 2π ( √ 1− c2 + (π − arccos c)c) √ K(x, x)K(x′, x′)
V′relu(K)(x, x ′) =
1
2π (π − arccos c)
where c = K(x, x′)/ √ K(x, x)K(x′, x′).
Fact F.2 (Neal (1995)). For any kernel K,
Verf(K)(x, x ′) =
2 π arcsin K(x, x′)√ (K(x, x) + 0.5)(K(x′, x′) + 0.5)
V′erf(K)(x, x ′) =
4 π √ (1 + 2K(x, x))(1 + 2K(x′, x′))− 4K(x, x′)2 .
Fact F.3. Let φ(x) = exp(x/σ) for some σ > 0. For any kernel K,
Vφ(K)(x, x ′) = exp
( K(x, x) + 2K(x, x′) +K(x′, x′)
2σ2
) .

heading: F.3 LINEAR EVOLUTION OF NEURAL NETWORK UNDER GD
text: Remarkably, the NTK governs the evolution of the neural network function under gradient descent in the infinite-width limit. First, let’s consider how the parameters θ and the neural network function f evolve under continuous time gradient flow. Suppose f is only defined on a finite input space X = {x1, . . . , xk}. We will visualize
f(X ) = f(x1)
... f(xk) , ∇fL =
∂L ∂f(x1)
... ∂L
∂f(xk)
, θ = θ1 ... θn , ∇θf =
∂f(x1) ∂θ1 · · · ∂f(x k)
∂θ1 ... . . . ...
∂f(x1) ∂θn · · · ∂f(x k) ∂θn
(best viewed in color). Then under continuous time gradient descent with learning rate η,
∂t θt = −η∇θL(ft) = −η ∇θft · ∇fL(ft) ,
∂t ft = ∇θft > · ∂t θt = −η ∇θft > · ∇θft · ∇fL(ft) = −η Θt · ∇fL(ft) (12)
where Θt = ∇θf>t · ∇θft ∈ Rk×k is of course the (finite width) NTK. These equations can be visualized as
∂t = −η · , ∂t = · ∂t = −η · · = −η ·
Thus f undergoes kernel gradient descent with (functional) loss L(f) and kernel Θt. This kernel Θt of course changes as f evolves, but remarkably, it in fact stays constant for f being an infinitely wide MLP (Jacot et al., 2018):
∂tft = −ηΘ · ∇fL(ft), (Training All Layers)
where Θ is the infinite-width NTK corresponding to f . A similar equation holds for the CK Σ if we train only the last layer,
∂tft = −ηΣ · ∇fL(ft). (Training Last Layer)
If L is the square loss against a ground truth function f∗, then ∇fL(ft) = 12k∇f‖ft − f ∗‖2 = 1 k (ft−f
∗), and the equations above become linear differential equations. However, typically we only have a training set X train ⊆ X of size far less than |X |. In this case, the loss function is effectively
L(f) = 1 2|X train| ∑ x∈X train (f(x)− f∗(x))2,
with functional gradient
∇fL(f) = 1
|X train| Dtrain · (f − f∗),
where Dtrain is a diagonal matrix of size k × k whose diagonal is 1 on x ∈ X train and 0 else. Then our function still evolves linearly
∂tft = −η(K ·Dtrain) · (ft − f∗) (13)
where K is the CK or the NTK depending on which parameters are trained.

heading: F.4 RELATIONSHIP TO GAUSSIAN PROCESS INFERENCE.
text: Recall that the initial f0 in Eq. (13) is distributed as a Gaussian process N (0,Σ) in the infinite width limit. As Eq. (13) is a linear differential equation, the distribution of ft will remain a Gaussian process for all t, whether K is CK or NTK. Under suitable conditions, it can be shown that (Lee et al., 2019), in the limit as t→∞, if we train only the last layer, then the resulting function f∞ is distributed as a Gaussian process with mean f̄∞ given by
f̄∞(x) = Σ(x,X train)Σ(X train,X train)−1f∗(X train)
and kernel Var f∞ given by
Var f∞(x, x ′) = Σ(x, x′)− Σ(x,X train)Σ(X train,X train)−1Σ(X train, x′).
These formulas precisely described the posterior distribution of f given prior N (0,Σ) and data {(x, f∗(x))}x∈X train . If we train all layers, then similarly as t→∞, the function f∞ is distributed as a Gaussian process with mean f̄∞ given by (Lee et al., 2019)
f̄∞(x) = Θ(x,X train)Θ(X train,X train)−1f∗(X train).
This is, again, the mean of the Gaussian process posterior given prior N (0,Θ) and the training data {(x, f∗(x))}x∈X train . However, the kernel of f∞ is no longer the kernel of this posterior, but rather is an expression involving both the NTK Θ and the CK Σ; see Lee et al. (2019).
In any case, we can make the following informal statement in the limit of large width
Training the last layer (resp. all layers) of an MLP infinitely long, in expectation, yields the mean prediction of the GP inference given prior N (0,Σ) (resp. N (0,Θ)).

heading: G A BRIEF REVIEW OF HILBERT-SCHMIDT OPERATORS AND THEIR SPECTRAL THEORY
text: In this section, we briefly review the theory of Hilbert-Schmidt kernels, and more importantly, to properly define the notion of eigenvalues and eigenfunctions. A function K : X 2 → R is called a Hilbert-Schmidt operator if K ∈ L2(X × X ), i.e.
‖K‖2HS def = E x,y∼X K(x, y)2 <∞.
‖K‖2HS is known as the Hilbert-Schmidt norm of K. K is called symmetric if K(x, y) = K(y, x) and positive definite (resp. semidefinite) if
E x,y∼X f(x)K(x, y)f(y) > 0 (resp. ≥ 0) for all f ∈ L2(X ) not a.e. zero.
A spectral theorem (Mercer’s theorem) holds for Hilbert-Schmidt operators.
Fact G.1. If K is a symmetric positive semidefinite Hilbert-Schmidt kernel, then there is a sequence of scalars λi ≥ 0 (eigenvalues) and functions fi ∈ L2(X ) (eigenfunctions), for i ∈ N, such that
∀i, j, 〈fi, fj〉 = I(i = j), and K(x, y) = ∑ i∈N λifi(x)fi(y)
where the convergence is in L2(X × X ) norm.
This theorem allows us to speak of the eigenfunctions and eigenvalues, which are important for training and generalization considerations when K is a kernel used in machine learning, as discussed in the main text.
A sufficient condition for K to be a Hilbert-Schmidt kernel in our case (concerning only probability measure on X ) is just that K is bounded. All Ks in this paper satisfy this property.

heading: H EIGENDECOMPOSITION OF NEURAL KERNEL ON DIFFERENT DOMAINS

heading: H.1 BOOLEAN CUBE
text: From the Fourier Series Perspective. We continue from the discussion of the boolean cube in the main text. Recall that T∆ is the shift operator on functions that sends Φ(·) to Φ(· −∆). Notice that, if we let Φ(t) = eκt for some κ ∈ C, then T∆Φ(s) = e−κ∆ · eκt. Thus Φ is an “eigenfunction” of the operator T∆ with eigenvalue e−κ∆. In particular, this implies that Proposition H.1. Suppose Φ(t) = et/σ 2
, as in the case whenK is the CK or NTK of a 1-layer neural network with nonlinearity exp(·/σ), up to multiplicative constant (Fact F.3). Then the eigenvalue µk over the boolean cube d equals
µk = 2 −d(1− exp(−∆/σ2))k(1 + exp(−∆/σ2))d−k · exp(1/σ2)
where ∆ = 2/d.
It would be nice if we can express any Φ as a linear combination of exponentials, so that Eq. (5) simplifies in the fashion of Prop H.1 — this is precisely the idea of Fourier series.
We will use the theory of Fourier analysis on the circle, and for this we need to discuss periodic functions. Let Φ̃ : [−2, 2]→ R be defined as
Φ̃(x) =  Φ(x) if x ∈ [−1, 1] Φ(2− x) if x ∈ [1, 2] Φ(−2− x) if x ∈ [−2,−1].
See Fig. 14 for an example illustration. Note that if Φ is continuous on [−1, 1], then Φ̃ is continuous as a periodic function on [−2, 2].
The Fourier basis on functions over [−2, 2] is the collection {t 7→ e 12πist}s∈Z. Under generic conditions (for example if Ψ ∈ L2[−2, 2]), a function Ψ has an associated Fourier series ∑ s∈Z Ψ̂(s)e 1 2πist. We briefly review basic facts of Fourier analysis on the circle. Recall the following notion of functions of bounded variation. Definition H.2. A function f : [a, b]→ R is said to have bounded variation if
sup P nP−1∑ i=0 |f(xi+1)− f(xi)| <∞,
where the supremum is taken over all partitions P of the interval [a, b], P = {x0, . . . , xnP }, x0 ≤ x1 ≤ · · · ≤ xnP .
Intuitively, a function of bounded variation has a graph (in [a, b]× R) of finite length. Fact H.3 (Katznelson (2004)). A bounded variation function f : [−2, 2]→ R that is periodic (i.e. f(−2) = f(2)) has a pointwise-convergent Fourier series:
lim T→∞ ∑ s∈[−T,T ] Ψ̂(s)e 1 2πist → Ψ(t), ∀t ∈ [−2, 2].
From this fact easily follows the following lemma. Lemma H.4. Suppose Φ is continuous and has bounded variation on [−1, 1]. Then Φ̃ is also continuous and has bounded variation, and its Fourier Series (on [−2, 2]) converges pointwise to Φ̃.
Proof. Φ̃ is obviously continuous and has bounded variation as well, and from Fact H.3, we know a periodic continuous function with bounded variation has a pointwise-convergent Fourier Series.
Certainly, T∆ sends continuous bounded variation functions to continuous bounded variation functions. Because T∆e 1 2πist = e− 1 2πis∆e 1 2πist,
T∆ ∑ s∈Z Ψ̂(s)e 1 2πist = ∑ s∈Z Ψ̂(s)e− 1 2πis∆e 1 2πist
whenever both sides are well defined. If Ψ is continuous and has bounded variation then T∆Ψ is also continuous and has bounded variation, and thus its Fourier series, the RHS above, converges pointwise to T∆Ψ. Now, observe
(I − T∆)k(I + T∆)d−kΦ̃(x) = d∑ r=0 Cd−k,kr Φ̃ (x− r∆)
(I − T∆)k(I + T∆)d−kΦ̃(1) = d∑ r=0 Cd−k,kr Φ (( d 2 − r ) ∆ ) = µk
Expressing the LHS in Fourier basis, we obtain Theorem H.5.
µk = ∑ s∈Z is(1− e− 12πis∆)k(1 + e− 12πis∆)d−k ˆ̃Φ(s)
where ˆ̃Φ(s) = 1
4 ∫ 2 −2 Φ̃(t)e− 1 2πist dt
= 1
4 ∫ 1 −1 Φ(t)(e− 1 2πist + (−1)se 12πist) dt
=
{ 1 2 ∫ 1 −1 Φ(t) cos( 1 2πst) dt if s is even
− i2 ∫ 1 −1 Φ(t) sin( 1 2πst) dt if s is odd
denote the Fourier coefficients of Φ̃ on [−2, 2]. (Here i is the imaginary unit here, not an index).
Recovering the values of Φ given the eigenvalues µ0, . . . , µd. Conversely, given eigenvalues µ0, . . . , µd corresponding to each monomial degree, we can recover the entries of the matrix K. Theorem H.6. For any x, y ∈ d with Hamming distance r,
K(x, y) = Φ
(( d 2 − r ) ∆ ) = d∑ k=0 Cd−r,rk µk,
where Cd−r,rk = ∑ j=0(−1)k+j ( d−r j )( r k−j ) as in Eq. (7).
Proof. Recall that for any S ⊆ [d], χS(x) = xS is the Fourier basis corresponding to S (see Eq. (3)). Then by converting from the Fourier basis to the regular basis, we get
Φ
(( d 2 − r ) ∆ ) = K(x, y) for any x, y ∈ d with Hamming distance r
= d∑ k=0 µk ∑ |S|=k χS(x)χS(y).
If x and y differ on a set T ⊆ [d], then we can simplify the inner sum
Φ
(( d 2 − r ) ∆ ) = d∑ k=0 µk ∑ |S|=k (−1)|S∩T | = d∑ k=0 µkC d−r,r k .
Remark H.7. If we let T be the operator that sends µ• 7→ µ•+1, then we have the following operator expression
Φ
(( d 2 − r ) ∆ ) = [(1 + T )d−r(1− T )rµ]0
Remark H.8. The above shows that the matrix C = {Cd−r,rk }dk,r=0 satisfies
C2 = 2dI.

heading: H.2 SPHERE
text: Now let’s consider the case when X = √ dSd−1 is the radius√ d sphere in Rd equipped with the uniform measure. Again, because x ∈ X all have the same norm, we will consider Φ as a univariate function with K(x, y) = Φ(〈x, y〉/‖x‖‖y‖) = Φ(〈x, y〉/d). As is long known (Schoenberg, 1942; Gneiting, 2013; Xu and Cheney, 1992; Smola et al., 2001), K is diagonalized by spherical harmonics. We review these results briefly below, as we will build on them to deduce spectral information of K on isotropic Gaussian distributions.
Review: spherical harmonics and Gegenbauer polynomials. Spherical harmonics are L2 functions on Sd−1 that are eigenfunctions of the Laplace-Beltrami operator ∆Sd−1 of Sd−1. They can be described as the restriction of certain homogeneous polynomials in Rd to Sd−1. Denote byHd−1,(l) the space of spherical harmonics of degree l on sphere Sd−1. Then we have the orthogonal decomposition L2(Sd−1) ∼= ⊕∞ l=0Hd−1,(l). It is a standard fact that dimHd−1,(l) = ( d−1+l d−1 ) − ( d−3+l d−1 ) .
There is a special class of spherical harmonics called zonal harmonics that can be represented as x 7→ p(〈x, y〉) for specific polynomials p : R→ R, and that possess a special reproducing property which we will describe shortly. Intuitively, the value of any zonal harmonics only depends on the “height” of x along some fixed axis y, so a typical zonal harmonics looks like Fig. 16. The polynomials pmust be one of the Gegenbauer polynomials. Gegenbauer polynomials {C(α)l (t)}∞l=0 are orthogonal polynomials with respect to the measure (1− t2)α− 12 on [−1, 1] (see Fig. 15 for examples), and here we adopt the convention that∫ 1
−1 C(α)n (t)C (α) l (t)(1− t 2)α− 1 2 dt =
π21−2αΓ(n+ 2α) n!(n+ α)[Γ(α)]2 I(n = l). (14)
Then for each (oriented) axis y ∈ Sd−1 and degree l, there is a unique zonal harmonic
Zd−1,(l)y ∈ Hd−1,(l), Zd−1,(l)y (x) def = c−1d,lC
( d−22 ) l (〈x, y〉)
for any x, y ∈ Sd−1, where cd,l = d−2d+2l−2 . Very importantly, they satisfy the following
Fact H.9 (Reproducing property (Suetin)). For any f ∈ Hd−1,(m),
E z∼Sd−1 Zd−1,(l)y (z)f(z) = f(y)I(l = m)
E z∼Sd−1 Zd−1,(l)y (z)Z d−1,(m) x (z) = Z d−1,(l) y (x)I(l = m) = c −1 d,lC
( d−22 ) l (〈x, y〉)I(l = m)
We also record a useful fact about Gegenbauer polynomials. Fact H.10 (Suetin).
C (α) l (±1) = (±1) l
( l + 2α− 1
l ) By a result of Schoenberg (1942), we have the following eigendecomposition of K on the sphere.
Theorem H.11 (Schoenberg). Suppose Φ : [−1, 1]→ R is in L2((1− t2) d−12 −1), so that it has the Gegenbauer expansion
Φ(t) a.e. = ∞∑ l=0 alc −1 d,lC ( d−22 ) l (t).
Then K has eigenspaces Hd−1,(l)√ d def = {f(x/
√ d) : f ∈ Hd−1,(l)} with corresponding eigenval-
ues al. Since ⊕∞
l=0H d−1,(l)√ d
is an orthogonal decomposition of L2( √ dSd−1), this describes all
eigenfunctions of K considered as an operator on L2( √ dSd−1).
For completeness, we include the proof of this theorem in Appendix I.
By Bezubik et al. (2008), we can express the Gegenbauer coefficients, and thus equivalently the eigenvalues, via derivatives of Φ: Theorem H.12 (Bezubik et al. (2008)). If the Taylor expansion of Φ at 0,
Φ(t) = ∞∑ n=0 Φ(n)(0) n! tn,
is absolutely convergent on the closed interval [−1, 1], then the Gegenbauer coefficients al in Thm H.11 in dimension d is equal to the absolute convergent series
al = Γ
( d
2 ) ∞∑ k=0
Φ(l+2k)(0) 2l+2kk!Γ ( d 2 + l + k ) . (15)
As the dimension d of the sphere tends to∞, the eigenvalues in fact simplify to the derivatives of Φ: Theorem H.13. Let K be the CK or NTK of an MLP on the sphere √ dSd−1. Then K can be expressed as K(x, y) = Φ(〈x, y〉/d) for some smooth Φ : [−1, 1] → R. Let a` denote K’s eigenvalues on the sphere (as in Thm H.11). If we fix ` and let d→∞, then
lim d→∞
d`a` = Φ (`)(0),
where Φ(`) denotes the `th derivative of Φ.
This theorem is the same as Thm I.6 except that it concerns the sphere rather than the boolean cube.
Proof. By Thm I.3, Φ’s Taylor expansion around 0 is absolutely convergent on [−1, 1], so that the condition of Thm H.12 is satisfied. Therefore, Eq. (15) holds and is absolutely convergent. By dominated convergence theorem, we can exchange the limit and the summation, and get
lim d→∞ dlal = lim d→∞ dlΓ
( d
2 ) ∞∑ k=0
Φ(l+2k)(0) 2l+2kk!Γ ( d 2 + l + k ) =
∞∑ k=0 Φ(l+2k)(0) lim d→∞
dlΓ ( d 2 ) 2l+2kk!Γ ( d 2 + l + k
) =
∞∑ k=0 Φ(l+2k)(0) lim d→∞ ( d 2 )−k (k!)−12−2k
= Φ(l)(0)
as desired.
H.3 ISOTROPIC GAUSSIAN
Now let’s consider X = Rd equipped with standard isotropic Gaussian N (0, I), so that K behaves like
Kf(x) = E y∼N (0,I) K(x, y)f(y) = E y∼N (0,I) Φ ( 〈x, y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) f(y)
for any f ∈ L2(N (0, I)). In contrast to the previous two sections, K will essentially depend on the effect of the norms ‖x‖ and ‖y‖ on Φ. Note that an isotropic Gaussian vector z ∼ N (0, I) can be sampled by independently sampling its direction v uniformly from the sphere Sd−1 and sampling its magnitude r from a chi distribution χd with d degrees of freedom. Proceeding along this line of logic yields the following spectral theorem: Theorem H.14. A function K : (Rd)2 → R of the form
K(x, y) = Φ ( 〈x, y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) forms a positive semidefinite Hilbert-Schmidt operator on L2(N (0, I)) iff Φ can be decomposed as
Φ(t, q, q′) = ∞∑ l=0 Al(q, q ′)c−1d,lC ( d−22 ) l (t) (16)
satisfying ∞∑ l=0 ‖Al‖2c−2d,l ‖C ( d−22 ) l ‖ 2 = ∞∑ l=0 ‖Al‖2c−2d,l π23−dΓ(l + d− 2) l!(l + d−22 )Γ ( d−2 2
)2 <∞, (17) where
• cd,l = d−2d+2l−2
• C( d−2 2 ) l (t) are Gegenbauer polynomials as in Appendix H.2, with ‖C ( d−22 ) l ‖ =√∫ 1 −1 C ( d−22 ) l (t) 2(1− t2) d−12 −1 dt denoting the norm of C( d−2 2 ) l in L 2((1− t2) d−12 −1).
• and Al are positive semidefinite Hilbert-Schmidt kernels on L2( 1dχ 2 d), the L 2 space over the probability measure of a χ2d-variable divided by d, and with ‖Al‖ denoting the HilbertSchmidt norm of Al.

heading: In addition, K is positive definite iff all Al are.
text: See Appendix I for a proof. As a consequence, K has an eigendecomposition as follows under the standard Gaussian measure in d dimensions. Corollary H.15. Suppose K and Φ are as in Thm H.14 and K is a positive semidefinite HilbertSchmidt operator, so that Eq. (16) holds, with Hilbert-Schmidt kernels Al. Let Al have eigendecomposition
Al(q, q ′) = ∞∑ i=0 aliuli(q)uli(q ′) (18)
for eigenvalues ali ≥ 0 and eigenfunctions uli ∈ L2( 1dχ 2 d) with Eq∼ 1dχ2d uli(q) 2 = 1. Then K has eigenvalues {ali : l, i ∈ [0,∞)}, and each eigenvalue ali corresponds to the eigenspace
uli ⊗Hd−1,(l) def = { uli ( ‖x‖2
d
) f ( x
‖x‖
) : f ∈ Hd−1,(l) } ,
where Hd−1,(l) is the space of degree l spherical harmonics on the unit sphere Sd−1 of ambient dimension d.
For certain simple F , we can obtain {Al}l≥0 explicitly. For example, suppose K is degree-s positive-homogeneous, in the sense that, for a, b > 0,
K(ax, by) = (ab)sK(x, y).
This happens when K is the CK or NTK of an MLP with degree-s positive-homogeneous. Then it’s easy to see that Φ(t, q, q′) = (qq′)sΦ̄(t) for some Φ̄ : [−1, 1]→ R, and
Φ(t, q, q′) = ∞∑ l=0 (qq′)salc −1 d,lC ( d−22 ) l (t)
where {al}l are the Gegenbauer coefficients of Φ̄,
Φ̄(t) = ∞∑ l=0 alc −1 d,lC ( d−22 ) l (t).
We can then conclude with the following theorem. Theorem H.16. Suppose K : (Rd)2 → R is a kernel given by
K(x, y) = R(‖x‖/d)R(‖x‖/d)Φ̄(〈x, y〉/‖x‖‖y‖)
for some functions R : [0,∞)→ R, Φ̄ : [−1, 1]→ R. Let
Φ̄(t) = ∞∑ l=0 alc −1 d,lC ( d−22 ) l (t)
be the Gegenbauer expansion of Φ̄. Also define λ def = √
E q∼ 1dχ 2 d R(q)2.
Then over the standard Gaussian in Rd, K has the following eigendecomposition
• For each l ≥ 0, λ−1R⊗Hd−1,(l) = {λ−1R(‖x‖2/d)f(x/‖x‖) : f ∈ Hd−1,(l)}
is an eigenspace with eigenvalue λ2al.
• For any S ∈ L2( 1dχ 2 d) that is orthogonal to R, i.e.
E q∼ 1dχ 2 d S(q)R(q) = 0,
the function S(‖x‖2/d)f(x/‖x‖)
for any f ∈ L2(Sd−1) is in the null space of K.
Proof. The Al in Eq. (16) for K are all equal to
Al(q, q ′) = λ2
R(q)
λ
R(q′)
λ .
This is a rank 1 kernel (on L2( 1dχ 2 d)), with eigenfunction R/λ and eigenvalue λ 2. The rest then follows straightforwardly from Thm H.11.
A common example where Thm H.16 applies is when K is the CK of an MLP with relu, or more generally degree-s positive homogeneous activation functions, so that the R in Thm H.16 is a polynomial.
In general, we cannot expect K can be exactly diagonalized in a natural basis, as {Al}l≥0 cannot even be simultaneously diagonalizable. We can, however, investigate the “variance due to each degree of spherical harmonics” by computing
al def = E
q∼d−1χ2d Al(q, q) (19)
which is the coefficient of Gegenbauer polynomials in
Φ̂d(t) def = E
q∼d−1χ2d Φ(t, q, q) = ∞∑ l=0 alc −1 d,lC ( d−22 ) l (t). (20)
Proposition H.17. Assume that Φ(t, q, q′) is continuous in q and q′. Suppose for any d and any t ∈ [−1, 1], the random variable Φ(t, q, q′) with q, q′ ∼ d−1χ2d has a bound |Φ(t, q, q′)| ≤ Y for some random variable Y with E |Y | <∞. Then for every t ∈ [−1, 1],
lim d→∞
|Φ̂d(t)− Φ(t, 1, 1)| = 0.
Proof. By the strong law of large number, d−1χ2d converges to 1 almost surely. Because Φ(t, q, q ′) is continuous in q and q′, almost surely we have Φ(t, q, q′) → Φ(t, 1, 1) almost surely. Since Φ is bounded by Y , by dominated convergence, we have
Φ̂d(t)− Φ(t, 1, 1) = E q,q′∼d−1χ2d
Φ(t, q, q′)− Φ(t, 1, 1)→ 0
as desired.
In the final part of this section, we show that in the limit of large input dimension d, the top eigenvalues of K over the standard Gaussian can be easily described (Thm H.19). First, we need to specify some conditions on Φ. Definition H.18. Φ : [−1, 1]× R+ × R+ → R is called reasonable if
• There is a Taylor expansion in t,
Φ(t, q, q′) = ∞∑ `=0 B`(q, q ′) `! t` (21)
that is absolutely convergent on (t, q, q′) ∈ [−1, 1]× R+ × R+, and such that each B` is smooth (C∞) on (q, q′) ∈ R+ × R+.
• for all (t, q, q′) ∈ [−1, 1] × R+ × R+, and for any l ∈ [0,∞), we have |Bl(t, q, q′)| ≤ C(1 + |q|r + |q′|r) for some constants C, r > 0 that may depend on l but not on t, q, q′.
Theorem H.19. Suppose K is a the CK or NTK of an MLP with polynomially bounded activation function. For every degree l, K over N (0, Id) has an eigenvalue al0 at spherical harmonics degree l (in the sense of Cor H.15) with
al0 ∈ [ E
q,q′∼ 1dχ 2 d
Al(q, q ′), E
q∼ 1dχ 2 d
Al(q, q)
] ,
where Al is as in Eq. (18). Furthermore,
lim d→∞ dl E q,q′∼ 1dχ 2 d Al(q, q ′) = lim d→∞ dl E q∼ 1dχ 2 d Al(q, q) = Φ (l)(0, 1, 1).
Here Φ(l) is the lth derivative of Φ(t, q, q′) against t.
Proof. By Lem H.20 below, Φ is reasonable. Let al = Eq∼ 1dχ2d Al(q, q) as in Eq. (19), and let bl def = Eq′∼ 1dχ2d Al(q, q
′). Let ali be the ith largest eigenvalue of Al as in Eq. (18), with al0 being the largest. Note that all of these quantities al, bl, ali depend on d, but we suppress this notationally. We seek to prove the following claims:
1. al0 ≤ al
2. al0 ≥ bl
3. dlbl → Φ(l)(0, 1, 1)
4. dlal → Φ(l)(0, 1, 1)
Claim 1: al0 ≤ al. First note that al is the trace of the operator Al, so that al = ∑∞ i=0 ali. Thus, any eigenvalue ali is at most al.
Claim 2: al0 ≥ bl. Now, by Min-Max theorem, the largest eigenvalue al0 of Al is equal to
al0 = sup f Eq,q′ f(q)f(q′)Al(q, q′) Eq f(q)2
where 0 6= f ∈ L2( 1dχ 2 d) and q, q ′ ∼ 1dχ 2 d. If we set f(q) = 1 identically, then we get
al0 ≥ E q,q′ Al(q, q ′) = bl,
as desired.
Claims 3 and 4. Now note that we have the following equalities in the Hilbert spaceL2((1−t2) d−12 −1), because of the absolute convergence in Eq. (17):
Φ̂(t) def = E
q∼ 1dχ 2 d
Φ(t, q, q) = ∞∑ `=0 E q∼ 1dχ 2 d A`(q, q)c −1 d,`C ( d−22 ) ` (t) = ∞∑ `=0 a`c −1 d,`C ( d−22 ) ` (t)
Φ̌(t) def = E
q,q′∼ 1dχ 2 d
Φ(t, q, q′) = ∞∑ `=0 E q,q′∼ 1dχ 2 d A`(q, q ′)c−1d,`C ( d−22 ) ` (t) = ∞∑ `=0 b`c −1 d,`C ( d−22 ) ` (t).
Thus, by Eq. (15),
al = Γ
( d
2 ) ∞∑ k=0
Φ̂(l+2k)(0) 2l+2kk!Γ ( d 2 + l + k ) . By the absolute convergence of Eq. (21), differentiation commutes with expectation:
Φ̂(l+2k)(0) = ∂l+2kt E q∼ 1dχ 2 d Φ(0, q, q) = E q∼ 1dχ 2 d Bl+2k(q, q) = E q∼ 1dχ 2 d Φ(l+2k)(0, q, q),
where Bl+2k is as in Eq. (21). Furthermore, as d → ∞, because Φ(l+2k)(t, q, q′) is smooth and polynomially bounded in q and q′, while 1dχ 2 d has exponential tails and converges a.s. to 1, we have
E q∼ 1dχ 2 d Φ(l+2k)(0, q, q)→ Φ(l+2k)(0, 1, 1).
Finally, by dominated convergence, we get
dlal = Γ
( d
2 ) ∞∑ k=0
dlΦ̂(l+2k)(0) 2l+2kk!Γ ( d 2 + l + k ) → Φ(l)(0, 1, 1) as desired.
The proof of lim dlbl = Φ(l)(0, 1, 1) is similar.
Lemma H.20. Let Φ be the function corresponding to the CK or NTK of an MLP with polynomially bounded activation functions. Then Φ is reasonable (Defn H.18).
Proof. The Taylor expansion Eq. (21), its absolute convergence, and the smoothness of Bl follows from the proof of Thm I.3. The polynomial-boundedness of Bl follows trivially from the polynomialboundedness of activation.
H.4 IN HIGH DIMENSION, BOOLEAN CUBE ≈ SPHERE ≈ STANDARD GAUSSIAN
For the same neural kernel K with K(x, y) = Φ(〈x, y〉/‖x‖‖y‖, ‖x‖2/d, ‖y‖2/d), let Φ̂d be as defined in Eq. (20), and let ΦS def = Φ(t, 1, 1). Thus ΦS is the univariate Φ that we studied in Appendix H.2 in the context of the sphere.
Empirical verification of the spectral closeness As d → ∞, 1dχ 2 d converges almost surely to 1, and each Al in Eq. (16) converges weakly to the simple operator that multiplies the input function by Al(1, 1). We verify in Fig. 17 that, for different erf kernels, Φ̂d approximates ΦS when d is large, which would then imply that their spectra become identical for large d. Note that this is tautologically true for relu kernels with no bias σb = 0 because they are positive-homogeneous.
Next, we compute the eigenvalues ofK on the boolean cube and the sphere, as well as the eigenvalues of the kernel K̂(x, y) def= Φ̂d(〈x, y〉/‖x‖‖y‖) on the sphere, which “summarizes” the eigenvalues of K on the standard Gaussian distribution. In Fig. 18, we compare them for different dimensions d up to degree 5 (there, “Gaussian” means K̂ on the sphere). We see that by dimension 128, all eigenvalues shown are very close to each other for each data distribution.
Theoretical verification of the spectral closeness Thms H.13, H.19 and I.6 show that when the dimension d is sufficiently large, the top eigenvalues of K as an operator on d is very close to the top eigenvalues of K as an operator on √ dSd−1 and those of K as an operator on N (0, Id):
Corollary H.21. Let K be the CK or NTK of an MLP with polynomially bounded activation function. ThenK can be expressed asK(x, y) = Φ(〈x, y〉/d) for some smooth Φ : [−1, 1]→ R. Let µk denote the eigenvalue of K as an operator on d corresponding to the eigenspace of degree k polynomials (see Thm 3.1). Let ak denote the eigenvalue of K as an operator on √ dSd−1 corresponding to the eigenspace of degree k spherical harmonics (see Thm H.11). Let ak0 denote the largest eigenvalue of K as an operator on N (0, Id) restricted to degree k spherical harmonics (see Thm H.19). Then for any fixed k ≥ 0, for sufficiently large d, we have
lim d→∞ dkµk = lim d→∞ dkak = lim d→∞ dkak0 = Φ (k)(0)
where Φ(k)(0) is the kth derivative of Φ at 0. If Φ(k)(0) 6= 0, then
lim d→∞ µk/ak = lim d→∞ µk/ak0 = 1.
If we fix k and let the input dimension d→∞, then the fractional variance of degree k converges to
(k!)−1Φ(k)(0)/Φ(1) = (k!)−1Φ(k)(0)∑ j≥0(j!) −1Φ(j)(0)
for all three input distributions d, √ dSd−1, and N (0, Id).
Practically speaking, only the top eigenvalues matter. Observe that, in the empirical and theoretical results above, we only verify that the top eigenvalues (µk, ak, or ak0 for k small compared to d) are close when d is large. While this result may seem very weak at face value, in practice, the closeness of these top eigenvalues is the only thing that matters. Indeed, in machine learning, we will only ever have a finite number, say N , of training samples to work with. Thus, we can only use a finite N ×N submatrix of the kernel K. This submatrix, of course, has only N eigenvalues. Furthermore, if these samples are collected in an iid fashion (as is typically assumed), then these eigenvalues approximate the largest N eigenvalues (top N counting multiplicity) of the kernel K itself (Tropp, 2015). As such, the smaller eigenvalues of K can hardly be detected in the training sample, and cannot affect the machine learning process very much.
Let’s discuss a more concrete example: Fig. 18 shows that the boolean cube eigenvalues µk are very close to the sphere eigenvalues ak for all k ≤ 5. Over the boolean cube, µ0, . . . , µ5 cover eigenspaces of total dimension ( d 0 ) + · · ·+ ( d 5 ) , which is 275,584,033 when d = 128. We need at least that many samples to be able to even detect the eigenvalue µ6 and the possible difference between it and the sphere eigenvalue a6. But note in comparison, Imagenet, one of the most common large datasets in use today, has only about 15 million samples, 10 times less than the number above.
Additionally, in this same comparison, d = 128 dramatically pales compared to Imagenet’s input dimension 3× 2562 = 196608, and even to the those of the smaller common datasets like CIFAR10 (d = 3× 322 = 3072) and MNIST (d = 242 = 576) — if we were to even use the input dimension of MNIST above, then µ0, . . . , µ5 would cover eigenspaces of 523 billion total dimensions! Hence, it is quite practically relevant to consider the effect of large d on the eigenvalues, while keeping k small. Again, we remark that even when one fixes k and increases d, the dimension of eigenspaces affected by our limit theorems Thms H.13, H.19 and I.6 increases like Θ(dk), which implies one needs an increasing number Θ(dk) of training samples to see the difference of eigenvalues in higher degrees k.
Finally, from the perspective of fractional variance, we also can see that only the top k spectral closeness matters: By Cor H.21, for any > 0, there is a k such that the total fractional variance of degree 0 to degree k (corresponding to eigenspaces of total dimension Θ(dk)) sums up to more than 1− , for the cube, the sphere, and the standard Gaussian simultaneously, when d is sufficiently large. This is because the asymptotic fractional variance is completely determined by the derivatives of Φ at t = 0.

heading: I OMITTED PROOFS
text: Theorem 3.1. On the d-dimensional boolean cube d, for every S ⊆ [d], χS is an eigenfunction of K with eigenvalue
µ|S| def = E x∈ d xSK(x,1) = E x∈ d xSΦ (∑ i xi/d ) , (4)
where 1 = (1, . . . , 1) ∈ d. This definition of µ|S| does not depend on the choice S, only on the cardinality of S. These are all of the eigenfunctions of K by dimensionality considerations.8
Proof. We directly verify KχS = µ|S|χS . Notice first that K(x, y) = Φ(〈x, y〉) = Φ(〈x y,1〉) = K(x y,1)
where is Hadamard product. We then calculate KχS(y) = E
x K(y, x)xS
= E x K(1, x y)(x y)SyS .
8Readers familiar with boolean Fourier analysis may be reminded of the noise operator Tρ, ρ ≤ 1 (O’Donnell, 2014, Defn 2.46). In the language of this work, Tρ is a neural kernel with eigenvalues µk = ρk.
Here we are using the fact that x and y are boolean to get xS = (x y)SyS . Changing variable z def = x y, we get
KχS(y) = y S E z K(1, z)zS = µ|S|χS(y)
as desired. Finally, note that µ|S| is invariant under permutation of [d], so indeed it depends only on the size of S.
Theorem H.11 (Schoenberg). Suppose Φ : [−1, 1]→ R is in L2((1− t2) d−12 −1), so that it has the Gegenbauer expansion
Φ(t) a.e. = ∞∑ l=0 alc −1 d,lC ( d−22 ) l (t).
Then K has eigenspaces Hd−1,(l)√ d def = {f(x/
√ d) : f ∈ Hd−1,(l)} with corresponding eigenval-
ues al. Since ⊕∞
l=0H d−1,(l)√ d
is an orthogonal decomposition of L2( √ dSd−1), this describes all
eigenfunctions of K considered as an operator on L2( √ dSd−1).
Proof. Let f ∈ Hd−1,(n). Then for any x ∈ √ dSd−1,
E y∈ √ dSd−1
K(x, y)f(y/ √ d) = E
y∈ √ dSd−1
Φ(〈x, y〉/d)f(y/ √ d)
= E ỹ∈Sd−1
Φ(〈x/ √ d, ỹ〉)f(ỹ)
= E ỹ∈Sd−1 ( ∞∑ l=0 al 1 cd,l C ( d−22 ) l (〈x/ √ d, ỹ〉) ) f(ỹ)
= E ỹ∈Sd−1 ( ∞∑ l=0 alZ d−1,(l) x/ √ d (ỹ) ) f(ỹ)
= an E ỹ∈Sd−1 Z d−1,(n) x/ √ d (ỹ)f(ỹ)
by orthogonality
= anf(x/ √ d)
by reproducing property.
Lemma 3.2. With µk as in Thm 3.1, µk = 2 −d(I − T∆)k(I + T∆)d−kΦ(1) (5)
= 2−d d∑ r=0 Cd−k,kr Φ (( d 2 − r ) ∆ ) (6)
where
Cd−k,kr def = ∑ j=0 (−1)r+j ( d− k j )( k r − j ) . (7)
Proof. Because ∑ i xi/d only takes on values {− d 2 ∆, (− d 2 + 1)∆, . . . , ( d 2 − 1)∆, d 2 ∆}, where ∆ = 2d , we can collect like terms in Eq. (4) and obtain
µk = 2 −d d∑ r=0 ∑ x has r ‘−1’s ( k∏ i=1 xi ) Φ (( d 2 − r ) ∆ ) which can easily be shown to be equal to
µk = 2 −d d∑ r=0 Cd−k,kr Φ (( d 2 − r ) ∆ ) ,
proving Eq. (6) in the claim. Finally, observe that Cd−k,kr is also the coefficient of x r in the polynomial (1− x)k(1 + x)d−k. Some operator arithmetic then yields Eq. (5).
Theorem H.14. A function K : (Rd)2 → R of the form
K(x, y) = Φ ( 〈x, y〉 ‖x‖‖y‖ , ‖x‖2 d , ‖y‖2 d ) forms a positive semidefinite Hilbert-Schmidt operator on L2(N (0, I)) iff Φ can be decomposed as
Φ(t, q, q′) = ∞∑ l=0 Al(q, q ′)c−1d,lC ( d−22 ) l (t) (16)
satisfying ∞∑ l=0 ‖Al‖2c−2d,l ‖C ( d−22 ) l ‖ 2 = ∞∑ l=0 ‖Al‖2c−2d,l π23−dΓ(l + d− 2) l!(l + d−22 )Γ ( d−2 2
)2 <∞, (17) where
• cd,l = d−2d+2l−2
• C( d−2 2 ) l (t) are Gegenbauer polynomials as in Appendix H.2, with ‖C ( d−22 ) l ‖ =√∫ 1 −1 C ( d−22 ) l (t) 2(1− t2) d−12 −1 dt denoting the norm of C( d−2 2 ) l in L 2((1− t2) d−12 −1).
• and Al are positive semidefinite Hilbert-Schmidt kernels on L2( 1dχ 2 d), the L 2 space over the probability measure of a χ2d-variable divided by d, and with ‖Al‖ denoting the HilbertSchmidt norm of Al.

heading: In addition, K is positive definite iff all Al are.
text: Proof. Note that an isotropic Gaussian vector z ∼ N (0, I) can be sampled by independently sampling its direction v uniformly from the sphere Sd−1 and sampling its magnitude r from a chi distribution χd with d degrees of freedom. In the following, we adopt the following notations
x, y ∈ Rd, q = ‖x‖2/d, q′ = ‖y‖2/d, v = x/‖x‖, v′ = y/‖y‖ Then, by the reasoning above,
Kf(x/ √ d) = E
v′∼Sd−1 q′∼d−1χ2d
Φ (〈v, v′〉, q, q′) f( √ q′v).
Here f ∈ L2(N (0, I/d)) (so that f(·/ √ d) ∈ L2(N (0, I))) and d−1χ2d is the distribution of a χ2d random variable divided by d. Because of this decomposition of N (0, I) into a product distribution of direction and magnitude, its space of L2 functions naturally decomposes into a tensor product of corresponding L2 spaces
L2(N (0, I)) = L2(Sd−1)⊗ L2(d−1χ2d). where every f ∈ L2(N (0, I)) can be written as a sum
f(x) = ∞∑ l=0 Rl(q)Pl(v)
where q = ‖x‖2/d, v = x/‖x‖ as above, Pl ∈ Hd−1,(l) is a spherical harmonic of degree l, and Rl ∈ L2(d−1χ2d) (i.e. Rl(·/d) ∈ L2(χ2d)). Here, the equality is understood as a convergence of the RHS partial sums in the Hilbert space L2(Sd−1)⊗ L2(d−1χ2d). Likewise, assuming K is Hilbert-Schmidt, then its Hilbert-Schmidt norm (measured against the uniform distribution over the sphere) is bounded:
‖K‖2HS = E v,v′∼Sd−1 q,q′∼d−1χ2d Φ(〈v, v′〉, q, q′)2
= E q,q′∼d−1χ2d ∫ 1 −1 Φ(t, q, q′)2(1− t2) d−1 2 −1 dt.
Thus Φ resides in the tensor product space
Φ ∈ F def= L2((1− t2) d−1 2 −1)⊗ L2(d−1χ2d)⊗ L2(d−1χ2d)
and therefore can be expanded as
Φ(t, q, q′) = ∞∑ l=0 Al(q, q ′)c−1d,lC ( d−22 ) l (t)
where C( d−2 2 )
l (t) are Gegenbauer polynomials as in Appendix H.2 and Al ∈ L2(d−1χ2d)⊗2. By Lem I.1 below, K being Hilbert-Schmidt implies that each Al is, as well. Furthermore, since ‖K‖2HS = ‖Φ‖2F is finite, we have
‖Φ‖2F = ∞∑ l=0 ‖Al‖2c−2d,l ‖C ( d−22 ) l ‖ 2 <∞.
Here, ‖Al‖ = √ Eq,q′∼ 1dχ2d Al(q, q ′)2 is the Hilbert-Schmidt norm of Al (i.e. its norm in
L2(d−1χ2d) ⊗2), and ‖C(
d−2 2 ) l ‖ = √∫ 1 −1 C ( d−22 ) l (t) 2(1− t2) d−12 −1 dt is the norm of C( d−2 2 ) l in
L2((1− t2) d−12 −1). Simplifying according to Eq. (14) yields the equality in the claim. Conversely, if each Al is Hilbert-Schmidt satisfying Eq. (17), then K is obviously a Hilbert-Schmidt kernel as it has finite Hilbert-Schmidt norm.
Lemma I.1. For each l, Al is a positive semidefinite Hilbert-Schmidt kernel. It is positive definite if K is.
Proof. With q = ‖x‖2/d, v = x/‖x‖ as above, let f(x/ √ d) = f( √ qv) = Rm(q)Pm(v) for some degree m nonzero spherical harmonics Pm ∈ Hd−1,(m) and some scalar function Rm ∈ L2(d−1χ2d) that is not a.e. zero. We have
Kφ(x/ √ d) = E
v′∼Sd−1 q′∼d−1χ2d
Φ(〈v, v′〉, q, q′)f( √ q′v)
= E v′∼Sd−1 q′∼d−1χ2d
Rm(q ′)Pm(v ′) ∞∑ l=0 Al(q, q ′)c−1d,lC ( d−22 ) l (〈v, v ′〉)
= E v′∼Sd−1 q′∼d−1χ2d
Rm(q ′)Pm(v ′) ∞∑ l=0 Al(q, q ′)Zd−1,(l)v (v ′)
= Pm(v) E q′∼d−1χ2d
Am(q, q ′)Rm(q ′).
Therefore,
E x∼N (0,I)
f(x)Kf(x/ √ d) = ( E
v∼Sd−1 Pm(v)
2 )( E
q,q′∼d−1χ2d Rm(q)Am(q, q
′)Rm(q ′) ) is nonnegative, and is positive if K is positive definite, by the assumption above that f is not a.e. zero. Since EPm(v)2 > 0 as Pm 6= 0, we must have
E q,q′∼d−1χ2d
Rm(q)Am(q, q ′)Rm(q ′) ≥ 0 (or > 0 if K is positive definite).
This argument holds for any Rm ∈ L2(d−1χ2d) that is not a.e. zero, so this implies that Al is positive semidefinite, and positive definite if K is so.

heading: Weak Spectral Simplicity Bias
text: Theorem 4.1 (Weak Spectral Simplicity Bias). Let K be the CK or NTK of an MLP on a boolean cube d. Then the eigenvalues µk, k = 0, . . . , d, satisfy
µ0 ≥ µ2 ≥ · · · ≥ µ2k ≥ · · · , µ1 ≥ µ3 ≥ · · · ≥ µ2k+1 ≥ · · · . (9)
Proof. Again, there is a function Φ such that K(x, y) = Φ(〈x, y〉/‖x‖‖y‖). By Thm I.3, Φ has a Taylor expansion with only nonnegative coefficients.
Φ(c) = a0 + a1c+ a2c 2 + · · · .
By Lem I.2, Eq. (9) is true for polynomials Φ(c) = cr,
µ0(c r) ≥ µ2(cr) ≥ · · · , µ1(cr) ≥ µ3(cr) ≥ · · · . (22) Then since each µk = µk(Φ) is a linear function of Φ, µk(Φ) = ∑∞ r=0 arµk(c
r) also follows the same ordering.
Lemma I.2. Let T∆ be the shift operator with step ∆ that sends a function Φ(·) to Φ(· −∆). Let Φ(c) = ct for some t. Let µdk be the eigenvalue of K(x, y) = Φ(〈x, y〉/‖x‖‖y‖) on the boolean cube d. Then for any 0 ≤ k ≤ d,
µdk = 0 if t+ k is odd (23)
1 ≥ µdk−2 ≥ µdk ≥ 0 if t+ k is even. (24)
We furthermore have the identity
µdk−2 − µdk = ( d− 2 d )t µd−2k−2
for any 2 ≤ k ≤ d.
Proof. As in Eq. (5),
2−d(I − T∆)k(I + T∆)d−kΦ(1) = 2−d d∑ r=0 Cd−k,kr Φ (( d 2 − r ) ∆ )
where Cd−k,kr def = ∑ j=0(−1)r+j ( d−k j )( k r−j ) . It’s easy to observe that
Cd−k,kr = C d−k,k d−r if k is even Cd−k,kr = −C d−k,k d−r if k is odd.
In the first case, when Φ(c) = ct with t odd, then by symmetry µk = 0. Similarly for the second case when t is even. This finishes the proof for Eq. (23).
Note that it’s clear from the form of Eq. (5) that µdk ≤ Φ(1) = 1 always. So we show via induction on d that the rest of Eq. (24) is true for any t ∈ N. The induction will reduce the case of d to the case of d− 2. So for the base case, we need both d = 0 and d = 1. Both can be shown by some simple calculations.
Now for the inductive step, assume Eq. (24) for d− 2, and we observe that
µdk−2 − µdk = 2−d(I − T2/d)k−2(I + T2/d)d−k((I + T2/d)2 − (I − T2/d)2)Φ(1) = 2−d(I − T2/d)k−2(I + T2/d)d−k(4T2/d)Φ(1) = 2−d+2(I − T2/d)k−2(I + T2/d)d−kΦ(1− 2/d).
By Eq. (5), we can expand this into
2−d+2 d−2∑ r=0 Cd−k,k−2r Φ ( 1− 2 d − r 2 d )
= 2−d+2 d−2∑ r=0 Cd−k,k−2r ( 1− 2 d − r 2 d )t
= ( d− 2 d )t 2−d+2 d−2∑ r=0 Cd−k,k−2r ( 1− r 2 d− 2 )t = ( d− 2 d )t µd−2k−2.
By induction, µd−2k−2 ≥ 0, so we have
µdk−2 − µdk = ( d− 2 d )t µd−2k−2 ≥ 0
as desired.
Theorem I.3. Let K be the CK or NTK of an MLP with domain √ dSd−1 ⊆ Rd. Then K(x, y) = Φ ( 〈x,y〉 ‖x‖‖y‖ ) where Φ : [−1, 1]→ R has a Taylor series expansion around 0
Φ(c) = a0 + a1c+ a2c 2 + · · · (25) with the properties that ai ≥ 0 for all i and ∑ i ai = Φ(1) ≤ ∞, so that Eq. (25) is absolutely convergent on c ∈ [−1, 1].
Proof. We first prove this statement for the CK of an L layer MLP with nonlinearity φ. Let Φl be the corresponding Φ-function for the CK Σl. It suffices to prove this by induction on depth L. For L = 1, Σ1(x, x′) = σ2w(n
0)−1〈x, x′〉+ σ2b by Eq. (CK), so clearly Φ1(c) = a0 + a1c where a0 and a1 are nonnegative.
Now for the inductive step, suppose Σl−1 satisfies the property that Φl−1 has Taylor expansion of the desired form. We seek to show the same for Σl and Φl. First notice that it suffices to show this for Vφ(Σ l−1), since multiplication by σ2w and addition of σ 2 b preserves the property. But
Vφ(Σ l−1)(x, x′) = Eφ(z)φ(z′)
where (z, z′) ∼ N ( 0, ( Σl−1(x, x) Σl−1(x, x′) Σl−1(x, x′) Σl−1(x′, x′) )) = N ( 0, ( Φl−1(1) Φl−1(c) Φl−1(c) Φl−1(1) )) .
Using the notation of Lem I.4, we can express this as
Vφ(Σ l−1)(x, x′) = Φl(1)φ̂
( Φl−1(c)
Φl−1(1)
) . (26)
Substituting the Taylor expansion of Φ l−1(c)
Φl−1(1) into the Taylor expansion of φ̂ given by Lem I.4 gives us
a Taylor series whose coefficients are all nonnegative, since those of Φ l−1(c)
Φl−1(1) and φ̂ are nonnegative
as well. In addition, plugging in 1 for c in this Taylor series shows that the sum of the coefficients equal φ̂(1) = 1, so the series is absolutely convergent for c ∈ [−1, 1]. This proves the inductive step. For the NTK, the proof is similar, except we now also have a product step where we multiply Vφ′(Σ
l−1) with Θl−1, and we simply just need to use the fact that product of two Taylor series with nonnegative coefficients is another Taylor series with nonnegative coefficients, and that the resulting radius of convergence is the minimum of the original two radii of convergence.
Lemma I.4. Consider any σ > 0 and any φ : R→ R square-integrable against N (0, σ2). Then the function
φ̂ : [−1, 1]→ [−1, 1], φ̂(c) def= Eφ(x)φ(y) Eφ(x)2
, where (x, y) ∼ N ( 0, ( σ2 σ2c σ2c σ2 )) ,
has a Taylor series expansion around 0
φ̂(c) = a0 + a1c+ a2c 2 + · · · (27) with the properties that ai ≥ 0 for all i and ∑ i ai = φ̂(1) = 1, so that Eq. (27) is absolutely convergent on c ∈ [−1, 1].
Proof. Since the denominator of φ̂ doesn’t depend on c, it suffices to prove this Taylor expansion exists for its numerator. Let φ̃(c) def= φ(σc). Then
E [ φ(x)φ(y) : (x, y) ∼ N ( 0, ( σ2 σ2c σ2c σ2 ))] = E [ φ̃(x)φ̃(y) : (x, y) ∼ N ( 0, ( 1 c c 1 ))] By Fact I.5, this is just
b20 + b 2 1c+ b 2 2c 2 + · · · where bi are the Hermite coefficients of φ̃. Since φ ∈ L2(N (0, σ2)),∑
i
b2i <∞
as desired.
Fact I.5 (O’Donnell (2014)). For any φ, ψ ∈ L2(N (0, 1)), E [ φ(x)ψ(y) : (x, y) ∼ N ( 0, ( 1 c c 1 ))] = a0b0 + a1b1c+ a2b2c 2 + · · ·
where ai and bi are respectively the Hermite coefficients of φ and ψ.
Fixing k, taking the d→∞ limit In this subsection, we prove that µk ∼ d−k as d→∞ with k fixed. More precisely,
Theorem I.6. Let K be the CK or NTK of an MLP on a boolean cube d. Then K can be expressed as K(x, y) = Φ(〈x, y〉/d) for some Φ : [−1, 1]→ R. If we fix k and let d→∞, then
lim d→∞
dkµk = Φ (k)(0),
where Φ(k) denotes the kth derivative of Φ.
An immediate consequence, using the fact that the degree k eigenspace has dimension ( d k ) ∼ dk/k! for large d, is Theorem 5.1 (Asymptotic Fractional Variance). Let K be the CK or NTK of an MLP on a boolean cube d. ThenK can be expressed asK(x, y) = Φ(〈x, y〉/d) for some analytic function Φ : R→ R. If we fix k and let the input dimension d→∞, then the fractional variance of degree k converges to
(k!)−1Φ(k)(0)/Φ(1) = (k!)−1Φ(k)(0)∑ j≥0(j!) −1Φ(j)(0)
where Φ(k) denotes the kth derivative of Φ.
Let’s first give some intuition for why Thm I.6 should be true. By Eqs. (8) and (5),
µk =
( 1 + T∆
2
)d−k ( 1− T∆
2
)k Φ(1)
= 1
2d d−k∑ r=0 ( d− k r ) Φ[k](1− r∆)
where Φ[k] = (1− T∆)k Φ is the kth backward finite difference with step ∆ = 2/d. Approximating the finite difference with derivative, we thus would expect that, as suggested by Lem I.8,
Φ[k](x) ≈ ∆kΦ(k)(x) = (2/d)kΦ(k)(x), when d is large,
where Φ(k) is the kth derivative of Φ. Since 1 2d−k ( d−k r ) is the probability mass at 1 − r∆ of the
binomial variable B def= 1d−k ∑d−k i=1 Xi, where Xi is the Bernoulli variable taking ±1 value with half chance each. This random variable converges almost surely to the delta distribution at 0, by law of large numbers. Therefore, we would expect
µk ∼ (1/d)kΦ(k)(0) as d→∞. There is a single difficulty when trying to formalize this intuition. One is the possibility that the kth derivative Φ(k) does not exist at the endpoints 1 and −1 (note that it must exist in the interior, because by Thm I.3, Φ is analytic on (−1, 1)). In this case, we must show that the portion of the interval [−1, 1] near the endpoints contributes exponentially little mass to the probability distribution of the binomial variable B, that it suppresses any blowup that can happen at the edge.
We formalize this reasoning in the proof below.
Proof. We first prevent any blowup that may happen at the edge. Let Φ[k] be the kth backward difference of Φ with step size ∆ = 2/d, as in Lem I.8. First note the trivial bound
Φ[k](x) = k∑ s=0 (−1)s ( k s ) Φ(x− s∆) ≤ k∑ s=0 ( k s ) Φ(x− s∆)
≤ 2kΦ(x) ≤ 2kΦ(1) since Φ is nondecreasing by Lem I.7. Then for any R ∈ [0, (d− k)/2− 1],
dkµk = dk
2d d−k∑ r=0 ( d− k r ) Φ[k](1− r∆)
= dk
2d d−k−R−1∑ r=R+1 ( d− k r ) Φ[k](1− r∆) + d k 2d ( R∑ r=0 + d−k∑ r=d−k−R )( d− k r ) Φ[k](1− r∆).
Now, in the second term, Φ[k](1 − r∆) ≤ 2kΦ(1) as above, so that, by the binomial coefficient entropy bound (Fact I.9),
0 ≤ dkµk − dk
2d d−k−R−1∑ r=R+1 ( d− k r ) Φ[k](1− r∆) ≤ d k 2d ( R∑ r=0 + d−k∑ r=d−k−R )( d− k r ) 2kΦ(1)
≤ d k
2d−k 2 · 2H(p)(d−k)Φ(1)
= dkΦ(1)
2(1−H(p))d−(1−H(p))k−1 , (28)
where p = R/(d− k). If we choose R = bp̃(d− k)c for a fixed p̃ < 1/10, then H(p) < 1, and (28) is decreasing exponentially fast to 0 with d.
Now we formalize the law of large number intuition. It then suffices to show that
dk
2d d−k−R−1∑ r=R+1 ( d− k r ) Φ[k](1− r∆)→ Φ(k)(0).
We will do so by presenting an upper and a lower bound which both converge to this limit.
We first discuss the upper bound. By Lem I.8,
dk
2d d−k−R−1∑ r=R+1 ( d− k r ) Φ[k](1− r∆) ≤ 1 2d−k d−k−R−1∑ r=R+1 ( d− k r ) Φ(k)(1− r∆).
The RHS can be upper bounded by an expectation over a binomial random variable: Let Xi be ±1 Bernoulli variables, and let Φ(k)1
2 p̃ be the function
Φ (k) 1 2 p̃ (x) =
{ Φ(k)(x) if x ∈ [−1 + 12 p̃, 1− 1 2 p̃]
0 otherwise
Note that Φ(k)1 2 p̃ (x) is a bounded function. Then with the binomial variableB = 1d ∑d−k i=1 (Xi+k/(d− k)), we have
1
2d−k d−k−R−1∑ r=R+1 ( d− k r ) Φ(k)(1− r∆) ≤ EΦ(k)1 2 p̃ (B).
By strong law of large numbers, B converges almost surely to 0 as d→∞ with k fixed. Thus the RHS converges to
Φ (k) 1 2 p̃ (0) = Φ(k)(0)
as desired.
A similar argument proceeds for the lower bound, after noting that, by Lem I.8,
1
2d−k d−k−R−1∑ r=R+1 ( d− k r ) Φ(k)(1− (r + k)∆) ≤ d k 2d d−k−R−1∑ r=R+1 ( d− k r ) Φ[k](1− r∆).
Lemma I.7. Let K be the CK or NTK of an MLP with domain √ dSd−1 ⊆ Rd. Then K(x, y) = Φ ( 〈x,y〉 ‖x‖‖y‖ ) where Φ : [−1, 1]→ R is analytic on (−1, 1), and all derivatives of Φ are nonnegative on (−1, 1).
Proof. Note that, by Thm I.3, Φ has a convergent Taylor expansion on (−1, 1), i.e Φ is analytic on (−1, 1). Thus, all derivatives of Φ exist (and are finite) on the interval (−1, 1) and have the obvious Taylor series derived from Eq. (25). For example,
Φ′(c) = a1c+ 2a2c+ 3a3c 2 + · · · .
Such Taylor series also converge on the open interval (−1, 1) (but could diverge on the endpoints −1 and 1). We get the desired result by noting that all ai are nonnegative.
By expressing finite difference as an integral, it’s easy to see that Lemma I.8. With the same setting as in Lem I.7, let Φ[k] be the kth backward finite difference with step size ∆ = 2/d,
Φ[k](x) def = (1− T∆)kΦ(x).
Then Φ[k] is analytic on (−1, 1), has all derivatives nonnegative, and 0 ≤ ∆kΦ(k)(x− k∆) ≤ Φ[k](x) ≤ ∆kΦ(k)(x), (29)
where Φ(k) is the kth derivative of Φ.
Proof. Everything follows immediately from
Φ[k](x) = ∫ ∆ 0 · · · ∫ ∆ 0 Φ(k)(x− h1 − · · · − hk) dh1 · · · dhk
and Lem I.7.
Fact I.9 (Entropy bound on sum of binomial coefficients). For any k ≤ d and R ≤ d/2, R∑ r=0 ( d− k r ) ≤ 2H(p)(d−k)
where p = R/(d− k), and H is the binary entropy H(p) = −p log2 p− (1− p) log2(1− p).
J THE {0, 1}d VS THE {±1}d BOOLEAN CUBE
Valle-Pérez et al. (2018) actually did their experiments on the {0, 1}d boolean cube, whereas here, we have focused on the {±1}d boolean cube. As datasets are typically centered before feeding into a neural network (for example, using Pytorch’s torchvision.transform.Normalize), {±1}d is much more natural. In comparison, using the {0, 1}d cube is equivalent to adding a bias in the input of a network and reducing the weight variance in the input layer, since any x ∈ {±1}d corresponds to 12 (x + 1) ∈ {0, 1}
d. As such, one would expect there is more bias toward low frequency components with inputs from {0, 1}d.
Nevertheless, here we verify that our observations of Section 4 above still holds over the {0, 1}d cube by repeating the same experiments as Fig. 1 in this setting (Fig. 19). Just like over the {±1}d cube, the relu network biases significantly toward certain functions, but with erf, and with increasing σ2w, this lessens. With depth 32 and σ2w, the boolean functions obtained from erf network see no bias at all.

Review summaries:
Aiming to resolve the question whether and why deep networks are biased towards simple functions, this paper gives a spectral analysis on neural networks' conjugate kernel(CK) and neural tangent kernel(NTK) on boolean cube. The eigenfunctions are identified and the eigenvalues are shown computable in polynomial time. Another main contribution of this paper is showing that the simplicity bias exists at least in a weak sense.