LEARNING CURVES FOR DEEP NEURAL NETWORKS: A FIELD THEORY PERSPECTIVE
Year: 2019

Abstract:
A series of recent works established a rigorous correspondence between very wide deep neural networks (DNNs), trained in a particular manner, and noiseless Bayesian Inference with a certain Gaussian Process (GP) known as the Neural Tangent Kernel (NTK). Here we leverage this correspondence to provide explicit analytical predictions for the learning curves of DNNs on any fixed target function. Focusing on datasets whose input measure is uniform on a hypersphere, our predictions match experimental curves very well and reveal a strong implicit bias towards functions which are low order polynomials of the input. Two novel analytical tools underlay our approach. First a renormalization-group approach is used to show that noiseless GP inference using NTK, which lacks a good analytical handle, can be well approximated by noisy GP inference on a related kernel we call the renormalized NTK. Second a field theory reformulation of the problem allows a controlled perturbative expansion of the test error in the inverse dataset size. We believe our approach lays the foundations for a more precise quantitative understanding of DNNs.

heading: 1 INTRODUCTION
text: Several pleasant features underlay the success of deep learning: The scarcity of bad minima encountered in their optimization (Draxler et al., 2018; Choromanska et al., 2014), their ability to generalize well despite being heavily over-parameterized (Neyshabur et al., 2018; 2014) and expressive (Zhang et al., 2016), and their ability to generate internal representations which generalize across different domains and tasks (Yosinski et al., 2014; Sermanet et al., 2013).
Due to the complexity of DNNs our current understanding of these features is still largely empirical. Notwithstanding, progress has been made recently in the highly over-parametrized regime (Daniely et al., 2016; Jacot et al., 2018) due to the fact that the networks’ parameters, in all non-linear layers, change in a minor yet important manner during training. This facilitated the derivation of various bounds (Allen-Zhu et al., 2018; Cao & Gu, 2019b;a) on generalization and, more relevant for this work, the following correspondence with GPs: Considering finite-depth DNNs which are much wider than the dataset-size, trained with MSE loss, no weight decay, and at vanishing learning rate (the NTK-regime) one finds that the initialization-averaged predictions are the same as those of Gaussian Processes Regression (GPR) with a kernel known as the NTK. Several subsequent works corroborated these results empirically (Lee et al., 2018; Lee et al., 2019; Arora et al., 2019) and extended them (Arora et al., 2019). For fully-connected DNNs, the NTK-regime (and GPs associated with DNNs in general (Lee et al., 2018; Novak et al., 2018)) seems to faithfully capture the generalization power of DNNs trained with MSE loss (Lee et al., 2019).
One of the most detailed objects quantifying generalization are learning-curves: graphs of how the test error diminishes with the number of datapoints (N ). There are currently no analytical predictions or bounds we are aware of for DNN learning-curves which are tight even just in terms of their scaling with N , let alone tight in an absolute sense. In contrast, the theory of GPR is rich with analytical tools which have yielded in the past high accuracy predictions for learning curves. One of the most transparent ones is the equivalence kernel (EK) (Rasmussen & Williams, 2005) which will be introduced in section 3. In short, EK gives an intuitive functional prediction of the expected GPR predictor of a fixed target function when averaged over datasets of size N , and consequently of the learning curves.
Clearly such a detailed understanding of generalization in DNNs is desirable. However, several technical issues prohibit the application of the EK and related results (Rasmussen & Williams, 2005; Malzahn & Opper, 2001) to DNNs trained in the NTK-regime. First, the NTK-regime corresponds to noiseless GPR where the DNN and corresponding GP both fit the training dataset exactly, whereas EK assumes Gaussian measurement noise on the training labels. In the case of noiseless GPR, EK and related approximations break down and predict perfect generalization with very small datasets. Second, the eigenfunctions (features) and eigenvalues of the NTK are needed so that the EK can be interpreted. Third, as EK results can be misleading, it is important to estimate the validity range of these approaches and, in a related manner, derive sub-leading corrections.
In this work we make the following contributions:
I We establish that noiseless GP inference using NTK can be well approximated by noisy GP inference on a related set of kernels we dub the renormalized NTKs.
II We obtain closed expression for the leading and sub-leading asymptotics of learning curves for any fixed target function (fixed-teacher learning curves) using an extension of the field-theory formalism of Malzahn & Opper (2001) for GPR.
III For uniform datasets, where the input is distributed uniformly on a hypersphere, these expression simplify considerably and, together with the features and eigenvalues we obtain for the renormalized NTKs, lead to clear relations between deep fully-connected networks and polynomial regression.
IV Again for uniform datasets, we provide analytical predictions for learning curves along with estimates for their range of validity. We believe our learning curves estimates stand-out in terms of accuracy as they get to within 3% of experimental generalization values. We emphasize that our predictions can be applied to fully-connected DNNs of any depth, trained in the NTK regime, and that our approach provides a systematic way of making further accuracy improvements.

heading: 2 PRIOR WORKS
text: Learning curves for GPs have been analyzed using a variety of techniques (see Rasmussen & Williams (2005) for a review) most of which focus on a GP-teacher averaged case where the target/teacher is drawn from the same GP used for inference (matched priors) and is furthermore averaged over. Fixedteacher or fixed-target learning curves have been analyzed using a similar grand-canonical/Poissonaveraged approach (Malzahn & Opper, 2001) as our, however, the treatment of the resulting partition function was variational whereas we take a perturbation-theory approach. In addition previous cited results for MSE-loss breakdown in the noiseless limit (Malzahn & Opper, 2001). To the best of our knowledge, noiseless GPs learning-curves have been analyzed analytically only in the teacheraveraged case and in the following settings: For matched priors, exact results are known for one dimensional data (Williams & Vivarelli, 2000; Rasmussen & Williams, 2005) and two dimensional data with some limitations of how one samples the inputs (in the context of optimal design) (Ritter, 2007; 1996). In addition Micchelli & Wahba (1979) derived a lower bound on generalization. For noiseless inference with partially mismatched-priors (matching features, mismatching eigenvalues) and at large input dimension the teacher and dataset averaging involved in obtained learning curves has been performed analytically and the resulting matrix traces analyzed numerically Sollich (2001). Notably none of these cited results apply in any straightforward manner in the NTK-regime.
Considering kernel eigenvalues, explicit expression for the features and eigenvalues of dot-product kernels (K = K(x ·x′)) were given in Azevedo & Menegatto (2015). The fact that the l-th eigenvalue of such kernels scales as d−l (d being the input dimension), which we used in our derivation of the bound, has been noticed in Sollich (2001). Kernels with a trimmed spectrum where the spectrum is trimmed after the first r’s leading eigenvalues, has previously been suggested as a way of reducing the computational cost of GP inference (Ferrari-Trecate et al., 1998). In contrast we trim the Taylor expansion of the kernel function rather than the spectrum (which has a very different effect) and show that an effective observation noise compensates for our trimming/renormalization procedure.
Several interesting recent works give bounds on generalization (Allen-Zhu et al., 2018; Cao & Gu, 2019b;a) which show O(1/ √ N) asymptotic decay of the learning-curve (at best). In contrast our predictions are typically well below this bound.

heading: 3 FIELD THEORY FORMULATION OF GP LEARNING-CURVES

heading: 3.1 GPR DEFINITION AND PROBLEM STATEMENT
text: We begin with standard definitions of GPs and Bayesian Inference on GPs. A GP is defined as a stochastic process of which any finite subset of random variables follow a multivariate normal distribution. In a similar fashion to multivariate normal variables, GPs are also determined by their first and second moments. The first is typically taken to be zero, and second is known as the covariance function or the kernel Kxx′ = E[f(x)f(x′)], where E[·] here denotes expectation under the GP distribution. The main appeal of GPs is that Bayesian Inference with GP priors is tractable. In GPR we use the posterior mean as the predictor g?, and it is given by:
g?(x?) = N∑ n,m=1 Kx?,xn [K(D) + σ 2I]−1nmgm (1)
where x? is a new datapoint, gm are the training targets, xn are the training data-points, [K(D)]nm = Kxn,xm is the covariance-matrix (the covariance-function projected on the training dataset D), and σ2 is the variance of the assumed Gaussian noise of the labels, which also acts as a regulator of the prediction. Some intuition for this formula can be gained by verifying that in the noiseless case (σ2 = 0) the prediction at some training point x? = xq coincides with that point’s label g? = gq .
The quantity of interest in this paper will be the dataset averaged generalization error which we define now. Throughout this paper we will assume that both train and test points are i.i.d. random variables drawn from a probability measure µ(x). With this in mind we define the generalization error of a prediction g? or MSE loss as
‖g(x?)− g?(x?)‖2 = ∫ dµx? (g(x?)− g?(x?)) 2 (2)
Note that g? is itself a function of N draws from µ which make up the training set DN . The dataset averaged generalization error is (2) averaged over the ensemble of all possible N sized training sets. We denote this average as 〈·〉DN , so for instance the dataset averaged generalization error would be 〈‖g? − g‖2〉DN . We see that in order to calculate learning-curves, one needs to average quantities like g? and g?2. These averages turn out to be difficult to handle analytically, so to facilitate their computation we adopt the approach of Malzahn & Opper (2001) and instead consider a related quantity given by the Poisson averaging of the former one
〈...〉η = e−η ∞∑ n=0 ηn n! 〈...〉Dn (3)
where ... can be any quantity, in particular g? and g?2. Borrowing jargon from physics we refer to the original data ensemble as the canonical ensemble and to the above as the grand-canonical ensemble. Taking η = N , means we are essentially averaging over values of N in an √ N vicinity of N . This means that as far as the leading asymptotic behavior is concerned, one can safely exchange N and η as the differences would be sub-leading. In App. A we compare learning curves as a function of N and η and show that they match very well.
Equation (1) determines the predictions, and therefore the learning-curves, but it is not very convenient for analytic exploration of the expected predictions. This fact is due to the (potentially very) large matrix inversion involved, and the additional averaging over DN required. Nonetheless there are some results for the expected prediction g?, the most famous of which is the equivalence kernel (EK) result for the prediction of a fixed target function as a function of the training set size N (Rasmussen & Williams, 2005):
g?EK,N (x) = ∑ n λn λn + σ2/N gnφn(x) (4)
Where λn and φn(x) here are the eigenvalues and eigenfunctions of the kernel w.r.t the input probability measure (µ) and g(x) = ∑ n gnφn(x) is the target. One notices immediately that the prediction breaks down completely in the noiseless case where (4) implies perfect estimation of the target with just one datapoint. To gain some intuition as to why having σ2 = 0 hinders predictions of 〈g?〉DN one can view it as a hard constraint (f(xn) = g(xn)), and hard constraints are typically less tractable than soft ones. Indeed this can also be seen as a motivation for considering the above Grand-canonical dataset ensemble. In a related view, finite σ2 can be seen as a form of averaging which smooths and regulates analytical expressions making them more tractable. Another limitation of the EK result is that (to the best of our knowledge) there is no systematic way to extend it in orders of 1/N and get a more detailed picture of GPR generalization.

heading: 3.2 FIELD THEORY FORMULATION
text: To facilitate the analysis of Eq. (1) we formulate the problem from a statistical-field-theory/pathintegral viewpoint (Schulman, 1996). These are well-studied, powerful approaches for performing integrations over a space of functions (the jargon is ”paths” when x in one dimensional and ”fields” when x is higher dimensional). This section is somewhat technical and the reader who is not interested in derivations is invited to skip to section 3.3. To get some familiarity with this formalism, consider first averages over the (centered) GP itself with no dataset. Using the path-integral formalism we define the probability density functional of a GP with kernel K:
P0[f ] = exp
( − 12 ∫ dxdx′f(x)K−1(x, x′)f(x′) )∫ Df̃ exp ( − 12 ∫ dxdx′f̃(x)K−1(x, x′)f̃(x′)
) (5) where ∫ Df denotes integration over the space of (well behaved) functions, for concreteness we
limit ∫ dx′ to some compact domain such as the hyper-sphere, K−1(x, x′) is the inverse covariance
function ( ∫ dx′K(x, x′)K−1(x′, x′′) = δ(x− x′′)). P0 is analogous to a PDF with the sample space being the set of real functions. To define the path-integrals one first chooses an orthonormal basis of functions φi(x) (with respect to ∫ dx) arranged in order of likeliness P0[φi] ≥ P0[φj ] for i > j (note that this comparison doesn’t require calculating the path integral in (5)). Second, one expands f = ∑ i fiφi(x), and defines the path-integral as a series of simple integrals
∫ DfF [f ] = ∫ df1 ∫ df2 . . .F [∑ i fiφi ] (6)
where F is some functional of f (like P0 in (5)). Finally, one makes this last expression well-defined by taking a limit procedure where the number of integrals is gradually taken to infinity (Schulman, 1996).
Performing the above procedure we show in App. E that a stochastic process with functional probability density P0 is equivalent to the original GP. Following a similar procedure, and denoting ‖f‖2K = ∫ dxdx′f(x)K−1(x, x′)f(x′) one can show an alternative representation of (1)
g?(x?) = Z −1 ∫ Df · f(x?) · exp ( −1
2 ‖f‖2K −
1
2σ2 N∑ n=1
(f(xn)− gn)2 )
(7)
Z = ∫ Df exp ( −1
2 ‖f‖2K −
1
2σ2 N∑ n=1
(f(xn)− gn)2 )
where Z is known as the partition function (see Rasmussen & Williams (2005) for equivalence of (1) and (7)). It is useful to define the partition function with a ”source term”:
Z[α(x)] = ∫ Df exp ( −1
2 ‖f‖2K +
∫ dxα(x)f(x)− 1
2σ2 N∑ n=1
(f(xn)− gn)2 )
(8)
Then the GPR prediction is simply
g?(x?) = ∂ log(Z[α(x?)])
∂α(x?)
∣∣∣∣ α(x?)=0
(9)
where ∂α(x)∂α(x?) = δ(x− x?). This form makes it clear that we can get quantities like dataset averaged predictions by finding the average of logZ. By using the grand canonical ensemble, averaging over draws from the dataset can be carried using the ”replica trick” (see for instance (Gardner & Derrida, 1988)), which aids in averaging over expressions like log(Z) and their derivatives via the equality log(Z) =M→0 ZM−1 M . Specifically:
〈g?(x?)〉η = 〈 lim M→0
∂
∂α(x?)
ZM − 1 M ∣∣∣∣ α(x?)=0 〉 η
(10)
〈ZM 〉η = ∫ M∏
m=1
Dfm exp
( −1
2 M∑ m=1 ‖fm‖2K + η ∫ dµxe − 1 2σ2 ∑M m=1(fm(x)−g(x)) 2 )
where, as standard in the replica formalism, the computation should be carried at positive integer M and the analytical result extrapolated to zero at the end. Neatly, a Taylor expansion in η of the above r.h.s. yields the 〈...〉η averaging appearing on the l.h.s. The main benefit of (10) over (1) is that it allows for a controlled expansion in 1/η. At large η (or similarly large N ) we expect the fluctuations in fm(x) to be small and centered around g(x). Indeed such a behavior is encouraged by the term multiplied by η in the exponent. We can therefore systematically Taylor expand the inner exponent
e− ∑M m=1(fm(x)−g(x)) 2 2σ2 = 1− ∑M m=1(fm(x)− g(x))2
2σ2 +
1
2
[∑M m=1(fm(x)− g(x))2
2σ2
]2 + ... (11)
and each term will yield a higher order of 〈g?(x?)〉η in 1/η. The first order of the systematic expansion we defined can be dealt with in an exact manner, in which case we recovers the aforementioned EK result (4) with the slight difference that N gets replaced with η (see App. F for the computation). The second order term and further terms render the theory non-Gaussian and cannot be dealt with exactly but rather through the use of perturbation-theory/Feynman-diagrams.

heading: 3.3 FIELD THEORY FORMULATION PREDICTIONS
text: Performing the calculation described in the previous section we obtain the following closed expression for large N or η behavior of noisy GPR using any kernel, in particular kernels associated with deep fully-connected networks and convolutional networks.
〈g?(x?)〉η = g?EK,η(x?) + g?SL,η(x?) +O(1/η3) (12)
g?SL,η(x?) = η
σ4 ∑ i,j,k σ2 η λi + σ2 η ( 1 λj + η σ2 )−1( 1 λk + η σ2 )−1 giφj (x?) ∫ dµxφi (x)φj (x)φ 2 k (x)
The first term is the familiar EK results but with N replaced by η. The second term (g?SL,η) is a sub-leading correction which captures the behavior at lower N ’s.
As shown App. F similar expressions for 〈g?2〉η are obtained using two replica indices. Interestingly we find that 〈g?2〉η = 〈g?〉2η + O(1/η3). Hence the averaged MSE error is simply (〈g?(x?)〉η − g(x?))
2 integrated over x?. Since the variance of g? came out to be O(1/η3) one finds that g? − g, which isO(1/η), is asymptotically much larger than its standard deviation. This implies self averaging at large η, or equivalently that our dataset-averaged results capture the behavior of a single fixed dataset.
Equation (12) and the average MSE error are one of our key results. They provides us with closed expressions for the dataset-averaged MSE loss as a function of η namely, the fixed-teacher learning curve. They hold without any limitations on the dataset or the kernel and yield a variant of the EK result along with its sub-leading correction. From an analytic perspective, once λi and φi(x) are known, the above expressions provide clear insights to how well the GP learns each feature and what unwanted cross-talk is generated between features due to the second sub-leading term. Notably for the renormalized NTK introduced below, the number of non-zero λi’s is finite, and so the above infinite summations reduce to finite ones. This makes these expressions computationally superior to directly performing the matrix-inversion in (1) along with an N−dimensional integral involved in dataset-averaging. In addition having the sub-leading correction allows us to estimate the range of validity of our approximation by comparing the sub-leading and leading contributions, as we shall do for the uniform case below.

heading: 4 UNIFORM DATASETS
text: To make the result (12) interpretable, φi(x) and λi are required. This can be done most readily for the case of datasets normalized to the hypersphere (‖xn‖ = 1) with a uniform probability measure and rotation-symmetric kernel functions. By the latter we mean Kx,x′ = KOx,Ox′ for any O, where O is an orthogonal matrix over the space of inputs. Although beyond the scope of the current work obvious extensions to consider are datasets which are uniform only in a sub-space of x and/or small perturbations to uniformity.
Importantly, the NTK associated with any DNN with a fully connected first layer and weights initialized from a normal distribution, has the above symmetry under rotations. This follows from the recursion relations defining the NTK (Jacot et al., 2018) along with fact that the kernel of the first fully-connected layer is only a function of x · x′. It follows that the NTK can be expanded as Kx,x′ = ∑ n bn(x · x′)n. An additional corollary (Azevedo & Menegatto, 2015) is that its features are hyperspherical harmonics (Ylm(x)) as these are the features of all dot product kernels. Hyperspherical harmonics are a complete (and orthonormal w.r.t a uniform measure) basis for functions on the hypersphere. For each l these can be written as a sum of polynomials in the input coordinates of degree l. The extra index m enumerates an orthogonal set of such polynomials (of size deg(l)). 1 For a kernel of the above form the eigenvalues are independent of m and given by (Azevedo & Menegatto, 2015)
λl = Γ ( d 2 ) √ π · 2l ∞∑ s=0 b2s+l (2s+ l)! (2s)! Γ ( s+ 12 ) Γ ( s+ l + d2
) (13) For ReLU and erf activations, the bn’s, can be obtained analytically up to any desirable order. Thus one can semi-analytically obtain the NTK eigenvalues up to any desired accuracy. For the particular case of depth 2 ReLU networks, we report in the App. H closed expression where the above summation can be carried out analytically. However as we shall argue soon, it is in fact desirable to trim the NTK in the sense of cutting-off its Taylor expansion at some order m, resulting in what we call the renormalized NTK. For such kernels, which would be our main focus next, the above result can be seen as a closed analytical expression for the eigenvalues.
Interestingly, for any fully-connected network and uniform datasets of dimension d on the hypersphere, there is a universal bound given by λl ≤ K/deg(l) ≈ O(d−l), where K is Kx,x which is a constant in x. Indeed note that Kx,x is finite and therefore its integral over the hypersphere is also finite and given by ∫ dµxKx,x = Kx,x = ∑ lm λl = ∑ l deg(l)λl. The degeneracy (deg(l)) is fixed from
properties of hyper spherical harmonics, and equals deg(l) = 2l+d−2l+d−2 ( l+d−2 l ) (Frye & Efthimiou, 2012) which goes as O(dl) for l d. This combined with the positivity of the λl’s implies the above bound.
1Note that usually the hyperspherical harmonics are normalized w.r.t Lebesgue measure on the hypersphere, but in this context the normalization is w.r.t a probability measure on the hypersphere.
Expressing our target in this feature basis g(x) = ∑ l,m glmYlm(x) (12) simplifies to
g? = g?EK,η(x?)− ∑ l,m [ η−1CK,σ2/η λl + σ2/η λl λl + σ2/η ] glmYlm(x?) (14)
where CK,σ2/η = ∑ lm(λ −1 l + η/σ
2)−1 and notably cross-talk between features has been eliminated at this order since ∑ m φlm(x) 2 is constant yielding ∑ m̃ ∫ dµxφlm(x)φl′m′(x)φ 2 l̃m̃
(x) ∝ δll′δmm′ . By splitting the sum, CK,σ2/η , to cases in which λl < σ2/η and their complement one has the bound CK,σ2/η < #Fσ 2/η + ∑ lm|λl>σ2/η λl, where #F is the number of non-zero kernel eigenvalues. Thus for kernels with a finite number of non-zero λi’s, as the renormalized NTK introduced below, CK,σ2/η has a η−1 asymptotic. This illustrates the fact the above terms are arranged by their orders in η.
Taking the leading order term one obtains the aforementioned EK result with N replaced by η. Equating the two contributions provides an estimate of when perturbation theory breaks down. Focusing on λl > σ2/η, the perturbation theory appears valid when CK,σ2/η σ2. In the limit σ2 → 0, and for trimmed kernels, this yield #F η. Notably it means that the original non-trimmed NTK cannot be analyzed perturbatively in the noiseless limit. In the next section we tackle this issue.

heading: 5 GENERALIZATION IN THE NOISELESS CASE AND THE RENORMALIZED NTK
text: The correspondence between DNNs trained in the NTK regime and GPR using NTK implies noiseless GPR (σ2 = 0) for which the perturbative analysis carried in previous sections fails. Here we show that the fluctuations of fm(x) associated with low λ’s can be traded for noise on the fluctuations of fm(x) associated with high λ’s thereby making our perturbative analysis applicable. As shown in the previous section, the lower λ’s correspond to higher spherical harmonics and hence have higher Fourier components. We argue that these higher Fourier modes can be marginalized over in a controlled manner to generate both noise and corrections to the high λ’s. This is very much in spirit of the renormalization group technique common in physics wherein high Fourier modes are integrated over to generate changes (renormalization) of some parameters in the probability distribution of the low Fourier modes.
We begin by defining a set of renormalized NTKs. As argued, an NTK of any fully-connected DNN can be expanded as Kx,x′ = ∑∞ q=0 bq(x · x′)q. The renormalized NTK at scale r is then simply K (r) x,x′ = ∑r q=0 bq(x · x′)q. Harmoniously with this notation we denote the prediction of GPR with the original kernel as g?∞. Our claim is that GPR with K and a noise of σ 2 can be well approximated by GPR with K(r) and noise σ2 + σ2r (where σ 2 r = ∑∞ q=r+1 bq), for sufficiently large r. Specifically that the discrepancy between GPR predictions scales as O( √ Nd−(r+1)/2/Kx1,x1), where d is the effective data-input dimension. As can be inferred from (13), the renormalized NTK at scale r has zero eigenvalues for all spherical Harmonics with l > r. Thus, as advertised, these high Fourier modes have been removed from the problem. In a related manner trimming the Taylor expansion after (x · x′)r effectively reduces our angular resolution and coarse grains the fine angular features captured by these spherical Harmonics with l > r.
To justify this approximation we consider the difference matrix K(r)xn,xm −Kxn,xm , first for n 6= m and xn’s drawn for a uniform distribution of a hypersphere of dimension d. The terms bq(xn · xm)q scale roughly as dq/2 (see App. H for a more accurate expression) due to the tendency of random vectors in high dimensions to be orthogonal. Consequently the above difference diminishes very quickly with r. Notably this also applies for K(r)x∗,xm − Kx∗,xn provided x∗ is a test point and not a train point. In contrast the diagonal part of the different matrix is simply σ2rδn,m and may diminish more slowly depending on the coefficients bq>r. Upon neglecting K (r) x∗,xm −Kx∗,xn and the off-diagonal elements of the difference matrix, one finds that (1) with these two GPRs yields identical predictions. As shown in App. H, these neglected off-diagonal elements yield a discrepancy which scales as √ Nd−(r+1)/2 (since they sum incoherently). Consequently the MSE error between the two GPRs should scale as N . This scaling with N should saturate when the accuracy is nearly perfect since then the predictions remain largely constant as N is increased.
Focusing back on the question of how to tackle noiseless GPR, we thus find that as long as the bq’s decays slowly enough with q, then at any finite N we can choose a large enough r such that two desirable properties are maintained: A. The discrepancy between the GPRs is small and B. σ2r is large enough to ensure convergence to our perturbative analysis. The required slow decay of bq is harmonious with the intuition that DNNs should be initialized at the edge of chaos (Schoenholz et al., 2016) where the output of the network has a fine and multi-scale sensitivity to small changes in the input. As Kx,x′ is the correlation of two outputs with inputs x and x′, having a power law decaying bq implies such fine and multi-scale sensitivity. Establishing relations between good initialization and effectiveness of our renormalized NTK is left for future work.
We have tested the accuracy of approximating noiseless NTK GPR with renormalized NTK GPR with the appropriate σ2r , both on artificial datasets (see next section) and on real world dataset such as CIFAR10 (see App B.). In both cases we found an excellent agreement between the two GPRs for r’s as small as 3 and 4.

heading: 6 GENERALIZATION IN THE NTK REGIME
text: Collecting the results of all the preceding sections, we can obtain a detailed and clear picture of generalization in fully connected DNNs trained in the NTK-regime on datasets with a uniform distribution normalized to some hypersphere in input space. We begin with a qualitative discussion and consider some renormalized NTK at scale r. From Sec. 4, we have that the features of this kernel are hyperspherical harmonics and that λl scales as d−l. We also recall that l is the maximal degree of the polynomial appearing in the Harmonic and all Harmonics up to the degree l span all polynomial on the hypersphere with degree up to l. Examining (14) we find that features with λl σ2/η are learnable and via the above scaling we find that we learn polynomials of degree O(log(η/σ2)/ log(d)) or less. In particular a function like parity, which is a polynomial of degree d is very hard to learn whereas a linear function is the easiest to learn. Thus despite having infinitely more parameters than data-points (due to infinite width) and despite being able to span almost any function (due to the richness of the kernel’s features), deep neural networks avoid over fitting by having a strong bias towards low degree polynomials.
To make more quantitative statements we now focus on a specific setting. We consider input data in dimension d = 50 and a scalar target function g(x) = ∑ l=1,2;m glmYlm(x) such that∑
l=1,m g 2 lm = ∑ l=2,m g 2 lm = 1/2, but otherwise iid glm’s. We generate several toy datasets DN consisting of N data points (xn) uniformly distributed on the hypersphere Sd−1 and their corresponding targets (g(xn)). We consider training a fully-connected DNN consisting of 4 layer with ReLU activations and width W which we initialize with variance (σ2w = σ 2 b = 1/d) for the input layer and (σ2w = σ 2 b = 1/W ) for the hidden layers (see for instance Lee et al. (2019) App. C and App. E for how to compute the kernel). To be in the NTK correspondence regime we consider training such a network at vanishing learning-rate, MSE loss, and with W N . One then has that the predictions of the DNN are given by GPR with σ2 = 0 and the K given by the NTK kernel (Jacot et al., 2018) 2.
For each such DNN we obtained the expected MSE loss ‖g?∞ − g‖ 2 of GPR with the NTK kernel by numerical integration over x? Repeating this process multiple times we obtained the dataset averaged loss 〈‖g?∞ − g‖
2〉DN for N = 1, 2, . . . , Nmax with a relative standard error of less then 5% (this typically required averaging over 10 datasets). For direct comparison with our prediction of the learning curve, we computed the Poisson averaged learning curve 〈‖g?∞ − g‖
2〉η in accordance with (3), neglecting the terms n > Nmax. We restricted ourselves to ηmax ≤ Nmax − 5 √ Nmax to make tail effects negligible. Notably the Poisson averaging makes the final statistical error negligible relative to the discrepancies coming from our large η approximations (see A).
Since the target function involves l = 0, 1, 2 the minimal scale for the renormalized NTK is r ≥ 2. To have some headroom we start from r = 3 which implies reasonable discrepancies of the MSE between the two GPRs of the order ofN/d4 = 5.6e−4, for the maximalN we have (N = 3500). Our
2To be more precise, Jacot et al. (2018) predict correspondence with GPR up to a random initialization factor, so to get exact match with GPR one would also need to average over initialization seeds. Recent research (Lee et al., 2019) suggests this caveat can be avoided under some conditions.
analytical expressions following (13) combined with known results (Jacot et al., 2018; Cho & Saul, 2009) about the Taylor coefficients (bn) yield λ0, ..., λ3 = {3.19, 7.27e− 3, 5.98e− 6, 1.62e− 7} and σ2r = 0.018. Since λ0, λ1 σ2/η λ2, λ3 for 50 < η < 3500, CKr,σ2/ησ−2 < [deg(0) + deg(1)]σ2/η +O(deg(2)10−6), thus CKr,σ2/ησ
−2 ≈ 51/η. Thus we expect perturbation theory to be valid for η 50. At η = 100 the l = 1 features are learned well since σ2/η = 1.8e− 4 λ1 and the l = 2 features neglected, at η = 1000 they are learned but suppressed by a factor of about 3. Had the target contained l = 3 features, they would have been entirely neglected at these η scale. Experimental learning curves along with our leading and sub-leading estimates are shown in Fig. 1. left panel showing an excellent agreement between theory and experiment.
While no actual DNNs were trained in the above experiments, the NTK correspondence means that this would be the exact behavior of a DNN trained in the NTK regime (Jacot et al., 2018; Lee et al., 2019; Arora et al., 2019). Furthermore since our aim was to estimate what the DNNs would predict rather than reach SOTA predictions, we focus on reasonable hyper-parameters but did not perform any hyper-parameter optimization.
Lastly we argue that the asymptotic behavior of learning-curve we predict is more accurate than the recent PAC based bounds (Allen-Zhu et al., 2018; Cao & Gu, 2019b;a). In App. C we show a log-log plot of the learning-curves contrasted with a 1/ √ η which is the most rapidly decaying bound appearing in those works. It can be seen that such an asymptotic cannot be made to fit the experimental learning-curve with any precision close to ours.

heading: 7 DISCUSSION AND OUTLOOK
text: In this work we laid out a formalism based on field theory tools for predicting learning-curves in the NTK regime. Well within the validly regime of our perturbative analysis we find excellent agreement to within 3% relative mismatch between our best estimate and the experimental curves. Central to our analysis was the renormalization-group transformation on the NTK leading to effective observation noise on the target. Our analysis could be readily extend in several ways: Going beyond the uniform dataset case should be possible for multi-variate Gaussian input distribution with a set of similar finite variances and a set of nearly zero variances. Adding weak randomness to Kx,x′ to study the difference between empirical and averaged NTKs. It would also be interesting to extend our analysis to simple CNNs. The renormalized kernel can also be used for spectral analysis of the NTK and other kernels associated with DNNs.

heading: A POISSON AVERAGING DEMONSTRATION
text: Here we demonstrate that Poisson averaging has no substantial effect on the learning curve. To this end we show the experimental learning curve from the main text pre- and post-averaging. It is evident that other than the unintended consequence of eliminating the experimental noise, the averaged learning curve is equivalent to the original for all practical intents.

heading: B COMPARISON OF NTK AND RENORMALIZED NTK PREDICTIONS ON NON-UNIFORM DATASET
text: While our lack of knowledge of the NTK eigenvalues and eigenfunctions with respect to a non-uniform measure prevents us from predicting learning curves, we would like to show that the renormalized NTK is still a valid approximation in this setting. To this end we compare the prediction of the NTK and renormalized NTK on the one-hot encoding of the cifar-10 dataset.

heading: C COMPARISON WITH RECENT BOUNDS
text: As mentioned in the main text, various recent bounds, relevant to the NTK regime, have been derived recently. Notwithstanding importance and rigor of these works, their bounds have at best a 1/ √ N asymptotic scaling. Here we show that given a functional behavior of the experimental learning curves such a bound cannot be nearly as tight as our predictions.

heading: D NOTATIONS FOR THE FIELD THEORY DERIVATION.
text: For completeness, here we re-state the notations used in the main-text.
Ω - Space pace of inputs.
x, x′, x∗ - Inputs (in Ω).
µx - Measure on Ω.
K (x, x′) - Kernel function (covariance) of a Gaussian process. Assumed to be symmetric and positive-semi-definite.
φi (x) - i’th eigenfunction of K (x, x′). By the spectral theorem, the set {φi}∞i=1 can be assumed to be orthonormal: ∫
x∈Ω
dµxφi (x)φj (x) = δij
λi - i’th eigenvalue of K (x, x′).∫ x′∈Ω dµx′K (x, x ′)φi (x ′) = λiφi (x) || · ||HK - RKHS norm. If f (x) = ∑ i fiφi (x) then ||f ||HK = ∑ i f2i λi
(where φi is an orthonormal set). Note that this norm is independent of µx.
g (x) - The target function.
σ2 - Noise variance. The noise is assumed to be Gaussian.
N - Number of inputs in the data-set.
DN - Data-set of size N , DN = {x1, ..., xN}. f∗DN ,σ2 (x) - The prediction function.

heading: E EXPLICIT PATH INTEGRAL COMPUTATIONS
text: Here we wish to prove the probability function defined in (5) of the main text yields the GP defined by a given kernel using explicit computation of the path integral. Denoting∫ dµxdµx′f(x)K −1 (x, x′) f (x′) as ‖f‖2K and noting that ‖f‖ 2 K = ∑ i f2i λi :
Notably, all other higher correlation functions split into products of the above correlation function due to standard properties of Gaussian integrals (Wick’s/Isserlis’ theorem).
∫ Df · f (x) · f (y) · exp ( − 12 ‖f‖ 2 K ) ∫ Df exp ( − 12 ‖f‖ 2 K ) = ∫ ∏i dfi ·∑i fiφi (x) ·∑j fjφj (y) · exp ( − 12 ∑ l f2l λl ) ∫ ∏ i dfi exp ( − 12 ∑ l f2l λl
) = = ∑ i ∫ df · f2 · exp ( − f 2 2λi ) ∫ df exp ( − f22λi ) ︸ ︷︷ ︸
λi
φi (x)φi (y) + ∑ i 6=j
∫ df · f · exp ( − f 2
2λi ) ∫ df exp ( − f22λi ) ︸ ︷︷ ︸
0
·
∫ df · f · exp ( − f 2
2λj ) ∫ df exp ( − f22λj ) ︸ ︷︷ ︸
0
=
= ∑ i λiφi (x)φi (y) = K (x, y)

heading: F GAUSSIAN PROCESS PREDICTION AS A FIELD THEORY
text: Let us assume a Gaussian process (GP) with mean 0 and co-variance function K (x, x′). For a data-set DN of size N and noisy targets (σ2 6= 0) {g (xi)}Ni=1, it is known that the posterior mean obtained by Bayesian inference is
f∗DN ,σ2 = arg minf
[ 1
2 ‖f‖2HK + ∑ xi∈DN (f (xi)− g (xi))2 2σ2
]
For a data-set DN of size N and noisy targets {g (xi)}Ni=1 , we present the GP canonical partition function:
ZDN ,σ2 [α (x)] def =
∫ Df exp ( −1
2 ‖f‖2HK +
∫ α (x) f (x) dx− ∑ xi∈DN (f (xi)− g (xi))2 2σ2 )
Where the Df notation stands for path integral. Notice that the functional derivative of log ( ZDN ,σ2 [α (x)] ) w.r.t α (x∗) at α (x) = 0 yields:
∂
∂α (x∗) ∣∣∣∣ α(x)=0 log ( ZDN ,σ2 [α (x)] ) =
1 ZDN ,σ2 [α (x) = 0] · ∂ ∂α (x∗) ∣∣∣∣ α(x)=0 ( ZDN ,σ2 [α (x)] ) =
=
∫ Df · f (x∗) exp ( − 12 ‖f‖ 2 HK − ∑ xi∈DN (f(xi)−g(xi))2 2σ2 ) ∫ Df exp ( − 12 ‖f‖ 2 HK − ∑ xi∈DN (f(xi)−g(xi))2 2σ2 ) = arg min f |x∗ [ 1 2 ‖f‖2HK + N∑ i=1 (f (xi)− g (xi))2 2σ2 ]
where the last equality is due to the fact that for Gaussian distributions, the expected value coincides with the most probable value. Therefore, the exact baysian inference mean:
f∗DN ,σ2 (x ∗) =
∂
∂α (x∗) ∣∣∣∣ α(x)=0 log ( ZDN ,σ2 [α (x)] ) F.1 CANONICAL ENSEMBLE FORMALISM
for evaluating the quality of a certain GP, we’re interested in the average prediction for all the data-sets of size N , meaning:
fCN,σ2 (x ∗) def = 〈 f∗DN ,σ2 (x ∗) 〉 DN∼µNx = ∫ dµx1 ∫ dµx2 . . . ∫ dµxN f ∗ DN={x1,...,xN},σ2 (x ∗) .
Using the replica trick we obtain:
fCN,σ2 (x ∗) = lim
M→0
∂
∂α (x∗) ∣∣∣∣ α(x)=0
〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx − 1
M
Now, let us calculate 〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx for an integer M :
ZMDN ,σ2 [α (x)] =
∫ ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx− M∑ j=1 ∑ xi∈DN (fj (xi)− g (xi))2 2σ2  ZN,M,σ2 = 〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx = ∫ dµx1 · · · ∫ dµxNZ M DN ,σ2 [α (x)] =
= ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx 〈exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉N x∼µx
so
fCN,σ2 (x ∗) = lim
M→0
∂ ∂α(x∗) ∣∣∣ α(x)=0 〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx
M
F.2 GRAND CANONICAL ENSEMBLE FORMALISM
We now wish to allow fluctuations in the value of N , meaning averaging over fCN,σ2 (x ∗) for different values of N . The motivation is to simplify the calculations, while averaging around a relatively confined region of Ns. Let us average the canonical prediction while weighting N according to Poisson distribution with expected value η:
fGCη,σ2 (x ∗) def = ∞∑ N=0 e−ηηN N ! fCN,σ2 (x ∗) .
and defining:
Zη,M,σ2 [α (x)] = ∞∑ N=0 e−ηηN N ! 〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx
we get:
fGCη,σ2 (x ∗) =
∂
∂α (x∗) ∣∣∣∣ α(x)=0 lim M→0 Zη,M,σ2 [α (x)] M
That is, the functional derivative w.r.t α (x∗) at α (x) = 0 yields the average prediction, averaged over different data-set sizes (the canonical averaging) and different data-sets for each size (the grand canonical averaging).
For a given η, the standard deviation of N is η, so the relative error is 1√η , decreases with η. Substituting 〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx in the expression for Zη,M,σ2 [α (x)] we obtain:
〈 ZMDN ,σ2 [α (x)] 〉 DN∼µNx =
= ∞∑ N=0 e−ηηN N ! ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx  〈
exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉N x∼µx =
= ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx  ∞∑ N=0 e−ηηN N ! 〈 exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉N x∼µx =
= e−η ∫ ... ∫ ︸ ︷︷ ︸
M
Df1...DfM
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx+ η 〈 exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉 x∼µx 
F.3 DERIVING THE EQUIVALENCE KERNEL USING THE GRAND CANONICAL FORMALISM
We wish to get rid of the exponent inside the exponent. Expending it using (first order) Taylor series:
Zη,M,σ2 [α (x)] = = e−η ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx+ η 〈 exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉 x∼µx  ≈ ≈ e−η ∫ ...
∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx+ η 〈 1− M∑ j=1 (fj (x)− g (x))2 2σ2 〉 x∼µx  = = ∫ Df exp −1
2 ‖f‖2HK +
∫ α (x) f (x) dx− η 〈 (f (x)− g (x))2
2σ2
〉 x∼µx M def= (ZEKη,σ2 [α (x)])M
ZEKη,σ2 [α (x)] = ∫ Df exp ( −1
2 ‖f‖2HK +
∫ α (x) f (x) dµx − η
2σ2
∫ dµx (f (x)− g (x))2 )
and without any additional approximation:
fGCη,σ2 (x ∗) =
∂
∂α (x∗) ∣∣∣∣ α(x)=0 lim M→0 Zη,M,σ2 [α (x)]− 1 M ≈ ∂ ∂α (x∗) ∣∣∣∣ α(x)=0 lim M→0
( ZEKη,σ2 [α (x)] )M − 1
M =
= ∂
∂α (x∗) ∣∣∣∣ α(x)=0 log ( ZEKη,σ2 [α (x)] ) =
1 ZEKη,σ2 [α (x) = 0] · ∂ ∂α (x∗) ∣∣∣∣ α(x)=0 ( ZEKη,σ2 [α (x)] ) =
= 1 ZEKη,σ2 [α (x) = 0] · ∫ Df · f (x∗)
exp ( −1
2 ‖f‖2HK +
∫ α (x) f (x) dx− η
2σ2
∫ dµx (f (x)− g (x))2 ) α(x)=0 =
=
∫ Df · f (x∗) exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2 ) ∫ Df exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2
) = = arg min
f |x∗
[ 1
2 ‖f‖2HK +
η
2σ2
∫ dµx (f (x)− g (x))2 ] def = fEKη,σ2 (x ∗)
and that is exactly the result for the equivalence kernel, where η is the data-set size (we regarded it as the mean of the data-set size).
Let us derive it explicitly. For f (x) = ∑ i fiφi (x) and g (x) = ∑ i giφi (x):
fEKη,σ2 (x ∗) =
∫ Df · f (x∗) exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2 ) ∫ Df exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2
) = = ∫ ∏ i dfi · ∑ i fiφi (x ∗) · exp ( − 12 ∑ i ( f2i λi + ησ2 (fi − gi) 2 ))
∫ ∏ i dfi exp ( − 12 ∑ i ( f2i λi + ησ2 (fi − gi) 2 )) =
= ∑ i φi (x ∗)
∫ dfi · fi · exp ( − f 2 i
2λi − η2σ2 (fi − gi)
2 )
∫ dfi exp ( − f 2 i
2λi − η2σ2 (fi − gi)
2 ) = ∑
i
λi
λi + σ2
η
giφi (x ∗)
F.4 EQUIVALENCE KERNEL AS FREE FIELD THEORY
Regarding the Equivalence Kernel as the free (quadratic) theory, we can denote fEKη,σ2 (x ∗) = 〈f (x∗)〉f∼EK = 〈f (x∗)〉0 = ∑ i
λi λi+ σ2
η
giφi (x ∗). ‘
Let us calculate the correlations in the free theory:
〈f (x) f (y)〉0 =
∫ Df · f (x) f (y) exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2 ) ∫ Df exp ( − 12 ‖f‖ 2 HK − η 2σ2 ∫ dµx (f (x)− g (x))2
) = = ∫ ∏ i dfi · ∑ i,j fifjφi (x)φj (y) · exp ( − 12 ∑ i ( f2i λi + ησ2 (fi − gi) 2 ))
∫ ∏ i dfi exp ( − 12 ∑ i ( f2i λi + ησ2 (fi − gi) 2 )) =
= ∑ i
∫ dfi · f2i · exp ( − f 2 i 2λi − η2σ2 (fi − gi) 2 )
∫ dfi exp ( − f 2 i
2λi − η2σ2 (fi − gi)
2 )
︸ ︷︷ ︸ λ2 i g2 i
(λi+σ 2 η )
2 + (
1 λi + η σ2
)−1 φi (x)φi (y) +
∑ i 6=j λigiλjgj( λi + σ2
η
)( λj + σ2
η
)φi (x)φj (y) =
= ∑ i ( 1 λi + η σ2 )−1 φi (x)φi (y) + ∑ i,j λigiλjgj( λi + σ2
η
)( λj + σ2
η )φi (x)φj (y) ︸ ︷︷ ︸
〈f(x)〉0〈f(y)〉0
Therefore:
Cov [f (x) , f (y)] = ∑ i ( 1 λi + η σ2 )−1 φi (x)φi (y)
and we see that the correlations are O (
1 η
) .
For rotationally invariant kernel, we get that
Var [f (x)] = ∞∑ l=0 deg(l)∑ m=0 ( 1 λl + η σ2 )−1 Y 2l,m (x) def = CK,η,σ2
is independent of x since deg(l)∑ m=0 Y 2l,m (x) = deg (l)
so
CK,η,σ2 = ∞∑ l=0 deg(l)∑ m=0
1
λ−1l + η/σ 2
F.5 PERTUBATIVE CORRECTION FOR THE EQUIVALENCE KERNEL
F.5.1 AVERAGING f
Going to the next order in the expansion:
Zη,M,σ2 [α (x)] = = e−η ∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx+ η 〈 exp − M∑ j=1 (fj (x)− g (x))2 2σ2 〉 x∼µx  ≈ ≈ e−η ∫ ...
∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK + M∑ j=1 ∫ α (x) fj (x) dx  exp
η〈1− M∑ j=1 (fj (x)− g (x))2 2σ2 + 1 2  M∑ j=1 (fj (x)− g (x))2 2σ2 2〉 x∼µx  = = ∫ ...
∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp  M∑ j=1 ( −1 2 ‖fj‖2HK + ∫ α (x) fj (x) dx− η 2σ2 ∫ dµx (fj (x)− g (x))2 ) exp  η 8σ4 M∑ j=1 M∑ l=1 ∫ dµx (fj (x)− g (x))2 · (fl (x)− g (x))2
 Note that:
fGCη,σ2 (x ∗) =
∂
∂α (x∗) ∣∣∣∣ α(x)=0 lim M→0 Zη,M,σ2 [α (x)]− 1 M =
= lim M→0
(ZEKη,σ2 [α (x) = 0])M · ∂ ∂α(x∗) ∣∣∣ α(x)=0
Zη,M,σ2 [α (x)]( ZEKη,σ2 [α (x) = 0]
)M · 1M  =
= lim M→0
( ZEKη,σ2 [α (x) = 0] )M︸ ︷︷ ︸ 1 · lim M→0
 ∂∂α(x∗) ∣∣∣ α(x)=0
Zη,M,σ2 [α (x)]( ZEKη,σ2 [α (x) = 0]
)M · 1M 
Calculating the first order pertubative corrections:
fGCη,σ2 (x ∗) = lim
M→0
∂ ∂α(x∗) ∣∣∣ α(x)=0 Zη,M,σ2 [α (x)]
M · ( ZEKη,σ2 [α (x) = 0] )M = = lim M→0 1
MZEKη,σ2 [0] M
∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj
exp −1 2 M∑ j=1 ‖fj‖2HK − η 2σ2 M∑ j=1 ∫ dµx (fj (x)− g (x))2  exp  η 8σ4 M∑ j=1 M∑ l=1 ∫ dµx (fj (x)− g (x))2 · (fl (x)− g (x))2  · M∑ i=1 fi (x ∗) =
= lim M→0
1
MZEKη,σ2 [0] M
∫ ... ∫ ︸ ︷︷ ︸
M
M∏ j=1 Dfj exp −1 2 M∑ j=1 ‖fj‖2HK − η 2σ2 M∑ j=1 ∫ dµx (fj (x)− g (x))2  1 + η
8σ4 M∑ j=1 M∑ l=1 ∫ dµx (fj (x)− g (x))2 · (fl (x)− g (x))2  · M∑ i=1 fi (x ∗) +O ( 1 η3 ) =
= lim M→0
1
M 〈1 + η 8σ4 M∑ j=1 M∑ l=1 ∫ dµx (fj (x)− g (x))2 · (fl (x)− g (x))2  · M∑ i=1 fi (x ∗) 〉 f1,...,fM∼EK
+O
( 1
η3
) =
= 〈f (x∗)〉0 + limM→0 1 M 〈 η 8σ4 M∑ j=1 M∑ l=1 ∫ dµx (fj (x)− g (x))2 · (fl (x)− g (x))2  · M∑ i=1 fi (x ∗) 〉 f1,...,fM∼EK
+O
( 1
η3
) =
= 〈f (x∗)〉0
+ lim M→0
1
M
η
8σ4
∫ dµx 〈 M∑ j=1 M∑ l=1 M∑ i=1 (fj (x)− g (x))2 · (fl (x)− g (x))2 fi (x∗) 〉 f1,...,fM∼EK +O ( 1 η3 ) .
Calculating the correction:
〈 M∑ j=1 M∑ l=1 M∑ i=1 (fj (x)− g (x))2 · (fl (x)− g (x))2 fi (x∗) 〉 f1,...,fM∼EK =
= M 〈 (f (x)− g (x))4 f (x∗) 〉
0 +M (M − 1) [ 2 〈 (f (x)− g (x))2 〉
0
〈 (f (x)− g (x))2 f (x∗) 〉 0 + 〈 (f (x)− g (x))4 〉 0 〈f (x∗)〉0 ] +M (M − 1) (M − 2) 〈 (f (x)− g (x))2 〉2 0 〈f (x∗)〉0
Note that we eventually divide by M and take the limit M → 0, so we only care about O (M) terms:
fGCη,σ2 (x ∗) = fEKη,σ2 (x ∗)
+ η
8σ4
∫ dµx [〈 (f (x)− g (x))4 f (x∗) 〉 0 − 2 〈 (f (x)− g (x))2 〉 0 〈 (f (x)− g (x))2 f (x∗) 〉 0
− 〈 (f (x)− g (x))4 〉
0 〈f (x∗)〉0 + 2
〈 (f (x)− g (x))2 〉2 0 〈f (x∗)〉0 ] +O ( 1 η3 ) These correlations can be calculated using Feynman diagrams, since the free theory (EK) is quadratic (Gaussian):
〈 (f (x)− g (x))4 f (x∗) 〉 0 =
= 3fEKη,σ2 (x ∗) Var [f (x)] 2 + 6fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )2 Var [f (x)]
+fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )4 + 4 ( fEKη,σ2 (x)− g (x) )3 Cov [f (x) , f (x∗)]
+12 ( fEKη,σ2 (x)− g (x) ) Var [f (x)] Cov [f (x) , f (x∗)]
〈 (f (x)− g (x))4 〉 0 〈f (x∗)〉0 =
= 3fEKη,σ2 (x ∗) Var [f (x)] 2 +6fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )2 Var [f (x)] + fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )4
〈 (f (x)− g (x))2 〉 0 〈 (f (x)− g (x))2 f (x∗) 〉 0 =
= 2Var [f (x)] Cov [f (x) , f (x∗)] ( fEKη,σ2 (x)− g (x) ) + fEKη,σ2 (x ∗) Var [f (x)] 2
+2Cov [f (x) , f (x∗)] ( fEKη,σ2 (x)− g (x) )3 +2fEKη,σ2 (x ∗) Var [f (x)] ( fEKη,σ2 (x)− g (x) )2 + fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )4
〈 (f (x)− g (x))2 〉2 0 〈f (x∗)〉0 =
= fEKη,σ2 (x ∗) Var [f (x)] 2
+2fEKη,σ2 (x ∗) Var [f (x)] · ( fEKη,σ2 (x)− g (x) )2 + fEKη,σ2 (x ∗) ( fEKη,σ2 (x)− g (x) )4
Summing everything up:
〈 (f (x)− g (x))4 f (x∗) 〉 0 − 〈 (f (x)− g (x))4 〉 0 〈f (x∗)〉0
−2 〈 (f (x)− g (x))2 〉
0
〈 (f (x)− g (x))2 f (x∗) 〉 0 + 2 〈 (f (x)− g (x))2 〉2 0 〈f (x∗)〉0 =
= 8 ( fEKη,σ2 (x)− g (x) ) Var [f (x)] Cov [f (x) , f (x∗)]
and all the bubble diagrams cancel as expected.
So we get:
fGCη,σ2 (x ∗) = fEKη,σ2 (x
∗) + η
σ4
∫ dµx ( fEKη,σ2 (x)− g (x) ) Var [f (x)] Cov [f (x) , f (x∗)] +O ( 1
η3
)
Substituting the expressions for the variance and covariance:
fGCη,σ2 (x ∗) =
= fEKη,σ2 (x ∗)− η
σ4 ∑ i,j,k σ2 η λi + σ2 η ( 1 λj + η σ2 )−1( 1 λk + η σ2 )−1 giφj (x ∗) ∫ dµxφi (x)φj (x)φ 2 k (x) +O ( 1 η3 )
F.5.2 AVERAGING f2
This time we must use two different replica indices:
〈[ f∗DN ,σ2 (x ∗) ]2〉
DN∼µNx = 〈 ∂ log (ZDN ,σ2 [α (x)]) ∂α (x∗) ∣∣∣∣∣ α(x)=0 2〉 DN∼µNx =
= 〈 lim M→0 1 M · ∂ZMDN ,σ2 [α (x)] ∂α (x∗) ∣∣∣∣∣ α(x)=0 2〉 DN∼µNx =
= lim M→0 lim M̃→0
1 MM̃ · ∫ Df1 . . . ∫ DfM ∫ Df̃1 . . . ∫ Df̃M̃ exp −1 2 M∑ m=1 ‖fm‖2HK − 1 2 M̃∑ m̃=1 ∥∥∥f̃m̃∥∥∥2 HK  M∑ m=1 fm (x ∗) M̃∑ m̃=1 f̃m̃ (x ∗) 〈 exp − M∑ m=1 (fm (x)− g (x))2 2σ2 − M̃∑ m̃=1 ( f̃m̃ (x)− g (x) )2 2σ2 〉 N
x∼µx
Averaging w.r.t poisson distribution:
〈〈[ f∗DN ,σ2 (x ∗) ]2〉
DN∼µNx 〉 N∼Poi(η) = ∞∑ N=0 e−ηηN N ! 〈[ f∗DN ,σ2 (x ∗) ]2〉 DN∼µNx =
= lim M→0 lim M̃→0
1 MM̃ · ∫ Df1 . . . ∫ DfM ∫ Df̃1 . . . ∫ Df̃M̃
exp −η − 1 2 M∑ m=1 ‖fm‖2HK − 1 2 M̃∑ m̃=1 ∥∥∥f̃m̃∥∥∥2 HK  exp η〈exp − M∑
m=1
(fm (x)− g (x))2
2σ2 − M̃∑ m̃=1
( f̃m̃ (x)− g (x) )2 2σ2 〉 x∼µx  M∑ m=1 fm (x ∗) M̃∑ m̃=1 f̃m̃ (x ∗) ≈
≈ lim M→0 lim M̃→0
1 MM̃ · ∫ Df1 . . . ∫ DfM ∫ Df̃1 . . . ∫ Df̃M̃
exp −1 2 M∑ m=1 ‖fm‖2HK − 1 2 M̃∑ m̃=1 ∥∥∥f̃m̃∥∥∥2 HK + η 〈− M∑ m=1 (fm (x)− g (x))2 2σ2 − M̃∑ m̃=1 ( f̃m̃ (x)− g (x) )2 2σ2 〉 x∼µx
+ η
2 〈− M∑ m=1 (fm (x)− g (x))2 2σ2 − M̃∑ m̃=1 ( f̃m̃ (x)− g (x) )2 2σ2  2〉
x∼µx
 M∑ m=1 fm (x ∗) M̃∑ m̃=1 f̃m̃ (x ∗) =
≈ lim M→0 lim M̃→0
1 MM̃ · ∫ Df1 . . . ∫ DfM ∫ Df̃1 . . . ∫ Df̃M̃
exp −1 2 M∑ m=1 ‖fm‖2HK − 1 2 M̃∑ m̃=1 ∥∥∥f̃m̃∥∥∥2 HK + η 〈− M∑ m=1 (fm (x)− g (x))2 2σ2 − M̃∑ m̃=1 ( f̃m̃ (x)− g (x) )2 2σ2 〉 x∼µx  1 + η2 〈 M∑ m=1 (fm (x)− g (x))2 2σ2 + M̃∑ m̃=1 ( f̃m̃ (x)− g (x) )2 2σ2  2〉
x∼µx
 M∑ m=1 fm (x ∗) M̃∑ m̃=1 f̃m̃ (x ∗) =
= ( fEKη,σ2 (x ∗) )2
+ lim M→0 lim M̃→0
1 MM̃ · η 8σ4
∫ dµx 〈 M∑ a=1 (fa (x)− g (x))2 + M̃∑ b=1 ( f̃b (x)− g (x) )22 M∑ c=1 fc (x ∗) M̃∑ d=1 f̃d (x ∗) 〉 0 =
= ( fEKη,σ2 (x ∗) )2
+ lim M→0 lim M̃→0
1 MM̃ · η 4σ4
∫ dµx 〈 M∑ a=1 (fa (x)− g (x))2 M∑ b=1 (fb (x)− g (x))2 M∑ c=1 fc (x ∗) M̃∑ d=1 f̃d (x ∗) 〉 0
+ 〈 M∑ a=1 (fa (x)− g (x))2 M̃∑ b=1 ( f̃b (x)− g (x) )2 M∑ c=1 fc (x ∗) M̃∑ d=1 f̃d (x ∗) 〉 0  = = ( fEKη,σ2 (x ∗) )2 + η
4σ4
∫ dµx lim
M→0
1
M 〈 M∑ a=1 (fa (x)− g (x))2 M∑ b=1 (fb (x)− g (x))2 M∑ c=1 fc (x ∗) 〉 0︸ ︷︷ ︸
8 ( fEK η,σ2 (x)−g(x) ) Var[f(x)]Cov[f(x),f(x∗)] as we saw in 〈f〉
fEKη,σ2 (x ∗)
+ η
4σ4
∫ dµx ( lim M→0 1 M 〈 M∑ a=1 (fa (x)− g (x))2 M∑ b=1 fb (x ∗) 〉 0 )2
and we’re left with:
lim M→0
1
M 〈 M∑ a=1 (fa (x)− g (x))2 M∑ b=1 fb (x ∗) 〉 0 =
lim M→0
1
M
[ M 〈 (f (x)− g (x))2 f (x∗) 〉
0 +M (M − 1) 〈f (x∗)〉
〈 (f (x)− g (x))2 〉 0 ] =
= 〈 (f (x)− g (x))2 f (x∗) 〉
0 − 〈f (x∗)〉0
〈 (f (x)− g (x))2 〉 0 = 2 ( fEKη,σ2 (x)− g (x) ) Cov [f (x) , f (x∗)]
but this correction gives O (
1 η3
) so:
〈〈[ f∗DN ,σ2 (x ∗) ]2〉
DN∼µNx 〉 N∼Poi(η) =
( fEKη,σ2 (x ∗) )2 + 2η
σ4 fEKη,σ2 (x ∗)
∫ dµx ( fEKη,σ2 (x)− g (x) ) Var [f (x)] Cov [f (x) , f (x∗)] +O ( 1
η3 ) and notably
〈 f2 〉 = 〈f〉2 +O ( 1
η3 ) F.6 PERUBATIVE CORRECTION FOR ROTATIONALLY INVARIANT KERNEL
We now wish to evaluate this expression for a rotationally invariant kernel and a uniform measure on the hypersphere. This simplyfies the expression for 〈f〉 to:
fGCη,σ2 (x ∗) =
= fEKη,σ2 (x ∗) +
η
σ4
∫ dµx ( fEKη,σ2 (x)− g (x) ) Var [f (x)] Cov [f (x) , f (x∗)] +O ( 1
η3
) =
= fEKη,σ2 (x ∗) +
η
σ4 CK,η,σ2
∫ dµx ( fEKη,σ2 (x)− g (x) ) Cov [f (x) , f (x∗)] =
= fEKη,σ2 (x ∗)− η
σ4 CK,η,σ2 ∑ i,j σ2 η λi + σ2 η ( 1 λj + η σ2 )−1 giφj (x ∗) ∫ dµxφi (x)φj (x)︸ ︷︷ ︸
δij
+O
( 1
η3
) =
= fEKη,σ2 (x ∗)− η
σ4 CK,η,σ2 ∑ i σ2 η λi + σ2 η ( 1 λi + η σ2 )−1 giφi (x ∗) +O ( 1 η3 ) =
= fEKη,σ2 (x ∗)− CK,η,σ2 ∑ l,m
gl.m σ2 (
2 + λlησ2 + σ2
λlη
)Yl,m (x∗) +O( 1 η3 ) .
The expression for 〈 f2 〉 is:
〈〈[ f∗DN ,σ2 (x ∗) ]2〉
DN∼µNx 〉 N∼Poi(η) =
= ( fEKη,σ2 (x ∗) )2 − 2fEKη,σ2 (x∗)CK,η,σ2 ∑
l,m
gl.m σ2 (
2 + λlησ2 + σ2
λlη
)Yl,m (x∗) +O( 1 η3 )
G VARIOUS INSIGHTS
G.1 CORRECTION MEANS WORSE GENERALIZATION
The correction always means worse generalization than what the EK suggests. Indeed
fGCη,σ2 (x ∗) = fEKη,σ2 (x ∗)− CK,η,σ2 ∑ l,m
gl,m σ2 (
2 + λlησ2 + σ2
λlη
)Yl,m (x∗) +O( 1 η3 ) =
= ∑ l,m λl λl + σ2 η gl,mYl,m (x ∗)− CK,η,σ2 ∑ l,m
gl,m σ2 (
2 + λlησ2 + σ2
λlη
)Yl,m (x∗) +O( 1 η3 ) =
= ∑ l,m  λl λl + σ2 η − CK,η,σ2 σ2 ( 2 + λlησ2 + σ2 λlη ) ︸ ︷︷ ︸
positive  ︸ ︷︷ ︸
< λl
λl+ σ2 η
<1
gl,mYl,m (x ∗)
G.2 EXACT EIGENVALUES FOR 2-LAYER RELU NTK WITH σ2b = 0
For the NTK associated with a 2-layer ReLU NTK without bias we were able to fined an exact expression for the eigenvalues for all l:
λ2k = σ2w1σ 2 w2 2π · d(1 + 2k) + (1− 2k) 2 8π
( Γ ( k − 12 ) Γ ( d 2 ) Γ ( k + d+12 ) )2 , λ2k+1 = σ2w1σ2w2 2π · π d δk,0
It is interesting to note that for all odd l > 1 λl = 0 so the expressive power of the kernel (and hence the neural network) is greatly reduced.

heading: H ACCURACY OF THE RENORMALIZED NTK
text: For two normalized datapoints x and x′, drawn from a uniform dataset on a hypersphere of radius 1, and at large d the random variable (x · x′) is approximately Gaussian with variance O(d−1). Since (x · x′) is bounded to [−1, 1], the random variable (x · x′)r must have a standard deviation which is decaying function of r. For r d and large d one can estimate the magnitude this standard deviation from exact known expressions and a saddle-point approximation yielding O((d/r)−r/2) ≈ O(d−r/2) 3. Considering next the tail of the Taylor expansion ∑ q>r bq(x · x′)q, projected on the dataset
( ∑ q>r bq(xn · xm)q). The resulting N by N matrix is ∑ q>r bq on the diagonal but O(d
−(r+1)/2) in all other entries. As we justified in the main text, our renormalization transformation amounts to keeping only the diagonal piece of this matrix and interpreting it as noise.
Consider then (1) for g? in two scenarios: (I) g?∞ with the full NTK (K(x, x ′)) and no noise and (II) g?r with the NTK trimmed after the r’th power (Kr(x, x′)) but with σ2r = ∑ q>r bq . The first K(x?, xn) piece, for x? drawn from the dataset distribution, obeys K(x?, xn)−Kr(x?, xn) = O(d−(r+1)/2). Next we compare Kr(xn, xm) + Inmσ2r and K(xn, xn). On their diagonal they agree exactly but their off-diagonal terms agree only up to a O(d−(r+1)/2) discrepancy. Denoting by δK the difference between these two matrices, we may expand K−1 = [Kr + σ2mI + δK] −1 = [Kr + σ 2 r I] −1[1 − δK[Kr + σ 2 r I] −1 + δK[Kr + σ 2 r I] −1δK[Kr + σ 2 r I] −1 + ...].
3A more accurate estimate is (
r r+d−3
)r/2 ( d−3
r+d−3
)d/4
We next argue that δK[Kr + σ2r I] −1 multiplied by target vector (g(xn)) is negligible compared to the identity for large enough r thereby establishing the equivalence of the two scenarios. Indeed consider the eigenvalues of δK[Kr + σ2r I]
−1. As δKnm is O(d−(r+1)/2) its typical eigenvalues are O( √ Nd−(r+1)/2) and bounded by O(Nd−(r+1)/2). The typical eigenvalues of [Kr + σ2mI]
−1 are of the same order as K(xn, xn) = K and bounded from below by σ2r . Thus typical eigenvalues of δK[Kr + σ 2 r I] −1 are O( √ Nd−(r+1)/2/K) and bounded from above by O(Nd−(r+1)/2/σ2r). The NTK has the desirable property that σ2r decays very slowly. Thus certainly in the typical case but even in the worse case scenario we expect good agreement at large r. In Fig. 1, right panel, we provide supporting numerical evidence.
We refer to Kr(x, x′) as the renormalized NTKs at the scale r. As follows from (13), λl’s with l ≥ r are zero. Therefore, as advertised, the high-energy-sector has been removed and compensated by noise on the target and a change of the remaining l < r (low-energy) eigenvalues. A proper choice of r involve two considerations. Requiring perturbation theory to hold well (CKr,σ2r/η < σ 2 r ) which puts an η-depended upper bound on r and requiring small discrepancy in predictions puts another η dependent lower bound on r (typically √ Nd−(r+1)/2 1).
Lastly we comment that our renormalization NTK approach is not limited to uniform datasets. The entire logic relies on having a rapidly decaying ratio of off-diagonal moments ((xn · xm)2r) and diagonal moments (xn · xn)2r as one increases r. We expect this to hold in real-world distributions. For instance for a multi-dimension Gaussian data distribution the input dimension (d) traded by an effective dimension (deff ) defined by the variance of (xm · xn).

Review summaries:
This paper explores how tools from perturbative field theory can be used to shed light on properties of the generalization error of Gaussian process/kernel regression, particularly on how the error depends on the number of samples N. For uniform data on the sphere, a controlled expansion is obtained in terms of the eigendecomposition of the kernel. Although the expansion breaks when the noise term goes to zero, a renormalized kernel is introduced for which an accurate perturbative expansion is possible. A variety of empirical results confirm the theoretical analysis.