Year: 2021

Abstract:
To understand how deep learning works, it is crucial to understand the training dynamics of neural networks. Several interesting hypotheses about these dynamics have been made based on empirically observed phenomena, but there exists a limited theoretical understanding of when and why such phenomena occur. In this paper, we consider the training dynamics of gradient flow on kernel leastsquares objectives, which is a limiting dynamics of SGD trained neural networks. Using precise high-dimensional asymptotics, we characterize the dynamics of the fitted model in two “worlds”: in the Oracle World the model is trained on the population distribution and in the Empirical World the model is trained on a sampled dataset. We show that under mild conditions on the kernel andL target regression function the training dynamics undergo three stages characterized by the behaviors of the models in the two worlds. Our theoretical results also mathematically formalize some interesting deep learning phenomena. Specifically, in our setting we show that SGD progressively learns more complex functions and that there is a “deep bootstrap” phenomenon: during the second stage, the test error of both worlds remain close despite the empirical training error being much smaller. Finally, we give a concrete example comparing the dynamics of two different kernels which shows that faster training is not necessary for better generalization.

heading: 1 INTRODUCTION
text: In order to fundamentally understand how and why deep learning works, there has been much effort to understand the dynamics of neural networks trained by gradient descent based algorithms. This effort has led to the discovery of many intriguing empirical phenomena (e.g. Frankle et al. (2020); Fort et al. (2020); Nakkiran et al. (2019a;b; 2020)) that help shape our conceptual framework for understanding the learning process in neural networks. Nakkiran et al. (2019b) provides evidence that SGD starts by first learning a linear classifier and over time learns increasingly complex functions. Nakkiran et al. (2020) introduces the “deep bootstrap” phenomenon: for some deep learning tasks the empirical world test error remains close to the oracle world error1 for many SGD iterations, even if the empirical training and test errors display a large gap. To better understand such phenomena, it is useful to study training dynamics in related but mathematically tractable settings.
One approach for theoretical investigation is to study kernel methods, which were recently shown to have a tight connection with over-parameterized neural networks (Jacot et al., 2018; Du et al., 2018). Indeed, consider a sequence of neural networks (fN (x;θ))N∈N with the widths of the layers going to infinity as N → ∞. Assuming proper parametrization and initialization, for large N the SGD dynamics on fN is known to be well approximated by the corresponding dynamics on the first-order Taylor expansion of fN around its initialization θ0,
fN,lin(x;θ) = fN (x;θ 0) + 〈∇θfN (x;θ0),θ − θ0〉.
Thus, in the large width limit it suffices to study the dynamics on the linearization fN,lin. When using the squared loss, these dynamics correspond to optimizing a kernel least-squares objective with the neural tangent kernel KN (x,x′) = 〈∇θfN (x;θ0),∇θfN (x′;θ0)〉.
1Their paper uses “Ideal World” for “Oracle World” and “Real World” for “Empirical World”.
Over the past few years, researchers have used kernel machines as a tractable model to investigate many neural network phenomena including benign overfitting, i.e., generalization despite the interpolation of noisy data (Bartlett et al., 2020; Liang & Rakhlin, 2020) and double-descent, i.e, risk curves that are not classically U-shaped (Belkin et al., 2020; Liu et al., 2021). Kernels have also been studied to better understand certain aspects of neural network architectures such as invariance and stability (Bietti & Mairal, 2017; Mei et al., 2021b). Although kernel methods cannot be used to explain some phenomena such as feature learning, they can still be conceptually useful for understanding other neural networks properties.

heading: 1.1 THREE STAGES OF KERNEL DYNAMICS
text: Despite much classical work in the study of gradient descent training of kernel machines (e.g. Yao et al. (2007); Raskutti et al. (2014)) there has been limited work understanding the high-dimensional setting, which is the setting of interest in this paper. Although solving the linear dynamics of gradient flow is simple, the statistical analysis of the fitted model requires involved random matrix theory arguments. In our analysis we study the dynamics of the Oracle World, where training is done on the (usually inaccessible) population risk, and the Empirical World, where training is done on the empirical risk (as is done in practice). Associated with the oracle world model fort and the empirical world model f̂t are the following quantities of interest: the empirical training error R̂n(f̂t), the empirical test error R(f̂t), and the oracle error R(fort ) defined in Eqs. (1), (2), (3) for which we derive expressions that are accurate in high dimensions.
Informally, our main results show that under reasonable conditions on the regression function and the kernel the training dynamics undergo the following three stages:
• Stage one: the empirical training error, the empirical test error, and the oracle error are all close. • Stage two: the empirical training error decays to zero, but the empirical test error and the oracle
error stay close and keep approximately constant. • Stage three: the empirical training error is still zero, the empirical test error stays approximately
constant, but the oracle test error decays to the approximation error.
We conceptually illustrate the error curves of the oracle and empirical world in Fig. 1 and provide intuition for the evolution of the learned models in Fig. 2. The existence of the first and third stages are not unexpected: at the beginning of training the model has not fit the dataset enough to distinguish the oracle and empirical world and at the end of training an expressive enough model with infinite samples will outperform one with finitely many. The most interesting stage is the second one where the empirical model begins to “overfit” the training set while still remaining close to the non-interpolating oracle model in the L2 sense (see Fig. 2).
In Section 2 we discuss some related work. In Section 3 we elaborate our description of the three stages and give a mathematical characterization for two particular settings in Theorems 1 and 2.
Although the three stages arise fairly generally, we remark that certain stages will vanish if the problem parameters are chosen in a special way (c.f. Remark 1). We connect our theoretical results to related empirical deep learning phenomena in Remark 3 and discuss the relation to deep learning in practice in Remark 4. In Section 4 we provide numerical simulations to illustrate the theory more concretely and in Section 5 we end with a summary and discussion of the results.

heading: 2 RELATED LITERATURE
text: The generalization error of the kernel ridge regression (KRR) solution has been well-studied in both the fixed dimension regime (Wainwright, 2019, Chap. 13), (Caponnetto & De Vito, 2007) and the high-dimensional regime (El Karoui, 2010; Liang & Rakhlin, 2020; Liu et al., 2021; Ghorbani et al., 2020; 2021; Mei et al., 2021a;b). Most closely related to our results is the setting of (Ghorbani et al., 2021; Mei et al., 2021a;b). Analysis of the entire KRR training trajectory has also been done (Yao et al., 2007; Raskutti et al., 2014; Cao et al., 2019) but only for the fixed dimensional setting. Classical non-parametric rates are often obtained by specifying a strong regularity assumption on the target function (e.g. the source condition in Fischer & Steinwart (2020)), whereas in our work the assumption on the target function is mild.
Another line of work directly studies the dynamics of learning in linear neural networks (Saxe et al., 2013; Li et al., 2018; Arora et al., 2019; Vaskevicius et al., 2019). Similar to us, these works show that some notion of complexity (typically effective rank or sparsity) increases in the linear network over the course of optimization.
The relationship between the speed of iterative optimization and gap between population and empirical quantities has been studied before in the context of algorithmic stability (Bousquet & Elisseeff, 2002; Hardt et al., 2016; Chen et al., 2018). These analyses certify good empirical generalization by using stability in the first few iterations to upper bound the gap between train and test error. In contrast, our analysis directly computes the errors at an arbitrary time t (c.f. Remark 2). The relationship between oracle and empirical training dynamics has been considered before in Bottou & LeCun (2004) and Pillaud-Vivien et al. (2018).

heading: 3 RESULTS
text: In this section we introduce the problem and present a specialization of our results to two concrete settings: dot product and group invariant kernels on the sphere (Theorems 1 and 2 respectively). The more general version of our results is described in Appendix A.3.

heading: 3.1 PROBLEM SETUP
text: We consider the supervised learning problem where we are given i.i.d. data (xi, yi)i≤n. The covariate vectors (xi)i≤n ∼iid Unif(Sd−1( √ d)) and the real-valued noisy responses yi = fd(xi) + εi for some unknown target function fd ∈ L2(Sd−1( √ d)) and (εi)i≤n ∼iid N (0, σ2ε). Given a function
f ∈ L2(Sd−1( √ d)), we define its test error R(f) and its training error R̂n(f) as
R(f) ≡ E(xnew,ynew){(ynew − f(xnew)) 2}, R̂n(f) ≡
1
n n∑ i=1 (yi − f(xi))2 (1)
where (xnew, ynew) is i.i.d. with (xi, yi)i≤n. The test error R(f) measures the fit of f on the population distribution and the training error R̂n(f) measures the fit of f to the training set. For a kernel function Hd : Sd−1( √ d)× Sd−1( √ d)→ R, we analyse the dynamics of the following two fitted models indexed by time t: the oracle model fort and the empirical model f̂t, which are given by the gradient flow on R and R̂n over the associated RKHSHd respectively
d dt fort (x) = −∇R(fort (x)) = E[Hd(x, z)(fd(z)− fort (z))], (2)
d dt f̂t(x) = −∇R̂n(f̂t(x)) = 1 n n∑ i=1 Hd(x,xi)(yi − f̂t(xi)), (3)
with zero initialization for0 ≡ f̂0 ≡ 0. These dynamics are motivated from the neural tangent kernel perspective of over-parameterized neural networks (Jacot et al., 2018; Du et al., 2018). A precise mathematical definition and derivation of these two dynamics are provided in Appendix E.1.
For our results we make some assumptions on the spectral properties of the kernels Hd similar to those in Mei et al. (2021a) that are discussed in detail in Section A.2. At a high-level we require that the diagonal elements of the kernel concentrate, that the kernel eigenvalues obey certain spectral gap conditions, and that the top eigenfunctions obey a hyperconctractivity condition which says they are “delocalized”. For the specific settings of Theorems 1 and 2 we give more specific conditions on the kernels that are more easily verified and imply the required spectral properties.

heading: 3.2 DOT PRODUCT KERNELS
text: In our first example, we consider dot product kernels Hd of the form
Hd(x1,x2) = hd(〈x1,x2〉/d), ∀x1,x2 ∈ Sd−1( √ d), (4)
for some function hd : [−1, 1] → R. Our results apply to general dot product kernels under weak conditions on hd given in Appendix C.2. In particular they apply to the random feature and neural tangent kernels associated to certain fully connected neural networks (Jacot et al., 2018).
Before presenting our results for this setting we introduce some notation. Denote by P≤` the orthogonal projection onto the subspace of L2(Sd−1( √ d)) spanned by polynomials of degree less than or
equal to `. The projectors P` and P>` are defined analogously (see Appendix G for details). We use od(·) for standard little-o relations, where the subscript d emphasizes the asymptotic variable. The statement f(d) = ωd(g(d)) is equivalent to g(d) = od(f(d)). We use od,P(·) in probability relations. Namely for two sequences of random variables Z1(d) and Z2(d), Z1(d) = od,P(Z2(d)) if for any ε, Cε > 0 there exists dε ∈ Z>0, such that P(|Z1(d)/Z2(d)| < Cε) ≤ ε for all d ≥ dε. The asymptotic notations Od, Od,P etc. are defined analogously. Theorem 1 (Dot Product Kernels). Let {fd ∈ L2(Sd−1( √ d))}d≥1 be a sequence of functions such that for some η > 0, ‖fd‖L2+η = Od(1), and let {Hd}d≥1 be a sequence of dot product kernels satisfying Assumption 4. Assume that for some fixed integers j, s ≥ 0 and some δ > 0 that
dj+δ ≤ t ≤ dj+1−δ, ds+δ ≤ n ≤ ds+1−δ. Then we have the following characterizations,
(a) (Oracle World) The oracle model learns every degree component of fd as time progresses R(fort ) = ∥∥P>jfd∥∥2L2 + σ2ε + od(1), ∥∥fort − P≤jfd∥∥2L2 = od(1).
(b) (Empirical World – Train) Empirical training error follows oracle error then goes to zero R̂n(f̂t) = ∥∥P>jfd∥∥2L2 + σ2ε + od,P(1) if t/n = od(1),
R̂n(f̂t) = od,P(1) if t/n = ωd(1).
(c) (Empirical World – Test) Empirical test error follows oracle error until the empirical model learns the degree-s component of fd
R(f̂t) = ∥∥P>min{j,s}fd∥∥2L2 + σ2ε + od,P(1), ∥∥∥f̂t − P≤min{j,s}fd∥∥∥2L2 = od,P(1).
The results are conceptually illustrated in an example in Fig. 3a which shows the stair-case phenomenon in high-dimensions and the three learning stages. We see that in both the oracle world and the empirical world, the prediction model increases in complexity over time. More precisely, the model learns the best polynomial fit to the target function (in an L2 sense) of increasingly higher degree. In the empirical world the maximum complexity is determined by the sample size n, which is in contrast to the oracle world where there are effectively infinite samples.
The results imply that generally (but not always c.f. Remark 1) there will be three stages of learning. In the first stage the oracle and empirical world models are close in L2 and fit a polynomial with degree determined by t. The first stage lasts from t = 0 to t = nd−ε n for some small ε > 0. As t approaches n, there is a phase transition and the empirical world training error goes to zero at t = ndε n. From time nd−ε till at least ds+1 is the second stage where the empirical and oracle models remain close in L2 but the gap between test and train error can be large. If the sample size n is not large enough for f̂t to learn the target function, then at some large enough t we will enter a third stage where fort improves in performance, outperforming f̂t which remains the same. On synthetic data in finite dimensions we can see a resemblance of the staircase shape which becomes sharper with increasing d (c.f. Appendix F).
Remark 1 (Degenerate stages of Learning). For special problem parameters we will not observe the second and/or third stages. The second stage will disappear if n dq for some q ∈ N (see Fig. 3b), or if fd is a degree-s polynomial. The third stage will not occur if P>sfd lies in the orthogonal complement of the RKHSHd.

heading: 3.3 GROUP INVARIANT KERNELS
text: We now consider our second setting which concerns the invariant function estimation problem introduced in Mei et al. (2021b). As before, we are given i.i.d. data (xi, yi)i≤n where the feature vectors (xi)i≤n ∼iid Unif(Sd−1( √ d)) and noisy responses yi = fd(xi) + εi. We now assume that the target function fd satisfies an invariant property. We consider a general type of invariance, defined by a group Gd that is represented as a subgroup of the orthogonal group in d dimensions. The group element g ∈ Gd acts on a vector x ∈ Rd via x 7→ g · x. We say that f? is Gd-invariant if f?(x) = f?(g · x) for all g ∈ Gd. We denote the space of square integrable Gd-invariant functions by L2(Sd−1( √ d),Gd). We focus on groups Gd that are groups of degeneracy α as defined below. As an example, we consider the cyclic group
Cycd = {g0, g1, . . . , gd−1} where for any x = (x1, . . . , xd)T ∈ Sd−1( √ d), the group action is defined by gi · x = (xi+1, xi+2, . . . , xd, x1, x2, . . . , xi)T. The cyclic group has degeneracy 1. Definition 1 (Groups of degeneracy α). Let Vd,k be the subspace of degree-k polynomials that are orthogonal to polynomials of degree at most (k − 1) in L2(Sd−1( √ d)), and denote by Vd,k(Gd) the subspace of Vd,k formed by polynomials that are Gd-invariant. We say that Gd has degeneracy α if for any integer k ≥ α we have dim(Vd,k/Vd,k(Gd)) dα (i.e., there exists 0 < ck ≤ Ck < +∞ such that ck ≤ dim(Vd,k/Vd,k(Gd)) ≤ Ck for any d ≥ 2).
To encode invariance in our kernel we consider Gd-invariant kernels Hd of the form
Hd,inv(x1,x2) = ∫ Gd h(〈x1, g · x2〉/d)πd(dg) (5)
where πd is the Haar measuare on Gd. Such kernels satisfy the following invariance property: for all g, g′ ∈ Gd and for Hd(x1,x2) = Hd(g · x1, g′ · x2) for every x1,x2. For the cyclic group, πd is the uniform measure. We now present our results for the group invariant setting.
Theorem 2 (Group Invariant Kernels). Let Gd be a group of degeneracy α ≤ 1 according to Definition 1. Let {fd ∈ L2(Sd−1( √ d),Gd)}d≥1 a sequence of Gd-invariant functions such that for some η > 0, ‖fd‖L2+η = Od(1), and let {Hd}d≥1 be a sequence of Gd-invariant kernels satisfying Assumption 5. Assume that for some fixed integers j ≥ 0, s ≥ 1 and some δ > 0 that
dj+δ ≤ t ≤ dj+1−δ, ds−α+δ ≤ n ≤ ds+1−α−δ.
Then we have the following characterizations,
(a) (Oracle World) The oracle model learns every degree component of fd as time progresses R(fort ) = ∥∥P>jfd∥∥2L2 + σ2ε + od(1), ∥∥fort − P≤jfd∥∥2L2 = od(1).
(b) (Empirical World – Train) Empirical training error follows oracle error then goes to zero R̂n(f̂t) = ∥∥P>jfd∥∥2L2 + σ2ε + od,P(1) if t/n = od(dα),
R̂n(f̂t) = od,P(1) if t/n = ωd(dα).
(c) (Empirical World – Test) Empirical test error follows oracle error until the empirical model learns the degree-s component of fd
R(f̂t) = ∥∥P>min{j,s}fd∥∥2L2 + σ2ε + od,P(1), ∥∥∥f̂t − P≤min{j,s}fd∥∥∥2L2 = od,P(1).
With respect to the dot product kernel setting (c.f. Theorem 1), in this setting the behavior of the oracle world is unchanged, but the empirical world behaves as if it has dα times as many samples. This is illustrated graphically in Fig. 3c. It can be shown that using an invariant kernel is equivalent to using a dot product kernel and augmenting the dataset to {(g · xi, yi) : g ∈ Gd, i ∈ [n]} (c.f. Appendix E.2). Hence for the cyclic group which has size dα = d, we arrive at the following intriguing conclusion: if the target function is invariant, then using a dot product kernel and augmenting n i.i.d samples to nd many samples is asymptotically equivalent to using nd i.i.d samples.
Remark 2 (Optimization Speed versus Generalization). Interestingly, training with an invariant kernel is slower than with a dot product kernel and takes longer to interpolate the dataset despite eventually generalizing better on invariant function estimation tasks (c.f. Fig. 3c). This conclusion is not an artifact of the continuous time analysis (c.f. Appendix E.3) and is observed empirically in Section 4.2 for discrete-time SGD. This example highlights the limitation of stability based analyses (e.g. Hardt et al. (2016)) which argue that faster SGD training leads to better generalization. While a faster rate leads to better generalization in the first stage when stability can control the gap between train and test error, our analysis shows that the duration length of the first stage also impacts the final generalization error.
Remark 3 (Connection with Deep Phenomena). The dynamics of high-dimensional kernel regression display behaviors that parallel some empirically observed phenomena in deep learning. For
kernel regression, we have shown that the complexity of the empirical model, measured as the number of learned eigenfunctions, depends on the time optimized when t n and the sample size when t n. At a high-level, we also expect a similar story for neural networks but for some other notion of complexity. It is believed that neural networks first learn simple functions and then progressively more complex ones, until the complexity saturates after interpolating at some time proportional to n (Nakkiran et al., 2019b). We have also shown that in kernel regression there is a non-trivial “deep boostrap” phenomenon (Nakkiran et al., 2020) during the second learning stage: the gap between the oracle world and empirical world test errors is negligible whereas the train and test errors exhibit a substantial gap. The gradient flow results for kernel regression can also provide insight into the deep bootstrap for random feature networks SGD as these results can approximately predict their behavior (see Section 4.2).
Remark 4 (Connection with Deep Learning Practice). Although we believe our results conceptually shed light on some of the interesting behaviors observed in the training dynamics of deep learning, due to our stylized setting we may not exactly see the predicted phenomena in practice. Accurately observing the three stages of kernel regression requires sufficiently high-dimensional data in order for the kernel eigenvalues to obey a staircase-like decay and for training to be sufficiently long as the time axis should be in log-scale. Our results hold for regression whereas for classification the empirical model may continue improving after classifying the train set correctly. Despite these caveats, certain conclusions can be observed in some realistic settings (c.f. Appendix E.4).

heading: 4 NUMERICAL SIMULATIONS
text: As mentioned previously, Fig. 3 is a “cartoon” of the conclusions stated in Theorems 1 and 2. In this section, we verify the qualitative predictions of our theorems using synthetic data. Concretely, throughout this section we take d = 400 and n = d1.5 = 8000, and following our theoretical setup (c.f. Section 3.1) generate covariates (xi)i≤n ∼iid Unif(Sd−1( √ d)) and responses yi = f?(xi)+εi with (εi)i≤n ∼iid N (0, σ2ε), for different choices of target function f?. All simulations in this section are for the noiseless case σ2ε = 0 but a noisy example is given in Appendix F.
In Section 4.1, we simulate the gradient flows of kernel least-squares with dot product kernels and cyclic kernels (Fig. 4) to reproduce the three stages as shown in Fig. 3. In Section 4.2 we show that SGD of (dot product and cyclic) random-feature models (Fig. 5) exhibit similar three stages phenomena, in which the second stage behaviors are consistent with the deep bootstrap phenomena observed in deep learning experiments (Nakkiran et al., 2020). Empirical quantities are averaged over 10 trials and the shaded regions indicate one standard deviation from the mean.

heading: 4.1 GRADIENT FLOW OF KERNEL LEAST-SQUARES
text: Under the synthetic data set-up mentioned earlier, we first simulate the oracle world and empirical world errors curves of gradient flow dynamics using dot product kernels of the form
H(x1,x2) = Ew∼Sd−1 [σ(〈w,x1〉)σ(〈w,x2〉)] (6)
for some activation function σ. We will examine a few different choices of f? and kernel H , which are specified in the descriptions of each figure.
The oracle world error is computed analytically but the empirical world errors require sampling train and test datasets. We compute empirical world errors by averaging over 10 trials. The results are visualized both in log-scale and linear-scale on the time axis. The log-scale plots allow for direct comparison with the cartoons in Fig. 3. The linear-scale plots are zoomed into the region 0 < t ≤ nd0.4 since: 1) plotting the full interval squeeze all curves to the left boundary which is uninformative 2) in practice one would not optimize for very long after interpolation.
In Figs. 4a and 4b we take the target function to be a polynomial of the form f?(x) = a0 He0(x1) + . . .+ ak Hek(x1), ∥∥Pjf?∥∥2L2 ≈ a2jj!, (7)
where Hei(t) is the ith Hermite polynomial (c.f. Appendix G.4) and the approximate equality in Eq. (7) holds in high-dimensions. In panel (4a) we consider use the ReLU activation function σ(t) = max(t, 0), and take f? to be a quadratic polynomial with (a0, a1, a2) = (1/2, 1/ √ 2, 1/ √ 8). With
such a choice of parameters, we can see the three stages phenomenon. In panel (4b) we choose f? to be a cubic polynomial with (a0, a1, a2, a3) = (1/2, 1/ √ 2, 0, 1/ √ 24) and σ(t) = max(t, 0) + 0.1 He3(t) (we need the third Hermite coefficient of σ to be non-zero for stage 3 to occur). This choice of coefficients for f? is such that ‖P>1f?‖2L2 ≈ ‖P>2f?‖2L2 , so that the second stage in (4b) is longer compared to (4a).
In Fig. 4c, we take the target function to be a cubic cyclic polynomial
f?(x) = 1√ 3d ( d∑ i=1 xi + d∑ i=1 xixi+1 + d∑ i=1 xixi+1xi+2 ) (8)
where the subindex addition in xi+k is understood to be taken modulo d. We compare the performance of the dot product kernel H and its invariant version Hinv (c.f. Eq. (5)) with activation function σ(t) = max(t, 0) + 0.1 He3(t). The kernel Hinv with n samples performs equivalently to H with nd samples (c.f. Remark 2), but is more computationally efficient since the size of the kernel matrix is still n× n. Using Hinv elongates the first stage by a factor d, delaying the later stages and ensuring that the empirical world model improves longer.
Although in the simulations, the dimension d is not yet high enough to see a totally sharp staircase phenomenon as in the illustrations of Fig. 3, even for this d we are still able to clearly see the three predicted learning stages and deep bootstrap phenomenon across a range of settings. To better understand the effect of dimension we show similar plots with varying d in Appendix F.

heading: 4.2 SGD FOR TWO-LAYER RANDOM-FEATURE MODELS
text: To more closely relate with deep learning practice and the deep bootstrap phenomenon (Nakkiran et al., 2020), we simulate the error curves of SGD training on random-feature (RF) models (i.e. twolayer networks with random first-layer weights and trainable second-layer weights), in the same synthetic data setup as before. In particular, we look at dot product RF models
f̂dot(x;a) = 1√ N N∑ i=1 aiσ(〈wi,x〉), wi ∼iid Unif(Sd−1),
and cyclic invariant RF models
f̂cyc(x;a) = 1√ N N∑ i=1 ai ∫ Gd σ(〈wi, g · x〉)πd(dg), wi ∼iid Unif(Sd−1).
For all following experiments we take the activation σ to be ReLU and N = 4× 105 ≈ n1.4. For a given data distribution and RF model we train two fitted functions, one on a finite dataset (empirical world) and the other on the data distribution (oracle world). More specifically, the training of the empirical model is done using multi-pass SGD on a finite training set of size n with learning rate η = 0.1 and batch size b = 50. The training of the oracle model is done using one-pass SGD with the same learning rate η and batch size b, but at each iteration a fresh batch is sampled from the population distribution. Both models f̂t, fort are initialized with ai = 0 for i ∈ [N ]. To speed up and stabilize optimization we use momentum β = 0.9. Note that if we took N → ∞, η → 0, and β = 0 we would be exactly in the dot product kernel gradient flow setting.
In Fig. 5, the data generating distributions of panels (5a), (5b), (5c) are respectively the same as that of panels (4a), (4b), and (4c) from Section 4.1. The top row of Fig. 5 shows SGD for {dot product, cyclic} RF models, and the bottom row shows the corresponding gradient flow for {dot product, cyclic} kernel least-squares. We see that the corresponding curves in these two rows exhibit qualitatively the same behaviors. Additionally, the results in panel (5c) show that as predicted, even for discrete SGD dynamics the dot product RF optimizes faster but fails to generalize, whereas the invariant RF optimizes slower but generalizes better.

heading: 5 SUMMARY AND DISCUSSION
text: In this paper, we used precise asymptotics to study the oracle world and empirical world dynamics of gradient flow on kernel least-squares objectives for high-dimensional regression problems. Under reasonable conditions on the target function and kernel, we showed that in this setting there are three learning stages based on the behaviors of the empirical and oracle models and also connected our results to some empirical deep learning phenomena.
Although our setting already captures some interesting aspects of deep learning training dynamics, there are some limitations which would be interesting to resolve in future work. We require very high-dimensional data in order for the asymptotics to be accurate, but real data distributions have low-dimensional structure. We work in a limiting regime of neural network training where the dynamics are linear and the step-size is infinitesimal. It is an important direction to extend this analysis to the non-linear feature learning regime and to consider discrete step-size minibatch SGD, as these are considered important aspects of network training. Our results hold for the square-loss in regression problems, but many deep learning problems involve classification using cross-entropy loss. Lastly, our analysis holds specifically for gradient flow, so it would be also interesting to consider other iterative learning algorithms such as boosting.

heading: A GENERAL SETTING
text: In this section we present our theory for training dynamics of kernel regression in an abstract setting similar to that of Mei et al. (2021a). We first introduce the setting of interest, then state the relevant assumptions, and finally we provide our theoretical results. We provide proofs of these results in Appendix B.

heading: A.1 PROBLEM SETUP
text: Consider a sequence of Polish probability spaces (Xd, νd), where νd is a probability measure on the configuration space Xd, indexed by an integer d. We denote by L2(Xd) = L2(Xd, νd) the space of square integrable functions on (Xd, νd). For p ≥ 1, we denote ‖f‖Lp(Xd) = Ex∼νd [|f(x)
p|]1/p the Lp norm of f . Let Dd ⊆ L2(Xd) be a closed linear subspace. In some simple applications we will consider Dd = L2(Xd), but the extra generality will be useful in certain applications. We are concerned with a supervised learning problem where we are given i.i.d data (yi,xi)i≤n. The feature vectors xi ∼iid νd are in Xd and the empirical-valued noisy responses yi are given by
yi = fd(xi) + εi
for some unknown target function fd ∈ Dd and εi ∼iid N (0, σ2ε). We consider a general RKHS defined on (Xd, νd) via the compact self-adjoint positive definite operator Hd : Dd → Dd which admits the representation
Hdg(x) = ∫ Xd Hd(x,x ′)g(x′)νd(dx ′),
where Hd ∈ L2(Xd ×Xd) with the property that ∫ Xd Hd(x,x ′)g(x′)νd(dx ′) = 0 for g ∈ D⊥d .
By the spectral theorem of compact operators, there exists an orthonormal basis (ψj)j≥1, span(ψj , j ≥ 1) = Dd ⊆ L2(Xd) and empirical eigenvalues (λd,j)j≥1 with nonincreasing absolute values |λd,1| ≥ |λd,2| ≥ · · · and ∑ j≥1 λ 2 d,j <∞ such that
Hd(x1,x2) = ∞∑ j=1 λ2d,jψj(x1)ψj(x2)
where convergence holds in L2(Xd ×Xd). For S ⊆ {1, 2, . . .} we denote PS to be the projection operator from L2(Xd) onto Dd,S := span(ψj , j ∈ S). We denote Hd,S to be the operator
Hd,S = ∞∑ j=1 λ2d,jψjψ ? j
and Hd,S the corresponding kernel.
If S = {j ∈ N : j ≤ `}we will write as short-hand Hd,≤` and analogously for S = {j ∈ N : j > `}. The trace of this operator is given by
Tr(Hd,S) ≡ ∑ j∈S λ2d,j = Ex∼νd [Hd,S(x,x)] <∞.
Define the test error R : L2(Xd)→ R and training error R̂n : L2(Xd)→ R
R(f) ≡ E(xnew,ynew){(ynew − f(xnew)) 2}, R̂n(f) ≡
1
n n∑ i=1 (yi − f(xi))2 (9)
where (x1, y1), . . . , (xn, yn), (xnew, ynew) are iid. For a kernelHd ∈ L2(Xd×Xd) we will consider the oracle model fort and the empirical model f̂t which satisfy the following gradient flows
d dt fort (x) = −∇R(fort (x)) = Ez∼νd [Hd(x, z)(fd(z)− fort (z))], (10)
d dt f̂t(x) = −∇R̂n(f̂t(x)) = 1 n n∑ i=1 Hd(x,xi)(yi − f̂t(xi)). (11)

heading: A.2 GENERAL ASSUMPTIONS
text: We now state our assumptions on the kernel and the sequence of probability spaces (Xd, νd). Assumption 1 ({n(d),m(d)}d≥1-Kernel Concentration Property). We say that the sequence of operators {Hd}d≥1 satisfies the Kernel Concentration Property (KCP) with respect to the sequence {n(d),m(d)}d≥1 if there exists a sequence of integers {r(d)}d≥1 with r(d) ≥ m(d), such that the following conditions hold.
(a) (Hypercontractivity of finite eigenspaces.) For any fixed q ≥ 1, there exists a constant C such that for any h ∈ Dd,≤r(d) = span(ψs, 1 ≤ s ≤ r(d)), we have
‖h‖L2q ≤ C‖h‖L2 .
(b) (Properly decaying eigenvalues) There exists fixed δ0 > 0, such that, for all d large enough, n(d)2+δ0 ≤ ( ∑∞ j=r(d)+1 λ 4 d,j)
2∑∞ j=r(d)+1 λ 8 d,j ,
n(d)2+δ0 ≤ ( ∑∞ j=r(d)+1 λ 2 d,j)
2∑∞ j=r(d)+1 λ 4 d,j .
(c) (Concentration of diagonal elements of kernel) For (xi)i∈[n(d)] ∼iid νd, we have:
max i∈[n(d)] ∣∣Ex∼νd [Hd,>m(d)(xi,x)2]− Ex,x′∼νd [Hd,>m(d)(x,x′)2]∣∣ = od,P(1) · Ex,x′∼νd [Hd,>m(d)(x,x′)2], max i∈[n(d)]
∣∣Hd,>m(d)(xi,xi)− Ex[Hd,>m(d)(x,x)]∣∣ = od,P(1) · Ex[Hd,>m(d)(x,x)]. Assumption 1(a) can be interpreted as requiring that the top eigenfunctions of Hd are delocalized. Assumption 1(b) concerns the tail of eigenvalues of Hd and is a mild assumption in high-dimensions. Lastly, 1(c) essentially requires that “most points” in Xd behave similarly in the sense of having similar values of the kernel diagonal Hd(x,x). Assumption 2 (Eigenvalue condition at level {(n(d),m(d))}d≥1). We say that the sequence of kernel operators {Hd}d≥1 satisfies the Eigenvalue Condition at level {(n(d),m(d))}d≥1 if the following conditions hold for all d large enough
(a) There exists a fixed δ0 > 0, such that
n(d)1+δ0 ≤ 1 λ4d,m(d)+1 ∞∑ k=m(d)+1 λ4d,k, (12)
n(d)1+δ0 ≤ 1 λ2d,m(d)+1 ∞∑ k=m(d)+1 λ2d,k. (13)
(b) There exists a fixed δ0 > 0, such that
n(d)1−δ0 ≥ 1 λ2d,m(d) ∞∑ k=m(d)+1 λ2d,k.
(c) There exists a fixed δ0 > 0, such that
m(d) ≤ n(d)1−δ0 .
Assumptions 2(a) and 2(b) can be seen as a spectral gap assumption. This ensures a clear separation between the eigenvalues in the subspace Dd,≤m(d) and the subspace Dd,>m(d). The technical requirement 2(c) is mild.
In the asymptotic setting, we will be interested in understanding the model learned at a time t = t(d) which scales with the dimension. The next set of assumptions give requirements for a valid scaling.
Assumption 3 (Admissible Time at {t(d), u(d), n(d),m(d)}d≥1). We say that the sequence {t(d), u(d), n(d),m(d)}d≥1 is an Admissible Time for the sequence of kernel operators {Hd}d≥1 if the following conditions hold
(a) There exists a fixed δ0 > 0, such that for d large enough 1
λ2d,u(d) ≤ t(d)1−δ0 ≤ t(d)1+δ0 ≤ 1 λ2d,u(d)+1 .
(b) If u(d) < m(d) for infinitely many d, then
t(d)
n(d) ∞∑ k=m(d)+1 λ2d,k = od(1).
(c) There exists a constant C such that m(d)∑ k=1 λ2d,k ≤ C ∞∑ k=m(d)+1 λ2d,k.
Assumption 3(a) is similar to the spectral gap condition Assumption 2(a), 2(b) for (t(d), u(d))d≥1. Assumption 3(b) relates the ordering of the indices u(d),m(d) to the relative growth of t(d), n(d). Assumption 3(c) requires that the eigenvalue tail does not decay too abruptly.

heading: A.3 MAIN RESULTS
text: In this section we give the main theoretical results. Recall the problem set-up and notation from Appendix A.1. We will characterize the gradient flow dynamics dynamics of the oracle model fort Eq. (10) and the empirical model f̂t Eq. (11) for a general kernel Hd. Theorem 3 (Oracle World). Let {fd ∈ Dd}d≥1 be a sequence of functions and {Hd}d≥1 be a sequence of kernel operators such that {(Hd, t(d), u(d))}d≥1 satisfies Assumption 3(a), then
R(fort ) = ∥∥P>u(d)fd∥∥2L2 + σ2ε + od(1) · ‖fd‖2L2 ,∥∥fort − P≤u(d)fd∥∥2L2 = od(1) · ‖fd‖2L2 ,
where P≤u(d) and P>u(d) are the projection operators onto the subspace spanned by the top u(d) kernel eigenfunctions and the orthogonal complement respectively, as defined in Appendix A.1.
The error of the oracle model is determined solely by optimization time t through u(d). Due to the spectral gap assumption 3(a) learning only occurs along the top u(d) eigenfunctions.
The next results describe the empirical model. First we characterize the training error. Theorem 4 (Empirical World - Train). Let {fd ∈ Dd}d≥1 be a sequence of functions, (xi)i∈[n(d)] ∼ νd independently, and {Hd}d≥1 be a sequence of kernel operators such that {(Hd, n(d),m(d), t(d), u(d))}d≥1 satisfies {(n(d),m(d))}d≥1-KPCP (Assumption 1), eigenvalue condition at level {(n(d),m(d))}d≥1 (Assumption 2), and {(n(d),m(d), t(d), u(d))}d≥1 is a valid time (Assumption 3). Define κH = Tr ( Hd,>m(d) ) and `(d) = min{u(d),m(d)}. Then for any η > 0 we have,
R̂n(f̂t) = ∥∥P>`(d)fd∥∥2L2 + σ2ε + od,P(1) · (‖fd‖2L2+η + σ2ε), if t = od(n/κH), R̂n(f̂t) = od,P(1) · (‖fd‖2L2+η + σ 2 ε), if t = ωd(n/κH).
In the early-time regime t n/κH the training error may be non-zero and matches the oracle world error if also u(d) ≤ m(d). In the late-time regime t n/κH the training error is negligible and the model interpolates the training set. The quantity κH arises since the empirical kernel matrix can be decomposed as H = H≤m +H>m and the second component is approximately a multiple of the identity: H>m ≈ Tr(H>m) · In = κH · In. This term acts as a self-induced ridge-regularizer. Our final result characterizes the test error of the empirical model.
Theorem 5 (Empirical World - Test). Let {fd ∈ Dd}d≥1 be a sequence of functions, (xi)i∈[n(d)] ∼ νd independently, and {Hd}d≥1 be a sequence of kernel operators such that {(Hd, n(d),m(d), t(d), u(d))}d≥1 satisfies {(n(d),m(d))}d≥1-KPCP (Assumption 1), eigenvalue condition at level {(n(d),m(d))}d≥1 (Assumption 2), and {(n(d),m(d), t(d), u(d))}d≥1 is a valid time (Assumption 3). Define `(d) = min{u(d),m(d)}. Then for any η > 0 we have,
R(f̂t) = ∥∥P>`(d)fd∥∥2L2 + σ2ε + od,P(1) · (‖fd‖2L2+η + σ2ε),∥∥∥f̂t − P≤`(d)fd∥∥∥2
L2 = od,P(1) · (‖fd‖2L2+η + σ 2 ε).
The result shows that the empirical model is essentially the projection of the regression function onto the first `(d) eigenfunctions. The quantity `(d) controls the complexity of the model which increases with time t up until u(d) ≥ m(d) after which the complexity is limited by n.

heading: B PROOF OF GENERAL SETTING

heading: B.1 ORACLE WORLD - PROOF OF THEOREM 3
text: The oracle model ODE Eq. (10) with initialization for0 ≡ 0 can be solved (c.f. Appendix E.1) to yield the solution fort = fd − e−tHdfd. Therefore the excess risk is given by
R(fort )−σ2ε = Ex∼νd(fd(x)−fort (x))2 = ∫ Xd fd(x)e −2tHdfd(x)νd(dx) = ∞∑ k=0 exp ( −2tλ2d,k ) f̂2k .
where λ2d,k are the kernel eigenvalues and f̂k := 〈fd, ψk〉L2 are the Fourier coefficients of fd in the kernel eigenbasis (c.f. Appendix A.1). We can control this quantity as follows,
|R(fort )− σ2ε − ‖P>ufd‖ 2 L2 |/ ‖fd‖ 2 L2 ≤ max { max k≤u exp ( −2tλ2d,k ) , max k≥u+1 1− exp ( −2tλ2d,k )} ≤ max { max k≤u exp ( −Ω(tλ2d,k) ) , max k≥u+1 O(tλ2d,k)
} ≤ max { exp ( −Ω(tδ0) ) , O(t−δ0) } = od(1)
where the last inequality follows from Assumption 3(a). This shows the first theorem statement,
‖fort − P≤ufd‖ 2 L2 = od(1) · ‖fd‖ 2 L2 .
Now observe that
P≤ufd − fort = u∑
k=0
e−tλ 2 d,k f̂kψk − ∑ k≥u+1 (1− e−tλ 2 d,k)f̂kψk
hence
‖fort − P≤ufd‖ 2 L2 / ‖fd‖ 2 L2 ≤ max { max k≤u exp ( −2tλ2d,k ) , max k≥u+1 (1− e−tλ 2 d,k)2 } ≤ max { max k≤u exp ( −2tλ2d,k ) , max k≥u+1 1− exp ( −2tλ2d,k )} = od(1)
where the second inequality follows from the fact that (1− e−x)2 ≤ 1− e−2x and the final equality is from the proof of the first part of the theorem. Thus,
‖fort − P≤ufd‖ 2 L2 = od(1) · ‖fd‖ 2 L2
completing the proof.

heading: B.2 EMPIRICAL WORLD – PRELIMINARIES
text: We will introduce some useful notations and the objects of analysis for studying the empirical world.

heading: B.2.1 TRAINING DYNAMICS
text: For the training of the empirical world c.f. Eq. (11), if u(t) = (f̂t(x1), . . . , f̂t(xn)) ∈ Rn then from Eq. (44) if u(0) = 0 then
u(t) = y − e−tH/ny = (In − e−tH/n)y, (14)
where H ∈ Rn×n with the (i, j)th element given by Hij = Hd(xi,xj) and y = f + ε ∈ Rn with f = (fd(x1), . . . , fd(xn)) T ∈ Rn and ε = (ε1, . . . , εn)T ∈ Rn.
For x ∈ Rd, define h(x) = (Hd(x1,x), . . . ,Hd(xn,x))T ∈ Rn. The training and test errors as defined in Eq. (1) can be written as
R̂n(f̂t) = 1
n ‖u(t)− y‖22 =
1 n yTe−2tH/ny,
R(f̂t) = Ex [ (fd(x)− u(t)TH−1h(x))2 ] .
Expanding R(f̂t) yields
R(f̂t) = Ex[fd(x)2]− 2u(t)TH−1E + u(t)TH−1MH−1u(t) (15)
where E = (E1, . . . , En)T ∈ Rn, M = (Mij)i,j∈[n] ∈ Rn×n, and H = (Hij)i,j∈[n] ∈ Rn×n with
Ei = Ex[fd(x)Hd(x,xi)], Mij = Ex[Hd(x,xi)Hd(x,xj)], Hij = Hd(xi,xj).

heading: B.2.2 DECOMPOSITIONS AND NOTATIONS
text: In this section we recall some useful decompositions of empirical quantities from Mei et al. (2021a). As mentioned earlier the eigendecomposition of Hd is given by
Hd(x,y) = ∞∑ k=1 λ2d,kψk(x)ψk(y).
We write the orthogonal decomposition of fd in the basis {ψk}k≥1 as
fd(x) = ∞∑ k=1 f̂d,kψk(x).
Define
ψk = (ψk(x1), . . . , ψk(xn)) T ∈ Rn,
D≤m = diag(λd,1, λd,2, . . . , λd,m) ∈ Rm×m, Ψ≤m = (ψk(xi))i∈[n],k∈[m] ∈ Rn×m,
f̂≤m = (f̂d,1, f̂d,2, . . . , f̂d,m) T ∈ Rm.
We have the following orthogonal basis decompositions of f ,H,E andM
f = f≤m + f>m, f≤m = Ψ≤mf̂≤m, f>m = ∞∑ k=m+1 f̂d,kψk,
H = H≤m +H>m, H≤m = Ψ≤mD 2 ≤mΨ T ≤m, H>m = ∞∑ k=m+1 λ2d,kψkψ T k ,
E = E≤m +E>m, E≤m = Ψ≤mD 2 ≤mf̂≤m, E>m = ∞∑ k=m+1 λ2d,kf̂d,kψk,
M = M≤m +M>m, M≤m = Ψ≤mD 4 ≤mΨ T ≤m, M>m = ∞∑ k=m+1 λ4d,kψkψ T k .
(16)
By Lemma 6 below, under Assumptions 1 and 2(a) the matricesH andM can be written as
H = Ψ≤mD 2 ≤mΨ T ≤m + κH(In + ∆H), (17)
M = Ψ≤mD 4 ≤mΨ T ≤m + κM (In + ∆M ), (18)
where
κH = Tr(Hd,>m) = ∞∑
k≥m+1
λ2d,k,
κM = Tr ( H2d,>m ) = ∞∑ k≥m+1 λ4d,k,
and max{‖∆H‖op , ‖∆M‖op} = od,P(1).
We will use α as shorthand for the scalar valued dimension dependent quantity e−(t/n)κH and take
K≤m := In − αe−(t/n)Ψ≤mD 2 ≤mΨ T ≤m . (19)
We also introduce the shrinkage matrix defined as
S≤m = ( Im +
κH n D−2≤m
)−1 = diag((sj)j∈[m]) ∈ Rm×m, where sj =
λ2d,j λ2d,j + κH n . (20)
If unspecified we will typically use ∆,∆′, etc. to denote matrices with operator norm od,P(1). For positive integers ` < m, define [`,m] = {`+ 1, . . . ,m}. We will use the following notation
S`m = diag((sj)j∈[`,m]) for sj defined in Eq. (20)
Ψ`m = (ψk(xi))i∈[n],k∈[`,m] ∈ Rn×(m−`) f`m = ∑
k∈[`,m]
f̂d,kψk

heading: B.2.3 AUXILIARY LEMMAS
text: Here we collect some lemmas which will be of use to us. Lemma 1 (Matrix Exponential Perturbation Inequality). For matrix operator norm ‖ · ‖, ifA,B ∈ Rn×n are symmetric then ∥∥eA − eB∥∥ ≤ ‖A−B‖max{∥∥eA∥∥,∥∥eB∥∥}. For general square matricesA,B ∈ Rn×n we have∥∥eA − eB∥∥ ≤ ‖A−B‖e‖A‖e‖B‖. Lemma 2. Let Z ∈ Rm×n,Ψ ∈ Rm×p,D ∈ Rp×p and t ∈ R. Denote A = ZTΨ ∈ Rn×p and B = ΨTΨ ∈ Rp×p. Then,
ZTetΨDΨ T Ψ = AetDB.
The notations in Lemmas 3-8 all follow the notations given in Section B.2.2. Lemma 3 (Lemma 12 from Mei et al. (2021a) with λ = 0). Let Assumptions 1 and 2 hold. Then,∥∥nH−1MH−1 −Ψ≤mS2≤mΨT≤m/n∥∥op = od,P(1). Lemma 4 (Theorem 6(b) from Mei et al. (2021a)). Let Assumptions 1(a), 2(c) hold. Then,∥∥ΨT≤mΨ≤m/n− Im∥∥op = od,P(1). Lemma 5 (Lemma 13 from Mei et al. (2021a) with λ = 0). Let Assumptions 1 and 2 hold. Then,∥∥ΨT≤mH−1Ψ≤mD2≤m − S≤m∥∥op = od,P(1). Lemma 6 (Theorem 6 from Mei et al. (2021a)). Let Assumptions 1 and 2(a) hold. Then we can decompose the kernel matrices as follows
H = Ψ≤mD 2 ≤mΨ T ≤m + κH(In + ∆H),
M = Ψ≤mD 4 ≤mΨ T ≤m + κM (In + ∆M ),
where
κH = Tr(Hd,>m) = ∞∑
k≥m+1
λ2d,k,
κM = Tr ( H2d,>m ) = ∞∑ k≥m+1 λ4d,k,
and max{‖∆H‖op , ‖∆M‖op} = od,P(1).
Lemma 7. Let Assumption 1(a) hold. Let S, T be disjoint subsets of N. Let D ∈ R|S|×|S| be a diagonal matrix. Then we have that for any η > 0, there exists C(η) (independent of d), such that
E[f̂TTΨTTΨSDΨTSΨT f̂T ]/n2 ≤ C(η) ‖PT fd‖ 2 L2+η Tr(D)/n,
where the expectation is with respect to the randomness in ΨS ,ΨT .
Proof. Let ι : S → [|S|] be the bijection such that ΨS = (ψι−1(k)(xi))i∈[n],k∈[|S|] then we have E[f̂TTΨTTΨSDΨTSΨT f̂T ]/n2 = ∑ u,v∈T ∑ s∈S ∑ i,j∈[n] Dι(s)ι(s)E[ψu(xi)ψs(xi)ψs(xj)ψv(xj)]f̂uf̂v/n2
= ∑ u,v∈T ∑ s∈S ∑ i∈[n] Dι(s)ι(s)E[ψu(xi)ψs(xi)ψs(xi)ψv(xi)]f̂uf̂v/n2
= 1
n ∑ s∈S Dι(s)ι(s)Ex[(PT fd(x))2ψs(x)2]
≤ 1 n ∑ s∈S Dι(s)ι(s) ‖PT fd‖ 2 L2+η ‖ψs‖ 2 L(4+2η)/η
≤ C(η) ‖PT fd‖2L2+η n∑ i=1 Dii/n,
where the second to last inequality is by Holder’s inequality and the last inequality used the hypercontractivity assumption as in Assumption 1(a).
Lemma 8 (Exponential Kernel Decomposition). Assume the conditions of Lemma 6 hold and assume there exists δ0 > 0, such that (t(d), u(d)) satisfies the condition
t(d)1+δ0 ≤ 1 λ2d,u(d)+1 . (21)
Let j(d) satisfy min{u(d),m(d)} ≤ j(d) ≤ m(d) then∥∥∥e−tHd/n − e−(t/n)(Ψ≤jD2≤jΨT≤j+κHIn)∥∥∥ op = od,P(1).
Proof. First consider the regime t = Od(n/κH). Recall the decomposition,
H = Ψ≤mD 2 ≤mΨ T ≤m + κH(In + ∆H)
where ‖∆H‖op = od,P(1). By Lemma 1 and Lemma 4,∥∥∥e−tH/n − e−(t/n)(Ψ≤jD2≤jΨT≤j+κHI)∥∥∥ op ≤ ∥∥∥∥∥∥ m∑
k≥j+1
tλ2d,kψkψ T k /n ∥∥∥∥∥∥ op + (t/n)κH ‖∆h‖op
≤ (t · max k≥j+1 ·λ2d,k) ∥∥ΨT≤mΨ≤m/n∥∥op + od,P(1)
≤ Od,P(t · max k≥j+1 λ2d,k) + od,P(1)
(a) = Od,P(t −δ0) + od,P(1) = od,P(1),
where equality (a) holds by assumption Eq. (21). Now if t = ωd(n/κH), then it is easy to see that
max {∥∥∥e−tH/n∥∥∥ op , ∥∥∥e−(t/n)(Ψ≤jD2≤jΨT≤j+κHI)∥∥∥ op } ≤ e−ωd,P(1) = od,P(1),
and therefore ∥∥∥e−tH/n − e−(t/n)(Ψ≤jD2≤jΨT≤j+κHI)∥∥∥ op
≤ 2 max {∥∥∥e−tH/n∥∥∥ op , ∥∥∥e−(t/n)(Ψ≤jD2≤jΨT≤j+κHI)∥∥∥ op } = od,P(1).

heading: B.3 EMPIRICAL WORLD - TRAIN
text: If t = ωd(n/κH), then by Eq. (17) it is easy to see that ∥∥e−2(t/n)H∥∥
op = od,P(1), hence
R̂n(f̂t) = od,P(1) · ‖fd‖2L2 .
From now on we focus on the case that t = od(n/κH). Let us decompose the training error as follows
R̂n(f̂t) = R1 +R2 +R3
where
R1 = 1
n fTe−2(t/n)Hf ,
R2 = 2
n εTe−2(t/n)Hf ,
R3 = 1
n εTe−2(t/n)Hε.
B.3.1 TERM R1
Let us start by analysing R1. We can write
R1 = 1
n (f≤` + f>`)
Te−2(t/n)H(f≤` + f>`)
= T1 + 2T2 + T3
where
T1 = 1
n fT≤`e −2(t/n)Hf≤`,
T2 = 1
n fT>`e −2(t/n)Hf≤`,
T3 = 1
n fT>`e −2(t/n)Hf>`.
Let us analyse the term T1,
T1 ≤ 1
n fT≤`e −2(t/n)Ψ≤`D2≤`Ψ T ≤`f≤`.
By Lemma 2, denotingB = ΨT≤`Ψ≤`/n,
1 n ΨT≤`e −2(t/n)Ψ≤`D2≤`Ψ T ≤`Ψ≤` = Be −2tD2≤`B.
We will show that ∥∥∥Be−2tD2≤`B∥∥∥ op = od,P(1). (22)
By Lemma 4, since ‖B‖op = Θd,P(1), for d largeB has a positive-definite square rootB1/2 w.h.p. Therefore we can write
Be−2tD 2 ≤`B = B1/2e−2tB 1/2D2≤`B 1/2 B1/2
and bound the operator norm∥∥∥Be−2tD2≤`B∥∥∥ op ≤ ‖B‖op ∥∥∥e−2tB1/2D2≤`B1/2∥∥∥ op
= ‖B‖op e −2tλmin(B1/2D2≤`B 1/2)
≤ ‖B‖op e −2tλmax(B)λmin(D2≤`) = od,P(1).
since by Assumption 3(a), tλmin(D2≤`) = Ω(t δ0). Therefore T1 = od,P(1) · ‖P≤`fd‖2L2 .
Now we analyse the term T3. Observe that by the inequality 1− x ≤ e−x ≤ 1, 1
n ‖f>`‖22 − (2t/n)
1 n fT>`(Ψ≤`D 2 ≤`Ψ T ≤` + κH(In + ∆H))f>` ≤ T3 ≤ 1 n ‖f>`‖22.
We have by Lemma 7,
tE[f̂T>`ΨT>`Ψ≤`D2≤`ΨT≤`Ψ>`f̂>`]/n2 ≤ C(η) ‖P>`fd‖ 2 L2+η
( t
n ∑̀ s=1 λ2d,s ) where the last quantity is od,P(1) · ‖P>`fd‖2L2+η by Assumption 3(c). Thus we see that
T3 = ‖P>`fd‖2L2 + od,P(1) · ‖P>`fd‖ 2 L2+η .
Observe that by the Cauchy-Schwarz inequality,
T2 ≤ (T1T3)1/2 ≤ od,P(1) ‖P≤`fd‖L2 ‖P>`fd‖L2+η . Putting everything together we see that,
R1 = ‖P>`fd‖2L2 + od,P(1) · (‖fd‖ 2 L2 + ‖P>`fd‖ 2 L2+η ).
B.3.2 TERM R2
Turning to R2, we take the second-moment with respect to ε 1
σ2ε Eε[R22] =
4
n2 Eε[εTe−2(t/n)HffTe−2(t/n)Hε]/σ2ε
= 4 n2 fTe−4(t/n)Hf ≤ 4 n (‖f‖22/n) = od,P(1) · ‖fd‖2L2 .
By Markov’s inequality R2 = od,P(1) · ‖fd‖2L2 σ2ε .
B.3.3 TERM R3
Now let us analyse R3. Recalling the definitionB = ΨT≤`Ψ≤`/n, we can compute the expectation
1
σ2ε (1− Eε[R3]) =
1 n Tr ( In − e−2(t/n)H ) (a) = 1
n e−κHt/n Tr
( In − e−2(t/n)(Ψ≤`D 2 ≤`Ψ T ≤`) ) + od,P(1)
(b) ≤ 1 n e−κHt/n Tr
( 2(t/n)(Ψ≤`D 2 ≤`Ψ T ≤`) ) + od,P(1)
(c) =
1 n e−κHt/n Tr
( 2tBD2≤` ) + od,P(1)
≤ 2t n e−κHt/n ‖B‖op Tr
( D2≤` ) + od,P(1) = od,P(1).
Equality (a) follows from Lemma 8. For (b) we used the inequality 1−e−x ≤ x and that trace is the sum of eigenvalues. Equality (c) uses the cyclic property of trace. In the last equality we used that by Lemma 4, ‖B‖op = Od,P(1) and Assumption 3(c) implies that (t/n) Tr ( D2≤` ) = Od((t/n)κH) = od(1) since by assumption we consider t = od(n/κH). Turning to the variance
Varε[R3] = Eε[R23]− Eε[R3]2
= 1
n2 Eε[(εTe−2(t/n)Hε)2]− Eε[R3]2
≤ 1 n2 Eε[(εTε)2]− σ4ε(1 + od,P(1))
(a) = O(1/n) · σ4ε + σ4ε − σ4ε(1 + od,P(1)) = od,P(1) · σ4ε ,
where the first inequality uses that ∥∥e−2(t/n)H∥∥
op ≤ 1 and (a) uses that for εi ∼ N (0, σ2ε), E[ε4i ] =
3σ4ε . Therefore by Chebyshev’s inequality, R3 = σ 2 ε(1 + od,P(1)).
Putting everything together yields
R̂n(f̂t) = R1 +R2 +R3 = ‖P>`‖2L2 + od,P(1) · (‖fd‖ 2 L2+η + σ 2 ε).

heading: B.4 EMPIRICAL WORLD - TEST
text: Recalling u(t) from Eq. (14), let u(t) = v(t) + ε(t)
where
v(t) = (In − e−tH/n)f , (23) ε(t) = (In − e−tH/n)ε. (24)
Recall the expansion of the test error from Eq. (15),
R(f̂t) = Ex[fd(x)2]− 2u(t)TH−1E + u(t)TH−1MH−1u(t) = ‖fd‖2L2 − 2T1 + T2 + T3 − 2T4 + 2T5,
where
T1 = v(t) TH−1E,
T2 = v(t) TH−1MH−1v(t),
T3 = ε(t) TH−1MH−1ε(t),
T4 = ε(t) TH−1E,
T5 = ε(t) TH−1MH−1v(t).
The proof for the test error is the most involved, but will follow a similar strategy of analysing each term in the expansion. We first begin by analysing T2 in Section B.4.1, then T1 in Section B.4.2, and finally terms T3, T4, T5, which are all simpler than the first two, in Section B.4.3. At the end of this Appendix section we present the proof of Theorem 5.
B.4.1 TERM T2
As before, we will analyze each term separately. We begin with term T2. Proposition 1 (Term T2).
T2 = v(t) TH−1MH−1v(t) = ∥∥∥S≤`f̂≤`∥∥∥2 2 + od,P(1) · ‖fd‖2L2+η
where we recall the shrinkage matrix S≤m defined in Eq. (20) and v(t) in Eq. (23).
Proof of Proposition 1. To analyze T2 we further decompose it into the following terms
T2 = (f≤` + f>`) T(In − e−tH/n)H−1MH−1(In − e−tH/n)(f≤` + f>`)
= T21 + T22 + T23,
where
T21 = f T ≤`(In − e−tH/n)H−1MH−1(In − e−tH/n)f≤`, (25)
T22 = 2f T ≤`(In − e−tH/n)H−1MH−1(In − e−tH/n)f>`, (26)
T23 = f T >`(In − e−tH/n)H−1MH−1(In − e−tH/n)f>`. (27)
Using Lemma 9 and 10 proven below, by the Cauchy-Schwarz inequality,
T22 ≤ (2T21T23)1/2 = od,P(1) · ‖P≤`fd‖L2 ‖fd‖L2+η and hence
T2 = T21 + T22 + T23 = ∥∥∥S≤`f̂≤`∥∥∥2
2 + od,P(1) · ‖fd‖2L2+η .
Lemma 9 (Term T21 Eq. (25)). T21 = ∥∥∥S≤`f̂≤`∥∥∥2
2 + od,P(1) · ‖P≤`fd‖2L2
Proof of Lemma 9. Recall the notation α andK≤` from Eq. (19). Define ∆ such that
∆ := (In − e−tH/n)−K≤` which by Lemma 8 satisfies ‖∆‖op = od,P(1). Then we can split T21 into
T21 = T211 + T212 + T213 + T214
where
T211 = f T ≤`K≤`Ψ≤`S 2 ≤`Ψ T ≤`K≤`f≤`/n 2,
T212 = f T ≤`K≤`Ψ`mS 2 `mΨ T `mK≤`f≤`/n 2,
T213 = f T ≤`K≤`(H −1MH−1 −Ψ≤mS2≤mΨT≤m/n2)K≤`f≤`, T214 = f T ≤`∆H −1MH−1K≤`f≤` + f T ≤`K≤`H −1MH−1∆f≤` + f T ≤`∆H −1MH−1∆f≤`. We will show that the dominant term is T211 and the others are of lower order. By Lemma 3,∥∥nH−1MH−1 −Ψ≤mS2≤mΨT≤m/n∥∥op = od,P(1), and sinceK≤` In, T213 = od,P(1) · ‖P≤`fd‖2L2 . By Lemma 3, Lemma 4 and since S≤m Im,∥∥nH−1MH−1∥∥ op ≤ ∥∥Ψ≤mS2≤mΨT≤m/n∥∥op + od,P(1) ≤ 1 + od,P(1)
hence it is easy to see that T214 = od,P(1) · ‖P≤`fd‖2L2 . Turning to T212 we have
T212 ≤ 1
n2 fT≤`(In − αe−(t/n)Ψ≤`D 2 ≤`Ψ T ≤`)Ψ`mΨ T `m(In − αe−(t/n)Ψ≤`D 2 ≤`Ψ T ≤`)f≤`
(a) =
1
n2 f̂T≤`(I− αe−(t/n)BD 2 ≤`)ΨT≤`Ψ`mΨ T `mΨ≤`(I− αe−(t/n)BD 2 ≤`)f̂T≤`
≤ ∥∥∥f̂≤`∥∥∥2
2 ∥∥ΨT≤`Ψ`m/n∥∥2op = od,P(1) · ‖P≤`fd‖2L2 , where the first inequality used S2`m Im−`, we used Lemma 2 in (a), and the last equality used Lemma 4 (note that ΨT≤`Ψ`m/n is an off-diagonal block of Ψ T ≤mΨ≤m/n which corresponds to an all zeroes submatrix of Im). Finally we look at the main term T211. Defining the matrices
A := 1
n ΨT≤`e −(t/n)Ψ≤`D2≤`Ψ T ≤`Ψ≤`,
B := 1
n ΨT≤`Ψ≤`,
we can write
T211 = 1
n2 f̂T≤`Ψ T ≤`(In − αe−(t/n)Ψ≤`D 2 ≤`Ψ T ≤`)Ψ≤`S 2 ≤`Ψ T ≤`(In − αe−(t/n)Ψ≤`D 2 ≤`Ψ T ≤`)Ψ≤`f̂≤`
= f̂T≤`(BS 2 ≤`B − 2αAS2≤`B + α2AS2≤`A)f̂≤`.
Now observe that by Lemma 2,
A = 1
n ΨT≤`e −(t/n)Ψ≤`D2≤`Ψ T ≤`Ψ≤` = Be −tD2≤`B
and by the same argument used to show Eq. (22), ‖A‖op = od,P(1). Since B = In + ∆′ and S≤` I`, we have
T211 = ∥∥∥S≤`f̂≤`∥∥∥2
2 + od,P(1) · ‖P≤`fd‖2L2 ,
hence combining terms T21 = T211 + T212 + T213 + T214 = ∥∥∥S≤`f̂≤`∥∥∥2
2 + od,P(1) · ‖P≤`fd‖2L2 .
Lemma 10 (Term T23 Eq. (27)).
T23 = f T >`(In − e−tH/n)H−1MH−1(In − e−tH/n)f>` = od,P(1) · ‖fd‖ 2 L2+η .
Proof of Lemma 10. Let us define the matrix
G = (In − e−tH/n)H−1MH−1(In − e−tH/n). Using similar reasoning from Lemma 9 we can write
G = 1
n2 K≤mΨ≤mS
2 ≤mΨ T ≤mK≤m + ∆
= 1
n2 K≤`Ψ≤mS
2 ≤mΨ T ≤mK≤` + ∆ ′.
for matrices ∆,∆′ satisfying max{‖∆‖op , ‖∆′‖op} = od,P(1). We can split T23 into
T23 = 1
n fT>`Gf>`
= 1
n (f>m + f`m)
TG(f>m + f`m)
= 1
n fT>mGf>m +
2 n fT`mGf>m + 1 n fT`mGf`m
= T231 + T232 + T233 + T234
where
T231 = 1
n2 fT>mK≤mΨ≤mS 2 ≤mΨ T ≤mK≤mf>m,
T232 = 2
n2 fT`mGf>m,
T233 = 1
n2 fT`mK≤`Ψ≤mS 2 ≤mΨ T ≤mK≤`f`m,
T234 = od,P(1) · ‖P>`fd‖2L2 .
Let us first start with analysing T231, definingB = ΨT≤mΨ≤m/n we have
T231 = 1
n2 fT>mΨ≤m(Im − αe−(t/n)D 2 ≤mB)S2≤m(Im − αe−(t/n)BD 2 ≤m)ΨT≤mf>m
(a) ≤ (1 + od,P(1))fT>mΨ≤mΨT≤mf>m/n2
= od,P(1) · ‖P>mfd‖2L2+η ,
where the first equality is by Lemma 2 and the last line follows from Lemma 7, Markov’s inequality, and by the fact that m/n = od(1) by Assumption 2(c). To see that inequality (a) holds note that∥∥∥Im − αe−(t/n)D2≤mB∥∥∥ op = ∥∥∥B−1/2(Im − αe−(t/n)B1/2D2≤mB1/2)B1/2∥∥∥ op
≤ ∥∥∥B−1/2∥∥∥
op ∥∥∥B1/2∥∥∥ op .
Since S2≤m Im,∥∥∥(Im − αe−(t/n)D2≤mB)S2≤m(Im − αe−(t/n)BD2≤m)∥∥∥ op ≤ ∥∥B−1∥∥ op ‖B‖op = 1 + od,P(1),
where the last equality is by Lemma 4. Now let us turn to T233, which we can further split into
T233 = T2331 + T2332 (28)
where
T2331 = 1
n2 fT`mK≤`Ψ≤`S 2 ≤`Ψ T ≤`K≤`f`m,
T2332 = 1
n2 fT`mK≤`Ψ`mS 2 `mΨ T `mK≤`f`m.
RedefiningB = ΨT≤`Ψ≤`/n and using a similar argument as for T231, the first term can be seen as
T2331 = 1
n2 f̂T`mΨ T `mΨ≤`(I− αe−(t/n)D 2 ≤`B)S2≤`(I− αe−(t/n)BD 2 ≤`)ΨT≤`Ψ`mf̂`m
≤ (1 + od,P(1)) ∥∥∥f̂`m∥∥∥2
2 ∥∥ΨT`mΨ≤`/n∥∥2op = od,P(1) · ‖P`mfd‖2L2 . (29) Turning to the second term T2332, let
A := 1
n ΨT`m
( In − e−(t/n)(Ψ≤`D 2 ≤`Ψ T ≤`+κHIn) ) Ψ`m,
then we can write this term as
T2332 = f̂ T `mAS 2 `mAf̂`m
≤ f̂T`mA 2 f̂`m ≤ f̂T`mAf̂`m + od,P(1) · ‖P`mfd‖ 2 L2 , (30)
where the last inequality holds sinceA ΨT`mΨ`m/n = Im−` + ∆ by Lemma 4 and so
A 2 = A 1/2 AA 1/2
A1/2(I + ∆)A1/2
= A+A 1/2 ∆A 1/2
= A+ ∆′.
By the inequality 1− e−x ≤ x, in the PSD order we see that
A t n2 ΨT`mΨ≤`D 2 ≤`Ψ T ≤`Ψ`m + κHt n
( 1
n ΨT`mΨ`m
) .
Therefore from Eq. (30) we have
T2332 ≤ t
n2 f̂T`mΨ T `mΨ≤`D 2 ≤`Ψ T ≤`Ψ`mf̂`m + f̂ T `m
κHt
n
( 1
n ΨT`mΨ`m
) f̂`m + od,P(1) · ‖P`mfd‖2L2 .
For the second term on the right, by Assumption 3(b) and by Lemma 4,
f̂T`m κHt
n
( 1
n ΨT`mΨ`m
) f̂`m = od,P(1) · ‖P`mfd‖2L2 .
For the first term by Lemma 7 and Assumptions 3(b), 3(c), t
n2 E[f̂T`mΨT`mΨ≤`D2≤`ΨT≤`Ψ`mf̂`m] ≤ (t/n)C(η) ‖P`mfd‖ 2 L2+η Tr
( D2≤` ) = od,P(1) · ‖P`mfd‖2L2+η ,
therefore by Markov’s inequality
T2332 = od,P(1) · ‖P`mfd‖2L2+η . Hence combining terms
T233 = T2331 + T2332 = od,P(1) · ‖P`mfd‖2L2+η . By the Cauchy-Schwarz inequality
T232 = 2
n2 fT`mGf>m
≤ 2 ( 1
n2 fT`mGf`m
1
n2 fT>mGf>m )1/2 = 2 ( T231 + od,P(1) · ‖P>`fd‖2L2 )1/2( T233 + od,P(1) · ‖P>`fd‖2L2
)1/2 ≤ od,P(1) · ‖P>`fd‖L2+η ‖P>mfd‖L2+η .
Putting everything together we see that
T23 = T231 + T232 + T233 + T234 = od,P(1) · ‖fd‖2L2+η .
B.4.2 TERM T1
Now we will analyse term T1. Proposition 2 (Term T1).
T1 = f T(In − e−tH/n)H−1E = ∥∥∥S1/2≤` f̂≤`∥∥∥2 2 + od,P(1) · ‖fd‖2L2+η .
Proof of Proposition 2. We break T1 into the following terms
T1 = T11 + T12 + T13
where
T11 = f T ≤`(I− e−tH/n)H−1E≤m, (31)
T12 = f T >`(I− e−tH/n)H−1E≤m, (32) T13 = f T(I− e−tH/n)H−1E>m. (33)
Using Lemma 11 and recalling Lemma 10, T12 ≤ (T23)1/2 ∥∥∥f̂≤m∥∥∥
2 = od,P(1) · ‖P≤mfd‖L2 ‖fd‖L2+η ,
where T23 is as given in Eq. (27). From the analysis of T11 in Lemma 12 and T13 in Lemma 13 we combine everything to get Proposition 2
T1 = ∥∥∥S1/2≤` f̂≤`∥∥∥2
2 + od,P(1) · ‖fd‖2L2+η .
Lemma 11. For a vector v ∈ Rn, vTH−1E≤m ≤ (vTH−1MH−1v)1/2 ∥∥∥f̂≤m∥∥∥
2
=
( 1
n2 vTΨ≤mS 2 ≤mΨ T ≤mv + od,P(1) ·
1 n ‖v‖2 )1/2∥∥∥f̂≤m∥∥∥ 2 .
Proof of Lemma 11. Recall that E≤m = Ψ≤mD2≤mf̂≤m. By the Cauchy-Schwarz inequality,
vTH−1E≤m = v TH−1Ψ≤mD 2 ≤mf̂≤m, ≤ (vTH−1Ψ≤mD4≤mΨT≤mH−1v)1/2 ∥∥∥f̂≤m∥∥∥
2 , ≤ (vTH−1MH−1v)1/2 ∥∥∥f̂≤m∥∥∥
2
=
( 1
n2 vTΨ≤mS 2 ≤mΨ T ≤mv + od,P(1) ·
1 n ‖v‖2 )1/2∥∥∥f̂≤m∥∥∥ 2 ,
where the last equality follows from Lemma 3.
Lemma 12 (Term T11 Eq. (31)). T11 = ∥∥∥S1/2≤` f̂≤`∥∥∥2
2 + od,P(1) · ‖P≤mfd‖2L2 .
Proof of Lemma 12. First note that by Lemma 8 for some ∆ such that ‖∆‖op = od,P(1), T11 = T111 + T112
where
T111 = f T ≤`(I− e−(t/n)(Ψ≤`D 2 ≤`Ψ T ≤`+κHIn))H−1E≤m,
T112 = f T ≤`∆H −1E≤m.
By Lemma 11, T112 ≤ ( ‖∆‖2op (‖f≤`‖ 2 2/n) ∥∥nH−1MH−1∥∥ op )1/2∥∥∥f̂≤m∥∥∥ 2 = od,P(1) · ‖P≤mfd‖2L2 .
We now consider T111 = T1111 − T1112
where
T1111 = f T ≤`H −1E≤m
T1112 = αf T ≤`e −(t/n)Ψ≤`D2≤`Ψ T ≤`H−1E≤m.
Note that by Lemma 5, ∥∥ΨT≤`H−1Ψ≤mD2≤m − [S≤`; 0]∥∥op = od,P(1), where 0 is a `× (m− `) matrix of zeros. Hence for the first term T1111,
T1111 = f̂ T ≤`Ψ T ≤`H −1Ψ≤mD 2 ≤mf̂≤m = ∥∥∥S1/2≤` f̂≤`∥∥∥2 2 + od,P(1) · ‖P≤mfd‖2L2 .
DefineH≤` := Ψ≤`D2≤`Ψ T ≤`. For the second term T1112, by Lemma 11 T1112 ≤ ( S + od,P(1) · ‖P≤`fd‖2L2 )1/2∥∥∥f̂≤m∥∥∥ 2 , (34)
where S = 1
n2 fT≤`e −(t/n)H≤`Ψ≤mS 2 ≤mΨ T ≤me −(t/n)H≤`f≤`.
DefineA := 1nΨ T ≤`e −(t/n)H≤`Ψ≤`. We have
S ≤ 1 n2 fT≤`e −(t/n)H≤`Ψ≤mΨ T ≤me −(t/n)H≤`f≤`
= 1
n2 fT≤`e −(t/n)H≤`Ψ≤`Ψ T ≤`e −(t/n)H≤`f≤` +
1
n2 fT≤`e −(t/n)H≤`Ψ`mΨ T `me −(t/n)H≤`f≤`
≤ 1 n2 fT≤`e −(t/n)H≤`Ψ≤`Ψ T ≤`e −(t/n)H≤`f≤` + ∥∥∥f̂≤`∥∥∥2 2 · ∥∥ΨT≤`Ψ`m/n∥∥2op = f̂T≤`A 2f̂≤` + od,P(1) · ‖P≤`fd‖2L2 = od,P(1) · ‖P≤`fd‖2L2 ,
since as noted before in Eq. (22), ‖A‖op = od,P(1). Therefore
T11 = T1111 + T1112 + T112 = ∥∥∥S1/2≤` f̂≤`∥∥∥2
2 + od,P(1) · ‖P≤mfd‖2L2 .
Lemma 13 (Term T13 Eq. (33)).
T13 = od,P(1) · ‖P>mfd‖L2 ‖fd‖L2 .
Proof of Lemma 13. We have
|T13| = |fT(I− e−tH/n)H−1E>m| ≤ ‖f‖2 ∥∥H−1∥∥ op ‖E>m‖2
Note that we have E[‖f‖22] = n ‖fd‖ 2 L2 . Further by Eq. (17), we have ∥∥H−1∥∥ op ≤ 2/κH with high probability. Finally, recalling the definition of E>m from Eq. (16), we have
E[‖E>m‖2] = n ∞∑
k=m+1
λ4d,kf̂ 2 k ≤ n [ max k≥m+1 λ4d,k ] ‖P>mfd‖2L2 .
As a result, we have
|T13| ≤ Od,P(1) · ‖P>mfd‖L2 ‖fd‖L2 [n 2 max k≥m+1 λ4d,k] 1/2/κH
= Od,P(1) · ‖P>mfd‖L2 ‖fd‖L2 [n max k≥m+1
λ2d,k] / ∑
k≥m+1
λ2d,k
= od,P(1) · ‖P>mfd‖L2 ‖fd‖L2 ,
where the last equality used Eq. (13) in Assumption 2(a).
B.4.3 TERMS T3, T4, T5
To analyse the terms T3, T4, T5 we can adapt the corresponding steps for the proof of Theorem 4 in Mei et al. (2021a). For the following analysis we recall the definition of ε(t) from Eq. (24).
Lemma 14 (Term T3).
T3 = ε(t) TH−1MH−1ε(t) = od,P(1) · σ2ε
Proof of Lemma 14.
1
σ2ε Eε[T3] = Tr
( (In − e−tH/n)2H−1MH−1 ) ≤ Tr ( H−1MH−1
) (a) = Tr ( Ψ≤mS 2 ≤mΨ T ≤m/n 2 ) + od,P(1)
(b) ≤ 1 n2
Tr ( Ψ≤mΨ T ≤m ) + od,P(1)
(c) =
1
n2 nm(1 + od,P(1)) + od,P(1) = od,P(1),
where (a) used Lemma 3, (b) used S≤m Im, and (c) used Lemma 4 and Assumption 2(c). The lemma then follows from Markov’s inequality.
Lemma 15 (Term T4).
T4 = ε(t) TH−1E = od,P(1) · (σ2ε + ‖fd‖ 2 L2).
Proof of Lemma 15. 1
σ2ε Eε[T 24 ] =
1
σ2ε Eε[εT(I− e−tH/n)H−1EETH−1(I− e−tH/n)ε]
= ETH−1(I− e−tH/n)2H−1E ≤ ETH−2E
Notice thatM Ψ≤LD4≤LΨT≤L for any L ∈ N, by the decomposition of Eq. (16). Therefore
sup L ∥∥D2≤LΨT≤LH−2Ψ≤LD2≤L∥∥op = sup L ∥∥H−1Ψ≤LD4≤LΨT≤LH−1∥∥op ≤ ∥∥H−1MH−1∥∥
op = od,P(1),
(35)
where the last inequality follows from Lemma 3. Hence,
ETH−2E = lim L→∞ ET≤LH −2ET≤L
(a) = lim
L→∞ f̂T≤L[D 2 ≤LΨ T ≤LH −2Ψ≤LD 2 ≤L]f̂≤L
(b) ≤ lim sup L→∞ ∥∥D2≤LΨT≤LH−2Ψ≤LD2≤L∥∥op · limL→∞∥∥∥f̂≤L∥∥∥22 (c)
≤ od,P(1) · ‖fd‖2L2 , where (a) follows from the definition ofE≤L, (b) follows from the definition of operator norm, and (c) follows from Eq. (35). Therefore we get
T4 = od,P(1) · σε · ‖fd‖L2 = od,P(1) · (σ 2 ε + ‖fd‖ 2 L2).
Lemma 16 (Term T5). T5 = ε(t) TH−1MH−1v(t) = od,P(1) · (σ2ε + ‖fd‖ 2 L2+η ).
Proof of Lemma 16. We can write term T5 as T5 = T51 + T52
where T51 = ε(t)
TH−1MH−1(In − e−tH/n)f≤`, T52 = ε(t)
TH−1MH−1(In − e−tH/n)f>`. Note as in Eq. (35), that by Lemma 3 and Lemma 4,∥∥∥M1/2H−2M1/2∥∥∥ op = ∥∥H−1MH−1∥∥ op = od,P(1),
and taking the second moment of T51 yields 1
σ2ε Eε[T 251] ≤
1
σ2ε Eε[εTH−1MH−1(In − e−tH/n)f≤`fT≤`(In − e−tH/n)H−1MH−1ε]
= fT≤`(In − e−tH/n)[H−1MH−1]2(In − e−tH/n)f≤` ≤ ∥∥∥M1/2H−2M1/2∥∥∥
op ∥∥∥M1/2H−1(In − e−tH/n)f≤`∥∥∥2 2
= od,P(1) · T21 = od,P(1) · ‖P≤`fd‖2L2 ,
where T21 is as given in Eq. (25). Similarly we get that 1
σ2ε Eε[T 252] = od,P(1) · T23 = od,P(1) · ‖fd‖ 2 L2+η ,
where T23 is as given in Eq. (27). By Markov’s inequality we deduce that
T5 = od,P(1) · σε · (‖P≤`fd‖L2 + ‖fd‖L2+η ) = od,P(1) · (σ 2 ε + ‖fd‖ 2 L2+η ).
Finally putting Propositions 1, 2 and Lemmas 14, 15, 16 together for terms T2, T1, T3, T4 and T5 respectively leads to the proof of Theorem 5
Proof of Theorem 5.
R(f̂t) = ‖fd‖2L2 − 2T1 + T2 + T3 − 2T4 + 2T5 = ‖P>`fd‖2L2 + ∥∥∥f̂`∥∥∥2 2 − 2 ∥∥∥S≤`f̂≤`∥∥∥2 + ∥∥∥S≤`f̂≤`∥∥∥2
+ od,P(1) · (‖fd‖2L2 + ‖fd‖ 2 L2+η + σ 2 ε) = ∥∥∥(I− S≤`)f̂≤`∥∥∥2
2 + ‖P>`fd‖2L2 + od,P(1) · (‖fd‖ 2 L2 + ‖fd‖ 2 L2+η + σ 2 ε).
By Assumption 2(b), κH/n = od(1) ·maxj≤` λ2d,j hence∥∥∥(I− S≤`)f̂≤`∥∥∥2 2
= od,P(1) · ‖fd‖2L2 and as a result we obtain the first part of the theorem
R(f̂t) = ‖P>`fd‖2L2 + od,P(1) · (‖fd‖ 2 L2 + ‖fd‖ 2 L2+η + σ 2 ε). (36)
Now observe that similar to Eq. (15) we have the following decomposition∥∥∥f̂t − P≤`fd∥∥∥2 L2 = ‖P≤`fd‖2L2 − 2u(t) TH−1E≤` + u(t) TH−1MH−1u(t), where u(t) is given in Eq. (14). Therefore we can write∥∥∥f̂t − P≤`fd∥∥∥2 L2 = R(f̂t)− ‖P>`fd‖2L2 + 2u(t) TH−1E>`. (37) We now focus on the term u(t)TH−1E>` = v(t) TH−1E>` + ε(t) TH−1E>`.
By choosing L = ` in the proof of Lemma 15, it follows that
ε(t)TH−1E≤` = od,P(1) · (σ2ε + ‖P≤`fd‖ 2 L2),
hence combining with the bound for T4 in Lemma 15 yields ε(t)TH−1E>` = T4 − ε(t)TH−1E≤`
= od,P(1) · (σ2ε + ‖fd‖ 2 L2)− od,P(1) · (σ 2 ε + ‖P≤`fd‖ 2 L2) = od,P(1) · (σ2ε + ‖fd‖ 2 L2).
We will now show that v(t)TH−1E>` = od,P(1) · ‖P>`fd‖L2 ‖fd‖L2 . (38) If `(d) = m(d), then Eq. (38) follows from Lemma 13. Otherwise, consider the case `(d) = u(d). Following similar logic to the proof of Lemma 13, because E[‖f‖22] = n ‖fd‖ 2 L2 and
E[‖E>u‖22] = n ∞∑
k=u+1
λ4d,kf̂ 2 k ≤ nλ4d,u+1 ‖P>ufd‖ 2 L2 ,
we also get Eq. (38) since
|v(t)TH−1E>`| = |fT(I− e−tH/nH−1)E>u| ≤ (t/n)‖f‖2 ∥∥∥(I− e−tH/n)(tH/n)−1∥∥∥
op ‖E>u‖2
(a) ≤ Od,P(1) · ‖P>ufd‖L2 ‖fd‖L2 tλ 2 d,u+1
(b) = od,P(1) · ‖P>ufd‖L2 ‖fd‖L2 ,
where (a) used the inequality (1 − e−x)/x ≤ 1 and (b) used Assumption 3(a). Therefore u(t)TH−1E>` = od,P(1) · (σ2ε + ‖fd‖ 2 L2), hence by Eq. (36) and Eq. (37) we obtain the final
part of the theorem ∥∥∥f̂t − P≤`fd∥∥∥2 L2 = od,P(1) · (‖fd‖2L2+η + σ 2 ε). (39)
C DOT PRODUCT KERNELS ON Sd−1( √ d)

heading: C.1 SETTING
text: We now apply our general theorems to the setting of dot product kernels on the sphere. Concretely we take Xd = Sd−1( √ d) and νd = Unif(Sd−1( √ d)) and consider dot product kernels Hd which take the form of Eq. (4). Note that by Eq. (61) any dot product kernel hd can be decomposed as
hd(〈x1,x2〉/d) = Ew∼Unif(Sd−1)[σd(〈w,x1〉)σd(〈w,x2〉)]
for some activation function σd. We state mild assumptions on σd and show that under these conditions we can apply the results in Section A.3.

heading: C.2 ASSUMPTIONS
text: We state our assumptions on σd after some definitions. See Appendix G for additional background. Denote by P≤` the orthogonal projection onto the subspace of L2(Sd−1( √ d)) spanned by polynomials of degree less than or equal to `. The projectors P` and P>` are defined analogously. Let us emphasize that the projectors P≤` are related but distinct from the P≤m: while P≤` projects onto the eigenspace of polynomials of degree at most `, P≤m projects onto the top m-eigenfunctions.
The assumptions given on the activations are the same as Assumption 3 of Mei et al. (2021a). Assumption 4 (Assumptions for Dot Product Kernels at level s ∈ N). Let {Hd}d≥1 be a sequence of dot product kernels with associated activation functions {σd}d≥1 as in Eq. (61). We assume the following hold
(a) There exists an integer k and constants c1 < 1 and c0 > 0, such that |σd(x)| ≤ c0 exp ( c1x 2/(4k) ) .
(b) We have
min k≤s ds−k ∥∥Pkσd(〈e, ·〉)∥∥2L2 = Ωd(1),∥∥P>2s+1σd(〈e, ·〉)∥∥2L2 = Ωd(1),
where e ∈ Sd−1 is a fixed vector (it is easy to see that these quantities do not depend on e).
Consider t(d), n(d) such that
dj+δ0 ≤ t ≤ dj+1−δ0 , ds+δ0 ≤ n ≤ ds+1−δ0
for some j, s ∈ N and δ0 > 0. We now verify that if {σd}d≥1 satisfies Assumption 4 at level s, then for an appropriate choice of (u(d),m(d)) the conditions in Appendix A.2 are satisfied and lead to Theorem 1. We set u(d) and m(d) to be the number of eigenvalues associated to spherical harmonics of degree less than or equal to j and s respectively
u = j∑ k=0 B(d, k) = Θd(d j), m = s∑ k=0 B(d, k) = Θd(d s).
The verification of Assumption 1 (Kernel Concentration Property) and Assumption 2 (Eigenvalue Condition) at level {(n(d),m(d)} is the same as the treatment in Theorem 2 of Mei et al. (2021a). We only need to verify Assumption 3. To see part 3(a), note that 1/λ2d,u(d) = Θd(d
j) and 1/λ2d,u(d) = Θd(d
j+1). For part 3(b), the condition holds because u(d) < m(d) for large d if and only if j < s. Assumption 3(c) is easily seen to hold since the trace of the kernel operator Tr(Hd) = Θd(1).
D GROUP INVARIANT KERNELS ON Sd−1( √ d)

heading: D.1 SETTING
text: We now apply our general theorems to the setting of group invariant kernels on the sphere. Concretely we take Xd = Sd−1( √ d) and νd = Unif(Sd−1( √ d)) and consider kernels Hd which take the form of Eq. (5) for some function h. By Eq. (61), for some activation function σd
Hd(x1,x2) = ∫ Gd Ex∼Unif(Sd−1)[σd(〈x1,w〉)σd(〈x2, g ·w〉)]πd(dg). (40)
We state mild assumptions on σd and show that under these conditions we can apply the results in Appendix A.3. For additional technical background refer to Appendix G.

heading: D.2 ASSUMPTIONS
text: We will assume that σd = σ for all d and make the following assumptions on σ which are the same as Assumption 1 in Mei et al. (2021b). Assumption 5 (Assumption on Group Invariant Kernel at level s). Let {Hd}d≥1 be a sequence of invariant kernels with associated activation functions σd = σ as in Eq. (40). We assume the following conditions hold
(a) For Gd = Cycd, we assume σ to be (s + 1) ∨ 3 differentiable and there exists constants c0 > 0 and c1 < 1 such that |σ(k)| ≤ c0ec1u
2/2 for any 2 ≤ k ≤ (s + 1) ∨ 3. For general Gd, we assume that σ is a (finite degree) polynomial function.
(b) The Hermite coefficients µk(σ) (c.f. Appendix) verify µk 6= 0 for any 0 ≤ k ≤ s.
(c) We assume that σ is not a polynomial with degrees less than or equal to s.
Consider t(d), n(d) such that
dj+δ0 ≤ t ≤ dj+1−δ0 , ds−α+δ0 ≤ n ≤ ds−α+1−δ0
for some j, s ∈ N and δ0 > 0. We now verify that if σ satisfies Assumption 5 at level s, then the conditions given in Appendix A.2 are satisfied for an appropriate choice of (u(d),m(d)) which leads to Theorem 2. We set u and m to be the number of eigenvalues invariant polynomials of degree less than or equal to j and s respectively
u = j∑ k=0 D(d, k) = Θd(d j−α), m = s∑ k=0 D(d, k) = Θd(d s−α),
where D(d, k) is the dimension of the subspace of invariant polynomials of degree k, c.f. Appendix G.6. The verification of Assumption 1 (Kernel Concentration Property) and Assumption 2 (Eigenvalue Condition) at level {(n(d),m(d)} is exactly the same as in Theorem 1 in Mei et al. (2021b). We must verify Assumption 3. To see part 3(a), note that 1/λ2d,u(d) = Θd(d
j) and 1/λ2d,u(d) = Θd(d
j+1). For part 3(b), the condition holds because u(d) < m(d) for large d if and only if j < s in which case
t(d) n(d) Tr ( Hd,m(d) ) ≤ d j+1−δ0 ds−α+δ0 Θ(d−α) = O(d−2δ0dj−s+1) = od(1).
Assumption 3(c) can be seen to hold from the fact that Tr ( Hd,>m(d) ) = Θ(ds−α) and
m∑ j=0 λ2d,j = s∑ k=0 ξ2d,kD(d, k) = Θ(sd −α) = Θ(d−α)
from which it follows that for some constant C m∑ j=0 λ2d,j ≤ C ∞∑ j>m λ2d,j .

heading: E AUXILIARY RESULTS

heading: E.1 SOLUTION TO KERNEL DYNAMICS
text: Recall that we are interested in the following dynamics given in Eqs. (2), (3),
d dt fort (x) = E[Hd(x, z)(fd(z)− fort (z))],
d dt f̂t(x) = 1 n n∑ i=1 Hd(x,xi)(yi − f̂t(xi))
with zero initialization for0 ≡ f̂0 ≡ 0. In this section we clarify the derivation, validity, and solution of these dynamics. Let us consider the maps R : Hd → R and R̂n : Hd → R defined in Eq. (1).
R(f) = ∫ Xd (f(x)− fd(x))2 dνd (x) + σ2ε ,
R̂n(f) = 1
n n∑ i=1 (f(xi)− yi)2.
First, we recall the definition of the Fréchet derivative of a functional V : Hd → R at f . The Fréchet derivative DV (f) is the linear functional such that for g ∈ Hd,
lim ‖g‖Hd→0 |V (f + g)− V (f)−DV (f)(g)| ‖g‖Hd = 0.
The gradient ∇V (f) ∈ Hd is defined such that 〈∇V (f), g〉Hd = DV (f)(g)
exists uniquely by the Riesz representation theorem. The gradients of the risk functionals are
∇R(f) = Hd(f − fd),
∇R̂n(f) = 1
n n∑ i=1 (f(xi)− yi)Hxi ,
where Hxi(x) := Hd(xi,x) and Hd is the kernel operator as in Section A.1. A proof of this fact is given in Proposition 2.1 of Yao et al. (2007). Taking fort (x) as shorthand for f
or(t, x) where for(t, ·) ∈ Hd is the oracle model at time t and similarly for f̂t(x), the following gradient flows with zero initialization are well-defined for t ≥ 0,
d dt fort = −∇R(fort ) = −Hd(fort − fd) = Ez[Hd(·, z)(fd(z)− fort (z)] (41)
d dt f̂t = −∇R̂n(f̂t) = − 1 n n∑ i=1 (f̂t(xi)− yi)Hxi = 1 n n∑ i=1 Hd(·,xi)(yi − f̂t(xi)). (42)
The oracle model ODE Eq. (41) is simply a linear differential equation which has the following solution involving the operator exponential exp(A) := ∞∑ k=0 Ak/k!
fort = fd + exp(−tHd)(for0 − fd) = fd − exp(−tHd)fd. (43) For the empirical model ODE Eq. (42) we first consider the system of scalar differential equations induced at the points {(xi, yi)}i∈[n]. Letting u(t) = (f̂t(x1), . . . , f̂t(xn))T, y = (y1, . . . , yn)T, andH = (Hd(xi,xj))i,j∈[n] we have
d dt u(t) = − 1 n H(u(t)− y),
with initial condition u(0) = 0. As this a linear ODE, the solution is given by
u(t) = y + e−tH/n(u(0)− y) = (In − e−tH/n)y. (44)
For x ∈ Rd, define h(x) = (Hd(x,x1), . . . ,Hd(x,xn))T ∈ Rn. Let a(t) = H−1u(t) ∈ Rn. We will show that the function f̂t(·) := 〈h(·),a(t)〉 ∈ Hd, which satisfies (f̂t(x1), . . . , f̂t(xn))T = u(t), satisfies the following equation
d dt f̂t(x) = 1 n 〈h(x),y − u(t)〉 = 1 n 〈h(x), e−tH/ny〉,
which is Eq. (42) at point x. Indeed, by the chain rule
d dt f̂t(x) = d dt 〈h(x),a(t)〉
= 〈h(x), d dt a(t)〉 = 〈h(x),H−1 1 n He−tH/ny〉 = 1
n 〈h(x), e−tH/ny〉
which is what we wanted to show.

heading: E.2 EQUIVALENCE BETWEEN INVARIANT KERNELS AND DATA AUGMENTATION
text: In this section we will show an equivalence between the (time rescaled) gradient flows for training invariant kernels and using an augmented dataset. Specifically consider a group G and a kernel H that is G-equivariant, that is
H(g · x1, g · x2) = H(x1,x2) ∀g ∈ G, ∀x1,x2 ∈ X .
Given a G-equivariant kernel H , we define a G-invariant kernel Hinv as the group averaged kernel
Hinv = ∫ G H(x1, g · x2)π(dg)
for the Haar measure π on G (c.f. Eq. (5)). Note that any dot product kernel is G-equivariant for G a subgroup the orthogonal group e.g. the cyclic group Cyc (c.f. Section 3.3).
Given a dataset (X,y) = {(xi, yi) : i ∈ [n]} consider the augmented dataset
(XG ,yG) = {(g · xi, yi) : g ∈ G, i ∈ [n]}.
We consider the (rescaled c.f. Remark 5) empirical dynamics Eq. (3) of the gradient flow on (X,y) using Hinv which we denote f̂t,inv
d dt f̂t,inv(x) = −m∇R̂n(f̂t,inv) = − m n n∑ i=1 (f̂t,inv(xi)− yi)Hinv(xi,x)
and the empirical dynamics of the gradient flow on (XG ,yG) using H which we denote f̂t,aug
d dt f̂t,aug(x) = − 1 n ∑ g∈G n∑ i=1 (f̂t,aug(g · xi)− yi)H(g · xi,x).
Proposition 3. Let G be a finite group with m elements. Given a G-equivariant kernel H , if π is the uniform measure on G then
f̂t,inv ≡ f̂t,aug, ∀t ≥ 0.
Proof. Let G = {g1, . . . , gm} where g1 is the identity. Define the output vectors
uinv(t) := (f̂t,inv(x1), . . . , f̂t,inv(xn)) T ∈ Rn,
ug(t) := (f̂t,aug(g · x1), . . . , f̂t,aug(g · xn))T ∈ Rn, uaug(t) := (ug1(t), . . . ,ugm(t)) T ∈ Rmn.
Furthermore, define the kernel matrices
Hg,g′ := [H(g · xi, g′ · xj)]i,j∈[n] ∈ Rn×n for g, g′ ∈ G, Haug := [Hg,g′ ]g,g′∈G ∈ Rmn×mn, Hinv = [Hinv(xi,xj)]i,j∈[n].
Note that by definitionHinv = 1m m∑ j=1 Hg1,gj . We will show that
uaug(t) = (uinv(t),uinv(t), . . . ,uinv(t)) for all t ≥ 0. (45)
From this the result follows by Theorem 4.1 in Li et al. (2019) since f̂t,inv, f̂t,aug are given by kernel regressions with targets uinv(t),uaug(t) and kernels Hinv, H respectively.
By Eq. (42), we can write
uinv(t) = (I− exp(−tmHinv/n))y, uaug(t) = (I− exp(−tHaug/n))yG ,
where yG = (y, . . . ,y) ∈ Rmn. By expanding the matrix exponential series and using linearity, it suffices to show that
HkaugyG = (m kHkinvy,m kHkinvy, . . . ,m kHkinvy) for all k ∈ N.
in order to show Eq. (45) holds. We prove the above by induction on k. For k = 1, observe that
HaugyG =  m∑ j=1 Hgi,gjy m i=1
=  m∑ j=1 Hg1,gjy m i=1 = (mHinvy) m i=1
where the second equality follows from G-equivariance of H . Assume the inductive hypothesis holds for k. Then
Hk+1aug yG = HaugH k augyG
= mk m∑ j=1 Hgi,gjH k invy m i=1
=  m∑ j=1 Hg1,gj  m∑ j′=1 Hg1,gj′ ky  m
i=1
=   m∑ j=1 Hg1,gj k+1y  m
i=1
= (mk+1Hk+1inv y) m i=1
where the second equality applies the induction hypothesis and the third equality uses equivariance. Thus the inductive claim is proved and the proof is complete.
Remark 5. The scaling factor m in the gradient flow for f̂t,inv, leads to a natural comparison with f̂t,aug as elaborated in Appendix E.3. As argued in that section, in the gradient descent discretization, it is natural to take a step-size inversely proportional to the maximum kernel eigenvalue. In the case of high-dimensional invariant kernels, note that
λmax(Haug) = mλmax(H) ∼ mλmax(Hinv) hence the step-size for the invariant kernel flow should be m times larger.

heading: E.3 DISCRETIZING TIME
text: Comparing different “speeds” of optimization algorithms only makes sense for discrete-time algorithms. Consider the following gradient descent dynamics with step-size η, obtained as the discretization of the empirical gradient flow Eq. (42)
f̂k+1 = f̂k − η∇R̂n(f̂k) = f̂k − η 1
n n∑ i=1 (f̂k(xi)− yi)H(·,xi), k = 0, 1, . . . (46)
We will argue that it is natural to take η ∼ n/λmax(H) whereH is the kernel matrix. Define the sampling operator S : Hd → Rn by S(f) = (f(xi))ni=1 ∈ Rn and let S∗ : Rn → Hd be its adjoint, defined by S∗(y) = 1n n∑ i=1 yiHxi (see Yao et al. (2007) Appendix B for more details). Then we can rewrite the gradient descent equation Eq. (46) as
f̂k+1 = f̂k − η(S∗Sf̂t − S∗y), k = 0, 1, . . . (47) Let T = S∗S and define H := SS∗ = 1nH to be the normalized kernel matrix. Let b := S∗H −1 y ∈ Hd and note that since Tb = S∗y, we can rewrite Eq. (47) as
f̂k+1 = f̂k − ηT (f̂k − b), k = 0, 1, . . . (48) Denote the eigenvalues of H as λ1 ≥ λ2 ≥ . . . ≥ λn > 0. By the spectral theorem, there exists a set of orthornomal eigenvectors φ1, . . . , φn ∈ Hd such that T =
n∑ i=1 λiφiφ ∗ i . Define αi(k) =
〈f̂k − b, φi〉 ∈ R. By taking Eq. (48) then subtracting b and taking the inner product with φi on both sides, we get the coordinate evolution equations for i = 1, . . . , n
αi(k + 1) = 〈(id− ηT )(f̂k − b), φi〉 = 〈(f̂k − b), (id− ηT )φi〉 = (1− ηλi)αi(k)
where the second equality holds since T is self-adjoint and the last equality is since φi is an eigenvector of T . It is easy to see that αi(k) = (1− ηλi)kαi(0). Therefore we see that gradient descent Eq. (46) is guaranteed to converge if η < 1/(2λ1) and may not otherwise. Therefore it is natural to choose the step-size η to scale asymptotically as η ∼ 1/λmax(H). For a dot product kernel H and its corresponding invariant kernel Hinv the kernel matrices have operator norms of the same order
λmax(H) ∼ λmax(Hinv) hence no time rescaling is need to compare the corresponding optimization speeds asymptotically.

heading: E.4 SIMILARITIES WITH EMPIRICAL PHENOMENA
text: In this section we elaborate upon Remark 4 and mention some connections with empirical observations in Nakkiran et al. (2020). Although the metric in our setting is the squared loss, we can still observe three stages in classification problems when measuring the soft error. In Fig. 6a taken from Nakkiran et al. (2020) we can observe stage 1 and stage 2. Either training has not continued long enough to observe stage 3 or n is large enough so that the models have converged to the approximation error of the neural network class (c.f. Remark 1). In Fig. 6b taken from Nakkiran et al. (2020), although the train errors are not plotted, by extrapolating from Fig. 6a, presumably for each n stage 1 and stage 2 occur. From the dimmer curves in Fig. 6b we can see that for n < 50000 stage 3 occurs as well.
In Fig. 7a, we see a parallel between the use of cyclic versus dot product kernels and the use data augmentation versus not for a Resnet-18 trained on CIFAR-5m (note that using a cyclic kernel is equivalent to using a dot product kernel with data-augmentation c.f. Appendix E.2). In both our theoretical results and in the empirical results of Nakkiran et al. (2020) we observe that the ideal world optimization speed of augmented and non-augmented training are the same, but for augmented training the real world training speed is slowed down, eventually leading to better generalization for long enough training.

heading: F ADDITIONAL FIGURES
text: To see the effects of varying the dimension d, in Fig. 8 we replicate the log-scale plots of kernel gradient flow with dot product kernels from Fig. 4. We take n = d1.5 and vary d ∈ {50, 100, 200, 400}. Each plot is averaged over 10 runs with the shaded region representing one standard deviation around the mean. We can see that as d increasing the standard deviation decreases and the curves approach the theoretical high-dimensional prediction.

heading: G TECHNICAL BACKGROUND

heading: G.1 NOTATIONS
text: For a positive integer, we denote by [n] the set {1, 2, . . . , n}. For vectors u,v ∈ Rd, we denote 〈u,v〉 = u1v1 + . . . + udvd their scalar product, and ‖u‖2 = 〈u,u〉1/2 the `2 norm. Given a matrix A ∈ Rn×m, we denote ‖A‖op = max‖u‖2=1 ‖Au‖2 its operator norm and by ‖A‖F =(∑
i,j A 2 ij
)1/2 its Frobenius norm. If A ∈ Rn×n is a square matrix, the trace of A is denoted by
Tr(A) = ∑ i∈[n]Aii.
We use Od( · ) (resp. od( · )) for the standard big-O (resp. little-o) relations, where the subscript d emphasizes the asymptotic variable. Furthermore, we write f = Ωd(g) if g(d) = Od(f(d)), and f = ωd(g) if g(d) = od(f(d)). Finally, f = Θd(g) if we have both f = Od(g) and f = Ωd(g).
We use Od,P( · ) (resp. od,P( · )) the big-O (resp. little-o) in probability relations. Namely, for h1(d) and h2(d) two sequences of random variables, h1(d) = Od,P(h2(d)) if for any ε > 0, there exists Cε > 0 and dε ∈ Z>0, such that P(|h1(d)/h2(d)| > Cε) ≤ ε, ∀d ≥ dε, and respectively: h1(d) = od,P(h2(d)), if h1(d)/h2(d) converges to 0 in probability. Similarly, we will denote h1(d) = Ωd,P(h2(d)) if h2(d) = Od,P(h1(d)), and h1(d) = ωd,P(h2(d)) if h2(d) = od,P(h1(d)). Finally, h1(d) = Θd,P(h2(d)) if we have both h1(d) = Od,P(h2(d)) and h1(d) = Ωd,P(h2(d)).

heading: G.2 FUNCTIONAL SPACES OVER THE SPHERE
text: For d ≥ 3, we let Sd−1(r) = {x ∈ Rd : ‖x‖2 = r} denote the sphere with radius r in Rd. We will mostly work with the sphere of radius √ d, Sd−1( √ d) and will denote by τd the uniform
probability measure on Sd−1( √ d). All functions in this section are assumed to be elements of L2(Sd−1( √ d), τd), with scalar product and norm denoted as 〈 · , · 〉L2 and ‖ · ‖L2 :
〈f, g〉L2 ≡ ∫ Sd−1( √ d) f(x) g(x) τd(dx) . (49)
For ` ∈ Z≥0, let Ṽd,` be the space of homogeneous harmonic polynomials of degree ` on Rd (i.e. homogeneous polynomials q(x) satisfying ∆q(x) = 0), and denote by Vd,` the linear space of functions obtained by restricting the polynomials in Ṽd,` to Sd−1( √ d). With these definitions, we have the following orthogonal decomposition
L2(Sd−1( √ d), τd) = ∞⊕ `=0 Vd,` . (50)
The dimension of each subspace is given by
dim(Vd,`) = B(d, `) = 2`+ d− 2 d− 2
( `+ d− 3
`
) . (51)
For each ` ∈ Z≥0, the spherical harmonics {Y (d)`,j }1≤j≤B(d,`) form an orthonormal basis of Vd,`:
〈Y (d)ki , Y (d) sj 〉L2 = δijδks.
Note that our convention is different from the more standard one, that defines the spherical harmonics as functions on Sd−1(1). It is immediate to pass from one convention to the other by a simple scaling. We will drop the superscript d and write Y`,j = Y (d) `,j whenever clear from the context. We denote by Pk the orthogonal projections to Vd,k in L2(Sd−1( √ d), τd). This can be written in terms of spherical harmonics as
Pkf(x) ≡ B(d,k)∑ l=1 〈f, Ykl〉L2Ykl(x). (52)
We also define P≤` ≡ ∑` k=0 Pk, P>` ≡ I−P≤` = ∑∞ k=`+1 Pk, and P<` ≡ P≤`−1, P≥` ≡ P>`−1.

heading: G.3 GEGENBAUER POLYNOMIALS
text: The `-th Gegenbauer polynomialQ(d)` is a polynomial of degree `. Consistently with our convention for spherical harmonics, we view Q(d)` as a function Q (d) ` : [−d, d] → R. The set {Q (d) ` }`≥0 forms
an orthogonal basis on L2([−d, d], τ̃1d ), where τ̃1d is the distribution of √ d〈x, e1〉 when x ∼ τd, satisfying the normalization condition:
〈Q(d)k ( √ d〈e1, ·〉), Q(d)j ( √ d〈e1, ·〉)〉L2(Sd−1(√d)) =
1
B(d, k) δjk . (53)
In particular, these polynomials are normalized so that Q(d)` (d) = 1. As above, we will omit the superscript (d) in Q(d)` when clear from the context. Gegenbauer polynomials are directly related to spherical harmonics as follows. Fix v ∈ Sd−1( √ d) and consider the subspace of V` formed by all functions that are invariant under rotations in Rd that keep v unchanged. It is not hard to see that this subspace has dimension one, and coincides with the span of the function Q(d)` (〈v, · 〉). We will use the following properties of Gegenbauer polynomials
1. For x,y ∈ Sd−1( √ d)
〈Q(d)j (〈x, ·〉), Q (d) k (〈y, ·〉)〉L2 =
1
B(d, k) δjkQ
(d) k (〈x,y〉). (54)
2. For x,y ∈ Sd−1( √ d)
Q (d) k (〈x,y〉) =
1
B(d, k) B(d,k)∑ i=1 Y (d) ki (x)Y (d) ki (y). (55)
These properties imply that, up to a constant, Q(d)k (〈x,y〉) is a representation of the projector onto the subspace of degree-k spherical harmonics
(Pkf)(x) = B(d, k) ∫ Sd−1( √ d) Q (d) k (〈x,y〉) f(y) τd(dy) . (56)
For a function σ ∈ L2([− √ d, √ d], τ1d ) (where τ 1 d is the distribution of 〈e1,x〉 when x ∼
Unif(Sd−1( √ d))), denoting its spherical harmonics coefficients ξd,k(σ) to be
ξd,k(σ) = ∫ [− √ d, √ d] σ(x)Q (d) k ( √ dx)τ1d (dx), (57)
then we have the following equation holds in L2([− √ d, √ d], τ1d ) sense
σ(x) = ∞∑ k=0 ξd,k(σ)B(d, k)Q (d) k ( √ dx).
For any rotationally invariant kernel Hd(x1,x2) = hd(〈x1,x2〉/d), with hd( √ d · ) ∈
L2([− √ d, √ d], τ1d ), we can associate a self adjoint operatorHd : L2(Sd−1( √ d))→ L2(Sd−1( √ d))
Hdf(x) ≡ ∫ Sd−1( √ d) hd(〈x,x1〉/d) f(x1) τd(dx1) . (58)
By rotational invariance, the space Vk of homogeneous polynomials of degree k is an eigenspace of Hd, and we will denote the corresponding eigenvalue by ξd,k(hd). In other words Hdf(x) ≡∑∞ k=0 ξd,k(hd)Pkf . The eigenvalues can be computed via
ξd,k(hd) = ∫ [− √ d, √ d] hd ( x/ √ d ) Q (d) k ( √ dx)τ1d (dx) . (59)
For a dot product kernel Hd(x,y) = hd(〈x,y〉/d) consider the Gegenbauer expansion of hd in L2([− √ d, √ d], τ1d )
hd(〈x,y〉/d) = ∞∑ k=0 ξk,d(hd)B(d, k)Q (d) k (〈x,y〉). (60)
Using Eq. (54) we can equivalently write the kernel as an expectation over random features for some activation σd hd(〈x,y〉/d) = Ew∼Unif(Sd−1)[σd(〈w,x〉)σd(〈w,y〉)] (61) by taking
σd(x) = ∞∑ k=0 ξd,k(hd) 1/2B(d, k)Q (d) k ( √ dx). (62)
Note that σd ∈ L2([− √ d, √ d], τ1d ) as long as h(1) <∞.

heading: G.4 HERMITE POLYNOMIALS
text: The Hermite polynomials {Hek}k≥0 form an orthogonal basis of L2(R, γ), where γ(dx) = e−x 2/2dx/ √
2π is the standard Gaussian measure, and Hek has degree k. We will follow the classical normalization (here and below, expectation is with respect to G ∼ N (0, 1)):
E { Hej(G) Hek(G) } = k! δjk . (63)
As a consequence, for any function g ∈ L2(R, γ), we have the decomposition
g(x) = ∞∑ k=0 µk(g) k! Hek(x) , µk(g) ≡ E { g(G) Hek(G)} . (64)
The Hermite polynomials can be obtained as high-dimensional limits of the Gegenbauer polynomials introduced in the previous section. Indeed, the Gegenbauer polynomials (up to a √ d scaling in domain) are constructed by Gram-Schmidt orthogonalization of the monomials {xk}k≥0 with respect to the measure τ̃1d , while Hermite polynomial are obtained by Gram-Schmidt orthogonalization with respect to γ. Since τ̃1d ⇒ γ (here⇒ denotes weak convergence), it is immediate to show that, for any fixed integer k,
lim d→∞
Coeff{Q(d)k ( √ dx)B(d, k)1/2} = Coeff
{ 1
(k!)1/2 Hek(x)
} . (65)
Here and below, for P a polynomial, Coeff{P (x)} is the vector of the coefficients of P . As a consequence, for any fixed integer k, we have
µk(σ) = lim d→∞
ξd,k(σ)(B(d, k)k!) 1/2, (66)
where µk(σ) and ξd,k(σ) are given in Eq. (64) and Eq. (57).

heading: G.5 THE INVARIANT FUNCTION CLASS AND THE SYMMETRIZATION OPERATOR
text: Let Gd be a group that is isomorphic to a subgroup of O(d), the orthogonal group in d dimension. That means, each element of Gd can be identified with a matrix in O(d) ⊆ Rd×d, and the group addition operation in Gd can be regarded as matrix multiplications in O(d). For any x ∈ Sd−1( √ d) and g ∈ Gd, we define group action g · x to be the multiplication of matrix representation of g with the vector x. We equip Gd with a probability measure πd, which is the uniform probability measure on Gd. More specifically, the Borel sigma algebra on Gd is defined as the Borel sigma algebra ofO(d) restricted on Gd. The uniform probability measure πd satisfies the property that, for any Borel-measurable set B ⊆ Gd and any g ∈ Gd, we have
πd(B) = πd(gB).
Let L2(Sd−1( √ d)) be the class of L2 functions on Sd−1( √ d) equipped with uniform probability measure Unif(Sd−1( √ d)). We define the invariant function class to be
L2(Sd−1( √ d),Gd) = { f ∈ L2(Sd−1( √ d)) : f(x) = f(g · x), ∀x ∈ Sd−1( √ d), ∀g ∈ Gd } .
We define the symmetrization operator S : L2(Sd−1( √ d))→ L2(Sd−1( √ d),Gd) to be
(Sf)(x) = ∫ Gd f(g · x)πd(dg).

heading: G.6 ORTHOGONAL POLYNOMIALS ON INVARIANT FUNCTION CLASS
text: We define Vd,≤k ⊆ L2(Sd−1( √ d)) to be the subspace spanned by all the degree ` polynomials, Vd,>k ≡ V ⊥d,≤k ⊆ L2(Sd−1( √ d)) to be the orthogonal complement of Vd,≤k, and Vd,k = Vd,≤k ∩ V ⊥d,≤k−1. In words, Vd,k contains all degree k polynomials that orthogonal to all polynomials of degree at most k − 1. We further define Vd,<k = Vd,≤k−1 and Vd,≥k = Vd,>k−1. Let P≤` to be the projection operator on L2(Sd−1( √ d),Unif) that project a function onto Vd,≤`, the space spanned by all the degree ` polynomials. Then it is easy to see that P≤` and S operator commute. This means, for any f ∈ L2(Sd−1( √ d)), we have
P≤`[S(f)] = S[P≤`(f)].
Similarly, we can define P`, P<`, P>`, P≥`, which commute with S. We denote Vd,`(Gd) ≡ P`(Sd−1( √ d),Gd) to be the space of polynomials in the images of P`S. Then we have
P`(Ad,Gd) = P`(L2(Sd−1( √ d),Gd)) = S[P`(L2(Ad))].
We denote D(d; k) = D(Sd−1( √ d);Gd; k) ≡ dim(Pk(Sd−1( √ d),Gd)) to be the dimension of Pk(Sd−1( √ d),Gd). We denote {Y (d)
kl }l∈[D(Sd−1(√d);k)] to be a set of orthonormal polynomial basis in Pk(Sd−1( √ d),Gd). That means
Ex∼Unif(Sd−1(√d))[Y (d) k1l1(x)Y (d) k2l2(x)] = 1{k1 = k2, l1 = l2},
and Y (d)
kl (x) = Y (d) kl (g · x), ∀x ∈ Sd−1( √ d), ∀g ∈ Gd.

Review summaries:
Studied the evolution of generalization error of the kernel gradient flow trajectory with respect to the training (empirical world) and population (ideal world) MSE loss. The analysis builds upon [Mei et al. 2021], which relates kernel ridge regression to projection onto low-degree polynomials. The authors showed that the estimator optimizing the empirical risk achieves vanishing training error, but the test error plateaus at certain value depending on the training set size, whereas the online (population) estimator can learn increasingly complex components of the target function as training proceeds.

This paper studies the learning dynamics of gradient flow for kernel ridge regression. The authors contrast to different setups: the "empirical world", where the model is trained on a finite data set, and the "oracle world", where the model is instead trained directly on the population loss, or in other words, by minimising the test error. Nakkiran et al. (2020) recently reported that the errors in both setups stay close to each throughout training for a range of deep neural networks. By studying the relation between empirical and oracle world in the setup of kernel ridge regression, the authors investigate this phenomenon in a setup under precise theoretical control.  The authors perform such a study by leveraging the analysis of Ghorbani et al. (2021) and Mei et al. (2021), who gave a detailed and precise analysis of the implicit bias of learning in random feature models, and in particular the learning of approximations of increasing complexity.   Contrasting the learning in the two worlds, the authors find that learning generally, but not always, proceeds in three stages. In the case of polynomial kernels, models in both the empirical and the oracle learn first the leading $\ell$-components of the target function, defined as those projections of the target function along the polynomials of small degree. After the training error of the empirical model reaches zero, a second phase ensues where the test errors of both models remain close, but the gap between training and test error in the empirical world can be large. Finally, either there is enough training data to learn the target function perfectly, or learning enters a third stage where the oracle model learns the target function perfectly, which the empirical model won't achieve.  The authors confirm this picture also for invariant kernels, say translation-invariant ones, using technical tools from Mei et al. (2021). Quantitatively, the presence of an invariance shows up in the training speed, which is slower than for the product kernel (Remark 2). Numerical experiments for kernel least-square confirm the theory, while the authors also report experiments for SGD for random feature regression.

The paper studies the training/learning dynamics of (inner product) kernel gradient descent in the high-dimensional setting. The main contribution is proving a three-stage learning dynamics which were observed in neural networks in existing work.  (1) Stage 1: train loss ~ test loss ~ oracle test loss (training set == whole (input, label) distribution)  (2) Stage 2: test loss ~ oracle loss; train loss ->0;  (3) Stage 3: test loss unchanged; oracle loss decays;   I think this is a nice application of the results developed in [1], showing that many seemingly surprising deep learning (DL) phenomena are indeed not unique in DL -- most of them are provable / observable in the kernel setting.   [1]Generalization error of random features and kernel methods: hypercontractivity and kernel matrix concentration.