Asymptotic learning curves of kernel methods: empirical data versus teacher–student paradigm

In supervised learning machines learn from a finite collection of n training data, and their generalization error is then evaluated on unseen data drawn from the same distribution. How many data are needed to learn a task is characterized by the learning curve relating generalization error to n. In various cases, the generalization error decays as a power law n^−β , with an exponent β that depends on both the data and the algorithm. In [1] β is reported for state-of-the-art (SOTA) deep neural networks for various tasks: in neural-machine translation β ≈ 0.3–0.36 (for fixed model size) or β ≈ 0.13 (for best-fit models at any n); language modeling shows β ≈ 0.06–0.09; in speech recognition β ≈ 0.3; SOTA models for image classification (on ImageNet) have exponents β ≈ 0.3–0.5. Currently there is no available theory of deep learning to rationalize these observations. Recently it was shown that for a proper initialization of the weights, deep learning in the infinite-width limit [2] converges to kernel learning. Moreover, it is nowadays part of the lore that there exist kernels whose performance is nearly comparable to deep networks [3, 4], at least for some tasks. It is thus of great interest to understand the learning curves of kernels. For regression, if the target function being learned is simply assumed to be Lipschitz, then the best guarantee is β = 1/d [5, 6] where d is the data dimension. Therefore, for large d, β is very small: learning is completely inefficient, a phenomenon referred to as the curse of dimensionality. As a result, various works on kernel regression make the much stronger assumption that the training points are sampled from a target function that belongs to the reproducing kernel Hilbert space (RKHS) of the kernel (for Gaussian r [7]). With this assumption β does not depend on d (for instance in [8] β = 1/2 is guaranteed). However, RKHS is a very strong assumption, which requires the smoothness of the target function to increase with d [6] (for Gaussian random fields see appendix H), which may not be realistic in large dimensions.

In section 3 we compute β empirically for kernel methods applied on MNIST and CIFAR10 datasets. We find β_MNIST ≈ 0.4 and β_CIFAR10 ≈ 0.1 respectively. Quite remarkably, we observe essentially the same exponents for regression and classification tasks, using either a Gaussian or a Laplace kernel. Thus, the exponents are not as small as 1/d (d = 784 for MNIST, d = 3072 for CIFAR10), but neither are they 1/2 as one would expect under the RKHS assumption. These facts call for frameworks in which assumptions on the smoothness of the data can be intermediary between Lipschitz and RKHS. Here we study such a framework for regression, in which the target function is assumed to be a Gaussian random field of zero mean with translation-invariant isotropic covariance ${K}_{\mathrm{T}}\left(\underline{x}\right)$. The data can equivalently be thought of as being synthesized by a 'teacher' kernel ${K}_{\mathrm{T}}\left(\underline{x}\right)$. Learning is performed with a 'student' kernel ${K}_{\mathrm{S}}\left(\underline{x}\right)$ that minimizes the mean-square error. In general ${K}_{\mathrm{T}}\left(\underline{x}\right)\ne {K}_{\mathrm{S}}\left(\underline{x}\right)$. In this set-up learning is very similar to a technique referred to as kriging, or Gaussian process regression, originally developed in the geostatistics community [9, 10].

To quantify learning, in section 4 we first perform numerical experiments for data points distributed uniformly at random on a hypersphere of varying dimension d, focusing on a Laplace kernel for the student, and considering a Laplace or Gaussian kernel for the teacher. We observe that in both cases β(d) is a decreasing function. In section 5, to derive β(d) we consider the simplified situation where the Gaussian random field is sampled at training points lying on a regular lattice. Building on the kriging literature [10], we show that β is controlled by the high-frequency scaling of both the teacher and student kernels: assuming that the Fourier transforms of the kernels decay as ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)={c}_{\mathrm{T}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)={c}_{\mathrm{S}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}\right)$, we obtain

$$\begin{equation}\beta =\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right).\end{equation} \tag{ 1 }$$

Importantly (i) equation (1) leads to a prediction for β(d) that accurately matches our numerical study for random training data points, leading to the conjecture that equation (1) holds in that case as well. We offer the following interpretation: ultimately, kernel methods are performing a local interpolation, the quality of which depends on the distance δ(n) between adjacent data points. δ(n) is asymptotically similar for random data or data sitting on a lattice. (ii) If the kernel K_S is not too sensitive to high frequencies, then learning is optimal as far as scaling is concerned and β = (α_T − d)/d. We will argue that the smoothness index s ≡ [(α_T − d)/2] characterizes the number of derivatives of the target function that are continuous. We thus recover the curse of dimensionality: s needs to be of order d to have non-vanishing β in large dimensions.

We show that in some regimes, the test error for Gaussian data is controlled by an exponent a describing how the coefficients of the true function in the eigenbasis of the kernel decay with rank. We estimate a by the kernel principal component analysis (kernel PCA) based on diagonalizing the Gram matrix. This measure yields a prediction for the learning curve exponent β that matches the numerical fit, with β_MNIST ≈ 0.36 and β_CIFAR10 ≈ 0.07. We show in appendix I using the recent formalism of [11], which does not assume Gaussianity but makes more technical assumptions, that the result of our theorem equation (1) is recovered, supporting further its validity for real data.

Finally, we discuss the following apparent paradox: β is significant for MNIST and CIFAR10, for which d is a priori very large, leading to a smoothness value s in the hundreds in both cases, which appears unrealistic. In section 7 the paradox is resolved by considering that real datasets actually live on lower-dimensional manifolds. As far as kernel learning is concerned, our findings support that the correct definition of dimension should be based on how the nearest-neighbors distance δ(n) scales with n: $\delta \left(n\right)\sim {n}^{-1/{d}_{\mathrm{e}\mathrm{f}\mathrm{f}}}$. Direct measurements of δ(n) support that MNIST and CIFAR10 live on manifolds of lower dimensions ${d}_{\mathrm{M}\mathrm{N}\mathrm{I}\mathrm{S}\mathrm{T}}^{\mathrm{e}\mathrm{ff}}\approx 15$ and ${d}_{\mathrm{C}\mathrm{I}\mathrm{F}\mathrm{A}\mathrm{R}\mathrm{10}}^{\mathrm{e}\mathrm{ff}}\approx 35$.

Part of the literature has investigated the problem of kernel regression from a different point of view, namely the optimal worst-case performance (see for instance [12–14]). The target function is not assumed to be generated by a Gaussian random field, but its regularity is controlled using a source condition that constrains the decay of its coefficients in the eigenbasis of the kernel. For uniform data distributions and isotropic kernels this is similar to controlling how the Fourier transform of the target function decays at high frequency. What we study in the present work is, on the contrary, the typical performance. Indeed, it turns out that both the worst-case and the typical learning curve decay as power laws, and the latter decays faster.

The teacher–student framework for kernel regression was previously introduced in [15, 16], where a formula for the learning curve was derived based on a few uncontrolled approximations. It is easy to show that their results match the predictions of our theorem, although in [16] the case where the performance is limited by the student (α_T − d > 2α_S) is ignored. More recently, [11] generalized this approach using similar approximations and extended it to kernel regression applied to any target function (or ensemble thereof). Kernel PCA on MNIST was used to support that these approximations hold well on real data. An asymptotic scaling relation between β and a was obtained; again the presence of other regimes was not noted. By contrast we perform an exact calculation for the asymptotic behavior of Gaussian data on a lattice. In appendix I we show that the two approaches are consistent and lead to the same asymptotic predictions for β.

Our set-up of teacher–student learning with kernels is also referred to as kriging, or Gaussian process regression, and it was originally developed in the geostatistics community [9]. In section 5 we present our theorem, that allows one to know the rate at which the test error decreases as we increase the number of training points, assumed to lie on a high-dimensional regular lattice. Similar results have been previously derived in the kriging literature [10] when sampling occurs on the regular lattice with the exception of the origin, where the inference is made. Here we propose an alternative derivation that some readers might find simpler. We also study a slightly different problem: instead of computing the test error when the inference is carried on at the origin, we compute the average error for a test point that lies at an arbitrary point, sampled uniformly at random and not necessarily on the lattice. Then, in what follows we show, via extensive numerical simulations, that such predictions are accurate even when the training points do not lie on a regular lattice, but are taken at random on a hypersphere. An exact proof of our result in such a general setting is difficult and cannot be found even in the kriging literature. To our knowledge the results that get closer to the point are those discussed in [17], where the author studies one-dimensional processes where the training data are not necessarily evenly spaced.

In this work the effective dimension of the data plays an import role, as it controls how the distance between nearest neighbors scales with the dataset size. Of course, there exists a vast literature [18–24] devoted to the study of effective dimensions, where other definitions are analyzed. The effective dimensions that we find are compatible with those obtained with more refined methods.

In what follows we apply kernel methods to the MNIST and CIFAR10 datasets, each consisting of a set of images ${\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}$. We simplify the problem by considering only two classes whose labels $Z\left({\underline{x}}_{\mu }\right)={\pm}1$ correspond to odd and even numbers for MNIST, and to two groups of five classes in CIFAR10. The goal is to infer the value of the label ${\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)$ of an image $\underline{x}$ that does not belong to the dataset. The S subscript reminds us that inference is performed using a positive definite kernel K_S. We perform inference in both a regression and a classification setting. The following algorithms and associated results can be found in [25].

Regression. Learning corresponds to minimizing a mean-square error:

$$\begin{equation}\enspace \mathrm{min}\sum _{\mu =1}^{n}{\left[{\hat{Z}}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right)-Z\left({\underline{x}}_{\mu }\right)\right]}^{2}.\end{equation} \tag{ 2 }$$

For algorithms seeking solutions of the form ${\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)={\sum }_{\mu }{a}_{\mu }{K}_{\mathrm{S}}\left({\underline{x}}_{\mu },\underline{x}\right)\equiv \underline{a}\cdot {\underline{k}}_{\mathrm{S}}\left(\underline{x}\right)$ by minimizing the mean-square loss over the vector $\underline{a}$, one obtains:

$$\begin{equation}{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)={\underline{k}}_{\mathrm{S}}\left(\underline{x}\right)\cdot {\mathbb{K}}_{\mathrm{S}}^{-1}\underline{Z},\end{equation} \tag{ 3 }$$

where the vector $\underline{Z}$ contains all the labels in the training set, $\underline{Z}\equiv {\left(Z\left({\underline{x}}_{\mu }\right)\right)}_{\mu =1}^{n}$, and ${\mathbb{K}}_{\mathrm{S},\mu \nu }\equiv {K}_{\mathrm{S}}\left({\underline{x}}_{\mu },{\underline{x}}_{\nu }\right)$ is the Gram matrix. The Gram matrix is always invertible if the kernel K_S is positive definite. The generalization error is then evaluated as the expected mean-square error on unseen data, estimated by averaging over a test set composed of n_test unseen data points:

$$\begin{align}\hfill \mathrm{M}\mathrm{S}\mathrm{E}=\frac{1}{{n}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}\sum _{\mu =1}^{{n}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}{\left[{\hat{Z}}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right)-Z\left({\underline{x}}_{\mu }\right)\right]}^{2}.\end{align} \tag{ 4 }$$

Classification. We perform kernel classification via the algorithm soft-margin SVM. The details can be found in appendix A. After learning from the training data with a student kernel K_S, performance is evaluated via the generalization error. It is estimated as the fraction of correctly predicted labels for data points belonging to a test set with n_test elements.

In figure 1 we present the learning curves for (binary) MNIST and CIFAR10, for regression and classification. Learning is performed both with a Gaussian kernel $K\left(\underline{x}\right)\propto \mathrm{exp}\left(-\vert \vert \underline{x}\vert {\vert }^{2}/\left(2{\sigma }^{2}\right)\right)$ and a Laplace one $K\left(\underline{x}\right)\propto \mathrm{exp}\left(-\vert \vert \underline{x}\vert \vert /\sigma \right)$. Remarkably, the power laws in the two tasks are essentially identical (although the estimated exponent appears to be slightly larger, in absolute value, for classification). Moreover, the two kernels display a very similar behavior, compatible with the same exponent: about −0.4 for MNIST and −0.1 for CIFAR10. The presented data are for σ = 1000; in appendix B we show that the same behavior is observed for different values.

Zoom In Zoom Out Reset image size

We study β in a simplified setting where the data is assumed to follow a Gaussian distribution with known covariance. It falls into the class of teacher–student problems, which are characterized by a machine (the teacher) that generates the data, and another machine (the student) that tries to learn from them. The teacher–student paradigm has been broadly used to study supervised learning [15, 16, 26–33]. He we restrict our attention to kernel methods: we assume that a target function is distributed according to a Gaussian random field $Z\sim \mathcal{N}\left(0,{K}_{\mathrm{T}}\right)$—the teacher—characterized by a translation-invariant isotropic covariance function ${K}_{\mathrm{T}}\left(\underline{x},{\underline{x}}^{\prime }\right)={K}_{\mathrm{T}}\left(\vert \vert \underline{x}-{\underline{x}}^{\prime }\vert \vert \right)$, and that the training dataset consists the finite set of n observations $\underline{Z}={\left(Z\left({\underline{x}}_{\mu }\right)\right)}_{\mu =1}^{n}$. This is equivalent to saying that the vector of training points follows a centered Gaussian distribution with a covariance matrix that depends on K_T and on the location of the points ${\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}$:

$$\begin{equation}\underline{Z}\sim \mathcal{N}\left(\underline{0},{\mathbb{K}}_{\mathrm{T}}\right),\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\quad {\mathbb{K}}_{\mathrm{T}}={\left({K}_{\mathrm{T}}\left({\underline{x}}_{\mu },{\underline{x}}_{\nu }\right)\right)}_{\mu ,\nu =1}^{n}.\end{equation} \tag{ 5 }$$

Once the teacher has generated the dataset, the rest follows as in the kernel regression described in the previous section. We use another translation-invariant isotropic kernel ${K}_{\mathrm{S}}\left(\underline{x},{\underline{x}}^{\prime }\right)$—the student—to infer the value of the field at another point, ${\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)$, with a regression task, i.e. minimizing the mean-square error in equation (2). The solution is therefore given again by equation (3).

Figures 2(a) and (b) shows the mean-square error obtained numerically. In the examples the student is always taken to be a Laplace kernel, and the teacher is either a Laplace kernel or a Gaussian kernel. The points ${\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}$ are taken uniformly at random on the unit d-dimensional hypersphere for several dimensions d and for several dataset sizes n. We take σ_S = σ_T = d as we observed that with this choice smaller datasets were enough to approach a limiting curve—in appendix C we show the plots for the case σ_S = σ_T = 10, which appears to converge to the same limit curve with increasing n, but at a smaller pace. The figure shows that when n is large enough, the mean-square error behaves as a power law (dashed lines) with an exponent that depends on the spatial dimension of the data, as well as on the kernels. The fitted exponents are plotted in figures 2(c) and (d) as a function of the spatial dimension d for different dataset sizes n. In the next section we will discuss the theoretical prediction, which in the figure is plotted as a thick black line. The figure shows that as the dataset gets bigger, the asymptotic exponent tends to our prediction. In appendix D we present the learning curves of Gaussian students with both a Laplace and a Gaussian kernel. When both kernels are Gaussian the test error decays exponentially fast, a result that matches our theoretical prediction. In appendix E we also provide further numerical results for the case where the teacher kernel is a Matérn kernel (as defined therein).

Zoom In Zoom Out Reset image size

In this section we compute analytically the exponent that describes the asymptotic decay of the generalization error when the number n of training data increases. In order to derive the result we assume that both the teacher Gaussian random fields live on a bounded hypercube, $\underline{x}\in \mathcal{V}\equiv {\left[0,L\right]}^{d}$, where L is a constant and d is the spatial dimension. The fields and the kernels can then be thought of as L-periodic along each dimension. Furthermore, to make the problem tractable we assume that the points ${\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}$ live on a regular lattice, covering all the hypercube $\mathcal{V}$. Therefore, the linear spacing between neighboring points is δ = Ln^−1/d . This is a different setting than the one used in the numerical simulations presented in the previous section, for which the data distribution is Gaussian, showing that our results below are robust to such differences.

Generalization error is then evaluated via the typical mean-square error

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}=\mathbb{E}{\left[Z\left(\underline{x}\right)-{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)\right]}^{2},\end{equation} \tag{ 6 }$$

where the expectation is taken over both the teacher process and the point $\underline{x}$ at which we estimate the field, assumed to be uniformly distributed in the hypercube $\mathcal{V}$. In appendix F we prove the following:

Theorem 1. Let ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)={c}_{\mathrm{T}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)={c}_{\mathrm{S}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}\right)$ as $\vert \vert \underline{w}\vert \vert \to \infty$, where ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)$ are the Fourier transforms of the kernels ${K}_{\mathrm{T}}\left(\underline{x}\right)$, ${K}_{\mathrm{S}}\left(\underline{x}\right)$ respectively, assumed to be positive definite. We assume that ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)$ have a finite limit as $\vert \vert \underline{w}\vert \vert \to 0$ and that $K\left(\underline{0}\right){< }\infty$. Then,

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}={n}^{-\beta }+o\left({n}^{-\beta }\right)\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad \beta =\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right).\end{equation} \tag{ 7 }$$

Moreover, in the case of a Gaussian kernel the result holds valid if we take the corresponding exponent to be α = ∞.

Apart from the specific value of the exponent in equation (7), theorem 1 implies that if the student kernel decays fast enough in the frequency domain, then β depends only on the data through the behavior of the teacher kernel at high frequencies. One then recovers β = (α_T − d)/d, also found for the Bayes-optimal setting where the student is identical to the teacher.

Consider the predictions of theorem 1 in the cases presented in figures 2(a) and (b) of Gaussian and Laplace kernels. If both kernels are Laplace kernels then α_T = α_S = d + 1 and $\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}\sim {n}^{-1/d}$, which scales very slowly with the dataset size in large dimensions. If the teacher is a Gaussian kernel (α_T = ∞) and the student is a Laplace kernel then β = 2(1 + 1/d), leading to β → 2 as d → ∞. In figures 2(c) and (d) we compare these predictions with the exponents extracted from figures 2(a) and (b). We plot $\mathrm{log}\enspace \mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}/\mathrm{log}\enspace n\equiv -\beta$, against the dimension d of the data, varying the dataset size n. The exponents extracted numerically tend to our analytical predictions when n is large enough.

Notice that, although the theory and the experiments do not assume the same distribution for the sampling points ${\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}$, this does not seem to yield any difference in the asymptotic behavior of the generalization error, leading to the conjecture that our predictions are exact even when the training set is random, and does not correspond to a lattice. The conjecture can be proven in one dimension following results of the kriging literature [17], but generalization to higher d is a much harder problem. Intuitively, for kernel learning performs an expansion, the quality of which is governed by the target function smoothness and the typical distance δ_min between a point and its nearest neighbors in the training set. Both for random points or on a lattice, one has δ_min ∼ n^−1/d when n is large enough, thus both situations lead to the same β. This is shown in figure 4 (left).

Theorem 1 underlines that kernel methods are subjected to the curse of dimensionality. Indeed, for appropriate students, one obtains β = (α_T − d)/d. Let us define the smoothness index $s\equiv \left[\left({\alpha }_{\mathrm{T}}-d\right)/2\right]=\beta d/2$, which must be $\mathcal{O}\left(d\right)$ to avoid β → 0 for large d. The two lemmas below, derived in the appendix, indicate that the target function is s time differentiable (in a mean-square sense). Thus learning with kernels in very large dimension can only occur if the target function is $\mathcal{O}\left(d\right)$ times differentiable, a condition that appears very restrictive in large d.

Lemma 1. Let $K\left(\underline{x},{\underline{x}}^{\prime }\right)$ be a translation-invariant isotropic kernel such that $\tilde {K}\left(\underline{w}\right)=c\vert \vert \underline{w}\vert {\vert }^{-\alpha }+o\left(\vert \vert \underline{w}\vert {\vert }^{-\alpha }\right)$ as ||w|| → ∞ and $\vert \vert \underline{w}\vert {\vert }^{d}\tilde {K}\left(\underline{w}\right)\to 0$ as ||w|| → 0. If α > d + n for some $n\in {\mathbb{Z}}^{+}$, then $K\left(\underline{x}\right)\in {C}^{n}$, that is, it is at least n-times differentiable. (Proof in appendix G).

Lemma 2. Let $Z\sim \mathcal{N}\left(0,K\right)$ be a d-dimensional Gaussian random field, with K ∈ C²ⁿ being a 2n-times differentiable kernel. Then Z is n-times mean-square differentiable in the sense that

• Derivatives of $Z\left(\underline{x}\right)$ are Gaussian random fields;
• $\mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)=0$;
• $\mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{1}}^{{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}^{\prime }}Z\left({\underline{x}}^{\prime }\right)={\partial }_{{x}_{1}}^{{n}_{1}+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }\right){< }\infty$ if the derivatives of K exist.

In particular, $\mathbb{E}{\partial }_{{x}_{i}}^{m}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{i}}^{m}Z\left({\underline{x}}^{\prime }\right)={\partial }_{{x}_{i}}^{2m}K\left(\underline{x}-{\underline{x}}^{\prime }\right){< }\infty \quad \forall \enspace m{\leqslant}n$. (Proof in appendix G).

Interpretation of theorem 1: when the student does not limit performance, i.e. when 2α_S > α_T − d and $\beta =\frac{{\alpha }_{\mathrm{T}}-d}{d}$, we can interpret the result as follows. An isotropic student kernel corresponds to a Gaussian prior on the Fourier coefficients of the target function being learned. The student puts large (low) power at low (high) frequencies, and it can then reconstruct a number of the order of n largest Fourier coefficients, which corresponds to frequencies $\underline{w}$ of norm $\vert \vert \underline{w}\vert \vert {\leqslant}1/\delta \sim {n}^{1/d}$. Fourier coefficients $\tilde {Z}\left(\underline{w}\right)$ at higher frequencies cannot be learned, and the mean square error is then simply of order of the sum of the squares of these coefficients:

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}\sim \sum _{\vert \vert \underline{w}\vert \vert {\geqslant}{n}^{1/d}}\vert \tilde {Z}\left(\underline{w}\right){\vert }^{2}\sim \sum _{\vert \vert \underline{w}\vert \vert {\geqslant}{n}^{1/d}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\sim {n}^{\frac{d-{\alpha }_{\mathrm{T}}}{d}}\sim {n}^{-\beta }.\end{equation} \tag{ 8 }$$

Equation (8) is not readily applicable to real data which are neither Gaussian nor uniformly distributed. However, it supports the following broader result: kernel methods can predict well of order n first coefficients of the true function in the eigenbasis of the kernel, but not the following ones. For any student kernel K_S, let λ₁ ⩾ λ_ρ ⩾ ... be its eigenvalues (positive and real, because of symmetry and positive definiteness) and ${\phi }_{\rho }\left(\underline{x}\right)$ the associated eigenfunctions:

$$\begin{equation}\int {\mathrm{d}}^{d}\underline{y}\enspace p\left(\underline{y}\right){K}_{\mathrm{S}}\left(\underline{x}-\underline{y}\right){\phi }_{\rho }\left(\underline{y}\right)={\lambda }_{\rho }{\phi }_{\rho }\left(\underline{x}\right),\end{equation} \tag{ 9 }$$

where $p\left(\underline{x}\right)$ is the density of the data points. Then the kernel can be decomposed in its eigenmodes:

$$\begin{equation}{K}_{\mathrm{S}}\left(\underline{x}-\underline{y}\right)=\sum _{\rho {\geqslant}1}{\lambda }_{\rho }{\phi }_{\rho }\left(\underline{x}\right){\phi }_{\rho }\left(\underline{y}\right).\end{equation} \tag{ 10 }$$

The eigenfunctions of K make a complete basis, and we can write any function $Z\left(\underline{x}\right)$ as

$$\begin{equation}Z\left(\underline{x}\right)=\sum _{\rho }{q}_{\rho }{\phi }_{\rho }\left(\underline{x}\right).\end{equation} \tag{ 11 }$$

The generalization of our result then simply reads:

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}\sim \sum _{\rho {\geqslant}n}{q}_{\rho }^{2}.\end{equation} \tag{ 12 }$$

Assuming a power-law behavior ${q}_{\rho }^{2}\sim {\rho }^{-a}$ then leads to β = a − 1.

To extract the exponent a and test this prediction for real data, we first approximate the eigenvalue equation equation (9) for the student kernel with the diagonalization of its finite-dimensional Gram matrix ${\mathbb{K}}_{\mathrm{S}}$ computed on a large dataset of size $\tilde {n}$:

$$\begin{equation}{\mathbb{K}}_{\mathrm{S}}{\underline{\phi }}_{\rho }\sim {\lambda }_{\rho }{\underline{\phi }}_{\rho },\end{equation} \tag{ 13 }$$

where now we have $\tilde {n}$ eigenvalues ${\lambda }_{1}{\geqslant}\cdots {\geqslant}{\lambda }_{\tilde {n}}$ and the eigenvectors ${\underline{\phi }}_{\rho }$ are $\tilde {n}$-dimensional eigenvectors. Computing this diagonalization for a given training set is referred to (uncentered) kernel PCA. This procedure is a discretized version of equation (10) and yields only an approximation to the largest $\tilde {n}$ eigenvalues of the kernel, that are exactly recovered as $\tilde {n}\to \infty$ (in figure J1 in appendix J we show that the eigenvalues of the Gram matrix converge when $\tilde {n}$ increases, and that their density displays the power-law behavior that one can extract from the kernel operator with a uniform distribution $p\left(\underline{x}\right)$). The coefficient q_ρ is then estimated by the scalar products $\left(\underline{Z}\cdot {\underline{\phi }}_{\rho }\right)$, where $\underline{Z}=\left(Z\left({\underline{x}}_{1}\right),\dots ,Z\left({\underline{x}}_{\tilde {n}}\right)\right)$ is the vector of the target function's values on the train set.

Finally, we approximate equation (12) as:

$$\begin{equation}\sum _{\rho =n}^{\tilde {n}}{\left(\underline{Z}\cdot {\underline{\phi }}_{\rho }\right)}^{2}.\end{equation} \tag{ 14 }$$

This quantity is plotted in figure 3, where we show that it correlates remarkably well with the true learning curve. Fitting these cumulative curves whose exponent is 1 − a for asymptotically large $\tilde {n}$ (here we plotted several curves for growing $\tilde {n}$) we extract an exponent a_MNIST = 1.36 leading to ${\hat{\beta }}_{\mathrm{M}\mathrm{N}\mathrm{I}\mathrm{S}\mathrm{T}}\approx 0.36$ and a_CIFAR10 = 1.07 leading to ${\hat{\beta }}_{\mathrm{C}\mathrm{I}\mathrm{F}\mathrm{A}\mathrm{R}\mathrm{10}}\approx 0.07$ that are very close to the exponents that we measured in section 3.

Zoom In Zoom Out Reset image size

Support for the genericity of equation (12) can be obtained from the recent paper [11], where the authors derived a formula for the generalization error based on the decomposition of the target function on the eigenbasis of the kernel. The formula is derived with uncontrolled approximations, but applies to a generic target function (or ensembles thereof) and a generic data point distribution $p\left(\underline{x}\right)$, and matches well their numerical experiments. In appendix I we show that the asymptotic limit (large n) of their formula yields equation (12). Furthermore, equation (7) is recovered if a power-law decay of the coefficient of the true function in the eigenbasis of the kernel is assumed, thus generalizing our result to non-Gaussian data.

Both our predictions and empirical observations support rather large values of β. From a Gaussian random process point of view, it is surprising: the exponent β avoids the curse of dimensionality only if the smoothness of the teacher is of the order of the data dimension. If these observations were to hold true also for real data, it would seem to imply that MNIST and CIFAR10 must be hundreds or thousands of times differentiable. However, there is a simple catch: real data actually live on a manifold of much lower dimensionality. Above we have argued that the quantity that governs the asymptotic learning curve is the typical distance δ_min between neighboring points in the training set. A simple way to measure the effective dimension d_eff of real data therefore consists of plotting the (asymptotic) dependence of δ_min on the number of points n in a random subset of the dataset, and fitting

$$\begin{align}\hfill {\delta }_{\mathrm{m}\mathrm{i}\mathrm{n}}\sim {n}^{-1/{d}_{\mathrm{e}\mathrm{f}\mathrm{f}}}.\end{align} \tag{ 15 }$$

In figure 4 (right) we show that for MNIST and CIFAR10 there is indeed a power-law relation linking δ_min to n, and that the effective dimension extracted this way is much smaller than the embedding dimension of the datasets:

$$\begin{align}\hfill {d}_{\mathrm{e}\mathrm{ff}}^{\mathrm{M}\mathrm{N}\mathrm{I}\mathrm{S}\mathrm{T}}\approx 15& \ll 784={d}^{\mathrm{M}\mathrm{N}\mathrm{I}\mathrm{S}\mathrm{T}},\hfill \end{align} \tag{ 16 }$$

$$\begin{align}\hfill {d}_{\mathrm{e}\mathrm{ff}}^{\mathrm{C}\mathrm{I}\mathrm{F}\mathrm{A}\mathrm{R}\mathrm{10}}\approx 35& \ll 3072={d}^{\mathrm{C}\mathrm{I}\mathrm{F}\mathrm{A}\mathrm{R}\mathrm{10}}.\hfill \end{align} \tag{ 17 }$$

Zoom In Zoom Out Reset image size

This measure is consistent with previous extrapolations of the intrinsic dimension of MNIST [19, 20, 22, 23].

In this work we have shown for CIFAR10 and MNIST respectively that kernel regression and classification display a power-law decay in the learning curves, quite remarkably with essentially the same exponent β regardless of task and kernel—a fact yet to be explained. These exponents are much larger than β = 1/d expected for Lipschitz target functions and smaller than β = 1/2 expected for RKHS target functions.

This observation led us to study a teacher–student framework for regression in which data are modeled as Gaussian random fields of varying smoothness, in which intermediary values of β are obtained. We find two regimes depending on the respective smoothness of the teacher and student kernels. If the student is smooth enough—i.e. it puts a sufficiently low prior on high frequency components—then β is entirely controlled by the teacher. We obtain that the smoothness index must scale with the dimension for β to be finite as d → ∞, recovering the curse of dimensionality.

In our calculations, the dimension enters as the parameter relating the number of points to the nearest-neighbor distance δ ∼ n^−1/d . Thus, in practice the parameter d considered should be the effective dimension d_eff of the data, which is much smaller than the number of pixels for MNIST and CIFAR data. It explains why β is not very small in these cases.

Finally, for Gaussian fields our result is equivalent to the statement that β is governed by the power of the true function past the first ∼n eigenvectors of the kernel. We test this more general idea both for CIFAR and MNIST and find that it correctly predicts β. Understanding what controls this power in a general setting (which include the effective dimension of the data and presumably a generalized quantity characterizing smoothness) thus appears necessary to understand how many data are required to learn a task.

We acknowledge G Biroli, C Hongler and F Gabriel for the discussions that stimulated this work, and we thank S d'Ascoli, A Jacot, C Pehlevan, L Sagun and M L Stein for discussions. This work was partially supported by the grant from the Simons Foundation (#454953 Matthieu Wyart). M W thanks the Swiss National Science Foundation for support under Grant No. 200021-165509.

The kernel classification task is performed via the algorithm known as soft-margin support vector machine.

We want to find a function ${\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)$ such that its sign correctly predicts the label of the data. In this context we model such a function as a linear prediction after projecting the data on a feature space via $\underline{x}\to \phi \left(\underline{x}\right)$:

$$\begin{equation}{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)=\underline{w}\cdot {\underline{\phi }}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right)+b,\end{equation} \tag{ 18 }$$

where $\underline{w},b$ are parameters to be learned from the training data. The kernel is related to the feature space via ${K}_{\mathrm{S}}\left(\underline{x},{\underline{x}}^{\prime }\right)={\underline{\phi }}_{\mathrm{S}}\left(\underline{x}\right)\cdot {\underline{\phi }}_{\mathrm{S}}\left({\underline{x}}^{\prime }\right)$. We require that $Z\left({\underline{x}}_{\mu }\right){\hat{Z}}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right){ >}1-{\xi }_{\mu }$ for all training points. Ideally, we want to have some large margins 1 − ξ_μ = 1, but we allow some of them to be smaller by introducing the slack variables ξ_μ and penalizing large values. To achieve this the following constrained minimization is performed:

$$\begin{equation}{\mathrm{min}}_{\underline{w},b,\underline{\xi }}\enspace \frac{1}{2}\vert \vert \underline{w}\vert {\vert }^{2}+C\sum _{\mu }{\xi }_{\mu }\enspace \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\enspace \mathrm{t}\mathrm{o}\enspace \forall \enspace \mu \enspace Z\left({\underline{x}}_{\mu }\right)\left[\underline{w}\cdot {\underline{\phi }}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right)+b\right]{\geqslant}1-{\xi }_{\mu },\quad {\xi }_{\mu }{\geqslant}0.\end{equation} \tag{ 19 }$$

This problem can be expressed in a dual formulation as

$$\begin{equation}{\mathrm{min}}_{\underline{a}}\enspace \frac{1}{2}\underline{a}\cdot {\mathbb{Q}}_{\mathrm{S}}\underline{a}-\sum _{\mu =1}^{n}{a}_{\mu }\enspace \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\enspace \mathrm{t}\mathrm{o}\enspace \underline{Z}\cdot \underline{a}=0,\quad 0{\leqslant}{a}_{\mu }{\leqslant}C,\end{equation} \tag{ 20 }$$

where ${\mathbb{Q}}_{\mathrm{S},\mu ,\nu }=Z\left({\underline{x}}_{\mu }\right)Z\left({\underline{x}}_{\nu }\right){K}_{\mathrm{S}}\left({\underline{x}}_{\mu },{\underline{x}}_{\nu }\right)$ and $\underline{Z}$ is the vector of the labels of the training points. Here C (=10⁴ in our simulations) controls the trade-off between minimizing the training error and maximizing the margins 1 − ξ_μ . For the details we refer to [25]. If ${\underline{a}}^{\star }$ is the solution to the minimization problem, then

$$\begin{align}\hfill & {\underline{w}}^{\star }=\sum _{\mu }{a}_{\mu }^{\star }{\underline{\phi }}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right),\hfill \end{align} \tag{ 21 }$$

$$\begin{align}\hfill & {b}^{\star }=Z\left({\underline{x}}_{\mu }\right)-\sum _{\nu }{a}_{\nu }{y}_{\nu }{K}_{\mathrm{S}}\left({\underline{x}}_{\mu },{\underline{x}}_{\nu }\right)\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{n}\mathrm{y}\enspace \mu \enspace \mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\enspace \mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}\enspace {a}_{\mu }{< }C.\hfill \end{align} \tag{ 22 }$$

The predicted label for unseen data points is then

$$\begin{equation}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)\right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\sum _{\mu }Z\left({\underline{x}}_{\mu }\right){a}_{\mu }{K}_{\mathrm{S}}\left({\underline{x}}_{\mu },\underline{x}\right)+{b}^{\star }\right)\end{equation} \tag{ 23 }$$

The generalization error is now defined as the probability that an unseen image has a predicted label different from the true one, and such a probability is again estimated as an average over a test set with n_test elements:

$$\begin{equation}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\frac{1}{{n}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}\sum _{\mu =1}^{{n}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}\theta \left[-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\hat{Z}}_{\mathrm{S}}\left({\underline{x}}_{\mu }\right)\right)Z\left({\underline{x}}_{\mu }\right)\right].\end{equation} \tag{ 24 }$$

In figure B1 we show the learning curves for kernel regression on the MNIST (parity) dataset—the same setting as in figure 1(a). Several Laplace kernels of varying variance σ are used. The variance ranges several orders of magnitude and the learning curves all decay with the same exponent, although for σ = 10 the algorithm achieves suboptimal performance and the test errors are increased by some factor.

Zoom In Zoom Out Reset image size

In figure C1 we show the learning curves for the teacher–student kernel regression problem, with a student kernel that is always Laplace and a teacher that can be either Gaussian or Laplace. We show how the test error decays with the size of the training dataset and how the asymptotic exponent depends on the spatial dimension. Every experiment is run with two different choices of the kernel variances: in one case σ_T = σ_S = d and in the other σ_T = σ_S = 10. We observed that scaling the variances with the spatial dimension leads more quickly to the results that we predicted in this paper, but overall the choice has little effect on the exponents (both tend towards the prediction as the dataset size is increased).

Zoom In Zoom Out Reset image size

In this appendix we present the learning curves of Gaussian students: the Fourier transform of these kernels decays faster than any power law and one can effectively consider α_S = ∞. If the teacher is Laplace (α_T = d + 1) then the predicted exponent is finite and takes the values $\beta =\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right)=\frac{1}{d}\enspace \mathrm{min}\left(1,\infty \right)=\frac{1}{d}$. Such a case is displayed in figure D1 (left) in dimension d = 6. However, if we consider the teacher to be Gaussian as well, then the predicted exponent would be $\beta =\frac{1}{d}\enspace \mathrm{min}\left(\infty ,\infty \right)=\infty$. This case corresponds to figure D1 (center): the test errors decay faster than a power law. In figure 7 (right) we compare the case where both kernels are Gaussian to the case where both kernels are Laplace: while the latter decays as a power law, the former decays much faster.

Zoom In Zoom Out Reset image size

To further test the applicability of our theory, we show here some numerical simulations for a teacher kernel that is a Matérn covariance function and a Laplace kernel as student. We ran the simulations in 1d: the data points are sampled uniformly on a one-dimensional circle embedded in ${\mathbb{R}}^{2}$. Matérn kernels are parametrized by a parameter ν > 0:

$$\begin{equation}{K}_{\mathrm{T}}\left(\underline{x}\right)=\frac{{2}^{1-\nu }}{{\Gamma}\left(\nu \right)}{z}^{\nu }{\mathcal{K}}_{\nu }\left(z\right),\end{equation} \tag{ 25 }$$

where $z=\sqrt{2\nu }\frac{\vert \vert \underline{x}\vert \vert }{\sigma }$ (σ being the kernel variance), Γ is the gamma function and ${\mathcal{K}}_{\nu }$ is the Bessel function of the second kind with parameter ν. Interestingly we recover the Laplace kernel for ν = 1/2 and the Gaussian kernel for ν = ∞. As one can find in e.g. [34], the exponent α_T that governs the decay at high frequency of this kernels is α_T = d + 2ν. Varying ν we can change the smoothness of the target function.

For d = 1 our prediction for the learning curve exponent β is

$$\begin{equation}\beta =\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right)=\mathrm{min}\left(2\nu ,4\right).\end{equation} \tag{ 26 }$$

In figure E1 we verify that our prediction matches the numerical results.

Zoom In Zoom Out Reset image size

We prove here theorem 1:

Theorem 1. Let ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)={c}_{\mathrm{T}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)={c}_{\mathrm{S}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}\right)$ as $\vert \vert \underline{w}\vert \vert \to \infty$, where ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)$ are the Fourier transforms of the kernels ${K}_{\mathrm{T}}\left(\underline{x}\right)$, ${K}_{\mathrm{S}}\left(\underline{x}\right)$ respectively, assumed to be positive definite. We assume that ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)$ have a finite limit as $\vert \vert \underline{w}\vert \vert \to 0$ and that $K\left(\underline{0}\right){< }\infty$. Then,

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}={n}^{-\beta }+o\left({n}^{-\beta }\right)\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad \beta =\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right).\end{equation} \tag{ 27 }$$

Moreover, in the case of a Gaussian kernel the result holds valid if we take the corresponding exponent to be α = ∞.

Proof. Our strategy is to compute how the mean-square test error scales with distance δ between two nearest neighbors on the d-dimensional regular lattice. At the end, we will use the fact that δ ∝ n^−1/d , where n is the number of sampled points on the lattice.

We denote by $\tilde {F}\left(\underline{w}\right)$ the Fourier transform of a function $F:\mathcal{V}\to \mathbb{R}$:

$$\begin{align}\hfill & \tilde {F}\left(\underline{w}\right)={L}^{-d/2}{\int }_{\mathcal{V}}\mathrm{d}\underline{x}\enspace {\text{e}}^{-\text{i}\underline{w}\cdot \underline{x}}\enspace F\left(\underline{x}\right),\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\enspace \underline{w}\in \mathbb{L}\equiv \frac{2\pi }{L}{\mathbb{Z}}^{d},\hfill \end{align} \tag{ 28 }$$

$$\begin{align}\hfill & F\left(\underline{x}\right)={L}^{-d/2}\sum _{\underline{w}\in \mathbb{L}}{\text{e}}^{\text{i}\underline{w}\cdot \underline{x}}\enspace \tilde {F}\left(\underline{w}\right).\hfill \end{align} \tag{ 29 }$$

If $Z\sim \mathcal{N}\left(0,K\right)$ is a Gaussian field with translation-invariant covariance K then by definition

$$\begin{equation}\mathbb{E}Z\left(\underline{x}\right)Z\left({\underline{x}}^{\prime }\right)=K\left(\underline{x}-{\underline{x}}^{\prime }\right).\end{equation} \tag{ 30 }$$

Properties of the Fourier transform of a Gaussian field:

$$\begin{align}\hfill & \tilde {K}\left(\underline{w}\right)=\tilde {K}\left(-\underline{w}\right)\in \mathbb{R},\hfill \end{align} \tag{ 31 }$$

$$\begin{align}\hfill & \mathbb{E}\tilde {Z}\left(\underline{w}\right)=0,\hfill \end{align} \tag{ 32 }$$

$$\begin{align}\hfill & \mathbb{E}\tilde {Z}\left(\underline{w}\right)\bar{\tilde {Z}\left({\underline{w}}^{\prime }\right)}={L}^{d/2}{\delta }_{\underline{w}{\underline{w}}^{\prime }}\tilde {K}\left(\underline{w}\right).\hfill \end{align} \tag{ 33 }$$

Equation (31) comes from the fact that $K\left(\underline{x}\right)$ is an even, real-valued function. The real and imaginary parts of $\tilde {Z}\left(\underline{w}\right)$ are Gaussian random variables. They are all independent except that $\tilde {Z}\left(-\underline{w}\right)=\bar{\tilde {Z}\left(\underline{w}\right)}$. Equation (33) follows from the fact that $Z\left(\underline{x}\right)$ and $K\left(\underline{x}\right)$ are L-periodic functions, and therefore ${\text{e}}^{\text{i}\underline{w}\cdot \underline{x}}\enspace \tilde {K}\left(\underline{w}\right)$ is the Fourier transform of $K\left(\cdot +\underline{x}\right)$ if $\underline{w}\in \frac{2\pi }{L}{\mathbb{Z}}^{d}$. □

The solution equation (3) for kernel regression has two interpretations. In section 4 we introduced it as the quantity that minimizes a quadratic error, but it can also be seen as the maximum-a-posteriori estimation of another formulation of the problem [34]. The field $Z\left(\underline{x}\right)$ is assumed to be drawn from a Gaussian distribution with covariance function ${K}_{\mathrm{S}}\left(\underline{x}\right)$: K_S therefore plays a role in the prior distribution of the data $\underline{Z}=\left(Z{\left({\underline{x}}_{\mu }\right)}_{\mu =1}^{n}\right)$. Inference about the value of the field ${\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)$ at another location is then performed by maximizing its posterior distribution,

$$\begin{equation}{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)\equiv \mathrm{a}\mathrm{r}\mathrm{g}\;\;\mathrm{m}\mathrm{a}\mathrm{x}\enspace \mathcal{P}\left(Z\left(\underline{x}\right)\vert \underline{Z}\right).\end{equation} \tag{ 34 }$$

Such a posterior distribution is Gaussian, and its mean—and therefore also the value that maximizes the probability—is exactly equation (3):

$$\begin{equation}{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)={\underline{k}}_{\mathrm{S}}\left(\underline{x}\right)\cdot {\mathbb{K}}_{\mathrm{S}}^{-1}\underline{Z},\end{equation} \tag{ 35 }$$

where $\underline{Z}={\left(Z\left({\underline{x}}_{\mu }\right)\right)}_{\mu =1}^{n}$ are the training data, ${\underline{k}}_{\mathrm{S}}\left(\underline{x}\right)={\left({K}_{\mathrm{S}}\left({\underline{x}}_{\mu },\underline{x}\right)\right)}_{\mu =1}^{n}$ and ${\mathbb{K}}_{\mathrm{S}}={\left({K}_{\mathrm{S}}\left({\underline{x}}_{\mu },{\underline{x}}_{\nu }\right)\right)}_{\mu ,\nu =1}^{n}$ is the Gram matrix, that is invertible since the kernel K_S is assumed to be positive definite. By Fourier transforming this relation we find

$$\begin{equation}{\tilde {Z}}_{\mathrm{S}}\left(\underline{w}\right)={\tilde {Z}}^{\star }\left(\underline{w}\right)\frac{{\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)},\end{equation} \tag{ 36 }$$

where we have defined ${F}^{\star }\left(\underline{w}\right)\equiv {\sum }_{\underline{n}\in {\mathbb{Z}}^{d}}F\left(\underline{w}+\frac{2\pi \underline{n}}{\delta }\right)$ for a generic function F.

Another way to reach equation (36) is to consider that we are observing the quantities

$$\begin{equation}{\tilde {Z}}^{\star }\left(\underline{w}\right)\equiv {\delta }^{d}{L}^{-d/2}\sum _{\underline{x}\in \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}}{\text{e}}^{-\text{i}\underline{w}\cdot \underline{x}}\enspace Z\left(\underline{x}\right)\equiv \sum _{\underline{n}\in {\mathbb{Z}}^{d}}\tilde {Z}\left(\underline{w}+\frac{2\pi \underline{n}}{\delta }\right).\end{equation} \tag{ 37 }$$

Given that we know the prior distribution of the Fourier components on the right-hand side in equation (37), we can infer their posterior distribution once their sums are constrained by the value of ${\tilde {Z}}^{\star }\left(\underline{w}\right)$, and it is straightforward to see that we recover equation (36).

The mean-square error can then be written using the Parseval–Plancherel identity,

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}={L}^{-d}\mathbb{E}{\int }_{\mathcal{V}}\mathrm{d}\underline{x}\enspace {\left[Z\left(\underline{x}\right)-{\hat{Z}}_{\mathrm{S}}\left(\underline{x}\right)\right]}^{2}={L}^{-d}\mathbb{E}\sum _{\underline{w}\in \mathbb{L}}{\left\vert \tilde {Z}\left(\underline{w}\right)-{\tilde {Z}}^{\star }\left(\underline{w}\right)\frac{{\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)}\right\vert }^{2}.\end{equation} \tag{ 38 }$$

By taking the expectation value with respect to the teacher and using equations (31)–(33) we can write the mean-square error as

$$\begin{align}\hfill \mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}& ={L}^{-d}\mathbb{E}\sum _{\underline{w}\in \mathbb{L}}\left[\tilde {Z}\left(\underline{w}\right)\bar{\tilde {Z}\left(\underline{w}\right)}-2\tilde {Z}\left(\underline{w}\right)\bar{{\tilde {Z}}^{\star }\left(\underline{w}\right)}\frac{{\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)}+{\tilde {Z}}^{\star }\left(\underline{w}\right)\bar{{\tilde {Z}}^{\star }\left(\underline{w}\right)}\frac{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}\right]\hfill \\ \hfill & ={L}^{-d}\mathbb{E}\sum _{\underline{w}\in \mathbb{L}}\tilde {Z}\left(\underline{w}\right)\bar{\tilde {Z}\left(\underline{w}\right)}-2\frac{{\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)}\sum _{\underline{n}\in {\mathbb{Z}}^{d}}\tilde {Z}\left(\underline{w}\right)\bar{\tilde {Z}\left(\underline{w}+\frac{2\pi \underline{n}}{\delta }\right)}\hfill \\ \hfill & \quad +\enspace \frac{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}\sum _{\underline{n},{\underline{n}}^{\prime }\in {\mathbb{Z}}^{d}}\tilde {Z}\left(\underline{w}+\frac{2\pi \underline{n}}{\delta }\right)\bar{\tilde {Z}\left(\underline{w}+\frac{2\pi {\underline{n}}^{\prime }}{\delta }\right)}\hfill \\ \hfill & ={L}^{-d/2}\sum _{\underline{w}\in \mathbb{L}}{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)-2\frac{{\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)}{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)+\frac{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}{{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\left(\underline{w}\right)}\sum _{\underline{n}\in {\mathbb{Z}}^{d}}{\tilde {K}}_{\mathrm{T}}\left(\underline{w}+\frac{2\pi \underline{n}}{\delta }\right)\hfill \\ \hfill & ={L}^{-d/2}\sum _{\underline{w}\in \mathbb{L}\cap \mathcal{B}}{\tilde {K}}_{\mathrm{T}}^{\star }\left(\underline{w}\right)-2\frac{{{\left[{{\tilde {K}}_{\mathrm{T}}}{\tilde {K}}_{\mathrm{S}}\right]}^{\star }}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }\left(\underline{w}\right)}+\frac{{\tilde {K}}_{\mathrm{T}}^{\star }\left(\underline{w}\right){\left[{\tilde{K}}_{{{\mathrm{S}}^{\star}}^{2}}\right]}^{\star }\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{\star }{{\left(\underline{w}\right)}^{2}}},\hfill \end{align} \tag{ 39 }$$

where $\mathcal{B}={\left[-\frac{\pi }{\delta },\frac{\pi }{\delta }\right]}^{d}$ is the Brillouin zone.

At high frequencies, ${\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)={c}_{\mathrm{T}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)$ and ${\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right)={c}_{\mathrm{S}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{S}}}\right)$. Therefore:

$$\begin{align}\hfill {\tilde {K}}_{\mathrm{T}}^{\star }\left(\underline{w}\right)& ={\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)+{\delta }^{{\alpha }_{\mathrm{T}}}{c}_{\mathrm{T}} \sum _{\underline{n}\in {\mathbb{Z}}^{d}{\backslash}\left\{\underline{0}\right\}} \vert \vert \underline{w}\delta +2\pi \underline{n}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)\hfill \\ \hfill & \equiv {\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)+{\delta }^{{\alpha }_{\mathrm{T}}}{c}_{\mathrm{T}}\enspace {\psi }_{\mathrm{T}}\left(\underline{w}\delta \right)+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right).\hfill \end{align} \tag{ 40 }$$

This equation defines the function ψ_T, and a similar equation holds for the student as well. The hypothesis ${K}_{\mathrm{T}}\left(\underline{0}\right)\propto \int \mathrm{d}\underline{w}\enspace {\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right){< }\infty$ implies α_T > d and therefore ${\sum }_{\underline{n}\in {\mathbb{Z}}^{d}}\vert \vert \underline{n}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}{< }\infty$ (and likewise for the student). Then, ${\psi }_{{\alpha }_{\mathrm{T}}}\left(\underline{0}\right),{\psi }_{{\alpha }_{\mathrm{S}}}\left(\underline{0}\right)$ are finite; furthermore, the $\underline{w}$'s in the sum equation (39) are at most of order $\mathcal{O}\left({\delta }^{-1}\right)$, therefore the terms ${\psi }_{\alpha }\left(\underline{w}\delta \right)$ are $\mathcal{O}\left({\delta }^{0}\right)$ and do not influence how equation (39) scales with δ. Applying equation (40), expanding for δ ≪ 1 and keeping only the leading orders, we find

$$\begin{align}\hfill \mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}=& ={L}^{-d/2}\left[\sum _{\underline{w}\in \mathbb{L}\cap \mathcal{B}}2{c}_{\mathrm{T}}{\psi }_{{\alpha }_{\mathrm{T}}}\left(\underline{w}\delta \right){\delta }^{{\alpha }_{\mathrm{T}}}+{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}{\delta }^{2{\alpha }_{\mathrm{S}}}\right.\hfill \\ \hfill & \quad \left.+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}}\right)+o\left(\vert \vert \underline{w}\vert {\vert }^{-2{\alpha }_{\mathrm{S}}}\right)\right]\hfill \\ \hfill & ={L}^{-d/2}\left[\sum _{\underline{w}\in \mathbb{L}\cap \mathcal{B}}2{c}_{\mathrm{T}}{\psi }_{{\alpha }_{\mathrm{T}}}\left(\underline{w}\delta \right){\delta }^{{\alpha }_{\mathrm{T}}}+{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}{\delta }^{2{\alpha }_{\mathrm{S}}}\right]\hfill \\ \hfill & \quad +o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}-d}\right)+o\left(\vert \vert \underline{w}\vert {\vert }^{-2{\alpha }_{\mathrm{S}}-d}\right).\hfill \end{align} \tag{ 41 }$$

We have neglected terms proportional to, for instance, ${\delta }^{{\alpha }_{\mathrm{T}}+{\alpha }_{\mathrm{S}}}$, since they are subleading with respect to ${\delta }^{{\alpha }_{\mathrm{T}}}$, but we must keep both ${\delta }^{{\alpha }_{\mathrm{T}}}$ and ${\delta }^{{\alpha }_{\mathrm{S}}}$ since we do not know a priori which one is dominant. The additional term δ^−d in the subleading terms comes from the fact that $\vert \mathbb{L}\cap \mathcal{B}\vert =\mathcal{O}\left({\delta }^{-d}\right)$.

The first term in equation (41) is the simplest to deal with: since $\vert \vert \underline{w}\delta \vert \vert$ is smaller than some constant for all $\underline{w}\in \mathbb{L}\cap \mathcal{B}$ and the function ${\psi }_{{\alpha }_{\mathrm{T}}}\left(\underline{w}\delta \right)$ has a finite limit, we have

$$\begin{equation}{\delta }^{{\alpha }_{\mathrm{T}}}\sum _{\underline{w}\in \mathbb{L}\cap \mathcal{B}}2{c}_{\mathrm{T}}{\psi }_{{\alpha }_{\mathrm{T}}}\left(\underline{w}\delta \right)=\mathcal{O}\left({\delta }^{{\alpha }_{\mathrm{T}}}\vert \mathbb{L}\cap \mathcal{B}\vert \right)=\mathcal{O}\left({\delta }^{{\alpha }_{\mathrm{T}}-d}\right).\end{equation} \tag{ 42 }$$

We then split the second term in equation (41) in two contributions:

Small $\vert \vert \underline{w}\vert \vert$. We consider 'small' all the terms $\underline{w}\in \mathbb{L}\cap \mathcal{B}$ such that $\vert \vert \underline{w}\vert \vert {< }{\Gamma}$, where Γ ≫ 1 is $\mathcal{O}\left({\delta }^{0}\right)$ but large. As δ → 0, ${\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\to {\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{0}\right)$ which is finite because $K\left(\underline{0}\right){< }\infty$. Therefore

$$\begin{equation}{\delta }^{2{\alpha }_{\mathrm{S}}}\sum _{\begin{subarray}{c}\underline{w}\in \mathbb{L}\cap \mathcal{B}\\ \vert \vert \underline{w}\vert \vert {< }{\Gamma}\end{subarray}}{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}\to {\delta }^{2{\alpha }_{\mathrm{S}}}{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{0}\right)\sum _{\begin{subarray}{c}\underline{w}\in \mathbb{L}\cap \mathcal{B}\\ \vert \vert \underline{w}\vert \vert {< }{\Gamma}\end{subarray}}\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}.\end{equation} \tag{ 43 }$$

The summand is real and strictly positive because the positive definiteness of the kernels implies that their Fourier transforms are strictly positive. Moreover, as δ → 0, $\mathbb{L}\cap \mathcal{B}\cap \left\{\vert \vert \underline{w}\vert \vert {< }{\Gamma}\right\}\to \mathbb{L}\cap \left\{\vert \vert \underline{w}\vert \vert {< }{\Gamma}\right\}$, which contains a finite number of elements, independent of δ. Therefore

$$\begin{equation}{\delta }^{2{\alpha }_{\mathrm{S}}}\sum _{\begin{subarray}{c}\underline{w}\in \mathbb{L}\cap \mathcal{B}\\ \vert \vert \underline{w}\vert \vert {< }{\Gamma}\end{subarray}}{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}=\mathcal{O}\left({\delta }^{2{\alpha }_{\mathrm{S}}}\right).\end{equation} \tag{ 44 }$$

Large $\vert \vert \underline{w}\vert \vert$. 'Large' $\underline{w}$ are those with $\vert \vert \underline{w}\vert \vert { >}{\Gamma}$: we recall that Γ ≫ 1 is $\mathcal{O}\left({\delta }^{0}\right)$ but large. This allows us to approximate ${\tilde {K}}_{\mathrm{T}}$, ${\tilde {K}}_{\mathrm{S}}$ in the sum with their asymptotic behavior:

$$\begin{align}\hfill & {\delta }^{2{\alpha }_{\mathrm{S}}}\sum _{\begin{subarray}{c}\underline{w}\in \mathbb{L}\cap \mathcal{B}\\ \vert \vert \underline{w}\vert \vert { >}{\Gamma}\end{subarray}}{c}_{\mathrm{S}}^{2}{\psi }_{2{\alpha }_{\mathrm{S}}}\left(\underline{w}\delta \right)\frac{{\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right)}{{\tilde {K}}_{\mathrm{S}}^{2}\left(\underline{w}\right)}\propto {\delta }^{2{\alpha }_{\mathrm{S}}}\sum _{\begin{subarray}{c}\underline{w}\in \mathbb{L}\cap \mathcal{B}\\ \vert \vert \underline{w}\vert \vert { >}{\Gamma}\end{subarray}}\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}+2{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}+2{\alpha }_{\mathrm{S}}}\right)\hfill \\ \hfill & \qquad \approx {\delta }^{2{\alpha }_{\mathrm{S}}}{\int }_{{\Gamma}}^{1/\delta }\mathrm{d}w\enspace {w}^{d-1-{\alpha }_{\mathrm{T}}+2{\alpha }_{\mathrm{S}}}+o\left(\vert \vert \underline{w}\vert {\vert }^{-{\alpha }_{\mathrm{T}}+2{\alpha }_{\mathrm{S}}}\right)=\mathcal{O}\left({\delta }^{\mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right)}\right).\hfill \end{align} \tag{ 45 }$$

Finally, putting equations (42), (44) and (45) together,

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}=\mathcal{O}\left({\delta }^{\mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right)}\right).\end{equation} \tag{ 46 }$$

The proof is concluded by considering that $\delta =\mathcal{O}\left({n}^{-1/d}\right)$.

In the case of a Gaussian kernel $K\left(\underline{x}\right)\propto \mathrm{exp}\left(-\vert \vert \underline{x}\vert {\vert }^{2}/\left(2{\sigma }^{2}\right)\right)$—and therefore $\tilde {K}\left(\underline{w}\right)\propto \mathrm{exp}\left(-{\sigma }^{2}\vert \vert \underline{w}\vert {\vert }^{2}/2\right)$—one has to redo the calculations starting from equation (39), but the final result can be easily recovered by taking the limit α → +∞ (Gaussian kernels decay faster than any power law). □

Lemma 1. Let $K\left(\underline{x},{\underline{x}}^{\prime }\right)$ be a translation-invariant isotropic kernel such that $\tilde {K}\left(\underline{w}\right)=c\vert \vert \underline{w}\vert {\vert }^{-\alpha }+o\left(\vert \vert \underline{w}\vert {\vert }^{-\alpha }\right)$ as ||w|| → ∞ and $\vert \vert \underline{w}\vert {\vert }^{d}\tilde {K}\left(\underline{w}\right)\to 0$ as ||w|| → 0. If α > d + n for some $n\in {\mathbb{Z}}^{+}$, then $K\left(\underline{x}\right)\in {C}^{n}$, that is, it is at least n-times differentiable.

Proof. The kernel is rotational invariant in real space ($K\left(\underline{x}\right)=K\left(\vert \vert \underline{x}\vert \vert \right)$) and therefore also in the frequency domain. Then, calling ${\hat{{\epsilon}}}_{1}=\left(1,0,\dots \enspace \right)$ the unitary vector along the first dimension x₁,

$$\begin{equation}K\left(x\right)\propto \int \mathrm{d}\underline{w}\enspace {\text{e}}^{\text{i}\underline{w}\cdot {\hat{{\epsilon}}}_{1}x}\enspace \tilde {K}\left(\vert \vert \underline{w}\vert \vert \right).\end{equation} \tag{ 47 }$$

It follows that

$$\begin{align}\hfill \left\vert {\partial }^{m}K\left(x\right)\right\vert \propto & \left\vert \int \mathrm{d}\underline{w}\enspace {\left(\underline{w}\cdot {\hat{{\epsilon}}}_{1}\right)}^{m}{\text{e}}^{\text{i}\underline{w}\cdot {\hat{{\epsilon}}}_{1}x}\enspace \tilde {K}\left(\vert \vert \underline{w}\vert \vert \right)\right\vert {< }\int \mathrm{d}\underline{w}\enspace \vert \underline{w}\cdot {\hat{{\epsilon}}}_{1}{\vert }^{m}\vert \tilde {K}\left(\vert \vert \underline{w}\vert \vert \right)\vert \hfill \\ \hfill \propto & {\int }_{0}^{\infty }\mathrm{d}w\enspace {w}^{d-1+m}\vert \tilde {K}\left(w\right)\vert {\int }_{0}^{\pi }\mathrm{d}{\phi }_{1}\vert \mathrm{cos}\left({\phi }_{1}\right){\vert }^{m}\propto {\int }_{0}^{\infty }\mathrm{d}w\enspace {w}^{d-1+m}\vert \tilde {K}\left(w\right)\vert .\hfill \end{align} \tag{ 48 }$$

We want to claim that this quantity is finite if m ⩽ n. Convergence at infinity requires m < α − d, that is always smaller than or equal to n because of the hypothesis of the lemma. Convergence in zero requires that ${w}^{d+m}\vert \tilde {K}\left(w\right)\vert \to 0$, and we want this to hold for all 0 ⩽ m < α − d, the most constraining one being the condition with m = 0. □

Lemma 2. Let $Z\sim \mathcal{N}\left(0,K\right)$ be a d-dimensional Gaussian random field, with K ∈ C²ⁿ being a 2n-times differentiable kernel. Then Z is n-times differentiable in the sense that

• Derivatives of $Z\left(\underline{x}\right)$ are a Gaussian random fields;
• $\mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)=0$;
• $\mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{1}}^{{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}^{\prime }}Z\left({\underline{x}}^{\prime }\right)={\partial }_{{x}_{1}}^{{n}_{1}+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }\right){< }\infty$ if the derivatives of K exist.

In particular, $\mathbb{E}{\partial }_{{x}_{i}}^{m}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{i}}^{m}Z\left({\underline{x}}^{\prime }\right)={\partial }_{{x}_{i}}^{2m}K\left(\underline{x}-{\underline{x}}^{\prime }\right){< }\infty \quad \forall \enspace m{\leqslant}n$.

Proof. Derivatives of $Z\left(\underline{x}\right)$ are defined as limits of sums and differences of the field Z evaluated at different points; therefore they are Gaussian random fields too, and furthermore it is straightforward to see that their expected value is always 0 if the field itself is zero centered.

The correlation can be computed via induction. Assume that $\mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{1}}^{{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}^{\prime }}Z\left({\underline{x}}^{\prime }\right)={\partial }_{{x}_{1}}^{{n}_{1}+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }\right)$ holds true. Then, if we increment n₁:

$$\begin{align}\hfill & \mathbb{E}{\partial }_{{x}_{1}}^{{n}_{1}+1}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)\cdot {\partial }_{{x}_{1}}^{{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}^{\prime }}Z\left({\underline{x}}^{\prime }\right)\hfill \\ \hfill & \qquad ={\mathrm{lim}}_{h\to 0}\enspace {h}^{-1}\mathbb{E}\left[{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}+h{\hat{{\epsilon}}}_{1}\right)-{\partial }_{{x}_{1}}^{{n}_{1}}\dots {\partial }_{{x}_{d}}^{{n}_{d}}Z\left(\underline{x}\right)\right]\cdot {\partial }_{{x}_{1}}^{{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}^{\prime }}Z\left({\underline{x}}^{\prime }\right)\hfill \\ \hfill & \qquad ={\mathrm{lim}}_{h\to 0}\enspace {h}^{-1}\left[{\partial }_{{x}_{1}}^{{n}_{1}+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }+h{\hat{{\epsilon}}}_{1}\right)\right.\hfill \\ \hfill & \qquad \quad \left.-\enspace {\partial }_{{x}_{1}}^{{n}_{1}+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }\right)\right]\hfill \\ \hfill & \qquad ={\partial }_{{x}_{1}}^{{n}_{1}+1+{n}_{1}^{\prime }}\dots {\partial }_{{x}_{d}}^{{n}_{d}+{n}_{d}^{\prime }}K\left(\underline{x}-{\underline{x}}^{\prime }\right).\hfill \end{align} \tag{ 49 }$$

Of course, by symmetry the same can be said about the increase of any other exponent. To conclude the induction proof we simply recall that by definition $\mathbb{E}Z\left(\underline{x}\right)Z\left({\underline{x}}^{\prime }\right)=K\left(\underline{x}-{\underline{x}}^{\prime }\right)$. □

It is important to note the high degree of smoothness underlying the RKHS hypothesis. Consider for instance realizations $Z\left(\underline{x}\right)$ of a teacher Gaussian process with covariance K_T and assume that they lie in the RKHS of the student kernel K_S (notice that they never belong to the RKHS of the same kernel K_T), namely

$$\begin{equation}\mathbb{E}\vert \vert Z\vert {\vert }_{{K}_{\mathrm{S}}}^{2}=\mathbb{E}\int \mathrm{d}\underline{x}\enspace \mathrm{d}\underline{\enspace }yZ\left(\underline{x}\right){K}_{\mathrm{S}}^{-1}\left(\underline{x}-\underline{y}\right)Z\left(\underline{y}\right)=\int \mathrm{d}\underline{w}\enspace {\tilde {K}}_{\mathrm{T}}\left(\underline{w}\right){\tilde {K}}_{\mathrm{S}}^{-1}\left(\underline{w}\right){< }\infty .\end{equation} \tag{ 50 }$$

If the teacher and student kernels decay in the frequency domain with exponents α_T and α_S respectively, convergence requires α_T > α_S + d, and ${K}_{\mathrm{S}}\left(\underline{0}\right)\propto \int \mathrm{d}\underline{w}\enspace {\tilde {K}}_{\mathrm{S}}\left(\underline{w}\right){< }\infty$ (true for many commonly used kernels) implies α_S > d. Then using lemmas 5.1 and 5.2 we can conclude that the realizations $Z\left(\underline{x}\right)$ must be at least ⌊d/2⌋-times mean-square differentiable to be RKHS.

In this appendix we show how to recover the prediction of theorem 1 using the approach presented in [11], that can be applied to a generic target function (not necessarily Gaussian nor evaluated only on a regular lattice). They derive their main formula, that we write below, by computing the amount of generalization error due to each eigenmode of the (student) kernel. In order to carry out the calculations they introduce a partial differential equation that they solve with two different approximations.

In the following we denote by λ₁ ⩾ ⋯ ⩾ λ_ρ ⩾ ⋯ the eigenvalues of the kernel, and by ${\phi }_{\rho }\left(\underline{x}\right)$ the corresponding eigenfunctions. In [11] they show that the generalization error can be written as:

$$\begin{align}\hfill & \mathbb{E}\mathrm{M}\mathrm{S}\mathrm{E}=\sum _{\rho }{E}_{\rho }\left(n\right),\hfill \end{align} \tag{ 51 }$$

$$\begin{align}\hfill & {E}_{\rho }\left(n\right)=\sum _{\rho }\frac{\mathbb{E}{w}_{\rho }^{2}}{{\lambda }_{\rho }}{\left(\frac{1}{{\lambda }_{\rho }}+\frac{n}{t\left(n\right)}\right)}^{-2}{\left(1-\frac{n\gamma \left(n\right)}{t{\left(n\right)}^{2}}\right)}^{-1},\hfill \end{align} \tag{ 52 }$$

$$\begin{align}\hfill & \gamma \left(n\right)=\sum _{\rho }\frac{{\lambda }_{\rho }^{2}}{{\left(1+{\lambda }_{\rho }\frac{n}{t\left(n\right)}\right)}^{2}},\hfill \end{align} \tag{ 53 }$$

$$\begin{align}\hfill & t\left(n\right)=\sum _{\rho }\frac{{\lambda }_{\rho }}{1+{\lambda }_{\rho }\frac{n}{t\left(n\right)}}.\hfill \end{align} \tag{ 54 }$$

The term $\mathbb{E}{w}_{\rho }^{2}$ is the variance of the coefficients of the target function in the kernel eigenbasis, defined as:

$$\begin{equation}{w}_{\rho }={\lambda }_{\rho }^{-1/2}\left\langle Z,{\phi }_{\rho }\right\rangle ,\end{equation} \tag{ 55 }$$

(the factor in front of the scalar product is to keep our notation consistent with that of [11]), to help the reader compare the two works. Notice that the variance can be computed with respect to an ensemble of target functions, but this ensemble may contain one deterministic function only.

In order to compute sums over the eigenmodes we will always replace them with integrals over eigenvalues. To do so, we must also introduce a density of eigenvalues $\mathcal{D}\left(\lambda \right)$: ${\sum }_{\rho }f\left({\lambda }_{\rho }\right)\to \int \mathrm{d}\lambda \enspace \mathcal{D}\left(\lambda \right)f\left(\lambda \right)$. The asymptotic behavior of this density for small eigenvalues can be derived as follows (for a given kernel whose Fourier transform decays with an exponent α):

$$\begin{equation}\mathcal{D}\left(\lambda \right)=\sum _{\rho }\delta \left(\lambda -{\lambda }_{\rho }\right)\sim \int {\mathrm{d}}^{d}\underline{w}\enspace \delta \left(\lambda -\vert \vert \underline{w}\vert {\vert }^{-\alpha }\right)\sim {\int }_{0}^{\infty }\mathrm{d}w\enspace {w}^{d-1}\delta \left(\lambda -{w}^{-\alpha }\right)={\lambda }^{-\theta },\end{equation} \tag{ 56 }$$

where we have defined the exponent $\theta \equiv 1+\frac{d}{\alpha }$. Notice that 1 < θ < 2, and that of course this exponent depends on the kernel. We can use this density also to derive a scaling behavior of small eigenvalues: indeed, the ρth (≫1) eigenvalue can be estimated by

$$\begin{equation}\rho \sim {\int }_{{\lambda }_{\rho }}^{{\lambda }_{1}}\mathrm{d}\lambda \enspace \mathcal{D}\left(\lambda \right)\sim {\int }_{{\lambda }_{\rho }}^{{\lambda }_{1}}\mathrm{d}\lambda \enspace {\lambda }^{-\theta }\sim {\lambda }_{\rho }^{-\left(\theta -1\right)}.\end{equation} \tag{ 57 }$$

The last equation follows from the fact that λ_ρ ≪ λ₁ and that θ > 1.

We now have to estimate the asymptotic behavior of the implicitly defined function t(n). It is easy to see that this function must go to 0 as n → ∞; therefore we can assume it is small. Splitting the integral according to whether the denominator in the definition of t(n) is dominated by the first or second term,

$$\begin{equation}t\left(n\right)\sim \int \mathrm{d}\lambda \enspace {\lambda }^{-\theta }\frac{\lambda }{1+\lambda \frac{n}{t\left(n\right)}}\sim {\int }_{0}^{\frac{t\left(n\right)}{n}}\mathrm{d}\lambda \enspace {\lambda }^{-\theta }+{\int }_{\frac{t\left(n\right)}{n}}^{{\lambda }_{1}}\mathrm{d}\lambda \enspace {\lambda }^{-\theta }\frac{t\left(n\right)}{n}\sim {\left(\frac{t\left(n\right)}{n}\right)}^{2-\theta }.\end{equation} \tag{ 58 }$$

Therefore, $t\left(n\right)\sim {n}^{-\frac{2-\theta }{\theta -1}}$, and with a similar approximation we can also deduce that $\gamma \left(n\right)\sim {n}^{-\frac{3-\theta }{\theta -1}}$. Injecting all we know in the formula for the generalization error and splitting the integral we find

$$\begin{align}\hfill & \mathbb{E}\mathrm{M}\mathrm{S}\mathrm{E}\sim \sum _{\rho }\frac{\mathbb{E}{w}_{\rho }^{2}}{{\lambda }_{\rho }}{\left(\frac{1}{{\lambda }_{\rho }}+{n}^{\frac{1}{\theta -1}}\right)}^{-2}{\left(1-{n}^{0}\right)}^{-1}\sim \sum _{\rho }\frac{\mathbb{E}{w}_{\rho }^{2}}{{\lambda }_{\rho }}{\left(\frac{1}{{\lambda }_{\rho }}+\frac{1}{{\lambda }_{n}}\right)}^{-2}\hfill \\ \hfill & \qquad \sim {\lambda }_{n}^{2}\sum _{\rho {\leqslant}n}\frac{\mathbb{E}{w}_{\rho }^{2}}{{\lambda }_{\rho }}+\sum _{\rho { >}n}\mathbb{E}{w}_{\rho }^{2}{\lambda }_{\rho }.\hfill \end{align} \tag{ 59 }$$

In the second equality we have used the fact that ${\lambda }_{n}\sim {n}^{-\frac{1}{\theta -1}}$ to introduce the nth eigenvalue into the formula. Then, we approximated the sum by splitting it in two sums, one over the first n eigenvalues (ρ ⩽ n, therefore λ_ρ ⩾ λ_n ) and one over the remaining ones (ρ > n). Notice that the second sum is indeed the sum in equation (12).

Next we assume that $\mathbb{E}{w}_{\rho }^{2}$ behaves asymptotically as a power law with respect to small eigenvalues, $\mathbb{E}{w}_{\rho }^{2}\sim {\lambda }_{\rho }^{q}$, with an exponent q that can be either positive or negative. We can now compute each of the integrals in the previous equation:

$$\begin{align}\hfill {\lambda }_{n}^{2}\sum _{\rho {\leqslant}n}\frac{\mathbb{E}{w}_{\rho }^{2}}{{\lambda }_{\rho }}& \sim {\lambda }_{n}^{2}{\int }_{{\lambda }_{n}}^{{\lambda }_{1}}\mathrm{d}\lambda \enspace {\lambda }^{-\theta }{\lambda }^{q-1},\hfill \end{align} \tag{ 60 }$$

$$\begin{align}\hfill \sum _{\rho { >}n}\mathbb{E}{w}_{\rho }^{2}{\lambda }_{\rho }& \sim {\int }_{0}^{{\lambda }_{n}}\mathrm{d}\lambda \enspace {\lambda }^{-\theta }{\lambda }^{q+1}\sim {\lambda }_{n}^{q-\theta +2}\sim {n}^{-\frac{q-\theta +2}{\theta -1}}.\hfill \end{align} \tag{ 61 }$$

For the second integral to converge we have assumed that the exponent q is larger than θ − 1. The first integral behaves differently according to whether q > θ or not: if q > θ, the integral scales as ${\lambda }_{n}^{2}\sim {n}^{-\frac{2}{\theta -1}}$; if q < θ, then it scales as ${\lambda }_{n}^{2-\theta +q}\sim {n}^{-\frac{q-\theta +2}{\theta -1}}$. Therefore,

$$\begin{equation}\mathbb{E}\enspace \mathrm{M}\mathrm{S}\mathrm{E}\sim {n}^{-\frac{\mathrm{min}\left(q-\theta ,0\right)+2}{\theta -1}}.\end{equation} \tag{ 62 }$$

A consequence of equations (60) and (61) is that if q < θ (which always occurs if the student is smooth enough, so that α characterizing the decay of the Fourier coefficient is small and θ is large), then the scaling of the generalization error is given by equation (61) alone, and we recover equation (12) from equation (59), justifying why this equation applies to real, non-Gaussian data.

Notice that if the target function is generated by a teacher Gaussian process, the exponent q takes the value $\frac{\theta -{\theta }_{\mathrm{T}}}{{\theta }_{\mathrm{T}}-1}$, where ${\theta }_{\mathrm{T}}=1+\frac{d}{{\alpha }_{\mathrm{T}}}$ and α_T is the exponent characterizing the decay of the Fourier transform of the teacher kernel. With some manipulations we then recover our theorem 1.

$$\begin{equation}\mathbb{E}\mathrm{M}\mathrm{S}\mathrm{E}\sim {n}^{-\frac{1}{d}\enspace \mathrm{min}\left({\alpha }_{\mathrm{T}}-d,2{\alpha }_{\mathrm{S}}\right)}.\end{equation} \tag{ 63 }$$

See figure J1.

Zoom In Zoom Out Reset image size