Spectrum Dependent Learning Curves in Kernel Regression and Wide Neural Networks

Introduction

Finding statistical patterns in data that generalize beyond a training set is a main goal of machine learning. Generalization performance depends on factors such as the number of training examples, the complexity of the learning task, and the nature of the learning machine. Identifying precisely how these factors impact the performance poses a theoretical challenge. Here, we present a theory of generalization in kernel machines (Schölkopf & Smola, 2001) and neural networks (LeCun et al., 2015) with wide hidden layers that addresses these questions.

The goal of our theory is not to provide worst case bounds on generalization performance in the sense of statistical learning theory (Vapnik, 1999), but to provide analytical expressions that explain the average or a typical performance in the spirit of statistical physics. The techniques we use are a continuous approximation to learning curves previously used in Gaussian processes (Sollich, 1999, 2002; Sollich & Halees, 2002) and the replica method of statistical physics (Sherrington & Kirkpatrick, 1975; Mézard et al., 1987).

We first develop an approximate theory of generalization in kernel regression that is applicable to any kernel. We then use our theory to gain insight into neural networks by using a correspondence between kernel regression and neural network training. When the hidden layers of a neural network are taken to infinite width with a certain initialization scheme, recent influential work (Jacot et al., 2018; Arora et al., 2019; Lee et al., 2019) showed that training a feedforward neural network with gradient descent to zero training loss is equivalent to kernel interpolation (or ridgeless kernel regression) with a kernel called the Neural Tangent Kernel (NTK) (Jacot et al., 2018). Our kernel regression theory contains kernel interpolation as a special limit (ridge parameter going to zero).

Our contributions and results are summarized below:

Using a continuous approximation to learning curves adapted from Gaussian process literature (Sollich, 1999, 2002), we derive analytical expressions for learning curves for each spectral component of a target function learned through kernel regression.

We present another way to arrive at the same analytical expressions using the replica method of statistical physics and a saddle-point approximation (Sherrington & Kirkpatrick, 1975; Mézard et al., 1987).

Analysis of our theoretical expressions show that different spectral modes of a target function are learned with different rates. Modes corresponding to higher kernel eigenvalues are learned faster, in the sense that a marginal training data point causes a greater percent reduction in generalization error for higher eigenvalue modes than for lower eigenvalue modes.

When data is sampled from a uniform distribution on a hypersphere, dot product kernels, which include NTK, admit a degenerate Mercer decomposition in spherical harmonics, $Y_{km}$. In this case, our theory predicts that generalization error of lower frequency modes of the target function decrease more quickly than higher frequency modes as the dataset size grows. Different learning stages are visible in the sense described below.

As the dimensions of data, $d$, go to infinity, learning curves exhibit different learning stages. For a training set of size $\ p\sim\mathcal{O}(d^{l})$, modes with $k<l$ are perfectly learned, $k=l$ are being learned, and $k>l$ are not learned.

We verify the predictions of our theory using numerical simulations for kernel regression and kernel interpolation with NTK, and wide and deep neural networks trained with gradient descent. Our theory fits experiments remarkably well on synthetic datasets and MNIST.

Our main approximation technique comes from the literature on Gaussian processes, which is related to kernel regression in a certain limit. Total learning curves for Gaussian processes, but not their spectral decomposition as we do here, have been studied in a limited teacher-student setting where both student and teacher were described by the same Gaussian process and the same noise in (Opper & Vivarelli, 1998; Sollich, 1999). We allow arbitrary teacher distributions. Sollich also considered mismatched models where teacher and student kernels had different eigenspectra and different noise levels (Sollich, 2002). The total learning curve from this model is consistent with our results when the teacher noise is sent to zero, but we also consider, provide expressions for, and analyze generalization in spectral modes. We use an analogue of the “lower-continuous” approximation scheme introduced in (Sollich & Halees, 2002), the results of which we reproduce through the replica method (Mézard et al., 1987).

Generalization bounds for kernel ridge regression have been obtained in many contexts (Schölkopf & Smola, 2001; Cucker & Smale, 2002; Vapnik, 1999; Gyorfi et al., 2003), but the rates of convergence often crucially depend on the explicit ridge parameter $\lambda$ and do not provide guarantees in the ridgeless case. Using a teacher-student setting, Spigler et al. (2019) showed that learning curves for kernel regression asymptotically decay with a power law determined by the decay rate of the teacher and the student. Such power law decays have been observed empirically on standard datasets (Hestness et al., 2017; Spigler et al., 2019). Recently, interest in explaining the phenomenon of interpolation has led to the study of generalization bounds on ridgeless regression (Belkin et al., 2018b, a, 2019b; Liang & Rakhlin, 2018). Here, our aim is to capture the average case performance of kernel regression, as opposed to a bound on it, that remains valid for the ridgeless case and finite sample sizes.

In statistical physics domain, Dietrich et al. (1999) calculated learning curves for support vector machines, but not kernel regression, in the limit of number of training samples going to infinity for dot product kernels with binary inputs using a replica method. Our theory applies to general kernels and finite size datasets. In the infinite training set limit, they observed several learning stages where each spectral mode is learned with a different rate. We observe similar phenomena in kernel regression. In a similar spirit, (Cohen et al., 2019) calculates learning curves for infinite-width neural networks using a path integral formulation and a replica analysis but does not discuss the spectral dependence of the generalization error.

In the infinite width limit, neural networks have many more parameters than training samples yet they do not overfit (Zhang et al., 2017). Some authors suggested that this is a consequence of the training procedure since stochastic gradient descent is implicitly biased towards choosing the simplest functions that interpolate the training data (Belkin et al., 2019a, 2018b; Xu et al., 2019a; Jacot et al., 2018). Other studies have shown that neural networks fit the low frequency components of the target before the high frequency components during training with gradient descent (Xu et al., 2019b; Rahaman et al., 2019; Zhang et al., 2019; Luo et al., 2019). In addition to training dynamics, recent works such as (Yang & Salman, 2019; Bietti & Mairal, 2019; Cao et al., 2019) have discussed how the spectrum of kernels impacts its smoothness and approximation properties. Here we explore similar ideas by explicitly calculating average case learning curves for kernel regression and studying its dependence on the kernel’s eigenspectrum.

Kernel Regression Learning Curves

We start with a general theory of kernel regression. Implications of our theory for dot product kernels including NTK and trained neural networks are described in Section 3.

We start by defining our notation and setting up our problem. Our initial goal is to derive a mathematical expression for generalization error in kernel regression, which we will analyze in the subsequent sections using techniques from the Gaussian process literature (Sollich, 1999, 2002; Sollich & Halees, 2002) and statistical physics (Sherrington & Kirkpatrick, 1975; Mézard et al., 1987).

The $\lambda\to 0$ limit is referred to as interpolating kernel regression, and, as we will discuss later, relevant to training wide neural networks. The unique minimum of the convex optimization problem is given by

Let $p(\mathbf{x})$ be the probability density function from which the input data are sampled. The generalization error is defined as the expected risk with expectation taken over new test points sampled from the same density $p(\mathbf{x})$. For a given dataset $\{\mathbf{x}_{i}\}$ and target function $f^{*}(\mathbf{x})$, let $f_{K}(\mathbf{x};\{\mathbf{x}_{i}\},f^{*})$ represent the function learned with kernel regression. The generalization error for this dataset and target function is

To calculate the average case performance of kernel regression, we average this generalization error over the possible datasets $\{\mathbf{x}_{i}\}$ and target functions $f^{*}$

Our aim is to calculate $E_{g}$ for a general kernel and a general distribution over teacher functions.

For our theory, we will find it convenient to work with the feature map defined by the Mercer decomposition. By Mercer’s theorem (Mercer, 1909; Rasmussen & Williams, 2005) , the kernel admits a representation in terms of its $M$ kernel eigenfunctions $\{\phi_{\rho}(\mathbf{x})\}$,

where $\psi_{\rho}(\mathbf{x})=\sqrt{\lambda_{\rho}}\phi(\mathbf{x})$ is the feature map we will work with. In our analysis, $M$ will be taken to be infinite, but for the derivation of the learning curves, we will first consider $M$ as a finite integer. The eigenfunctions and eigenvalues are defined with respect to the probability measure that generates the data $d\mu(\mathbf{x})=p(\mathbf{x})d\mathbf{x}$

We will also find it convenient to work with a vector representation of the RKHS functions in the feature space. Kernel eigenfunctions form a complete orthonormal basis, allowing the expansion of the target function $f^{*}$ and learned function $f$ in terms of features $\{\psi_{\rho}(\mathbf{x})\}$

Hence, $M$-dimensional vectors $\mathbf{w}$ and $\mathbf{\overline{w}}$ constitute a representation of $f$ and $f^{*}$ respectively in the feature space.

Another novelty of our theory is the decomposition of the generalization error into its contributions from different eigenmodes. The feature space expression of the generalization error after averaging over the data distribution can be written as:

where we identify $E_{\rho}$ as the generalization error in mode $\rho$.

We introduce a matrix notation for RKHS eigenvalues $\mathbf{\Lambda}_{\rho,\gamma}\equiv\delta_{\rho,\gamma}\lambda_{\rho}$ for convenience. Finally, with our notation set up, we can present our first result about generalization error.

For the $\mathbf{w}$ that minimizes the training error (eq. (8)), the generalization error (eq. (4)) is given by

which can be decomposed into modal generalization errors

We leave the proof to SI Section 2 but provide a few cursory observations of this result. First, note that all of the dependence on the teacher function comes in the matrix $\mathbf{D}$ whereas all of the dependence on the empirical samples is in $\mathbf{G}$. In the rest of the paper, we will develop multiple theoretical methods to calculate the generalization error given by expression (11).

Averaging over the target weights in the expression for $\mathbf{D}$ is easily done for generic weight distributions. The case of a fixed target is included by choosing a delta-function distribution over $\mathbf{\overline{w}}$.

We present two methods for computing the nontrivial average of the matrix $\mathbf{G}^{2}$ over the training samples $\{\mathbf{x}_{i}\}$. First, we consider the effect of adding a single new sample to $\mathbf{G}$ to derive a recurrence relation for $\mathbf{G}$ at different number of data points. This method generates a partial differential equation that must be solved to compute the generalization error. Second, we use a replica method and a saddle point approximation to calculate the matrix elements of $\mathbf{G}$. These approaches give identical predictions for the learning curves of kernel machines.

For notational simplicity, in the rest of the paper, we will use $\braket{\ldots}$ to mean $\braket{\ldots}_{{\{\mathbf{x}_{i}\},\mathbf{\overline{w}}}}$ unless stated otherwise. In all cases, the quantity inside the brackets will depend either on the data distribution or the distribution of target weights, but not both.

2 Continuous Approximation to Learning Curves

First, we adopt a method following Sollich (1999, 2002) and Sollich & Halees (2002) to calculate the generalization error. We generalize the definition of $\mathbf{G}$ by introducing an auxiliary parameter $v$, and make explicit its dataset size, $p$, dependence:

Note that the quantity we want to calculate is given by

By considering the effect of adding a single randomly sampled input $\mathbf{x^{\prime}}$, and treating $p$ as a continuous parameter, we can derive an approximate quasi-linear partial differential equation (PDE) for the average elements of $\mathbf{G}$ as a function of the number of data points $p$ (see below for a derivation):

Performing the average of the last term on the right hand side is difficult so we resort to an approximation, where the numerator and denominator are averaged separately.

where we used the fact that $\left<\phi_{\rho}(\mathbf{x^{\prime}})\phi_{\gamma}(\mathbf{x}^{\prime})\right>_{\mathbf{x}^{\prime}\sim p(\mathbf{x}^{\prime})}=\delta_{\rho,\gamma}$.

Treating $p$ as a continuous variable and taking a continuum limit of the finite differences given above, we arrive at (17). ∎

Next, we present the solution to the PDE (17) and the resulting generalization error.

This implicit solution is obtained from the method of characteristics which we provide in Section 3 of the SI.

Under the PDE approximation (17), the average error $E_{\rho}$ associated with mode $\rho$ is

where $t(p)\equiv\sum_{\rho}g_{\rho}(p,0)$ is the solution to the implicit equation

The full proof of this proposition is provided in Section 3 of the SI. We show the steps required to compute theoretical learning curves numerically in Algorithm 1.

In eq. (21), the target function sets the overall scale of $E_{\rho}$. That $E_{\rho}$ depends only on $\bar{w}_{\rho}$, but not other target modes, is an artifact of our approximation scheme, and in a full treatment may not necessarily hold. The spectrum of the kernel affects all modes in a nontrivial way. When we apply this theory to neural networks in Section 3, the information about the architecture of the network will be in the spectrum $\{\lambda_{\rho}\}$. The dependence on number of samples $p$ is also nontrivial, but we will consider various informative limits below.

We note that though the mode errors fall asymptotically like $p^{-2}$ (SI Section 4), the total generalization error $E_{g}$ can scale with $p$ in a nontrivial manner. For instance, if $\overline{w}^{2}_{\rho}\lambda_{\rho}\sim{\rho}^{-a}$ and $\lambda_{\rho}\sim{\rho}^{-b}$ then a simple computation (SI Section 4) shows that $E_{g}\sim p^{-\min\{a-1,2b\}}$ as $p\to\infty$ for ridgeless regression and $E_{g}\sim p^{-\min\{a-1,2b\}/b}$ for explicitly regularized regression. This is consistent with recent observations that total generalization error for neural networks and kernel regression falls in a power law $E_{g}\sim p^{-\beta}$ with $\beta$ dependent on kernel and target function (Hestness et al., 2017; Spigler et al., 2019).

3 Computing Learning Curves with Replica Method

The result of the continuous approximation can be obtained using another approximation method, which we outline here and detail in SI Section 5. We perform the average of matrix $\mathbf{G}(p,v)$ over the training data, using the replica method (Sherrington & Kirkpatrick, 1975; Mézard et al., 1987) from statistical physics and a finite size saddle-point approximation, and obtain identical learning curves to Proposition 3. Our starting point is a Gaussian integral representation of the matrix inverse

where $Z=\int d\mathbf{u}\ e^{-\frac{1}{2}\mathbf{u}^{\top}(\frac{1}{\lambda}\mathbf{\Phi}\mathbf{\Phi}^{\top}+\mathbf{\Lambda}^{-1}+v\mathbf{I})\mathbf{u}}$. Since $Z$ also depends on the dataset (quenched disorder) $\mathbf{\Phi}$, to make the average over $\mathbf{\Phi}$ tractable, we use the following limiting procedure: $Z^{-1}=\lim_{n\to 0}Z^{n-1}$. As is common in the physics of disordered systems (Mézard et al., 1987), we compute $R(p,v,\mathbf{h})$ for integer $n$ and analytically continue the expressions in the $n\to 0$ limit under a symmetry ansatz. This procedure produces the same average matrix elements as the continuous approximation discussed in Proposition 2, and therefore the same generalization error given in Proposition 3. Further detail is provided in SI Section 5.

4 Spectral Dependency of Learning Curves

We can get insight about the behavior of learning curves by considering ratios between errors in different modes:

For small $p$ this ratio approaches $\frac{E_{\rho}}{E_{\gamma}}\sim\frac{\lambda_{\rho}\braket{\overline{w}_{\rho}^{2}}}{\lambda_{\gamma}\braket{\overline{w}_{\gamma}^{2}}}$. For large $p$, $\frac{E_{\rho}}{E_{\gamma}}\sim\frac{\braket{\overline{w}_{\rho}^{2}}/\lambda_{\rho}}{\braket{\overline{w}_{\gamma}^{2}}/\lambda_{\gamma}}$, indicating that asymptotically ($p\to\infty$), the amount of relative error in mode $\rho$ grows with the ratio $\braket{\overline{w}_{\rho}^{2}}/\lambda_{\rho}$, showing that the asymptotic mode error is relatively large if the teacher function places large amounts of power in modes that have small RKHS eigenvalues $\lambda_{\rho}$.

We can also examine how the RKHS spectrum affects the evolution of the error ratios with $p$. Without loss of generality, we take $\lambda_{\gamma}>\lambda_{\rho}$ and show in SI Section 6 that

In this sense, the marginal training data point causes a greater percent reduction in generalization error for modes with larger RKHS eigenvalues.

5 Multiple Outputs

The solution to the learning problem depends on the same kernel but different targets for each function:

Our theory can be used to generate predictions for the generalization error of each of the $C$ learned functions, $f_{c}(\mathbf{x})$, and then summed to obtain the total error.

Here, $N(d,k)$ is the dimension of the subspace spanned by $d$-dimensional spherical harmonics of degree $k$. Rotation invariance renders the eigenspectrum degenerate since each of the $N(d,k)$ modes of frequency $k$ share the same eigenvalue $\lambda_{k}$. A review of these topics is given in SI Sections 7 and 8.

We briefly comment on another fact that will later be used in our numerical simulations. Dot product kernels admit an expansion in terms of Gegenbauer polynomials $\{Q_{k}\}$, which form a complete and orthonormal basis for the uniform measure on the sphere (Dai & Xu, 2013): $\kappa(z)=\sum_{k=0}^{\infty}\lambda_{k}N(d,k)Q_{k}(z)$. The Gegenbauer polynomials are related to spherical harmonics $\{Y_{km}\}$ through $Q_{k}(\mathbf{x}^{\top}\mathbf{x}^{\prime})=\frac{1}{N(d,k)}\sum_{m=1}^{N(d,k)}Y_{km}(\mathbf{x})Y_{km}(\mathbf{x}^{\prime})$ (Dai & Xu, 2013) (see SI Sections 7 and 8 for a review).

In the special case of dot product kernels with monotonically decaying spectra, results given in Section 2.4 indicate that the marginal training data point causes greater reduction in relative error for low frequency modes than for high frequency modes. Monotonic RKHS spectra represent an inductive bias that preferentially favors fitting lower frequencies as more data becomes available. More rapid decay in the spectrum yields a stronger bias to fit low frequencies first.

where the constant is given in SI Section 9. In other words, $k<l$ modes are perfectly learned, $k=l$ are being learned with an asymptotic $1/\alpha^{2}$ rate, and $k>l$ are not learned.

This simple calculation demonstrates that the lower modes are learned earlier with increasing sample complexity since the higher modes stays stationary until $p$ reaches to the degeneracy of that mode.

2 Neural Tangent Kernel and its Spectrum

For fully connected architectures, the NTK is a rotation invariant kernel that describes how the predictions of infinitely wide neural networks evolve under gradient flow (Jacot et al., 2018). Let $\theta_{i}$ index all of the parameters of the neural network and let $f_{\theta}(\mathbf{x})$ be the output of the network. Here, we focus on scalar network outputs for simplicity, but generalization to multiple outputs is straightforward, as discussed in Section 2.5. Then the neural tangent kernel is defined as

Note that this corresponds to ridgeless, interpolating regression where $\lambda=0$. We will use this correspondence and our kernel regression theory to explain neural network learning curves in the next section. For more information about NTK for fully connected architectures see SI Sections 10 and 11.

Experiments

In this section, we test our theoretical results for kernel regression, kernel interpolation and wide networks for various kernels and datasets.

We first test our theory in a kernel regression task with NTK demonstrating the spectral decomposition. In this experiment, the target function is a linear combination of a kernel evaluated at randomly sampled points $\{\mathbf{\overline{x}}_{i}\}$:

Figure 2 shows the errors for each frequency $k$ as a function of sample size $p$. In Figure 2(a), we show that the mode errors sequentially start falling when $p\sim N(d,k)$. Figure 2(b) shows the mode error corresponding to $k=1$ for kernel regression with 3-layer NTK across different dimensions. Higher input dimension causes the frequency modes to be learned at larger $p$. We observe an asymptotic $\sim 1/\alpha^{2}$ decay in modal errors. Finally, we show the effect of regularization on mode errors with a 10-layer NTK in Figure 2(c). With increasing $\lambda$, learning begins at larger $p$ values.

2 Learning Curves for Finite Width Neural Networks

Having established that our theory accurately predicts the generalization error of kernel regression with NTK, we now compare the generalization error of finite width neural networks trained on a quadratic loss with the theoretical learning curves for NTK. For these experiments, we use the Neural-Tangents Library (Novak et al., 2020) which supports training and inference for both finite and infinite width neural networks.

First, we use “pure mode” teacher functions, meaning the teacher is composed only of spherical harmonics of the same degree. For “pure mode” $k$, the teacher is constructed with the following rule:

where again $\overline{\alpha}_{i}\sim\mathcal{B}(1/2)$ and $\mathbf{\overline{x}}_{i}\sim p(\mathbf{x})$ are sampled randomly. Figure 3(a) shows the learning curve for a fully connected 2-layer ReLU network with width $N=10000$, input dimension $d=30$ and $p^{\prime}=10000$. As before, we see that the lower $k$ pure modes require less data to be fit. Experimental test errors for kernel regression with NTK on the same synthetic datasets are plotted as triangles. Our theory perfectly fits the experiments.

Results from a 4-layer NN simulation are provided in Figure 3(b). Each hidden layer had $N=500$ hidden units. We again see that the $k=2$ mode is only learned for $p>200$. $k=4$ mode is not learned at all in this range. Our theory again perfectly fits the experiments.

Lastly, we show that our theory also works for composite functions that contain many different degree spherical harmonics. In this setup, we randomly initialize a two layer teacher neural network and train a student neural network

3 Gaussian Kernel Regression and Interpolation

4 MNIST: Discrete Data Measure and Kernel PCA

Conclusion

In this paper, we presented an approximate theory of the average generalization performance for kernel regression. We studied our theory in the ridgeless limit to explain the behavior of trained neural networks in the infinite width limit (Jacot et al., 2018; Arora et al., 2019; Lee et al., 2019). We demonstrated how the RKHS eigenspectrum of NTK encodes a preferential bias to learn high spectral modes only after the sample size $p$ is sufficiently large. Our theory fits kernel regression experiments remarkably well. We further experimentally verified that the theoretical learning curves obtained in the infinite width limit provide a good approximation of the learning curves for wide but finite-width neural networks. Our MNIST result suggests that our theory can be applied to datasets with practical value.

Acknowledgements

We thank Matthieu Wyart and Stefano Spigler for comments and pointing to a recent version of their paper (Spigler et al., 2019) with an independent derivation of the generalization error scaling for power law kernel and target spectra (see (4)). C. Pehlevan thanks the Harvard Data Science Initiative, Google and Intel for support.

References

Background on Kernel Machines

If such a kernel exists for a Hilbert space, then it is unique and defined as the reproducing kernel for the RKHS (Evgeniou et al., 1999; Schölkopf & Smola, 2001).

2 Mercer’s Theorem

Let $\mathcal{H}$ be a RKHS with kernel $K$. Mercer’s theorem (Mercer, 1909; Rasmussen & Williams, 2005) allows the eigendecomposition of $K$

3 Representer Theorem

Let $\mathcal{H}$ be a RKHS with inner product $\left<.,.\right>_{\mathcal{H}}$. Consider the regularized learning problem

where $\hat{\mathcal{L}}[f]$ is an empirical cost defined on the discrete support of the dataset and $\lambda>0$. The optimal solution to the optimization problem above can always be written as (Schölkopf & Smola, 2001)

4 Solution to Least Squares

Specializing to the case of least squares regression, let

Using the representer theorem, we may reformulate the entire objective in terms of the $p$ coefficients $\alpha_{i}$

Optimizing this loss with respect to $\bm{\alpha}$ gives

Therefore the optimal function evaluated at a test point is

Derivation of the Generalization Error

Define the student’s eigenfunction expansion $f(\mathbf{x})=\sum_{\rho}w_{\rho}\psi_{\rho}(\mathbf{x})$ and decompose the risk in the basis of eigenfunctions:

Next, it suffices to calculate the weights $\mathbf{w}$ learned through kernel regression. Define a matrix with elements $\mathbf{\Psi}_{\rho,i}=\psi_{\rho}(\mathbf{x}_{i})$. The training error for kernel regression is

since $\left<\psi_{\rho}(\cdot),\psi_{\gamma}(\cdot)\right>_{\mathcal{H}}=\delta_{\rho,\gamma}$ (Bietti & Mairal, 2019). This fact can be verified by invoking the reproducing property of the kernel and it’s Mercer decomposition. Let $g(\cdot)=\sum_{\rho}a_{\rho}\psi_{\rho}(\cdot)$. By the reproducing property

Due to the arbitrariness of $a_{\rho}$, we must have $\left<\psi_{\rho}(\cdot),\psi_{\gamma}(\cdot)\right>_{\mathcal{H}}=\delta_{\rho,\gamma}$. We stress the difference between the action of the Hilbert inner product and averaging feature functions over a dataset $\left<\psi_{\rho}(\mathbf{x})\psi_{\gamma}(\mathbf{x})\right>_{\mathbf{x}}=\lambda_{\rho}\delta_{\rho,\gamma}$ which produce different results. We will always decorate angular brackets with $\mathcal{H}$ to denote Hilbert inner product.

where the target function is produced according to $\mathbf{y}=\mathbf{\Psi}^{\top}\mathbf{\overline{w}}$.

Plugging in the $\mathbf{w}$ that minimizes the training error into the formula for the generalization error, we find

and identifying the terms in (SI.19) with these definitions, we obtain the desired result. Then each component of the mode error is given by:

Solution of the PDE Using Method of Characteristics

Here we derive the solution to the PDE in equation 17 of the main text by adapting the method used by (Sollich, 1999). We will prove both Propositions 2 and 3.

It follows from equation 17 that $t$ obeys the PDE

with an initial condition $t(0,v)=\text{Tr}(\mathbf{\Lambda}^{-1}+v\mathbf{I})^{-1}$. The solution to first order PDEs of the form is given by the method of characteristics (Arfken, 1985), which we describe below, and prove Proposition 2.

One such normal vector is $(-1,\frac{\partial t}{\partial p},\frac{\partial t}{\partial v})$.

The PDE can be written as a dot product involving this normal vector,

The first of these equations indicate that $t$ is constant along each characteristic curve. Integrating along the parameter, $p=s+p_{0}$ and $v=-\frac{s}{\lambda+t}+v_{0}$ where $p_{0}$ is the value of $p$ when $s=0$ and $v_{0}$ is the value of $v$ at $s=0$. Without loss of generality, take $p_{0}=0$ so that $s=p$. At $s=0$, we have our initial condition

Since $t$ takes on the same value for each characteristic

which gives an implicit solution for $t(p,v)$. Now that we have solved for $t(p,v)$, remembering (SI.24), we may write

This equation proves Proposition 2 of the main text. ∎

Next, we compute the modal generalization errors $E_{\rho}$ and prove Proposition 3.

Computing generalization error of kernel regression requires the differentiation with respect to $v$ at $v=0$ (eq.s (11) and (16) of main text). Since $\left<\mathbf{G}^{2}\right>$ is diagonal, the mode errors only depend on the diagonals of $\mathbf{D}$ and on $\left<\mathbf{G}_{\rho,\rho}^{2}\right>=-\frac{\partial g_{\rho}}{\partial v}|_{v=0}$:

We proceed with calculating the derivative in the above equation.

We need to calculate $\frac{\partial t(p,v)}{\partial v}|_{v=0}$

so it suffices to numerically solve for $t(p,0)$ to recover predictions of the mode errors. Equations (SI.29) (evaluated at $v=0$), (SI.34) and (SI.37) collectively prove Proposition 3. ∎

Learning Curve for Power Law Spectra

For $\lambda>0$, the mode errors asymptotically satisfy $E_{\rho}\sim\mathcal{O}(p^{-2})$ since $\frac{p}{\lambda+t}\sim\frac{p}{\lambda}$ and $\frac{(\lambda+t)^{2}}{(\lambda+t)^{2}-\gamma p}\sim\mathcal{O}_{p}(1)$ (see below). Although each mode error decays asymptotically like $p^{-2}$, the total generalization error can have nontrivial scaling with $p$ that depends on both the kernel and the target function.

To illustrate the dependence of the learning curves on the choice of kernel and target function, we consider a case where both have power law spectra. Specifically, we assume that $\lambda_{\rho}=\rho^{-b}$ and $a_{\rho}^{2}\equiv\overline{w}_{\rho}^{2}\lambda_{\rho}=\rho^{-a}$ for $\rho=1,2,...$. We introduce the variable $z=t+\lambda$ to simplify the computations below. We further approximate the sums over modes with integrals

We use the same approximation technique to study the behavior of $z(p)$

where $F(b,p,z)=\int_{(z/p)^{1/b}}^{\infty}\frac{du}{1+u^{b}}$. If $p\gg\lambda^{-1/(b-1)}$ then $z\approx\lambda$, otherwise $z\approx p^{1-b}F(b,p,z)^{b}$. Further, the scaling $z\sim\mathcal{O}(p^{1-b})$ is self-consistent since the lower endpoint of integration $(z/p)^{1/b}\sim p^{-1}\to 0$ so $F(b,z,p)$ approaches a constant $F(b)$ for $p\to\infty$

We similarly find that $p\gamma(p)\sim\mathcal{O}(p^{2-2b})$ if $p\ll\lambda^{-1/(b-1)}$. The mode-independent prefactor is approximately constant $\frac{z^{2}}{z^{2}-\gamma p}\sim\mathcal{O}_{p}(1)$.

We can use all of these facts to identify scalings of $E_{g}$. We will first consider the case where $p\ll\lambda^{-1/(b-1)}$:

If $2b>a-1$ then the second term dominates, indicating that higher frequency modes $k>p$ provide a greater contribution to the error due to the slow decay in the target power. In this case $E_{g}\sim p^{-(a-1)}$. If, on the other hand, $2b<a-1$ then lower frequency modes $k<p$ dominate the error and $E_{g}\sim p^{-2b}$.

Now, suppose that $p>\lambda^{-1/(b-1)}$. In this regime

Here there are two possible scalings. If $2b>a-1$ then $E_{g}\sim p^{-(a-1)/b}$ while $2b<a-1$ implies $E_{g}\sim p^{-2}$.

A verification of this scaling is provided in Figure SI.1, which shows the behavior of $z$ and $E_{g}$ in these two regimes. When the explicit regularization is low (or zero) ($p<\lambda^{-1/(b-1)}$), our equations reproduce the power law scalings derived with Fourier analysis in (Spigler et al., 2019)We note that in a recent version of their paper, Spigler et al. (2019) used our formalism to independently derive the scalings in (4) for the ridgeless ($\lambda=0$) case. Our calculation in an earlier preprint had missed the possible $\sim p^{-2b}$ and $\sim p^{-2}$ scalings, which we corrected after their paper..

The slower asymptotic decays in generalization error when explicit regularization $\lambda$ is large relative to the sample size indicates that explicit regularization hurts performance. The decay exponents also indicate that the RKHS eigenspectrum should decay with exponent at least as large as $b^{*}>\frac{a-1}{2}$ for optimal asymptotics. Kernels with slow decays in their RKHS spectra induce larger errors.

Replica Calculation

In this section, we present the replica trick and the saddle-point approximation summarized in main text Section 2.3. Our goal is to show that the continuous approximation of the main paper and previous section can be interpreted as a finite size saddle-point approximation to the replicated system under a replica symmetry ansatz. We will present a detailed treatment of the thermodynamic limit and the replica symmetric ansatz in a different paper.

and make use of the identity $Z^{-1}=\lim_{n\to 0}Z^{n-1}$ to rewrite the entire average in the form

Following the replica method from the physics of disordered systems, we will first restrict ourselves to integer $n$ and then analytically continue the resulting expressions to take the limit of $n\to 0$.

Averaging over the quenched disorder (dataset) with the assumption that the residual error $(\mathbf{w}-\mathbf{\overline{w}})\cdot\mathbf{\Psi}(\mathbf{x}_{i})$ is a Gaussian process, we find

where order parameters $Q_{ab}=\mathbf{u}^{a}\cdot\mathbf{u}^{b}$ have been introduced.

To enforce the definition of these order parameters, Dirac delta functions are inserted into the expression for $R$. We then represent each delta function as a Fourier integral so that integrals over $\mathbf{u}^{a}$ can be computed

After inserting delta functions to enforce order parameter definitions, we are left with integrals over the thermal degrees of freedom

We now make a replica symmetric ansatz $Q_{ab}=q\delta_{ab}+q_{0}$ and $2i\hat{Q}_{ab}=\hat{q}\delta_{ab}+\hat{q}_{0}$. Under this ansatz $R(\mathbf{h})$ can be rewritten as

In the limit $p\to\infty$, $R(\mathbf{h})$ is dominated by the saddle point of the free energy where $\nabla\mathcal{F}(q,\hat{q},q_{0},\hat{q}_{0})=0$. The saddle point equations are

We see that $q^{*}$ is exactly equivalent to $t(p,v)$ defined in SI.29 for the continuous approximation. Under the saddle point approximation we find

Taking the $n\to 0$ limit as promised, we obtain the normalized average

Using our formula for the mode errors, we find

consistent with our result from the continuous approximation.

Spectral Dependence of Learning Curves

We want to calculate how different mode errors change as we add one more sample. We study:

where $E_{\rho}$ is given by eq. (21). Evaluating the derivative, we find:

where we identified the sum with $\gamma$. Inserting this, we obtain:

Finally, solving for $\partial t/\partial p$ from (6), we get:

proving that $\partial t/\partial p<0$. Taking $\lambda_{\gamma}>\lambda_{\rho}$ without loss of generality, it follows that

Spherical Harmonics

The Laplace Beltrami Operator can be decomposed into the radial and angular parts, allowing

Using this decomposition, the spherical harmonics are eigenfunctions of the surface Laplacian

Let $K(\mathbf{x},\mathbf{x}^{\prime})=\kappa(\mathbf{x}^{\top}\mathbf{x}^{\prime})$. The kernel’s orthogonal decomposition is

To numerically calculate the kernel eigenvalues of $\kappa$, we use Gauss-Gegenbauer quadrature (Abramowitz & Stegun, 1972) for the measure $d\tau(z)$ so that for a quadrature scheme of order $r$

where $z_{i}$ are the $r$ roots of $Q_{r}(z)$ and the weights $w_{i}$ are chosen with

Frequency Dependence of Learning Curves in d→∞→𝑑d\to\infty Limit

Here, we consider an informative limit where the number of input data dimension, $d$, goes to infinity.

Denoting the index $\rho=(k,m)$, we can write mode error (SI.37), after some rearranging, as:

where $t$ and $\gamma$, after performing the sum over degenerate indices, are:

In the limit $d\to\infty$, the degeneracy factor (SI.64) approaches to $N(d,k)\sim\mathcal{O}(d^{k})$. We note that for dot-product kernels $\lambda_{k}$ scales with $d$ as $\lambda_{k}\sim d^{-k}$ (Smola et al., 2001) (Figure 1), which leads us to define the $\mathcal{O}(1)$ parameter $\bar{\lambda}_{k}=d^{k}\lambda_{k}$. Plugging these in, we get:

where $g_{k}=p/d^{k}$ is the ratio of sample size to the degeneracy. Furthermore, we want to calculate the ratio $E_{km}(p)/E_{km}(0)$ to probe how much the mode errors move from their initial value:

Let us consider an integer $l$ such that the scaling $P=\alpha d^{l}$ holds. This leads to three different asymptotic behavior of $g_{k}$s:

If we assume $t\sim\mathcal{O}(1)$, we get an asymptotically consistent set of equations:

Then using (SI.75), (9) and (9), we find the errors associated to different modes as:

Neural Tangent Kernel

For a neural network, it is convenient to compute this recursively in terms of the Neural Network Gaussian Process (NNGP) kernel which corresponds to only training the readout weights from the final layer (Jacot et al., 2018; Arora et al., 2019). We will restrict our attention to networks with zero bias and nonlinear activation function $\sigma$. Then

If $\sigma$ is chosen to be the ReLU activation, then we can analytically simplify the expression. Defining the following function

where $f^{\circ(L-1)}(z)$ is the function $f$ composed into itself $L-1$ times.

This simplification gives an exact recursive formula to compute the kernel as a function of $z=\mathbf{x}^{\top}\mathbf{x}^{\prime}$, which is what we use to compute the eigenspectrum with the quadrature scheme described in the previous section.

Spectra of Fully Connected ReLU NTK

By the representer theorem, let $f(\mathbf{x})=\sum_{i=1}^{p}\alpha_{i}K(\mathbf{x},\mathbf{x}_{i})$. By Green’s theorem, the variance of the $n$-th derivative can be rewritten as

where $\alpha^{*}=\max_{j}|\alpha_{j}|$ and $|Q_{k}(z)|\leq CN(d,k)$ for a universal constant $C$. A sufficient condition for this sum to converge is that $\lambda_{k}^{2}k^{n}(k+d-2)^{n}N(d,k)^{2}\sim\mathcal{O}(k^{-1})$ which is equivalent to demanding $\lambda_{k}N(d,k)\sim\mathcal{O}(k^{-n-1/2})$ since $(k+d-2)^{n}\sim k^{n}$ as $k\to\infty$. ∎

Decomposition of Risk for Numerical Experiments

As we describe in Section 4.1 of the main text, the teacher functions for the kernel regression experiments are chosen as

where the coefficients $\overline{\alpha_{i}}\sim\mathcal{B}(1/2)$ are randomly sampled from a centered Bernoulli distribution on $\{\pm 1\}$ and the points $\mathbf{\overline{x}}_{i}\sim p(\mathbf{x})$ are drawn from the same distribution as the training data. In general $p^{\prime}$ is not the same as the number of samples $p$. Choosing a function of this form is very convenient for producing theoretical predictions of mode errors as we discuss below.

Since the matrix elements $\left<\mathbf{G}_{\rho\rho}^{2}\right>$ are determined completely by the kernel eigenvalues $\{\lambda_{\rho}\}$, it suffices to calculate the diagonal elements of $\mathbf{D}$ to find the generalization error. For the teacher function sampled in the way described above, there is a convenient expression for $\mathbf{D}_{\rho\rho}$.

The teacher function admits an expansion in the basis of kernel eigenfunctions

Using the Mercer decomposition of the kernel we can identify the coefficients

Comparing each term in these two expressions, we identify the coefficient of the $\rho$-th eigenfunction

We now need to compute the $D_{\rho\rho}$, by averaging $\overline{w}_{\rho}^{2}$ over all possible teachers

since $\left<\psi_{\rho}(\mathbf{x})\psi_{\rho}(\mathbf{x})\right>=\lambda_{\rho}$. Thus it suffices to calculate $\frac{\partial}{\partial v}g_{\rho}(p,v)$ for each mode and then compute mode errors with

where $\frac{\partial g_{\rho}}{\partial v}|_{v=0}$ is evaluated in terms of the numerical solution for $t(p,0)$.

2 Empirical Mode Errors

By the representer theorem, we may represent the student function as $f(\mathbf{x})=\sum_{i=1}^{P}\alpha_{i}K(\mathbf{x},\mathbf{x}_{i})$. Then, the generalization error is given by

On the $d$-sphere, by defining $E_{k}=\sum_{m=1}^{N(d,k)}E_{km}$ we arrive at the formula

We randomly sample the $\overline{\alpha}$ variables for the teacher and fit $\alpha=(\mathbf{K}+\lambda\mathbf{I})^{-1}\mathbf{y}$ to the training data. Once these coefficients are known, we can obtain empirical mode errors.

Neural Network Experiments

For the “pure mode” experiments with neural networks, the target function was

whereas, for the composite experiment, the target function was a randomly sampled two layer neural network with ReLU activations

This target model is a special case of eq. (SI.90) so the same technology can be used to compute the theoretical learning curves. We can use a similar trick as that shown in equation (12.1) to determine $\overline{w}_{\rho}$ for the NN teacher experiment. Let the Gegenbauer polynomial expansion of $\sigma(z)$ be $\sigma(z)=\sum_{k=0}^{\infty}a_{k}N(d,k)Q_{k}(z)$. Then the mode error for mode $k$ is $E_{k}=\frac{a_{k}^{2}}{\lambda_{k}^{2}}\left<g_{k}^{2}\right>$ where $\left<g_{k}^{2}\right>$ is computed with equation (SI.37).

A sample of some training error and generalization errors from pure mode experiments are provided below in Figures SI.3 and SI.4.

The choice of the number of hidden units $N$ was based primarily on computational considerations. For two layer neural networks, the total number of parameters scales linearly with $N$, so to approach the overparameterized regime, we aimed to have $N\approx 10p_{max}$ where $p_{max}$ is the largest sample size used in our experiment. For $p_{max}=500$, we chose $N=4000,10000$.

For the three and four layer networks, the number of parameters scales quadratically with $N$, making simulations with $N>10^{3}$ computationally expensive. We chose $N$ to give comparable training time for the 2 layer case which corresponded to $N=500$ after experimenting with $\{100,250,500,1000,5000\}$.

We found that the learning rate needed to be quite large for the training loss to be reduced by a factor of $\approx 10^{6}$. For the 2 layer networks, we tried learning rates $\{10^{-3},10^{-2},1,10,32\}$ and found that a learning rate of 32 gave the lowest training error. For the three and four layer networks, we found that lower learning rates worked better and used learning rates in the range from $[0.5,3]$.

Discrete Measure and Kernel PCA

For this measure, the integral eigenvalue equation becomes

Evaluating $\mathbf{x}^{\prime}$ at each of the points $\mathbf{x}_{i}$ in the dataset yields a matrix equation. Let $\mathbf{\Phi}_{\rho,i}=\phi_{\rho}(\mathbf{x}_{i})$ and $\mathbf{\Lambda}_{\rho,\gamma}=\delta_{\rho,\gamma}\lambda_{\rho}$