# Precise Learning Curves and Higher-Order Scaling Limits for Dot Product
Kernel Regression

Lechao Xiao
Google Research, Brain Team
<EMAIL_ADDRESS>Jeffrey Pennington
Google Research, Brain Team
<EMAIL_ADDRESS>
###### Abstract

As modern machine learning models continue to advance the computational
frontier, it has become increasingly important to develop precise estimates
for expected performance improvements under different model and data scaling
regimes. Currently, theoretical understanding of the learning curves that
characterize how the prediction error depends on the number of samples is
restricted to either large-sample asymptotics ($m\to\infty$) or, for certain
simple data distributions, to the high-dimensional asymptotics in which the
number of samples scales linearly with the dimension ($m\propto d$). There is
a wide gulf between these two regimes, including all higher-order scaling
relations $m\propto d^{r}$, which are the subject of the present paper. We
focus on the problem of kernel ridge regression for dot-product kernels and
present precise formulas for the test error, bias, and variance, for data
drawn uniformly from the sphere in the $r$th-order asymptotic scaling regime
$m\to\infty$ with $m/d^{r}$ held constant. We observe a peak in the learning
curve whenever $m\approx d^{r}/r!$ for any integer $r$, leading to multiple
sample-wise descent and nontrivial behavior at multiple scales.

##  1 Introduction

Modern machine learning has entered an era in which scaling is arguably the
most critical ingredient to improve performance. Recent breakthroughs such as
GPT-3 (Kaplan et al., 2020) and PaLM (Chowdhery et al., 2022) have
demonstrated that performance of various learning algorithms improves in a
predictable manner as the amount of data and computational resources used in
training increases. The functional relationships between performance and
resources are loosely referred to as learning curves. While extrapolation of
empirical learning curves is widely used to make predictions about how a model
might perform when extra resources become available, a rigorous theoretical
understanding is lacking. A fundamental obstacle in developing a detailed
theoretical model of such learning curves is that they depend on many moving
parts, e.g. the data distribution, the network architecture, the training
algorithm, among others. In addition, even in the simplest possible settings,
the learning curves can exhibit non-trivial structure that naive scaling laws
fail to model, e.g. the well-known double-descent phenomenon (Belkin et al.,
2019; Advani et al., 2020).

In the past couple years, a large amount of effort from the community has
improved our theoretical understanding of such phenomena and in some cases
precise characterizations of learning curves have been obtained (see e.g.
(Hastie et al., 2022; Adlam and Pennington, 2020a; Mei and Montanari, 2022)).
These results have helped clarify several puzzling empirical observations,
such as the origin of the double-descent peak (Adlam and Pennington, 2020b;
Lin and Dobriban, 2021; d’Ascoli et al., 2020) and linear trends between in-
and out-of-distribution generalization performance (Tripuraneni et al., 2021a,
b; Mania and Sra, 2020), among many others. However, the precise predictions
from many of these analyses have been possible only in the linear high-
dimensional scaling regime in which the number of training samples $m$ scales
linearly with the dimension $d$, i.e. $m\propto d$. In these asymptotics, the
model’s effective capacity is limited to linear functions of the features. In
contrast, many state-of-the-art models operate in a regime where the amount of
data is much larger than the data dimensionality; for example, large text
corpora can contain trillions of tokens, whereas the effective input
dimensionality of language models is at most millions . Therefore, going
beyond the linear scaling regime ($m\propto d$) to higher-order scaling
regimes ($m\propto d^{r}$) is essential in improving our understanding of
modern machine learning systems, and is the focus of the current paper.

Several works have investigated the behavior of the learning curves for
nonlinear scalings in the dot-product kernel or random features setting, but
they have done so only in the noncritical regime $m\not\propto d^{r}$
(Ghorbani et al., 2021; Mei et al., 2021; Liang et al., 2020). Extending these
result to include the critical regime $m\propto d^{r}$ is nontrivial, both
from the technical perspective and also from the phenomenological perspective,
as we shall see the critical behavior induces nonmonotonicity and multiple
sample-wise descent.

Figure 1: Precise Sample-wise Learning Curves for One-hidden Layer CNN
kernels. The theoretical predictions (Eq. (18), solid lines) agree with
finite-size simulations (markers) across several orders of magnitude and
captures cases in which the curves are relatively simple (monotonically
decreasing, small spectral gap) and complex ( multiple-descent, large spectral
gap). Simulations are obtained from kernel regression with one-layer CNN
kernels averaged over 50 runs. The input is of shape $d=d_{0}\times p$ with
size $d_{0}=14$ and number of patches $p=6$. We vary the kernels by varying
the ratio (aka spectral gap) between consecutive eigenspaces, where the ratio
$\mathrm{Gap}\in[{\bf\color[rgb]{0.9609375,0.39453125,0.39453125}2},{\bf\color[rgb]{0.28125,0.734375,0.46875}8},{\bf\color[rgb]{0.9296875,0.5390625,0.2109375}32},{\bf\color[rgb]{0.2578125,0.6015625,0.8828125}128}]$.

In this work, we obtain precise formulas for the sample-wise learning curves
in the kernel ridge regression setting for a family of dot-product kernels for
spherical input data in the polynomial scaling regimes $m\propto d^{r}$ for
all $r\in{\mathbb{N}}^{*}$. This family of kernels includes the neural network
Gaussian Process (NNGP) kernels and Neural Tangent Kernels (NTK) associated
with multi-layer fully-connected networks or convolutional networks. Both
kernels serve as important starting points towards a deeper understanding of
neural networks as they often capture the first order learning dynamics of
neural networks in certain scaling limits (Jacot et al., 2018; Lee et al.,
2019; Arora et al., 2019).

###  1.1 Contributions

Our primary contributions are to establish the following, for data drawn
uniformly from the sphere:

  1. 1.

The empirical spectral density of the Gram matrix induced by degree-$r$
spherical harmonics converges to a Marchenko-Pastur distribution when
$(d^{r}/r!)/m$ converges to a positive constant as $d\to\infty$ (Theorem 1);

  2. 2.

A precise closed-form formula for the sample-wise learning curves for dot-
product kernel regression when $m\propto d^{r}$ for all $r\in{\mathbb{N}}^{*}$
as $d\to\infty$ (Theorem 2);

  3. 3.

Empirically, the theoretical predictions agree with finite-size simulations
surprisingly well even in the strong finite-size correction regime (Fig. 1);

  4. 4.

An extension of the above results to convolutional kernels (Section 5).

While finalizing this manuscript, we became aware of two concurrent works also
studying the learning curves of dot-product kernel regression in the higher-
order scaling limits (Misiakiewicz, 2022; Hu and Lu, 2022). Their main results
are similar to our Theorems 1 and 2 and are obtained using similar techniques.
Relative to those works, our paper establishes slightly more general results
on convolutional kernels (Theorem 4), and provides stronger empirical support
through finite-size simulations (Fig. 1) and improved formulas (Corollaries 1
and 2).

##  2 Notation and Setup

Let $\bm{{\mathcal{X}}}={\mathbb{S}}_{d-1}$ denote the input space, where
${\mathbb{S}}_{d-1}$ is the unit sphere in $\mathbb{R}^{d}$ and
$\bm{{\mathcal{X}}}$ is equipped with the normalized uniform measure $\sigma$.
We use $\Delta_{d}$ to represent any quantity (a scalar, vector or a matrix)
with $|||\Delta_{d}|||\to 0$ as $d\to\infty$ (in probability if $\Delta_{d}$
is stochastic), where $|||\cdot|||$ can be the absolute value of a scalar, the
norm of a vector or the operator norm of a matrix.

Let ${\bm{X}}\in{\mathbb{R}}^{m\times d}$ be the training inputs where the
$i$-th row of ${\bm{X}}$ is ${\bm{x}}_{i}^{T}$. We assume
$\\{{\bm{x}}_{i}\\}_{i\in[m]}$ is sampled uniformly, iid from
$\bm{{\mathcal{X}}}$. The label function $f:{\mathbb{S}}_{d-1}\to{\mathbb{R}}$
will be defined in Section 4. Let $K=K^{(d)}$ be a dot-product kernel defined
on ${\mathbb{S}}_{d-1}\times{\mathbb{S}}_{d-1}$, i.e.,
$K({\bm{x}},{\bm{x}}^{\prime})=h({\bm{x}}^{T}{\bm{x}}^{\prime})$ for some
function $h\in[-1,1]\to\mathbb{R}$. We assume $h$ has the following
decomposition

$\displaystyle
h(t)=\sum_{k=1}^{\infty}\hat{h}_{k}^{2}P_{k}(t),\quad\mathrm{with}\quad\sum_{k=1}^{\infty}\hat{h}_{k}^{2}<\infty\,,$
(1)

where $P_{k}$ is the $k$-th order Legendre polynomials in $d$ dimensions. For
simplicity, we assume $\hat{\bm{h}}=(\hat{h}_{k})_{k\geq 1}$ is a sequence
that is independent of $d$ and $\hat{h}_{k}\neq 0$ for all $k\leq k_{0}$ where
$k_{0}$ is sufficiently large. As such, we can decompose the kernel function
using sperical harmonics,

$\displaystyle
K({\bm{x}},{\bm{x}}^{\prime})=\sum_{k=1}^{\infty}\sigma_{k}^{2}\sum_{l\in[N(d,k)]}Y_{k,l}({\bm{x}})Y_{k,l}({\bm{x}}^{\prime})=\sum_{k=1}^{\infty}\sigma_{k}^{2}Y_{k}({\bm{x}})^{T}Y_{k}({\bm{x}}^{\prime})\,,$
(2)

where $Y_{k,l}$ is the $l$-th spherical harmonic of degree $k$,
$N(d,k)=d^{k}/k!+O(d^{k-1})$ is the total number of degree $k$ spherical
harmonics in $d$ dimensions, $\sigma_{k}^{2}=\hat{h}_{k}^{2}/N(d,k)$ is the
eigenvalue of $Y_{kl}$, and $Y_{k}({\bm{x}})$ is the column vector
$[Y_{k,l}({\bm{x}})]_{l\in[N(d,k)]}^{T}$. We also denote by $Y_{k}({\bm{X}})$
the $m\times N(d,k)$ matrix whose $i$-th row is $Y_{k}({\bm{x}}_{i})^{T}$.

##  3 Structure of the Empirical Kernel and Marchenko-Pastur Distribution

The structure of the empirical kernel matrix $K({\bm{X}},{\bm{X}})$ plays a
critical role in characterizing the sample-wise test error for the kernel
ridge regressor associated to $K$. We assume the training set size scales
polynomially, i.e. $m\sim d^{r}$ for some positive integer
$r\in{\mathbb{N}}^{*}$. Decompose this kernel into low-, critical- and high-
frequency modes as follows,

$\displaystyle
K({\bm{X}},{\bm{X}})=\sum_{k<r}\sigma_{k}^{2}Y_{k}({\bm{X}})Y_{k}({\bm{X}})^{T}+\sigma_{r}^{2}Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}+\sum_{k>r}\sigma_{k}^{2}Y_{k}({\bm{X}})Y_{k}({\bm{X}})^{T}\,.$
(3)

The low- and high-frequency parts have simple structures since $N(d,k)/m$
either diverges to infinity or converges to zero with rate as least $d^{\pm
1}$, yielding concentration that results in significant simplification. To be
precise, for high-frequency modes $k>r$, $Y_{k}({\bm{X}})$ is a “fat" matrix
and $Y_{k}({\bm{X}})Y_{k}({\bm{X}})^{T}/N(d,k)={\bm{I}}_{m}+\Delta_{d}$ where
$\Delta_{d}$ vanishes as $d\to\infty$ (Mei et al., 2021). Thus, the high-
frequency parts behave like a regularizer in the following sense,

$\displaystyle\sum_{k>r}\sigma_{k}^{2}Y_{k}({\bm{X}})Y_{k}({\bm{X}})^{T}=\sum_{k>r}\sigma_{k}^{2}N(d,k){\bm{I}}_{m}+\Delta_{d}=\left(\sum_{k>r}\hat{h}_{k}^{2}\right){\bm{I}}_{m}+\Delta_{d}\,.$
(4)

On the other hand, when $k<r$, $Y_{k}({\bm{X}})$ is a $m\times N(d,k)$ “tall"
matrix with $N(d,k)/m=O(d^{k}/m)=O(d^{-(r-k)})\to 0$. Similarly, Mei et al.
(2021) show that
$Y_{k}({\bm{X}})^{T}Y_{k}({\bm{X}})/m={\bm{I}}_{N(d,k)}+\Delta_{d}$, implying
that when restricted to the subspace spanned by low-frequency functions
$\\{Y_{kl}\\}_{k<r}$, the regressor associated to the empirical kernel
$K({\bm{X}},{\bm{X}})$ acts like a pure multiplicative scaling.

It remains to understand the critical-frequency mode
$Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$. It turns out that if
$N(d,k)/m\to\alpha\in(0,\infty)$, then the empirical spectral measure of the
random matrix $Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})/m$ converges to the
Marchenko-Pastur distribution $\mu_{\alpha}$, whose density is given by

$\displaystyle\mu_{\alpha}(t)=\left(1-\frac{1}{\alpha}\right)^{+}\delta_{0}(t)+\frac{\sqrt{(\alpha_{+}-t)(t-\alpha_{-})}}{2\pi\alpha
t}\bm{1}_{[\alpha_{-},\alpha_{+}]}(t),\,\mathrm{where}\,\,\alpha_{\pm}=(1\pm\sqrt{\alpha})^{2}.$
(5)

See Fig. 2 for visualizations of $\mu_{\alpha}$. The $r=1$ case is obvious as
$Y_{1}({\bm{X}})=c_{d}{\bm{X}}$ for some normalizing constant $c_{d}$ and it
is clear
$\frac{1}{m}Y_{1}({\bm{X}})^{T}Y_{1}({\bm{X}})=\frac{c_{d}^{2}}{m}{\bm{X}}^{T}{\bm{X}}$
converges to the Marchenko-Pastur distribution $\mu_{\alpha}$ if
$d/m\to\alpha\in(0,\infty)$ as $d\to\infty$ (Tao, 2012). Our first result show
that this result continues to hold for all degrees.

Figure 2: Marchenko-Pastur Distribution of Spherical Harmonics. Top: the
empirical distribution of product kernels
$Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}/N(d,r)$ vs theory prediction from
$\mu_{\alpha}$ for various degrees $r$, input dimensions $d$ and number of
samples $m$ as indicated in the titles. Bottom: the empirical distribution of
the CNN kernel $Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}/pN(d,r)$ vs theoretical
prediction. We fix $r,d$ and $m$ but varying the number of patches
$p\in\\{6,10,20\\}$.

######  Theorem 1.

For fixed $r\in{\mathbb{N}}$ and $\alpha\in(0,\infty)$, if
$N(d,r)/m\to\alpha\in(0,\infty)$ as $d\to\infty$, then the empirical spectral
distribution of $\frac{1}{m}Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$ converges in
distribution to the Marchenko-Pastur distribution $\mu_{\alpha}$ .

In the top panel of Fig. 2, we generate the empirical spectra111In the plot,
we generate the spectra of the kernel matrix
$Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}$ instead of the covariate matrix
$Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$. Although both of them have the same set
of non-zero eigenvalues, the former can be easily implemented via Legendre
polynomials $P_{r}({\bm{x}}^{T}{\bm{x}}^{\prime})$. of
$\frac{1}{m}Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$ for various values of $r$,
$d$, and $\alpha$. The Marchenko-Pastur distribution $\mu_{\alpha}$ perfectly
captures the empirical measures of the random matrices
$\frac{1}{m}Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$ for all $r$ considered. We
sketch the main steps of the proof of the theorem below; see Sec. B for the
whole proof.

###### Sketch of Proof.

From Bai and Zhou (2008, Theorem 1.1), it suffices to prove concentration of
the following quadratic forms: for every sequence of $N(d,k)\times N(d,k)$
matrices $\\{{\bm{A}}_{d}\\}$ with operator norm
$\|{\bm{A}}_{d}\|_{\mathrm{op}}\leq 1$, the variance

$\displaystyle
N(d,r)^{-2}{\mathbb{V}}(Y_{r}({\bm{x}})^{T}{\bm{A}}_{d}Y_{r}({\bm{x}})-\mathrm{Tr}({\bm{A}}_{d}))\to
0\quad\mathrm{as}\quad d\to\infty\,.$ (6)

For the purpose of illustration, we assume ${\bm{A}}\equiv{\bm{A}}_{d}$ is a
diagonal matrix. Then we only need to show

$\displaystyle
N(d,k)^{-2}\sum_{l,l^{\prime}\in[N(d,k)]}A_{ll}A_{l^{\prime}l^{\prime}}(\mathbb{E}_{\bm{x}}Y_{k,l}^{2}({\bm{x}})Y_{k,l^{\prime}}^{2}({\bm{x}})-1)\to
0.$ (7)

By hypercontractivity of spherical harmonics (Beckner, 1992),

$\displaystyle\mathbb{E}_{\bm{x}}Y_{k,l}^{2}({\bm{x}})Y_{k,l^{\prime}}^{2}({\bm{x}})\leq\left(\mathbb{E}_{\bm{x}}Y_{k,l}({\bm{x}})^{4}\mathbb{E}_{\bm{x}}Y_{k,l^{\prime}}^{4}({\bm{x}})\right)^{1/2}\leq
C_{k}\mathbb{E}_{\bm{x}}Y_{k,l}({\bm{x}})^{2}\mathbb{E}_{\bm{x}}Y_{k,l^{\prime}}({\bm{x}})^{2}=C_{k},$
(8)

where $C_{k}$ is some absolute constant. Since
$|A_{ll}|\leq\|{\bm{A}}\|_{\mathrm{op}}\leq 1$, we can drop any
$o(N(d,l)^{2})$ pairs of $(l,l^{\prime})$ in Eq. (7). We show that for the
remaining pairs $(l,l^{\prime})$, the eigenfunctions are asymptotically
uncorrelated in the sense

$\displaystyle\mathbb{E}_{\bm{x}}Y_{k,l}^{2}({\bm{x}})Y_{k,l^{\prime}}^{2}({\bm{x}})=\mathbb{E}_{\bm{x}}Y_{k,l}^{2}({\bm{x}})\mathbb{E}_{\bm{x}}Y_{k,l^{\prime}}^{2}({\bm{x}})+O(d^{-1})=1+O(d^{-1})$
(9)

which implies Eq. (7). ∎

##  4 Generalization Error of Dot-Product Kernel Regression

In this section, we establish the average generalization error for the kernel
regression in the asymptotic regime $N(d,r)/m\to\alpha$, for some
$\alpha\in(0,\infty)$ and $r\geq 1$ fixed. We assume the label function $f\in
L^{2}(\mathbb{S}_{d-1})$ is given by

$\displaystyle f({\bm{x}})=\sum_{k\geq
1}\sum_{l\in[N(d,k)]}\hat{f}_{kl}Y_{kl}({\bm{x}})=\sum_{k\geq
1}\hat{\bm{f}}_{k}^{T}Y_{k}({\bm{x}})\,,$ (10)

where $\hat{f}_{kl}$ are the “Fourier" coefficients and
$\hat{\bm{f}}_{k}=[f_{kl}]_{l\in N(d,k)}^{T}$. We need to make a technical
assumption that for $k^{\prime},k\geq r$

$\displaystyle\mathbb{E}\hat{\bm{f}}_{k}=\bm{0},\,\quad\quad\mathbb{E}\hat{\bm{f}}_{k}\hat{\bm{f}}_{k}=\frac{\hat{F}_{k}^{2}}{N(d,k)}{\bm{I}}_{N(d,k)}\quad\mathrm{and}\quad\mathbb{E}\hat{\bm{f}}_{k}\hat{\bm{f}}_{k^{\prime}}=\bm{0}_{N(d,k)\times
N(d,k^{\prime})}$ (11)

i.e. $\hat{\bm{f}}_{k}$ is centered with isotropic covariance and
$\\{\hat{\bm{f}}_{k}\\}_{k\geq r}$ are mutually uncorrelated. Note that we
allow $\hat{f}_{kl}$ to be deterministic for $k<r$. We let
${\bm{F}}=(\hat{F}_{k})_{\geq 1}$ be a fixed sequence with $\sum_{k\geq
1}\hat{F}_{k}^{2}<\infty$, where
$\hat{F}_{k}^{2}=\sum_{l\in[N(d,k)]}\hat{f}_{kl}^{2}$ for $k<r$. For
convenience, set $\hat{F}_{>j}^{2}=\sum_{k>j}\hat{F}_{k}^{2}$ (similarly for
$\hat{F}^{2}_{\geq j}$, $\hat{F}^{2}_{\leq j}$, etc.) and use ${\bm{f}}$ to
denote the random vector $\\{\hat{f}_{kl}\\}_{kl}$. Given training inputs
${\bm{X}}$ and observed labels ${\bm{Y}}=f({\bm{X}})+{\bm{\epsilon}}$, where
$\bm{\epsilon}\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}{\bm{I}}_{m})$ is
the noise, the prediction using kernel function $K$ is given by

$\displaystyle
y({\bm{x}})=K({\bm{x}},{\bm{X}})(K({\bm{X}},{\bm{X}})+\lambda{\bm{I}}_{m})^{-1}(f({\bm{X}})+\bm{\epsilon})\,.$
(12)

Here $\lambda\geq 0$ is the regularization. As such, the mean test error over
the random labels is given by

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},\hat{\bm{h}})=\mathbb{E}_{{\bm{f}}}\mathrm{Err}({\bm{X}};\lambda,{\bm{f}},\hat{\bm{h}})\quad\mathrm{where}\quad\mathrm{Err}({\bm{X}};\lambda,{\bm{f}},\hat{\bm{h}})=\mathbb{E}_{{\bm{x}},\bm{\epsilon}}|y({\bm{x}})-f({\bm{x}})|^{2}\,.$
(13)

To state our results, we need to introduce two functions $\chi_{B}$ and
$\chi_{V}$ which are related to the bias and variance in the generalization
error,

$\displaystyle\chi_{B}(\alpha,\xi)=\int(1+\xi
t)^{-2}\mu_{\alpha}(t)dt\quad\mathrm{and}\quad\chi_{V}(\alpha,\xi)=\alpha\xi^{2}\int{t}(1+\xi
t)^{-2}\mu_{\alpha}(t)dt\,.$ (14)

Both $\chi_{B}$ and $\chi_{V}$ have closed-form representations; see Sec C.5.
Define the effective regularization associated to the $r$-th order scaling to
be

$\displaystyle\xi_{r}(\hat{\bm{h}},\lambda,\alpha)=\frac{\hat{h}_{r}^{2}}{\alpha(\lambda+\hat{h}_{>r}^{2})}$
(15)

Finally, we define the bias and variance associated to the $r$-th order
scaling to be

$\displaystyle\mathrm{B}_{r}(\alpha)=\mathrm{B}_{r}(\alpha;\lambda,{\bm{F}},\hat{\bm{h}})=$
$\displaystyle\chi_{B}(\alpha,\xi_{r}(\hat{\bm{h}},\lambda,\alpha))\hat{F}_{r}^{2}+\hat{F}_{>r}^{2}$
(16)
$\displaystyle\mathrm{V}_{r}(\alpha)=\mathrm{V}_{r}(\alpha;\lambda,{\bm{F}},\hat{\bm{h}})=$
$\displaystyle\chi_{V}\left(\alpha,\xi_{r}(\hat{\bm{h}},\lambda,\alpha)\right)\left(\hat{F}_{>r}^{2}+\sigma_{\epsilon}^{2}\right)$
(17)

The following is our main result, which characterizes the test error in the
asymptotic regime $m\propto d^{r}$.

######  Theorem 2.

Let $\alpha\in(0,\infty)$ and $r\geq 1$ be fixed. Assume $N(d,r)/m\to\alpha$
as $d\to\infty$. Then the average test error is given by

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},\hat{\bm{h}})=\mathrm{B}_{r}(\alpha;\lambda,{\bm{F}},\hat{\bm{h}})+\mathrm{V}_{r}(\alpha;\lambda,{\bm{F}},\hat{\bm{h}})+\Delta_{d},$
(18)

where $\Delta_{d}\to 0$ in probability.

###  4.1 Interpretations

We provide some high-level interpretations of the bias term $\mathrm{B_{r}}$
and variance term $\mathrm{V_{r}}$.

#### The Bias.

From Eq. (16), the regressor learns all low-frequency modes ($k<r$) but none
of the high-frequency modes ($k>r$) as the bias $\mathrm{B}_{r}$ contains no
low-frequency modes (i.e. $k<r$) but all high-frequency modes
$\hat{F}_{>r}^{2}$. Importantly, the regressor is progressively learning the
critical-frequency mode $Y_{r}$ as the training size
$m=\frac{1}{\alpha}N(d,r)$ increases, i.e. from $\alpha=\infty$ to
$\alpha=0^{+}$ since $\chi_{B}(\alpha,\xi_{r}(\hat{\bm{h}},\lambda,\alpha))\to
1$ if $\alpha\to\infty$ and
$\chi_{B}(\alpha,\xi_{r}(\hat{\bm{h}},\lambda,\alpha))\to 0$ if $\alpha\to
0^{+}$. See Fig.3 for the illustration.

Figure 3: Multi-scale Bias-Variance Decomposition. Theoretical predictions of
the bias and variance from Eq.16 and Eq.17. For each $r$, the variance is non-
monotonic and has a peak at $N(d,\leq r)=\sum_{k\leq r}N(d,k)$.

#### The Variance.

From Eq. (17), the variance term $\chi_{\mathrm{V}}$ treats all high-frequency
modes $\hat{F}_{>r}^{2}$ the same as the noise term $\bm{\epsilon}$. Moreover,
$\chi_{\mathrm{V}}\to 0$ as $\alpha\to 0$ or $\infty$ and is peaked at
$\alpha=1$. The height of the peak depends on the effective regularization
$\xi_{r}$ and it diverges to infinity with rate $\xi_{r}^{\frac{1}{2}}$ as
$\xi_{r}\to\infty$. Indeed, when $\alpha=1$ and $\xi^{-1}/2\leq
t\leq\xi^{-1}$, we have $t(1+\xi t)^{-2}\mu_{\alpha}(t)\propto\xi^{-1/2}$
which implies $\chi_{V}(1,\xi)\geq\xi^{2}\int t(1+\xi
t)^{-2}\mu_{\alpha}(t)\bm{1}_{\xi^{-1}/2\leq
t\leq\xi^{-1}}dt\propto\xi^{1/2}$.

Finally, Eq. (18) not only gives precise generalization formula (up to a
vanishing term $\Delta_{d}$) when $m\approx N(d,r)\sim d^{r}$ but also when
$d^{r-1+\delta}\lesssim m\lesssim d^{r-\delta}$ (i.e. when “$\alpha=\infty$")
and when $d^{r+\delta}\lesssim m\lesssim m^{r+1-\delta}$ (i.e. when
“$\alpha=0^{+}$") for any $\delta\in(0,1/2)$. Indeed, in the non-critical
scaling regime $d^{r-1+\delta}\lesssim m\lesssim d^{r-\delta}$,
$\alpha=N(d,k)/m\to\infty$ as $d\to\infty$ and

$\displaystyle\mathrm{B_{r}}(\alpha=\infty)+\mathrm{V_{r}}(\alpha=\infty)=\hat{F}_{\geq
r}^{2}+0=\hat{F}_{\geq r}^{2}\,.$ (19)

As such, the regressor learns all low-frequency modes but none of the
critical- and high-frequency modes ($k\geq r$), which is consistent with the
result in (Ghorbani et al., 2021; Mei et al., 2021). A similar argument also
shows
$\mathrm{B_{r}}(\alpha=0^{+})+\mathrm{V_{r}}(\alpha=0^{+})=\hat{F}_{>r}^{2}$,
namely, the regressor also learns the $r$-frequency mode. This observation
implies that we can glue together Eq. (18) for $r\geq 1$ and remove all
duplicate terms to generate a sample-wise learning curve (LC):

$\displaystyle\mathrm{LC}(m;\lambda,{\bm{f}},\hat{\bm{h}})=\sum_{r\geq
1}\left(\mathrm{B}_{r}\left(\frac{N(d,r)}{m};\lambda,{\bm{f}},\hat{\bm{h}}\right)-(r-1)\hat{F}_{r}^{2}\right)+\mathrm{V}_{r}\left(\frac{N(d,\leq
r)}{m};\lambda,{\bm{f}},\hat{\bm{h}}\right)$ (20)

where $N(d,\leq r)=\sum_{k=1}^{r}N(d,k)$. The “$-(r-1)\hat{F}_{r}^{2}$" term
in the above equation is due to the fact that $\hat{F}_{r}^{2}$ is over-
counted $(r-1)$ many times (one in each $\mathrm{B}_{k}$ for $k=1,...,(r-1)$.)
It is worth mentioning that using $\alpha=N(d,\leq r)/{m}$ rather than
$\alpha={N(d,r)}/{m}$ in the variance $\mathrm{V}_{r}$ captures the finite-
size correction more accurately. See Eq. (207).

######  Corollary 1.

If, for $1\leq r\in{\mathbb{N}}$, (1) $N(d,r)/m\to\alpha$ for some
$\alpha\in(0,\infty)$, or (2) $d^{r-1+\delta}\lesssim m\lesssim d^{r-\delta}$
for some $\delta\in(0,1/2)$, then

$\displaystyle\mathrm{Err}({\bm{X}};{\bm{f}},\hat{\bm{h}},\lambda)=\mathrm{LC}(m;\lambda,{\bm{f}},\hat{\bm{h}})+\Delta_{d}$
(21)

where $\Delta_{d}\to 0$ in probability as $d\to\infty$.

Recall that for each $r$, the variance term $\mathrm{V}_{r}$ could diverge to
infinity as $\xi_{r}\to\infty$ at $\alpha=1$. Thus we might expect a peak in
the learning curve for each $r$, yielding the multiple-descent phenomenon, as
shown in Fig.1. However, such phenomena can disappear by making the heights of
the peaks small via choosing $\xi_{r}$ small. We will discuss this point in
the experimental section.

###  4.2 Proof Sketch

The proof of this theorem is quite involved; see Sec.C. For simplicity, we
assume the observed labels are noiseless, i.e. $\sigma_{\epsilon}^{2}=0$. An
ingredient is to understand the structure of the operator

$\displaystyle{\bm{T}}_{K}f({\bm{x}})=K({\bm{x}},{\bm{X}})(K({\bm{X}},{\bm{X}})+\lambda{\bm{I}}_{m})^{-1}f({\bm{X}})\,.$
(22)

The high-level strategy is as follows. We decompose the function into low-,
critical- and high-frequency parts $f=f_{<r}+f_{r}+f_{>r}$. As such, the test
error is roughly

$\displaystyle\mathrm{Err}({\bm{X}})\approx\mathrm{Err}_{<r}({\bm{X}})+\mathrm{Err}_{r}({\bm{X}})+\mathrm{Err}_{>r}({\bm{X}})\,\,\mathrm{where}\,\,\mathrm{Err}_{r}({\bm{X}})=\mathbb{E}_{\bm{f}}\mathbb{E}_{\bm{x}}|{\bm{T}}_{K}f_{r}({\bm{x}})-f_{r}({\bm{x}})|^{2}\,,$

and similarly for $\mathrm{Err}_{<r}({\bm{X}})$ and
$\mathrm{Err}_{>r}({\bm{X}})$. The next step is to estimate each part
separately. Low-frequencies. Using the fact that the low-frequency parts of
the kernel function $K$ is almost an isometric operator on the column space of
$Y_{<r}({\bm{X}})$, one can show that
$\mathbb{E}_{\bm{x}}|T_{K}(f_{<r})({\bm{x}})-f_{<r}({\bm{x}})|^{2}=\Delta_{d}\to
0$ in probability, pointwisely. Critical-frequency. Up to a vanishing term,
one can remove all non-critical frequencies in the kernel function $K$ in
$T_{K}$ in the sense of making the following substitutions

$\displaystyle
K({\bm{x}},{\bm{X}})\to\sigma_{r}^{2}Y_{r}({\bm{x}})^{T}Y_{r}({\bm{X}})^{T}\quad\mathrm{and}\quad
K_{r}({\bm{X}},{\bm{X}}^{\prime})\to\sigma_{r}^{2}Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}+\hat{h}^{2}_{>r}{\bm{I}}_{m}\,.$
(23)

Thus, with $\gamma=(\lambda+\hat{h}_{>r}^{2})$,
$T_{K}f_{r}({\bm{x}})-f_{r}({\bm{x}})=Y_{r}({\bm{x}})^{T}{\bm{M}}_{r}({\bm{X}}){\bm{f}}_{r}+\Delta_{d}$,
where

$\displaystyle{\bm{M}}_{r}({\bm{X}})=\left({\bm{I}}_{N(d,r)}-\sigma_{r}^{2}Y_{r}({\bm{X}})^{T}(\sigma_{r}^{2}Y_{r}({\bm{X}})Y_{r}({\bm{X}})^{T}+\gamma{\bm{I}}_{m})^{-1}Y_{r}({\bm{X}})\right)\,.$
(24)

Taking expectation with respect to ${\bm{x}}$ (using orthogonality of
$Y_{r}({\bm{x}})$) and then with respect to ${\bm{f}}_{r}$,

$\displaystyle\mathbb{E}_{{\bm{f}}_{r}}\mathbb{E}_{\bm{x}}|Y_{r}({\bm{x}})^{T}{\bm{M}}_{r}({\bm{X}}){\bm{f}}_{r}|^{2}=\mathbb{E}_{{\bm{f}}_{r}}|{\bm{M}}_{r}({\bm{X}}){\bm{f}}_{r}|^{2}=\hat{F}_{r}^{2}\mathrm{Tr}({\bm{M}}_{r}^{2})/N(d,r)\,.$
(25)

Applying the Sherman–Morrison–Woodbury formula and then Theorem 1,

$\displaystyle\hat{F}_{r}^{2}\mathrm{Tr}\left(Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})/m+m\sigma_{r}^{2}/\gamma{\bm{I}}_{N(d,r)}\right)^{-2}/N(d,r)\to\hat{F}_{r}^{2}\int\frac{\mu_{\alpha}(t)}{(t+\xi_{r})^{2}}dt\,.$

High-frequencies. The cross term
$\mathbb{E}_{\bm{f}}\mathbb{E}_{\bm{x}}{\bm{T}}_{K}f_{>r}({\bm{x}})f_{>r}({\bm{x}})=\Delta_{d}$
and thus

$\displaystyle\mathbb{E}_{{\bm{f}},{\bm{x}}}|{\bm{T}}_{K}f_{>r}({\bm{x}})-f_{>r}({\bm{x}})|^{2}=\mathbb{E}_{{\bm{f}},{\bm{x}}}|{\bm{T}}_{K}f_{>r}({\bm{x}})|^{2}+\mathbb{E}_{{\bm{f}},{\bm{x}}}|f_{>r}({\bm{x}})|^{2}+\Delta_{d}\,.$
(26)

The second term is equal to $\hat{F}_{>}^{2}$. The calculation of the first
term is similar to that of the critical frequency above (namely, we remove all
high-/low-frequency components in $K$.)

##  5 Convolutional Kernels

###  5.1 One hidden layer

Our analysis can be extended to analyzing NNGP kernel and NT kernel for one-
layer convolution (Novak et al., 2018, 2019; Yang, 2020). In this case, we
assume the input space is $\bm{{\mathcal{X}}}={\mathbb{S}}_{d_{0}-1}^{p}$,
where $d_{0}$ is the dimension of a patch, $p$ is the number of patches, and
$d=pd_{0}$ is the total dimensions of the inputs. The measure associated to
$\bm{{\mathcal{X}}}$ is the product of the uniform measure on
${\mathbb{S}}_{d_{0}-1}$. We assume that both the filter size and stride of
the convolution are equal to $d_{0}$. As such, after the first convolutional
layer, the input is reduced to a vector of dimension $p$. We then apply a non-
linearity and a dense layer to map this $p$-dimensional vector to a scalar.
The NNGP and NT kernel have the following general form. Let
${\bm{x}}=({\bm{x}}_{i})_{i\in[p]}\in\bm{{\mathcal{X}}}$, where
${\bm{x}}_{i}\in{\mathbb{S}}_{d_{0}-1}$ is the $i$-th patch

$\displaystyle
K({\bm{x}},{\bm{x}}^{\prime})=\frac{1}{p}\sum_{i\in[p]}h({\bm{x}}_{i}^{T}{\bm{x}}_{i}^{\prime})=\frac{1}{p}\sum_{i\in[p]}\sum_{k\geq
1}\hat{h}_{k}^{2}P_{k}({\bm{x}}_{i}^{T}{\bm{x}}_{i}^{\prime})=\sum_{k\geq
1}\frac{\sigma_{k}^{2}}{p}\sum_{i\in[p]}\sum_{l\in[N(d_{0},k)]}Y_{kl}({\bm{x}}_{i})Y_{kl}({\bm{x}}_{i}^{\prime})\,.$
(27)

Denote $Y_{kl}^{(i)}({\bm{x}})=Y_{kl}({\bm{x}}_{i})$ and
$Y_{k}({\bm{x}})=[Y_{kl}^{(i)}({\bm{x}})^{T}]_{l\in[N(d_{0},k)],i\in[p]}^{T}$.
Then $Y_{k}({\bm{x}})$ is the degree $k$ spherical harmonics associated to
this kernel, which span a space of dimension $pN(d_{0},k)$.

######  Theorem 3.

Let $r\in{\mathbb{N}}^{*}$ and $\alpha\in(0,\infty)$ be fixed. If
$pN(d_{0},r)/m\to\alpha\in(0,\infty)$ as $d_{0}\to\infty$ and the rows of
${\bm{X}}$ are sampled uniformly, iid from ${\mathbb{S}}_{d_{0}-1}^{p}$, then
the empirical spectral distribution of
$\frac{1}{m}Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})$ tends to the Marchenko-Pastur
distribution $\mu_{\alpha}$ as $d_{0}\to\infty$.

The assumptions on the label function are similar to that of dot-product
kernel, e.g.222The Gaussian assumption is unessential. We use it here for
convenience.

$\displaystyle f({\bm{x}})=\sum_{k\geq
1}{\bm{f}}_{k}^{T}Y_{k}({\bm{x}}),\quad\mathrm{with}\quad{\bm{f}}_{k}\sim{\mathcal{N}}\left(0,\frac{\hat{F}_{k}^{2}}{pN(d_{0},k)}{\bm{I}}_{pN(d_{0},k)}\right)\,\,\mathrm{if}\,k\geq
r\,,$ (28)

otherwise ${\bm{f}}_{k}$ is deterministic with
$\|{\bm{f}}_{k}\|_{2}^{2}=\hat{F}_{k}^{2}$.

######  Theorem 4.

Let $\alpha\in(0,\infty)$ and $r\geq 1$ be fixed. Assume
$pN(d_{0},r)/m\to\alpha$ as $d_{0}\to\infty$. Then the average test error is
given by

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},\hat{\bm{h}})=\mathrm{B}_{r}(\alpha;\lambda,{\bm{f}},\hat{\bm{h}})+\mathrm{V}_{r}(\alpha;\lambda,{\bm{f}},\hat{\bm{h}})+\Delta_{d_{0}},$
(29)

where $\Delta_{d_{0}}\to 0$ in probability as $d_{0}\to\infty$.

######  Corollary 2.

If, for $1\leq r\in{\mathbb{N}}$, (1) $pN(d_{0},r)/m\to\alpha$ for some
$\alpha\in(0,\infty)$, or (2) $pd_{0}^{r-1+\delta}\lesssim m\lesssim
pd_{0}^{r-\delta}$ for some $\delta\in(0,1/2)$, then

$\displaystyle\mathrm{Err}({\bm{X}};{\bm{f}},\hat{\bm{h}},\lambda)=\mathrm{LC}(m;\lambda,{\bm{f}},\hat{\bm{h}})+\Delta_{d_{0}}$
(30)

where $\Delta_{d_{0}}\to 0$ in probability as $d_{0}\to\infty$.

###  5.2 Deep Convolutional Kernels

The eigenstructure of general CNN kernels are much more complicated as they
depend on both the frequencies (i.e. the order of the polynomials) and the
topologies of the networks (Xiao, 2021). To rigorously describe the
eigenstructure, a heavy dose of notation must be introduced, which is beyond
the scope of the paper. Nevertheless, the approach developed here is readily
extended to cover general CNN kernels. We briefly describe the main ideas.

Following (Xiao, 2021), we assume the input space is still
$\bm{{\mathcal{X}}}={\mathbb{S}}_{d_{0}-1}^{p}$, where $p$ is the number of
patches. For simplicity, we assume the network has $L$ convolutional layers
and in each layer, the filter size and the stride are all equal to $d_{0}$.
Thus the spatial dimension of the input is reduced to 1 after $L$
convolutional layers. We then add a non-linearality and a dense layer to
generate the logits. The kernel has the following form

$\displaystyle
K({\bm{x}},{\bm{x}}^{\prime})=\sum_{{\bm{k}}\in{\mathbb{N}}^{p}}\sum_{{\bm{l}}\in\prod_{i\in[p]}[N(d_{0},k_{i})]}\sigma^{2}_{{\bm{k}},{\bm{l}}}Y_{{\bm{k}},{\bm{l}}}({\bm{x}})Y_{{\bm{k}},{\bm{l}}}({\bm{x}}^{\prime}),\,\mathrm{where}\quad
Y_{{\bm{k}},{\bm{l}}}({\bm{x}})=\prod_{i\in[p]}Y_{k_{i},l_{i}}({\bm{x}}_{i})\,.$
(31)

Unlike dot-product kernels in which ${\bm{k}}$ is a scalar and the eigenvalues
depend only on $|{\bm{k}}|$ (i.e. the frequencies),
$\sigma^{2}_{{\bm{k}},{\bm{l}}}$ depends on both $|{\bm{k}}|$ and the spatial
structure of the vector ${\bm{k}}$ in a rather complicated manner.
Nevertheless, as $d_{0}\to\infty$, $\sigma^{2}_{{\bm{k}},{\bm{l}}}\sim
d_{0}^{-j_{{\bm{k}}}}=d^{-j_{\bm{k}}/L}$ for some $L\leq
j_{\bm{k}}\in{\mathbb{N}}$. We can then categorize the eigenvectors according
to the decay order of $\sigma^{2}_{{\bm{k}},{\bm{l}}}$. Unlike the case of
dot-product kernels or the one-hidden layer CNN kernels, in which eigenvectors
with same-order eigenvalues are in the same eigenspace (i.e. the eigenvalues
are the same), multiple-layer CNN kernels can have multiple eigenspaces with
the same-order eigenvalues. Although this results in extra challenges (see
below), our overall approach carries over. Consider the critical scaling
regime $m\sim d^{r}$, for $r=j/L$ for some $L\leq j\in{\mathbb{N}}$. Likewise,
we can decompose the kernel into low-, critical- and high-frequency parts
according to $j_{\bm{k}}<r$, $j_{\bm{k}}=r$ and $j_{\bm{k}}>r$, resp.
Following similar assumptions on the labels and eigenvalues, the bias and the
variance can be essentially reduced to computing

$\displaystyle\chi_{B}=$
$\displaystyle\frac{1}{N_{r}}\operatorname{Tr}{\bm{R}}_{r}^{2}({\bm{R}}_{r}+Y_{r}({\bm{X}})^{T}Y_{r}({\bm{X}})/m)^{-2}$
(32) $\displaystyle\chi_{V}=$ $\displaystyle\frac{N_{\leq
r}}{m}\frac{1}{N_{\leq r}}\operatorname{Tr}({\bm{R}}_{\leq r}+Y_{\leq
r}({\bm{X}})^{T}Y_{\leq r}({\bm{X}})/m)^{-2}Y_{\leq r}({\bm{X}})^{T}Y_{\leq
r}({\bm{X}})/m$ (33)

where ${\bm{R}}_{r}$ and ${\bm{R}}_{\leq r}$ are diagonal matrices whose
entries are determined by the eigenvalues of the critical-frequency modes. In
the dot-product kernels or one-hidden layer CNN kernels setting,
${\bm{R}}_{r}$/${\bm{R}}_{\leq r}$ is a scaled identity matrix and simple,
closed-form expressions for the above traces straightforwardly follow from the
Marcenko-Pastur distribution. However, for a general diagonal matrix
${\bm{R}}_{r}$ with bounded limiting spectra, ${\bm{R}}_{r}$ does not commute
with $Y_{r}^{T}({\bm{X}})Y_{r}({\bm{X}})$, and a more detailed random matrix
analysis is needed. See the supplementary material for more details.

##  6 Experiments

We provide experiments to show that our learning curves (Eq. (21)) accurately
capture empirical sample-wise learning curve even when the ambient dimensions
remains small. Even though our theoretical results require averaging the test
error over random labels (aka, mean test error), our experimental results
suggest this is unnecessary, i.e. the learning curve Eq. (21) can capture the
test error accurately for any given draw of label function.

#### Experimental setup.

We generate a polynomial kernel function
$h(t)=\sum_{k=1}^{7}\hat{h}_{k}^{2}P_{k}(t)$, where $P_{k}$ is the degree-$k$
Legendre polynomial in $d$ dimensions. The kernel function can be efficiently
computed via $K({\bm{x}},{\bm{x}}^{\prime})=h({\bm{x}}^{T}{\bm{x}}^{\prime})$.
We choose the label function to be
$f({\bm{x}})=\sum_{k=1}^{7}\hat{F}_{k}y_{k}({\bm{x}})$, where
$y_{k}({\bm{x}})=\sum_{j\in[d]}w_{k,j}\prod_{i=j}^{j+k-1}{x_{i}}$, and the
coefficients $w_{k,j}$ are randomly sampled from a Gaussian and then
normalized so that $\mathbb{E}_{\bm{x}}|y_{k}({\bm{x}})|^{2}=1$ for each $k$.
Therefore
$\mathbb{E}_{\bm{x}}|f({\bm{x}})|^{2}=\sum_{k=1}^{7}\hat{F}_{k}^{2}$. For
simplicity, we also set $\sigma_{\epsilon}^{2}=0$ (i.e. noiseless) and
$\lambda=0$ (i.e. ridgeless). Note that when $m\lesssim d^{r}$, the regressor
still contains “effective noise" $\hat{F}_{>r}^{2}$ from un-learnable high-
frequency modes and “effective regularization" $\hat{h}_{>r}^{2}$. Finally, in
our experiments, we choose $\hat{F}_{k}^{2}=k^{-2}$ and
$\hat{h}_{k}^{2}=\mathrm{Gap}^{-(k-1)}$, where we will vary the value of the
spectral gap: $\mathrm{Gap}=\hat{h}_{k}^{2}/\hat{h}_{k+1}^{2}$. Under this
setup, the predicted learning curve
$\mathrm{LC}(m)=\mathrm{LC}(m;\mathrm{Gap})$ depends only on the spectral gap
of the kernel.

To simulate higher-order scaling ($r\geq 3$), the dimension $d$ has to be very
small as we need to invert a sequence of matrices of size ranging from $m=1$
to $m\propto d^{r}/r!$. Due to the constraints in compute and memory, the
largest $m$ we can have is typically $m_{\mathrm{max}}\approx 25,000$ for one
single GPU and $d$ in our experiments is typically around $d=24$. As such we
are in a regime with strong finite-size corrections. Finally, all experiments
are run in a single A-100 using Google Cloud Colab Notebooks.

#### Learning Curves Accurately Capture Simulations.

In Fig. 4, we generate the empirical sample-wise learning curve by applying
kernel regression Eq. (13) with training set ${\bm{X}}$. We vary the training
set size $m$ densely in $[1,m_{\mathrm{max}}]$ and for each $m$ we sample 20
independent ${\bm{X}}$ to get the errorbar plot for the test error. The
closed-form learning curve is obtained from Eq. (21) and the calculation is
done in Sec.C.5. Even in the low-dimensional regime with $d=24$ for dot-
product kernel ($d_{0}=20$ and $p=6$ for one-hidden layer CNN kernel), the
predicted learning curve captures the empirical learning curve surprisingly
well, which has a highly non-trivial multiple-descent behavior. It is worth
mentioning that, from the simulation, the deviation of the test error from its
mean is relatively large when $m$ is small but vanishes quickly as $m$ becomes
larger. This suggests Theorem 2 and Corollary 1 should hold in a pointwise
fashion, i.e. without averaging the test error over random labels.

(a) Dot Product Kernel

(b) One-layer CNN Kernel

Figure 4: Simulation vs Prediction. We generate the learning curves obtain
from  kernel regression by densely varying $m$ from 1 to $24000$. For each
$m$, we average the MSE over 20 runs. The closed-form prediction from Eq. (18)
captures the simulations surprising well even for small $d$. Left: dot product
kernel with $d=24$. Right: one-hidden layer CNN kernel with $d_{0}=20$ and
$p=6$. The spectral gap is $\mathrm{Gap}=32$ in both plots.

#### Small Spectral Gap Eliminates Multiple-descents.

In Fig. 1, we plot both the predicted learning curves and simulations when
$\mathrm{Gap}$ ranging in $[2,8,32,128]$. For $1\leq r\leq 6$, we have
$\xi_{r}=\mathrm{Gap}^{-(r-1)}/\sum_{k=r}^{7}\mathrm{Gap}^{-(k-1)}$ and
$\xi_{r}\approx\mathrm{Gap}$ when $\mathrm{Gap}$ is large. Recall that the
variance term $\mathrm{V}_{r}$ peaks at $\alpha=1$ and the peak scales like
$\xi_{r}^{1/2}\approx\mathrm{Gap}^{1/2}$. When $\mathrm{Gap}$ is large, e.g.
$\mathrm{Gap}=32,128$, the variance is also large near $\alpha=1$, the
multiple-descent phenomena are more prominent. On the other hand, when
$\mathrm{Gap}$ is small, e.g. $\mathrm{Gap}=8,2$, such phenomena disappear and
learning curves become monotonic.

##  7 Conclusion

In this work, we establish precise asymptotic formulas for the sample-wise
learning curves in the kernel ridge regression setting for a family of dot-
product kernels in the polynomial scaling regimes $m\propto d^{r}$ for all
$r\in{\mathbb{N}}^{*}$. We demonstrate that these formulas can capture
empirical learning curves surprisingly well even in the regime where strong
finite-size corrections would be expected. We rigorously prove that the
learning curves can be non-monotonic near $m\propto d^{r}/r!$ for each
$r\in{\mathbb{N}}^{*}$. There are a couple limitations of our approach which
could be improved in future work. The first one is the strong assumption on
the distribution of the input data, namely, the uniform distribution on the
spherical type of data. In addition, the learning curves are obtained only in
the kernel regression setting and extending the results to the feature
learning setting would be a meaningful future direction.

## Acknowledgments and Disclosure of Funding

We thank Ben Adlam for providing valuable feedbacks on a draft. This work was
performed at and funded by Google.

## References

  * Adlam and Pennington [2020a] Ben Adlam and Jeffrey Pennington.  The neural tangent kernel in high dimensions: Triple descent and a multi-scale theory of generalization.  In _International Conference on Machine Learning_ , pages 74–84. PMLR, 2020a.
  * Adlam and Pennington [2020b] Ben Adlam and Jeffrey Pennington.  Understanding double descent requires a fine-grained bias-variance decomposition.  _Advances in neural information processing systems_ , 33:11022–11032, 2020b.
  * Advani et al. [2020] Madhu S Advani, Andrew M Saxe, and Haim Sompolinsky.  High-dimensional dynamics of generalization error in neural networks.  _Neural Networks_ , 132:428–446, 2020.
  * Arora et al. [2019] Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, Russ R Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Bai and Zhou [2008] Zhidong Bai and Wang Zhou.  Large sample covariance matrices without independence structures in columns.  _Statistica Sinica_ , pages 425–442, 2008.
  * Beckner [1992] William Beckner.  Sobolev inequalities, the poisson semigroup, and analysis on the sphere sn.  _Proceedings of the National Academy of Sciences_ , 89(11):4816–4819, 1992.
  * Belkin et al. [2019] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  _Proceedings of the National Academy of Sciences_ , 116(32):15849–15854, 2019.
  * Bryson et al. [2021] Jennifer Bryson, Roman Vershynin, and Hongkai Zhao.  Marchenko–pastur law with relaxed independence conditions.  _Random Matrices: Theory and Applications_ , page 2150040, 2021.
  * Chowdhery et al. [2022] Aakanksha Chowdhery, Sharan Narang, Jacob Devlin, Maarten Bosma, Gaurav Mishra, Adam Roberts, Paul Barham, Hyung Won Chung, Charles Sutton, Sebastian Gehrmann, et al.  Palm: Scaling language modeling with pathways.  _arXiv preprint arXiv:2204.02311_ , 2022.
  * d’Ascoli et al. [2020] Stéphane d’Ascoli, Maria Refinetti, Giulio Biroli, and Florent Krzakala.  Double trouble in double descent: Bias and variance (s) in the lazy regime.  In _International Conference on Machine Learning_ , pages 2280–2290. PMLR, 2020.
  * Far et al. [2006] Reza Rashidi Far, Tamer Oraby, Wlodzimierz Bryc, and Roland Speicher.  Spectra of large block matrices.  _arXiv preprint cs/0610045_ , 2006.
  * Frye and Efthimiou [2012] Christopher Frye and Costas J Efthimiou.  Spherical harmonics in p dimensions.  _arXiv preprint arXiv:1205.3548_ , 2012.
  * Ghorbani et al. [2021] Behrooz Ghorbani, Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Linearized two-layers neural networks in high dimension.  _The Annals of Statistics_ , 49(2):1029–1054, 2021.
  * Hastie et al. [2022] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  _The Annals of Statistics_ , 50(2):949–986, 2022\.
  * Hu and Lu [2022] Hong Hu and Yue M Lu.  Sharp asymptotics of kernel ridge regression beyond the linear regime.  _arXiv preprint arXiv:2205.06798_ , 2022.
  * Jacot et al. [2018] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Kaplan et al. [2020] Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models.  _arXiv preprint arXiv:2001.08361_ , 2020.
  * Lee et al. [2019] Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  _Advances in neural information processing systems_ , 32, 2019.
  * Liang et al. [2020] Tengyuan Liang, Alexander Rakhlin, and Xiyu Zhai.  On the multiple descent of minimum-norm interpolants and restricted lower isometry of kernels.  In _Conference on Learning Theory_ , pages 2683–2711. PMLR, 2020\.
  * Lin and Dobriban [2021] Licong Lin and Edgar Dobriban.  What causes the test error? going beyond bias-variance via anova.  _Journal of Machine Learning Research_ , 22(155):1–82, 2021.
  * Mania and Sra [2020] Horia Mania and Suvrit Sra.  Why do classifier accuracies show linear trends under distribution shift?  _arXiv preprint arXiv:2012.15483_ , 2020.
  * Mei and Montanari [2022] Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and the double descent curve.  _Communications on Pure and Applied Mathematics_ , 75(4):667–766, 2022.
  * Mei et al. [2021] Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Generalization error of random feature and kernel methods: hypercontractivity and kernel matrix concentration.  _Applied and Computational Harmonic Analysis_ , 2021.
  * Mingo and Speicher [2017] James A Mingo and Roland Speicher.  _Free probability and random matrices_ , volume 35.  Springer, 2017.
  * Misiakiewicz [2022] Theodor Misiakiewicz.  Spectrum of inner-product kernel matrices in the polynomial regime and multiple descent phenomenon in kernel ridge regression.  _arXiv preprint arXiv:2204.10425_ , 2022.
  * Novak et al. [2018] Roman Novak, Lechao Xiao, Jaehoon Lee, Yasaman Bahri, Greg Yang, Jiri Hron, Daniel A Abolafia, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Bayesian deep convolutional networks with many channels are gaussian processes.  _arXiv preprint arXiv:1810.05148_ , 2018.
  * Novak et al. [2019] Roman Novak, Lechao Xiao, Jiri Hron, Jaehoon Lee, Alexander A Alemi, Jascha Sohl-Dickstein, and Samuel S Schoenholz.  Neural tangents: Fast and easy infinite neural networks in python.  _arXiv preprint arXiv:1912.02803_ , 2019.
  * Tao [2012] Terence Tao.  _Topics in random matrix theory_ , volume 132.  American Mathematical Soc., 2012.
  * Tripuraneni et al. [2021a] Nilesh Tripuraneni, Ben Adlam, and Jeffrey Pennington.  Covariate shift in high-dimensional random feature regression.  _arXiv preprint arXiv:2111.08234_ , 2021a.
  * Tripuraneni et al. [2021b] Nilesh Tripuraneni, Ben Adlam, and Jeffrey Pennington.  Overparameterization improves robustness to covariate shift in high dimensions.  _Advances in Neural Information Processing Systems_ , 34, 2021b.
  * Vershynin [2010] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Xiao [2021] Lechao Xiao.  Eigenspace restructuring: a principle of space and frequency in neural networks.  _arXiv preprint arXiv:2112.05611_ , 2021.
  * Yang [2020] Greg Yang.  Tensor programs ii: Neural tangent kernel for any architecture.  _arXiv preprint arXiv:2006.14548_ , 2020.

##  Appendix A Appendix Guidelines

The appendix is organized as follows. We prove Theorem 1 and Theorem 3 in Sec.
B. In Sec. C, we prove the test errors for dot-product kernels, namely,
Theorem 2 and for one-hidden-layer convolutional kernels, namely, Theorem 4.
The proof also shows how to reduce the test error of multiple-layer
convolutional kernels to evaluating Eq. (32) and Eq. (33). Finally, in Sec.D,
we provide additional plots to empirically verify that the finite-size
correction becomes smaller as $d$ grows larger.

##  Appendix B Proof of Theorem 1

We begin with some notations. For positive numbers $a$ and $b$, we use
$a\lesssim b$ to mean there is a constant independent of $d$ such that $a\leq
Cb$. In addition, $a\sim b$, if $a\lesssim b$ and $b\lesssim a$.

The proof of Theorem 3 is similar. We only present the proof of Theorem 1. Our
proof is based on the following result from Bai and Zhou [2008].

######  Lemma 1 (Bai and Zhou [2008]).

Let ${\bm{x}}_{p}\in\mathbb{R}^{p}$ be random vectors and
${\bm{X}}=[{\bm{x}}_{p1},\dots\bm{x}_{pn}]$ be a $p\times m$ matrix with iid
columns. If for every $\\{{\bm{A}}_{p}\\}_{p}$, $p\times p$ matrix with
uniform operator norm,

$\displaystyle\frac{1}{p^{2}}\mathbb{E}|{\bm{x}}_{p}^{T}{\bm{A}}_{p}{\bm{x}}_{p}-{\rm
Tr}({\bm{A}}_{p}))|^{2}\to 0,$ (34)

then the empirical spectral distribution of $\frac{1}{m}{\bm{X}}{\bm{X}}^{T}$
converges to $\mu_{\alpha}$ weakly if $p/m\to\alpha\in(0,\infty)$.

We prove a slightly more general version.

######  Theorem 5.

Let $r\in{\mathbb{N}}$ and $\alpha\in(0,\infty)$ be fixed. Assume $m=m(d)$
with $N(d,r)/m\to\alpha\in(0,\infty)$ as $d\to\infty$. Let
${\bm{u}}={\bm{u}}({\bm{d}})$ be a sequence of functions defined on
${\mathbb{S}}_{d-1}$ such that

  1. (1.)

the cardinality of ${\bm{u}}$ satisfies $|{\bm{u}}|/d^{r}\to 0$ as
$d\to\infty$;

  2. (2.)

the functions in ${\bm{u}}$ and $Y_{r}$ are mutually orthogonal;

  3. (3.)

for any unit vector $\theta$, let
$\mathbb{E}_{\bm{x}}|\theta^{T}Z_{r}({\bm{x}})|^{4}\lesssim 1$ uniformly of
$d$ and $\theta$, where $Z_{r}({\bm{x}})$ is the concatenation of
${\bm{u}}({\bm{x}})$ and $Y_{r}({\bm{x}})$.

Let $Z_{r}({\bm{X}})$ be the concatenation of $Y_{r}({\bm{X}})$ and
${\bm{u}}({\bm{X}})$. Then the empirical spectral distribution of
$\frac{1}{m}Z_{r}({\bm{X}})^{T}Z_{r}({\bm{X}})$ converges in distribution to
the Marchenko-Pastur distribution $\mu_{\alpha}$.

We mainly use the case when ${\bm{u}}$ is the empty set, i.e. $Z_{r}=Y_{r}$
and the case when ${\bm{u}}=[Y_{kl}]_{k<r}^{T}$, i.e. $Z_{r}=Y_{\leq r}$.

###### Proof of Theorem 5.

We apply Lemma 1 to $Z_{r}^{T}({\bm{X}})$. We only need to show for matrices
${\bm{A}}={\bm{A}}^{(d)}$ with $\|{\bm{A}}\|_{op}\leq 1$,

$\displaystyle(|{\bm{u}}|+N(d,r))^{-2}\mathbb{E}_{x}|Z_{r}({\bm{x}})^{T}{\bm{A}}Z_{r}({\bm{x}})-{\rm
Tr}({\bm{A}})|^{2}\to 0\,\quad\mathrm{as}\quad d\to\infty\,,$ (35)

or, equivalently

$\displaystyle
N(d,r)^{-2}\mathbb{E}_{x}|Z_{r}({\bm{x}})^{T}{\bm{A}}Z_{r}({\bm{x}})-{\rm
Tr}({\bm{A}})|^{2}\to 0\,\quad\mathrm{as}\quad d\to\infty.$ (36)

since $|{\bm{u}}|=o(d^{r})$. The assumption $\|{\bm{A}}\|_{op}\leq 1$ implies
that the absolute values of all entries of ${\bm{A}}$ are bounded 1. A key
observation in proving the above estimate is that, up to a unitary
transformation, almost all functions in $\\{Y_{r,l}({\bm{x}})\\}_{l}$ are
monomials of the form

$\displaystyle g_{\bm{i}}({\bm{x}})=C_{d,r}\prod_{i\in{\bm{i}}}x_{i}$ (37)

where ${\bm{i}}\subseteq[d]$ with $|{\bm{i}}|=r$ and $C_{d,r}$ is a
normalizing factor such that

$\displaystyle
C_{d,r}^{2}\int_{{\mathbb{S}}_{d-1}}\prod_{i\in{\bm{i}}}|x_{i}|^{2}d{\bm{x}}=1\,.$
(38)

We prove later that for any finite integer $r\geq 1$,

$\displaystyle\int_{{\mathbb{S}}_{d-1}}\prod_{i\in{\bm{i}}}|x_{i}|^{2}d{\bm{x}}=\prod_{i\in{\bm{i}}}\int_{{\mathbb{S}}_{d-1}}|x_{i}|^{2}d{\bm{x}}+O(d^{-r-1})=d^{-r}+O(d^{-r-1})$
(39)

Now we proceed to prove Eq. (36). First note that
${\mathbb{G}}=\\{g_{\bm{i}}:{\bm{i}}\subseteq[d],|{\bm{i}}|=r\\}$ is an
orthonormal set. This can be proved by noticing that if
${\bm{i}}\neq{\bm{j}}$, then there is $i\in{\bm{i}}$ but $i\notin{\bm{j}}$.
Clearly, the symmetries of the measure on ${\mathbb{S}}_{d-1}$ implies

$\displaystyle\int_{{\mathbb{S}}_{d-1}}g_{\bm{i}}({\bm{x}})g_{\bm{j}}({\bm{x}})d{\bm{x}}=0\,.$
(40)

We choose ${\mathbb{B}}=\\{b_{j}\\}_{j\in[p]}$ so that
${\mathbb{G}}\sqcup{\mathbb{B}}$ forms an orthonormal basis of
$Z_{r}({\bm{x}})$. Note that the cardinality of ${\mathbb{B}}$ is $o(d^{r})$.
Indeed,

$\displaystyle
p=|{\bm{u}}|+N(d,r)-{{d}\choose{r}}=|{\bm{u}}|+{{d+r-2}\choose{r}}+{{d+r-3}\choose{r-1}}-{{d}\choose{r}}=|{\bm{u}}|+O(d^{r-1})=o(d^{r})$
(41)

Thus $\frac{p}{N(d,r)}\to 0$ as $d\to\infty$. After a change of basis, we can
assume $Z_{r}({\bm{x}})=[{\bm{g}}({\bm{x}})^{T},{\bm{b}}({\bm{x}})^{T}]^{T}$,
where ${\bm{g}}=[g_{\bm{i}}]_{\bm{i}}^{T}$ and
${\bm{b}}=[b_{j}]_{j\in[p]}^{T}$. Here we use the fact that Eq. (36) holds for
all ${\bm{A}}$ with uniform operator norms is equivalent to that it holds for
${\bm{Q}}^{T}{\bm{A}}{\bm{Q}}$ for all such ${\bm{A}}$ and any unitary matrix
${\bm{Q}}$. We write,

$\displaystyle{\bm{A}}=\begin{bmatrix}{\bm{A}}_{11}&{\bm{A}}_{12}\\\
{\bm{A}}_{21}&{\bm{A}}_{22}\end{bmatrix}$ (42)

where ${\bm{A}}_{11}$ is the upper left ${{d}\choose{r}}\times{{d}\choose{r}}$
block of ${\bm{A}}$, ${\bm{A}}_{22}$ is the lower right $p\times p$ block of
${\bm{A}}$ and the other two blocks are defined similarly. Note that
$\|{\bm{A}}_{ij}\|_{op}\leq\|{\bm{A}}\|_{op}\leq 1$ for $i,j\in\\{1,2\\}$. As
such, we have

$\displaystyle
N(d,r)^{-2}\mathbb{E}_{\bm{x}}|Z_{r}({\bm{x}})^{T}{\bm{A}}Z_{r}({\bm{x}})-{\rm
Tr}({\bm{A}})|^{2}\leq 5(I_{1}+I_{2}+I_{3}+I_{4}+I_{5})$ (43)

where

$\displaystyle I_{1}$
$\displaystyle=N(d,r)^{-2}\mathbb{E}_{\bm{x}}|{\bm{g}}({\bm{x}})^{T}{\bm{A}}_{11}{\bm{g}}({\bm{x}})-{\rm
Tr}({\bm{A}}_{11})|^{2}$ (44) $\displaystyle I_{2}$
$\displaystyle=N(d,r)^{-2}\mathbb{E}_{\bm{x}}|{\bm{b}}({\bm{x}})^{T}{\bm{A}}_{22}{\bm{b}}({\bm{x}})|^{2}$
(45) $\displaystyle I_{3}$
$\displaystyle=N(d,r)^{-2}\mathbb{E}_{\bm{x}}|{\bm{g}}({\bm{x}})^{T}{\bm{A}}_{12}{\bm{b}}({\bm{x}})|^{2}$
(46) $\displaystyle I_{4}$
$\displaystyle=N(d,r)^{-2}\mathbb{E}_{\bm{x}}|{\bm{b}}({\bm{x}})^{T}{\bm{A}}_{21}{\bm{g}}({\bm{x}})|^{2}$
(47) $\displaystyle I_{5}$ $\displaystyle=N(d,r)^{-2}|{\rm
Tr}({\bm{A}}_{22})|^{2}$ (48)

We prove $I_{i}\to 0$ for $1\leq i\leq 5$. The $i=1$ case is the most
difficult and the others are straightforward since $pN(d,r)^{-1}\to 0$. E.g.,
when $i=3$

$\displaystyle I_{3}$ $\displaystyle\leq
N(d,r)^{-2}\mathbb{E}_{\bm{x}}\|{\bm{g}}({\bm{x}})\|_{l_{2}}^{2}\|{\bm{A}}_{12}\|_{op}^{2}\|{\bm{b}}({\bm{x}})\|^{2}_{l_{2}}$
(49) $\displaystyle\leq
N(d,r)^{-2}\left(\mathbb{E}_{\bm{x}}\|{\bm{g}}({\bm{x}})\|_{l_{2}}^{4}\mathbb{E}_{\bm{x}}\|{\bm{b}}({\bm{x}})\|^{4}_{l_{2}}\right)^{1/2}$
(50) $\displaystyle\leq
N(d,r)^{-2}\max_{{\bm{i}},j}(\mathbb{E}_{\bm{x}}|g_{\bm{i}}({\bm{x}})|^{4}\mathbb{E}_{\bm{x}}|b_{j}({\bm{x}})|^{4})^{1/2}pN(d,r)$
(51) $\displaystyle=CpN(d,r)^{-1}\to 0\,$ (52)

where we have set
$C=(\mathbb{E}_{\bm{x}}|g_{\bm{i}}({\bm{x}})|^{4}\mathbb{E}_{\bm{x}}|b_{j}({\bm{x}})|^{4})^{1/2}$,
which is $O(1)$ due to assumption (3.) in Theorem 5. The bounds for $I_{2}$
and $I_{4}$ can be obtained similarly. For $I_{5}$, we simply use ${\rm
Tr}({\bm{A}}_{22})\leq p\|{\bm{A}}\|_{op}\leq p$.

It remains to control $I_{1}$. To ease the notation, denote
${\bm{B}}={\bm{A}}_{11}$. We split $I_{1}\leq I_{11}+I_{12}$, where $I_{11}$
and $I_{12}$ are the diagonal and the off-diagonal parts, resp.,

$\displaystyle I_{11}$
$\displaystyle=2N(d,r)^{-2}\mathbb{E}_{\bm{x}}|\sum_{{\bm{i}}}B_{{\bm{i}}{\bm{i}}}(g_{\bm{i}}^{2}({\bm{x}})-1)|^{2}$
(53) $\displaystyle I_{12}$
$\displaystyle=2N(d,r)^{-2}\mathbb{E}_{\bm{x}}|\sum_{{\bm{i}}\neq{\bm{j}}}B_{{\bm{i}}{\bm{j}}}g_{\bm{i}}g_{\bm{j}}({\bm{x}})|^{2}$
(54)

#### Bounding the diagonal part $I_{11}$.

Using $\mathbb{E}_{\bm{x}}g_{\bm{i}}({\bm{x}})^{2}=1$, we have

$\displaystyle
I_{11}=2N(d,r)^{-2}\mathbb{E}_{\bm{x}}\sum_{{\bm{i}},{\bm{j}}}B_{{\bm{i}}{\bm{i}}}B_{{\bm{j}}{\bm{j}}}(g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1)$
(55)

We spit the proof into two cases: ${\bm{i}}\cap{\bm{j}}=\emptyset$ and
${\bm{i}}\cap{\bm{j}}\neq\emptyset$. The following is the key estimate to
handle the first case.

######  Lemma 2.

If ${\bm{i}}\cap{\bm{j}}=\emptyset$,

$\displaystyle\max_{{\bm{i}},{\bm{j}}}|\mathbb{E}_{\bm{x}}g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1|\lesssim
d^{-1}\,.$ (56)

We prove this lemma later. We show how to use this lemma to handle the
${\bm{i}}\cap{\bm{j}}=\emptyset$ case. Recall that
$|B_{{\bm{i}},{\bm{j}}}|\leq 1$ and the number of tuples $({\bm{i}},{\bm{j}})$
is fewer than $N(d,r)^{2}$. We have

$\displaystyle
N(d,r)^{-2}|\mathbb{E}_{\bm{x}}\sum_{{\bm{i}},{\bm{j}},{\bm{i}}\cap{\bm{j}}=\emptyset}B_{{\bm{i}}{\bm{i}}}B_{{\bm{j}}{\bm{j}}}(g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1)|$
(57) $\displaystyle\leq$ $\displaystyle
N(d,r)^{-2}\sum_{{\bm{i}},{\bm{j}},{\bm{i}}\cap{\bm{j}}=\emptyset}\max_{{\bm{i}},{\bm{j}}}|\mathbb{E}_{\bm{x}}g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1|$
(58) $\displaystyle\leq$
$\displaystyle\max_{{\bm{i}},{\bm{j}}}|\mathbb{E}_{\bm{x}}g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1|\lesssim
d^{-1}\,.$ (59)

We turn to $|{\bm{i}}\cap{\bm{j}}|=t$, $1\leq t\leq r$. For each fixed
${\bm{i}}$, the number of choices of ${\bm{j}}$ is

$\displaystyle\sum_{1\leq t\leq r}{r\choose t}{d-r\choose
r-t}\lesssim\sum_{1\leq t\leq r}d^{r-t}\sim d^{r-1}$ (60)

As such,

$\displaystyle
2N(d,r)^{-2}\mathbb{E}_{\bm{x}}\sum_{{\bm{i}},{\bm{j}},{\bm{i}}\cap{\bm{j}}\neq\emptyset}|B_{{\bm{i}}{\bm{i}}}B_{{\bm{j}}{\bm{j}}}(g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1)|$
(61) $\displaystyle\lesssim$ $\displaystyle
N(d,r)^{-2}N(d,r)d^{r-1}\max_{{\bm{i}},{\bm{j}}}\mathbb{E}_{\bm{x}}|g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1|$
(62) $\displaystyle\lesssim$ $\displaystyle d^{-1}$ (63)

#### Off-diagonal terms $I_{12}$.

Bounding the off-diagonal terms can be reduced to a combinatorics problem,
which is similar to the random tensor model considered in Bryson et al.
[2021]. We need to estimate

$\displaystyle
I_{12}=2N(d,r)^{-2}\sum_{{\bm{i}}\neq{\bm{j}},{\bm{l}}\neq{\bm{k}}}B_{{\bm{i}}{\bm{j}}}B_{{\bm{l}}{\bm{k}}}\mathbb{E}_{\bm{x}}g_{\bm{i}}({\bm{x}})g_{\bm{j}}({\bm{x}})g_{\bm{l}}({\bm{x}})g_{\bm{k}}({\bm{x}})$
(64)

By symmetries of the uniform measure on the sphere, we can assume the monomial
$g_{\bm{i}}({\bm{x}})g_{\bm{j}}({\bm{x}})g_{\bm{l}}({\bm{x}})g_{\bm{k}}({\bm{x}})$
has no linear factor, that is the degree of any $x_{i}$ in this monomial must
be at least $2$ if not 0. In addition, for such monomials, the Holder
inequality and hypercontractivities yield,

$\displaystyle|\mathbb{E}_{\bm{x}}g_{\bm{i}}({\bm{x}})g_{\bm{j}}({\bm{x}})g_{\bm{l}}({\bm{x}})g_{\bm{k}}({\bm{x}})|\leq\max_{{\bm{i}}}\mathbb{E}_{\bm{x}}|g_{\bm{i}}({\bm{x}})|^{4}\lesssim\max_{{\bm{i}}}(\mathbb{E}_{\bm{x}}|g_{\bm{i}}({\bm{x}})|^{2})^{2}=1\,.$
(65)

As such, we only need to show that the growth of the number of such quadruples
$({\bm{i}},{\bm{j}},{\bm{k}},{\bm{l}})$, as a function of $d$, is slower than
$N(d,r)^{2}\sim d^{2r}$. We proceed to prove this claim. For each fixed
${\bm{i}}$, let $t=|{\bm{i}}\cap{\bm{j}}|$ where $0\leq t\leq r-1$ ($t\neq r$
since ${\bm{i}}\neq{\bm{j}}$). Let $J({\bm{i}};t)$ denote the set of such
${\bm{j}}$, whose cardinality is

$\displaystyle|J({\bm{i}};t)|={r\choose t}{d-r\choose
r-t}<r^{t}d^{r-t}\lesssim d^{r-t}\,.$ (66)

Next we estimate the number of tuples $({\bm{l}},{\bm{k}})$. Let
$w=|({\bm{l}}\cup{\bm{k}})\backslash({\bm{i}}\cup{\bm{j}})|$. Since
$|{\bm{l}}\cup{\bm{k}}\cup{\bm{i}}\cup{\bm{j}}|\leq 2r$ and
$|{\bm{i}}\cup{\bm{j}}|=2r-t$, we have $w\leq t$. The cardinality of choosing
such ${\bm{k}}\cup{\bm{l}}$ cannot exceed

$\displaystyle\sum_{0\leq w\leq t}{d-(2r-t)\choose
w}\sum_{v=0}^{2r-w}{2r\choose v}\lesssim d^{t}\,.$ (67)

With ${\bm{k}}\cup{\bm{l}}$ given, the pair of $({\bm{k}},{\bm{l}})$ cannot
exceed ${2r\choose r}^{2}\lesssim 1$. Thus, with ${\bm{i}},{\bm{j}}$ and $t$
fixed, the number of pairs of $({\bm{k}},{\bm{l}})$ is $\lesssim d^{t}$. Using
$|B_{{\bm{l}}{\bm{k}}}|\leq 1$ and
$\max_{{\bm{i}}}(\sum_{{\bm{j}}}B_{{\bm{i}}{\bm{j}}}^{2})^{1/2}\leq\|{\bm{B}}\|_{op}\leq
1$, we have

$\displaystyle N(d,r)^{2}I_{12}$
$\displaystyle\lesssim\sum_{{\bm{i}}}\sum_{0\leq t\leq r-1}\sum_{j\in
J({\bm{i}};t)}B_{{\bm{i}}{\bm{j}}}d^{t}$ (68)
$\displaystyle\leq\sum_{{\bm{i}}}\sum_{0\leq t\leq
r-1}|J({\bm{i}};t)|^{1/2}(\sum_{{\bm{j}}\in
J({\bm{i}};t)}B_{{\bm{i}}{\bm{j}}}^{2})^{1/2}d^{t}$ (69)
$\displaystyle\lesssim\|{\bm{A}}\|_{op}\sum_{{\bm{i}}}\sum_{0\leq t\leq
r-1}d^{(r-t)/2}d^{t}$ (70) $\displaystyle\lesssim N(d,r)d^{2r-1/2}$ (71)

which gives $I_{12}\lesssim d^{-1/2}$. ∎

###### Proof of Lemma 2.

It suffices to prove that for any finite integer $j>1$,

$\displaystyle\displaystyle{\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j}x_{t}^{2}d{\bm{x}}}=d^{-j}+O(d^{-j-1}).$ (72)

Indeed, assuming this estimate, we have $C_{d,j}^{-2}=d^{-j}+O(d^{-j-1})$ and
$C_{d,j}^{2}=d^{j}+O(d^{j-1})$. For any ${\bm{i}}$ and ${\bm{j}}$ with
${\bm{i}}\cap{\bm{j}}=\emptyset$,

$\displaystyle\mathbb{E}_{\bm{x}}g_{\bm{i}}^{2}({\bm{x}})g_{\bm{j}}^{2}({\bm{x}})-1=C_{d,r}^{4}\int_{{\mathbb{S}}_{d-1}}\prod_{t\in{\bm{i}}\cup{\bm{j}}}x_{t}^{2}d{\bm{x}}=C_{d,r}^{4}C_{d,2r}^{-2}-1=O(d^{-1})$
(73)

It remains to prove Eq. (72). By symmetries,

$\displaystyle\int_{{\mathbb{S}}_{d-1}}x_{t}^{2}d{\bm{x}}=\frac{1}{d}\int_{{\mathbb{S}}_{d-1}}\sum_{1\leq
t\leq
d}x_{t}^{2}d{\bm{x}}=\frac{1}{d}\int_{{\mathbb{S}}_{d-1}}1d{\bm{x}}=\frac{1}{d}\,.$
(74)

By symmetries again,

$\displaystyle\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j-1}x_{t}^{2}d{\bm{x}}$ (75) $\displaystyle=$
$\displaystyle\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j-1}x_{t}^{2}\left(\sum_{1\leq i\leq j-1}x_{i}^{2}+\sum_{j\leq i\leq
d}x_{i}^{2}\right)d{\bm{x}}$ (76) $\displaystyle=$
$\displaystyle(d-j+1)\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j}x_{t}^{2}d{\bm{x}}+(j-1)\int_{{\mathbb{S}}_{d-1}}x_{1}^{2}\prod_{1\leq t\leq
j-1}x_{t}^{2}d{\bm{x}}$ (77)

We use hypercontractivities to bound the error term, namely, the second term.
Recall that for any $q\geq 2$ and any polynomial defined on the sphere,

$\displaystyle\left(\int_{{\mathbb{S}}_{d-1}}|f({\bm{x}})|^{q}d{\bm{x}}\right)^{1/q}\leq(q-1)^{\deg(f)/2}\left(\int_{{\mathbb{S}}_{d-1}}|f({\bm{x}})|^{2}d{\bm{x}}\right)^{1/2}\,.$
(78)

Setting $f({\bm{x}})=x_{t}$ (with $\deg(f)=1$) gives

$\displaystyle\int_{{\mathbb{S}}_{d-1}}|x_{t}|^{q}d{\bm{x}}\leq(q-1)^{q/2}\left(\int_{{\mathbb{S}}_{d-1}}|x_{t}|^{2}d{\bm{x}}\right)^{q/2}=(q-1)^{q/2}d^{-\frac{q}{2}}\,.$
(79)

By Holder’s inequality and symmtries

$\displaystyle\int_{{\mathbb{S}}_{d-1}}x_{1}^{2}\prod_{1\leq t\leq
j-1}x_{t}^{2}d{\bm{x}}$
$\displaystyle\leq\left(\int_{{\mathbb{S}}_{d-1}}|x_{1}|^{2j}d{\bm{x}}\prod_{1\leq
t\leq j-1}\int_{{\mathbb{S}}_{d-1}}|x_{t}|^{2j}d{\bm{x}}\right)^{\frac{1}{j}}$
(80)
$\displaystyle=\int_{{\mathbb{S}}_{d-1}}|x_{1}|^{2j}d{\bm{x}}\leq(2j-1)^{j}d^{-j}$
(81)

Thus

$\displaystyle\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j-1}x_{t}^{2}d{\bm{x}}=(d-j+1)\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j}x_{t}^{2}d{\bm{x}}+O(d^{-j})$ (82)

and

$\displaystyle\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j}x_{t}^{2}d{\bm{x}}=$
$\displaystyle\frac{1}{d-(j-1)}\int_{{\mathbb{S}}_{d-1}}\prod_{1\leq t\leq
j-1}x_{t}^{2}d{\bm{x}}+O(d^{-j-1})\,.$ (83)

Finally, Eq. (72) is a consequence of this estimate and induction.

∎

##  Appendix C Generalization

We aim to obtain the asymptotic formulas for the test error in this section.
In the high-level, we decompose the empirical kernel $K({\bm{X}},{\bm{X}})$
into low-, critical- and high-frequency modes, where we have concentration in
the low- and high-frequency parts of the kernel. The test error associated to
these two parts are easier to handle. The critical-frequency part is more
difficult in which random matrix behaviors emerge, namely, the Marchenko-
Pastur distribution. As such, our first step is to remove the contribution in
the test error coming from the non-critical frequency parts. After that, the
remaining is essentially equivalent to computing the trace of certain
functional forms related to the Marchenko-Pastur distribution.

We consider a general setting that includes the dot-product kernels, the one-
hidden-layer and the multiple-layer convolutional kernels (NNGP and NT
kernels.) In what follows, we use
$\Delta_{d},\Delta_{d}^{\prime},\Delta_{d}^{\prime\prime}$, etc. to represent
quantities that converge to 0 in probability (the absolute value of a scalar,
the norm of a vector, the operator norm of a matrix, etc.), whose exact form
may change from line to line.

###  C.1 Setup

For $d\in{\mathbb{N}}^{*}$, let
$\bm{{\mathcal{X}}}^{(d)}\subseteq{\mathbb{R}}^{d}$ be the input space
associated with a probability measure $\sigma^{(d)}$ and a kernel function
$K^{(d)}$. Assume the kernel function has the following eigen-structure

$\displaystyle K^{(d)}({\bm{x}},{\bm{x}}^{\prime})=\sum_{k\geq
1}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{2}\sum_{l\in
N_{kn}^{(d)}}\phi^{(d)}_{knl}({\bm{x}})\phi^{(d)}_{knl}({\bm{x}}^{\prime})$
(84)

in the sense $K^{(d)}$, as the integral operator from
$L^{2}(\bm{{\mathcal{X}}}^{(d)},\sigma^{(d)})$ to itself,

$\displaystyle K^{(d)}\phi_{knl}^{(d)}({\bm{x}})=\int
K^{(d)}({\bm{x}},{\bm{x}}^{\prime})\phi^{(d)}_{knl}({\bm{x}}^{\prime})\sigma^{(d)}(d{\bm{x}}^{\prime})=(\sigma_{kn}^{(d)})^{2}\phi^{(d)}_{knl}({\bm{x}}).$
(85)

Here $\\{\phi^{(d)}_{knl}\\}_{knl}$ is an orthonormal basis of
$L^{2}(\bm{{\mathcal{X}}}^{(d)},\sigma^{(d)})$. We also assume $K^{(d)}$ is a
trace-class operator, i.e.,

$\displaystyle\sum_{knl}\langle
K^{(d)}\phi_{knl}^{(d)},\phi_{knl}^{(d)}\rangle=\sum_{k\geq
1}\sum_{n\in[E_{k}]}N_{kn}^{(d)}(\sigma_{kn}^{(d)})^{2}<\infty\,.$ (86)

In the above notations, we use the triplet $(k,n,l)$ to index the
eigenfunctions $\phi_{knl}^{(d)}$. The tuple $(k,n)$ determines the
eigenspace, whose eigenvalue is of the form
“$(\sigma_{kn}^{(d)})^{2}=C_{n}d^{-s_{k}}+\mathrm{Lower\,\,Order}$" and $l$
lists all eigenfunctions in the $kn$-eigenspace. We make the following
assumptions.

#### Kernel Assumptions.

  1. (1.)

Spectral Gap. There are $\delta_{0}>0$ and a sequence of strictly increasing
positive real numbers $\\{s_{k}\\}$ with $|s_{k}-s_{k-1}|\geq\delta_{0}$ for
all $k\geq 2$ such that

$\displaystyle(\sigma_{kn}^{(d)})^{2}\sim d^{-s_{k}}\sim(N_{kn}^{(d)})^{-1}$
(87)

Moreover, $\\{E_{k}\\}\subseteq{\mathbb{N}}^{*}$ is independent from $d$ which
grows at most exponentially. We also assume that there is a sequence of real
numbers $\\{\hat{h}_{kn}^{2}\\}_{kn}$ with $\hat{h}_{kn}^{2}\neq 0$ unless $k$
is sufficiently large and

$\displaystyle\sum_{k}\sum_{n\in[E_{k}]}\hat{h}_{kn}^{2}<\infty\quad\mathrm{and}\quad(\sigma_{kn}^{(d)})^{2}N_{kn}^{(d)}=\hat{h}_{kn}^{2}\quad\mathrm{as}\quad
d\to\infty$ (88)

  2. (2.)

Hypercontractivity Inequalities. For any $p\geq 2$ there are constant
$C_{p,k}$ such that for any function $f$ in the closure of
$\mathrm{Span}\\{\phi^{(d)}_{jnl}\\}_{j\leq k}$

$\displaystyle\|f\|_{p}\leq C_{p,k}\|f\|_{2}$ (89)

  3. (3.)

Concentration of Quadratic Forms. Let $\phi^{(d)}_{k}({\bm{x}})$ denote the
column vector consists of elements
$\\{\phi^{(d)}_{knl}({\bm{x}})\\}_{l\in[N_{kn}{(d)}n\in[E_{k}]}$. For every
sequence of matrices $\\{{\bm{A}}^{(d)}\\}$ with uniformly bounded operator
norm,

$\displaystyle\left(\sum_{n\in[E_{k}]}N_{kn}^{(d)}\right)^{-2}\mathbb{E}_{{\bm{x}}}|\phi^{(d)}_{k}({\bm{x}})^{T}{\bm{A}}^{(d)}\phi^{(d)}_{k}({\bm{x}})-\operatorname{Tr}{\bm{A}}^{(d)}|^{2}\to
0\quad\mathrm{as}\quad d\to\infty.$ (90)

  4. (4.)

Addition Theorem. For $k\in{\mathbb{N}}^{*}$ and $n\in[E_{k}]$ and
${\bm{x}}\in\bm{{\mathcal{X}}}^{(d)}$

$\displaystyle\sum_{l\in[N_{kn}^{(d)}]}\phi_{knl}^{(d)}({\bm{x}})^{2}=N_{kn}^{(d)}$
(91)

Let us briefly explain the assumptions. The Spectral Gap assumption basically
says, we can classify the eigenvectors into countably many categories indexed
by $k\in{\mathbb{N}}^{*}$. In the $k$-th category, it has exactly $E_{k}$ many
eigenspaces, each of them has dimensions $\sim d^{s_{k}}$ and eigenvalues
$d^{-s_{k}}$. It also implies the number of eigenfunctions with eigenvalues
$\lesssim d^{-s_{k}}$ is $\sim d^{s_{k}}$. Assumptions (1.), (2.) and (4.)
together are stronger than those in Theorem 6 in Mei et al. [2021] (and
slightly less technical), which allow us to apply kernel concentration from
Mei et al. [2021]. In particular, they imply concentration of the low- and
high-frequency parts of the empirical kernel $K^{(d)}({\bm{X}},{\bm{X}})$.
Finally, Assumption (3) is designed to meet the requirements in Lemma 1, which
allows us to claim Marchenko-Pastur type behavior of the gram matrix induced
by the feature map $\phi_{k}$. We provide a couple examples.

######  Example 1 (Dot-product Kernels).

When $\bm{{\mathcal{X}}}^{(d)}={\mathbb{S}}_{d-1}$ and $K^{(d)}$ is the dot-
product kernel, we have $E_{k}=1$, $s_{k}=k$, $N_{kn}=N(d,k)\sim d^{k}/k!$,
and $\phi_{knl}=Y_{kl}$ (note that $n=0$ since $E_{k}=1$.) Note that by the
Addition Theorem of spherical harmonics (Theorem 4.11 in Frye and Efthimiou
[2012]),

$\displaystyle\sum_{l\in[N(d,k)]}Y_{kl}({\bm{x}})^{2}=N(d,k)$ (92)

######  Example 2 (One-hidden-layer Convolutional Kernels).

Sightly more general setting is the one-layer convolutional kernel (NNGP or NT
kernels). In this case, $\bm{{\mathcal{X}}}^{(d)}={\mathbb{S}}_{d_{0}-1}^{p}$
where $p$ is the number of patches and the input dimension is $d=pd_{0}$. We
can set either $p=O(1)$ (i.e. independent of $d_{0}\to\infty$) or $p\sim
d^{\alpha_{p}}$ for some $\alpha_{p}>0$. This kernel is essentially the sum of
$p$ dot-product kernels. As such, $E_{k}=1$, $N_{kn}{(d)}=pN(d_{0},k)\sim
pd_{0}^{k}/k!$ and $(\sigma_{kn}^{u}pd)^{2}\sim(pd_{0}^{k})^{-1}$. If $p\sim
d^{\alpha_{p}}$ and $d_{0}\sim d^{\alpha_{d_{0}}}$ with
$\alpha_{d_{0}}+\alpha_{p}=1$, we have $s_{k}=\alpha_{p}+k\alpha_{d_{0}}$ and
$d^{-s_{k}}$ is the decay rate of the $k$-th order spherical harmonics.

######  Example 3 (Multiple-layer Convolutional Kernels).

General convolutional kernels are much more complicated [Xiao, 2021]. In this
case, $\bm{{\mathcal{X}}}^{(d)}={\mathbb{S}}_{d_{0}-1}^{p}$ where $p$ is the
number of patches and the input dimension is $d=pd_{0}$. We additionally
assume, $p=k_{0}^{L-1}$ for some $k_{0}\in\mathbb{N}^{*}$ and the network has
$L$ convolutional layers with filter size and strides being the same in each
layer (equal to $d_{0}$ in the first layer and to $k_{0}$ for the remaining
$(L-1)$ layers.) The eigenstructures of such kernels are studied in Xiao
[2021]. The eigenfunctions are tensor products of spherical harmonics defined
on copies of ${\mathbb{S}}_{d_{0}-1}$,

$\displaystyle
Y_{{\bm{k}},{\bm{l}}}({\bm{x}})=\prod_{i\in[p]}Y_{k_{i}l_{i}}({\bm{x}}_{i})$
(93)

The eigenvalues are more complicated to compute as they depend on both the
frequencies of $Y_{{\bm{k}},{\bm{l}}}$ and the topologies of the networks.
When $d_{0}\propto d^{\alpha_{d_{0}}}$ and $k_{0}\propto d^{\alpha_{k_{0}}}$
with $\alpha_{d_{0}}+(L-1)\alpha_{k_{0}}=1$ and
$\alpha_{k_{0}},\alpha_{d_{0}}>0$, the eigenvalue of $Y_{{\bm{k}}{\bm{l}}}$ is
$\propto d^{-({\mathscr{F}}({\bm{k}})+{\mathscr{S}}({\bm{k}}))}$, as
$d\to\infty$. Here ${\mathscr{F}}({\bm{k}})\equiv|{\bm{k}}|\alpha_{d_{0}}$ is
the frequency index of $Y_{{\bm{k}}{\bm{l}}}$ and
${\mathscr{S}}({\bm{k}})=J_{\bm{k}}\alpha_{k_{0}}$ is the spatial index, where
$J_{\bm{k}}$ is the number of edges in the sub-tree connecting all interacting
patches (i.e. $k_{i}\neq 0$) to the output; see Xiao [2021] for more details.

As there can possibly exist ${\bm{k}}$ and ${\bm{k}}^{\prime}$ with
${\mathscr{F}}({\bm{k}})+{\mathscr{S}}({\bm{k}})={\mathscr{F}}({\bm{k}}^{\prime})+{\mathscr{S}}({\bm{k}}^{\prime})$
(i.e., same order of decay) but
$({\mathscr{F}}({\bm{k}}),{\mathscr{S}}({\bm{k}}))\neq({\mathscr{F}}({\bm{k}}^{\prime}),{\mathscr{S}}({\bm{k}}^{\prime}))$
(i.e. different space-frequency combination), there can exist more than one
eigenspaces whose eigenvalues decay to zero with the same rate
$d^{-({\mathscr{F}}({\bm{k}})+{\mathscr{S}}({\bm{k}}))}$, but with different
leading coefficients. This is the main reason why we need to allow $|E_{k}|>1$
in Eq. (84).

Next we discuss the assumptions on the label function. Let ${\bm{X}}$ be the
training set with $m\sim d^{s_{r}}$ many training samples for some
$r\in{\mathbb{N}}^{*}$ fixed. Then let the ground true label function to be

$\displaystyle
f({\bm{x}})=\sum_{k\in{\mathbb{N}}^{*}}\sum_{n\in[E_{k}]}\sum_{l\in[N_{kn}^{(d)}]}\hat{f}_{knl}\phi_{knl}^{(d)}({\bm{x}}).$
(94)

Let $N_{k}^{(d)}=\sum_{n\in[E_{k}]}N_{kn}^{(d)}$. We assume, for $k\geq r$,
$\hat{\bm{f}}_{kn}=\\{\hat{f}_{knl}\\}_{l\in[N_{kn}^{(d)}]}$ is a random
vector with

$\displaystyle\mathbb{E}\hat{\bm{f}}_{kn}=\bm{0}\quad\mathrm{and}\quad\mathbb{E}\hat{\bm{f}}_{kn}\hat{\bm{f}}_{kn}^{T}=\frac{\hat{F}_{kn}^{2}}{N_{kn}^{(d)}}{\bm{I}}_{N_{kn}^{(d)}}.$
(95)

and $\\{\hat{\bm{f}}_{kn}\\}_{n\in[E_{k}],k\geq r}$ are mutually independent.
One concrete example is

$\displaystyle\hat{\bm{f}}_{kn}\sim\mathcal{N}\left(\bm{0},\frac{\hat{F}_{kn}^{2}}{N_{kn}^{(d)}}{\bm{I}}_{N_{kn}^{(d)}}\right)\,.$
(96)

For $k<r$, we assume the coefficients are deterministic with
$\sum_{l}\hat{f}_{knl}^{2}=\hat{F}_{kn}^{2}$,
$\sum_{n}\hat{F}_{kn}^{2}=\hat{F}_{k}^{2}$ and

$\displaystyle\sum_{k\in{\mathbb{N}}^{*}}\hat{F}_{k}^{2}<\infty$ (97)

Our goal is to compute the average test error over the random labels defined
above in the scaling limit $m\sim d^{s_{r}}$.

###  C.2 Structure of the Empirical Kernels

For convenience, denote

$\displaystyle\phi^{(d)}_{\leq k}({\bm{x}})$
$\displaystyle=[\phi^{(d)}_{jnl}({\bm{x}})]_{l\in[N_{jn}^{(d)}],1\leq j\leq
k,n\in[E_{j}],}^{T}$ (98) $\displaystyle N_{k}^{(d)}$
$\displaystyle=\sum_{n\in[E_{k}]}N_{kn}^{(d)}$ (99) $\displaystyle N_{\leq
k}^{(d)}$ $\displaystyle=\sum_{1\leq j\leq k}N_{\leq j}^{(d)}$ (100)

Let ${\bm{x}}_{i}^{T}$ be the $i$-th row of the training matrix ${\bm{X}}$.
Similarly,

$\displaystyle
Z_{k}({\bm{X}})=[\phi^{(d)}_{k}({\bm{x}}_{0}),\dots,\phi^{(d)}_{k}({\bm{x}}_{m-1})]^{T}\quad\quad$
$\displaystyle Z_{\leq k}({\bm{X}})=[\phi^{(d)}_{\leq
k}({\bm{x}}_{0}),\dots,\phi^{(d)}_{\leq k}({\bm{x}}_{m-1})]^{T}$ (101)
$\displaystyle{\bm{\Lambda}}_{k}=\mathrm{diag}\left([(\sigma_{kn}^{(d)})^{2}{\bm{I}}_{N_{kn}^{(d)}}]_{n\in[E_{k}]}\right)$
$\displaystyle{\bm{\Lambda}}_{\leq
k}=\mathrm{diag}\left([(\sigma_{jn}^{(d)})^{2}{\bm{I}}_{N_{jn}^{(d)}}]_{n\in[E_{j}],1\leq
j\leq k}\right)$ (102)

Note that $Z_{k}({\bm{X}})$ ($Z_{\leq k}({\bm{X}})$) is an $m\times
N_{k}^{(d)}$ ($m\times N_{\leq k}^{(d)}$) matrix and ${\bm{\Lambda}}_{k}$
(${\bm{\Lambda}}_{\leq k}$) is an $N_{k}^{(d)}\times N_{k}^{(d)}$ ($N_{\leq
k}^{(d)}\times N_{\leq k}^{(d)}$) diagonal matrix. The followings are defined
similarly,

$\displaystyle Z_{<k}({\bm{X}}),\quad
Z_{kn}({\bm{X}}),\quad{\bm{\Lambda}}_{<k},\quad{\bm{\Lambda}}_{kn},\quad
N_{<k}^{(d)},\quad N_{kn}^{(d)}\,.$ (103)

Next, we decompose the train-train kernel into two parts: the $\leq r$
frequency part and the $>r$ frequency parts,

$\displaystyle K^{(d)}({\bm{X}},{\bm{X}})$
$\displaystyle=\sum_{k\in{\mathbb{N}}^{*}}Z_{k}({\bm{X}}){\bm{\Lambda}}_{k}Z_{k}({\bm{X}})^{T}=Z_{\leq
r}({\bm{X}}){\bm{\Lambda}}_{\leq}Z_{\leq r}({\bm{X}})^{T}+\sum_{k\geq
r+1}Z_{k}({\bm{X}}){\bm{\Lambda}}_{k}Z_{k}({\bm{X}})^{T}$ (104)
$\displaystyle=Z_{\leq r}({\bm{X}}){\bm{\Lambda}}_{\leq}Z_{\leq
r}({\bm{X}})^{T}+\sum_{k\geq
r+1}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{2}Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}$
(105) $\displaystyle\equiv K_{\leq
r}^{(d)}({\bm{X}},{\bm{X}})+K^{(d)}_{>r}({\bm{X}},{\bm{X}})$ (106)

Assumptions (1.) (2.) allow us to apply kernel concentration [Ghorbani et al.,
2021, Mei et al., 2021], which implies that the low-frequency and high-
frequency parts of the empirical kernels are concentrated. By saying
concentration in the high-frequency part, we mean

######  Claim 1.

Let

$\displaystyle\Delta_{kn}^{(d)}\equiv\frac{1}{N_{kn}^{(d)}}Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}-{\bm{I}}_{m}\,.$
(107)

Then

$\displaystyle\mathbb{E}\sum_{k>r}\sum_{n\in[E_{k}]}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}\to
0$ (108)

The proof of this claim essentially follows from the arguments and results in
Theorem 6 of [Mei et al., 2021]; see Sec. C.7. Thus

$\displaystyle K^{(d)}_{>r}({\bm{X}},{\bm{X}})$ $\displaystyle=\sum_{k\geq
r+1}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{2}N_{kn}^{(d)}(\displaystyle{\bm{I}}_{m}+\Delta_{kn}^{(d)})$
(109) $\displaystyle=\sum_{k\geq r+1}\hat{h}_{k}^{2}{\bm{I}}_{m}+\sum_{k\geq
r+1}\sum_{n\in[E_{k}]}\hat{h}_{kn}^{2}\Delta_{kn}^{(d)}$ (110)

Denote

$\displaystyle\hat{h}_{>r}^{2}$ $\displaystyle=\sum_{k\geq
r+1}\sum_{n\in[E_{k}]}N_{kn}^{(d)}(\sigma_{kn}^{(d)})^{2}\sim 1$ (111)
$\displaystyle\Delta_{>r}^{(d)}$ $\displaystyle=\sum_{k\geq
r+1}\sum_{n\in[E_{k}]}\hat{h}_{kn}^{2}\Delta_{kn}^{(d)}\,,$ (112)

we can write

$\displaystyle
K^{(d)}_{>r}({\bm{X}},{\bm{X}})=\hat{h}_{>r}^{2}{\bm{I}}_{m}+\Delta_{>r}^{(d)},\quad\mathrm{where}\quad\mathbb{E}\|\Delta_{>r}^{(d)}\|_{\mathrm{op}}\to
0$ (113)

By saying the low-frequency part of the kernel concentrates, we mean (Theorem
6 [Mei et al., 2021])

$\displaystyle\frac{1}{m}Z_{<}({\bm{X}})^{T}Z_{<}({\bm{X}})={\bm{I}}_{N_{<}^{(d)}}+\Delta_{<r}^{(d)},\quad\mathrm{where}\quad\mathbb{E}\|\Delta_{<r}^{(d)}\|_{\mathrm{op}}\to
0$ (114)

Finally, Lemma 1 and Assumption (3.) imply that if
$N_{r}^{(d)}/m\to\alpha\in(0,\infty)$, then the empirical measure of the
critical part of the kernel matrix
$\frac{1}{m}Z_{r}({\bm{X}})^{T}Z_{r}({\bm{X}})$ , and the low-and-critical
frequency parts $\frac{1}{m}Z_{\leq r}({\bm{X}})^{T}Z_{\leq r}({\bm{X}})$
converge to the Marchenko-Pastur distribution $\mu_{\alpha}$ weakly by Lemma
1. In particular,
$\|\frac{1}{m}Z_{r}({\bm{X}})^{T}Z_{r}({\bm{X}})\|_{\mathrm{op}}+\|\frac{1}{m}Z_{\leq
r}({\bm{X}})^{T}Z_{\leq r}({\bm{X}})\|_{\mathrm{op}}=O(1)$ in probability as
$d\to\infty$.

For convenience, we summarize the structure of the empirical kernel in the
following.

######  Corollary 3.

Assume Assumptions (1.-4.). Let $r\in{\mathbb{N}}^{*}$ and $\alpha>0$ be fixed
and $m=m^{(d)}$ be such that $N_{r}^{(d)}/m\to\alpha\in(0,\infty)$ as
$d\to\infty$. Let ${\bm{X}}$, of shape $m\times d$, be the training set matrix
whose rows are drawn, uniformly, iid from $\bm{{\mathcal{X}}}^{(d)}$. Then we
have the following structure for the empirical kernel matrix

High-frequency Features: $\displaystyle
K^{(d)}_{>r}({\bm{X}},{\bm{X}})=\hat{h}_{>r}^{2}{\bm{I}}_{m}+\Delta_{>r}^{(d)},$
(115) Low-frequency Features:
$\displaystyle\frac{1}{m}Z_{<}({\bm{X}})^{T}Z_{<}({\bm{X}})={\bm{I}}_{N_{<}^{(d)}}+\Delta_{<r}^{(d)}$
(116) Critical-frequency Features: $\displaystyle\text{the empirical measure
of}\,\,\frac{1}{m}Z_{r}({\bm{X}})^{T}Z_{r}({\bm{X}})\to\mu_{\alpha}$ (117)
Low-and-critical-frequency Features: $\displaystyle\text{the empirical measure
of}\,\,\frac{1}{m}Z_{\leq r}({\bm{X}})^{T}Z_{\leq r}({\bm{X}})\to\mu_{\alpha}$
(118)

Let $0\leq\lambda=O(1)$ be the regularization and
$\gamma=\lambda+\hat{h}_{>r}^{2}$ be the effective regularization. To ease the
notations, denote

$\displaystyle\overline{\bm{Z}}_{<}=\frac{1}{\sqrt{m}}Z_{<r}({\bm{X}})\quad$
$\displaystyle\overline{\bm{Z}}_{\leq}=\frac{1}{\sqrt{m}}Z_{\leq
r}({\bm{X}})\quad$
$\displaystyle\overline{\bm{Z}}_{>}=\frac{1}{\sqrt{m}}Z_{>r}({\bm{X}})\quad$
$\displaystyle\overline{\bm{Z}}_{r}=\frac{1}{\sqrt{m}}Z_{r}({\bm{X}})$
$\displaystyle\overline{\bm{\Lambda}}_{<}=\gamma^{-1}m{\bm{\Lambda}}_{<}\quad$
$\displaystyle\overline{\bm{\Lambda}}_{\leq}=\gamma^{-1}m{\bm{\Lambda}}_{\leq}\quad$
$\displaystyle\overline{\bm{\Lambda}}_{>}=\gamma^{-1}m{\bm{\Lambda}}_{>}\quad$
$\displaystyle\overline{\bm{\Lambda}}_{r}=\gamma^{-1}m{\bm{\Lambda}}_{r}$

Clearly, Corollary 3 and the assumption on the spectra imply that in
probability as $d\to\infty$,

$\displaystyle\|\overline{\bm{Z}}_{<}\|_{\mathrm{op}}+\|\overline{\bm{Z}}_{\leq}\|_{\mathrm{op}}+\|\overline{\bm{Z}}_{r}\|_{\mathrm{op}}+\|\overline{\bm{\Lambda}}_{<}\|_{\mathrm{op}}+\|\overline{\bm{\Lambda}}_{\leq}\|_{\mathrm{op}}+\|\overline{\bm{\Lambda}}_{r}\|_{\mathrm{op}}\lesssim
1$ (119)

Then we can write $K^{(d)}$ and $K_{\lambda}^{(d)}$ as

$\displaystyle K^{(d)}_{\lambda}({\bm{X}},{\bm{X}})\equiv
K^{(d)}({\bm{X}},{\bm{X}})+\lambda{\bm{I}}_{m}=\gamma(\overline{\bm{Z}}_{\leq}\overline{\bm{\Lambda}}_{\leq}\overline{\bm{Z}}_{\leq}^{T}+{\bm{I}}_{m})+\Delta_{>r}^{(d)}\equiv{\bm{K}}+\Delta_{>r}^{(d)}$
(120)

where

$\displaystyle{\bm{K}}=\gamma(\overline{\bm{Z}}_{\leq}\overline{\bm{\Lambda}}_{\leq}\overline{\bm{Z}}_{\leq}^{T}+{\bm{I}}_{m})$
(121)

Then by Sherman–Morrison–Woodbury formula

$\displaystyle{\bm{K}}^{-1}$
$\displaystyle=\gamma^{-1}\left({\bm{I}}_{m}-\overline{\bm{Z}}_{\leq}\left(\overline{\bm{\Lambda}}_{\leq}^{-1}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)^{-1}\overline{\bm{Z}}_{\leq}^{T}\right)=\gamma^{-1}\left({\bm{I}}_{m}-\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\right)$
(122)

where

$\displaystyle\overline{\bm{D}}$
$\displaystyle=\overline{\bm{\Lambda}}_{\leq}^{-1}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\,.$
(123)

The matrix $\overline{\bm{D}}$ plays a critical role in the remaining
analysis. We have the following control regarding its eigenvalues, which says
the eigenvalues of $\overline{\bm{D}}$ are away from $0$ and $\infty$

######  Lemma 3.

There are constants $0<\lambda_{1}<\lambda_{2}$ independent of $d$ such that,
in probability as $d\to\infty$,

$\displaystyle\lambda_{1}{\bm{I}}\prec\overline{\bm{D}}\prec\lambda_{2}{\bm{I}}$
(124)

We will prove the lemma later in Sec.C.6.

Note that

$\displaystyle
K^{(d)}_{\lambda}({\bm{X}},{\bm{X}})^{-1}=(K^{(d)}_{\lambda}({\bm{X}},{\bm{X}})^{-1}{\bm{K}}){\bm{K}}^{-1}=({\bm{I}}_{m}+{\bm{K}}^{-1}\Delta_{>r}^{(d)})^{-1}{\bm{K}}^{-1}=({\bm{I}}_{m}+\Delta_{d}^{\prime}){\bm{K}}^{-1}$
(125)

where $\|\Delta_{d}^{\prime}\|_{\mathrm{op}}\to 0$ in probability since
$\|{\bm{K}}^{-1}\|_{\mathrm{op}}\leq\gamma^{-1}$ and
$\|\Delta_{>r}^{(d)}\|_{\mathrm{op}}\to 0$ in probability.

Similarly, we can write

$\displaystyle M({\bm{X}},{\bm{X}})\equiv$
$\displaystyle\mathbb{E}K({\bm{X}},x)K({\bm{X}},x)^{T}$ (126) $\displaystyle=$
$\displaystyle\sum_{1\leq k\leq
r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{4}Z_{k}({\bm{X}})Z_{k}({\bm{X}})^{T}+\sum_{k>r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{4}Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}$
(127) $\displaystyle=$ $\displaystyle
m^{-1}\gamma^{2}\overline{\bm{Z}}_{\leq}\overline{\bm{\Lambda}}_{\leq}^{2}\overline{\bm{Z}}_{\leq}^{T}+m^{-1}\Delta_{d}^{\prime\prime}$
(128)

where

$\displaystyle\Delta_{d}^{\prime\prime}\equiv
m\sum_{k>r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{4}Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}=m\sum_{k>r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{4}N_{kn}^{(d)}({\bm{I}}_{m}+\Delta_{kn}^{(d)})$
(129)

We have $\|\Delta_{d}^{\prime\prime}\|_{\mathrm{op}}\to 0$ in probability as
$d\to\infty$, since

$\displaystyle
m\sum_{k>r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{4}N_{kn}^{(d)}\sim$
$\displaystyle
m\sum_{k>r}d^{-s_{k}}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{2}N_{kn}^{(d)}$
(130) $\displaystyle\lesssim$ $\displaystyle
d^{s_{r}}d^{-s_{r+1}}\sum_{k>r}\sum_{n\in[E_{k}]}(\sigma_{kn}^{(d)})^{2}N_{kn}^{(d)}$
(131) $\displaystyle\leq$ $\displaystyle
d^{s_{r}}d^{-s_{r+1}}\sum_{k>r}\sum_{n\in[E_{k}]}\hat{h}_{kn}^{2}\lesssim
d^{-\delta_{0}}\to 0$ (132)

Finally, the above estimates imply

$\displaystyle{\bm{H}}$ $\displaystyle\equiv
K_{\gamma}^{(d)}({\bm{X}},{\bm{X}})^{-1}M({\bm{X}},{\bm{X}})K_{\gamma}^{(d)}({\bm{X}},{\bm{X}})^{-1}$
(133)
$\displaystyle=({\bm{I}}_{m}+\Delta_{d}^{\prime}){\bm{K}}^{-1}M({\bm{X}},{\bm{X}}){\bm{K}}^{-1}({\bm{I}}_{m}+\Delta_{d}^{\prime})$
(134)
$\displaystyle=m^{-1}({\bm{I}}_{m}+\Delta_{d}^{\prime}){\bm{K}}^{-1}(\gamma^{2}\overline{\bm{Z}}_{\leq}\overline{\bm{\Lambda}}_{\leq}^{2}\overline{\bm{Z}}_{\leq}^{T}+\Delta_{d}^{\prime\prime}){\bm{K}}^{-1}({\bm{I}}_{m}+\Delta_{d}^{\prime})$
(135)
$\displaystyle=m^{-1}\left(\overline{\bm{Z}}_{\leq}\left(\overline{\bm{\Lambda}}_{\leq}^{-1}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)^{-2}\overline{\bm{Z}}_{\leq}^{T}+\Delta_{d}^{\prime\prime\prime}\right)$
(136)
$\displaystyle=m^{-1}\left(\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}+\Delta_{d}^{\prime\prime\prime}\right)$
(137)

with the error term $\|\Delta_{d}^{\prime\prime\prime}\|_{\mathrm{op}}\to 0$
in probability.

###  C.3 Reduction I: Reducing the MSE to Traces

We will repeatedly use the following simple results.

######  Lemma 4.

Let ${\bm{u}}$ be a random vector with $\mathbb{E}{\bm{u}}=\bm{0}$ and
$\mathbb{E}{\bm{u}}{\bm{u}}^{T}=\sigma^{2}{\bm{I}}_{k}$. Then for any $k\times
k$ deterministic matrix ${\bm{A}}$,

$\displaystyle\mathbb{E}_{\bm{u}}{\bm{u}}^{T}{\bm{A}}{\bm{u}}$
$\displaystyle=\sigma^{2}\operatorname{Tr}({\bm{A}})$ (138)
$\displaystyle\mathbb{E}_{\bm{u}}\|{\bm{A}}{\bm{u}}\|_{2}^{2}$
$\displaystyle=\sigma^{2}\operatorname{Tr}({\bm{A}}^{T}{\bm{A}})\,.$ (139)

Next, we compute the loss by decomposing it as follows. Recall that the
observed labels is $f({\bm{X}})+\bm{\epsilon}$ where $\bm{\epsilon}$ is the
iid noise term, which is centered and has variance $\sigma_{\epsilon}^{2}$.
Thus the average test error is

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},{\bm{h}})$
$\displaystyle=$
$\displaystyle\mathbb{E}_{{\bm{f}},\bm{\epsilon},{\bm{x}}}\left|f({\bm{x}})-K^{(d)}({\bm{x}},{\bm{X}})K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}(f({\bm{X}})+\bm{\epsilon})\right|^{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\bm{f}}\mathbb{E}_{\bm{x}}f^{2}({\bm{x}})-2\mathbb{E}_{{\bm{f}}}\mathbb{E}_{{\bm{x}}}f({\bm{x}})K^{(d)}({\bm{x}},{\bm{X}})K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}f({\bm{X}})+\mathbb{E}_{{\bm{f}}}f^{T}({\bm{X}}){\bm{H}}f({\bm{X}})+\sigma_{\epsilon}^{2}\mathrm{Tr}({\bm{H}})$
$\displaystyle\equiv$
$\displaystyle\mathbb{E}_{{\bm{f}}}I_{1}+\mathbb{E}_{{\bm{f}}}I_{2}+\mathbb{E}_{{\bm{f}}}I_{3}+I_{4}.$

Here

$\displaystyle I_{1}$ $\displaystyle=\mathbb{E}_{\bm{x}}f^{2}({\bm{x}})$ (140)
$\displaystyle I_{2}$
$\displaystyle=-2\mathbb{E}_{{\bm{f}}}\mathbb{E}_{{\bm{x}}}f({\bm{x}})K^{(d)}({\bm{x}},{\bm{X}})K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}f({\bm{X}})$
(141) $\displaystyle I_{3}$
$\displaystyle=\mathbb{E}_{{\bm{f}}}f^{T}({\bm{X}}){\bm{H}}f({\bm{X}})$ (142)
$\displaystyle I_{4}$
$\displaystyle=\sigma_{\epsilon}^{2}\mathrm{Tr}({\bm{H}})$ (143)

We estimate each $I_{i}$ individually.

####  Estimate $I_{1}$.

We simply keep it unchanged at the moment.

####  Estimate $I_{2}$.

Note that

$\displaystyle\mathbb{E}_{\bm{x}}f({\bm{x}})K^{(d)}({\bm{x}},{\bm{X}})$
$\displaystyle=\sum_{k}\sum_{n,l}\hat{f}_{knl}(\sigma_{kn}^{(d)})^{2}\phi^{(d)}_{knl}({\bm{X}})^{T}=\sqrt{m}\hat{\bm{f}}_{\leq}^{T}{\bm{\Lambda}}_{\leq}\overline{\bm{Z}}_{\leq}^{T}+\sum_{k>r}\sqrt{m}\hat{\bm{f}}_{k}^{T}{\bm{\Lambda}}_{k}\overline{\bm{Z}}_{k}^{T}$
(144)
$\displaystyle=\gamma\frac{1}{\sqrt{m}}\hat{\bm{f}}_{\leq}^{T}\overline{\bm{\Lambda}}_{\leq}\overline{\bm{Z}}_{\leq}^{T}+\sum_{k>r}\gamma\frac{1}{\sqrt{m}}\hat{\bm{f}}_{k}^{T}\overline{\bm{\Lambda}}_{k}\overline{\bm{Z}}_{k}^{T}$
(145)

where $\hat{\bm{f}}_{\leq}$ is the column vector with elements
$\\{\hat{f}_{knl}\\}_{k\leq r}$ and $\hat{\bm{f}}_{<}$, $\hat{\bm{f}}_{>}$,
$\hat{\bm{f}}_{k}$, etc. are defined similarly. In addition,

$\displaystyle
f({\bm{X}})=\sqrt{m}\overline{\bm{Z}}_{\leq}\hat{\bm{f}}_{\leq}+\sqrt{m}\sum_{k>r}\overline{\bm{Z}}_{k}\hat{\bm{f}}_{k}$
(146)

We then use the fact that ${\bm{f}}_{>}$ is centered to eliminate all cross
terms between $\hat{\bm{f}}_{\leq}$ and $\hat{\bm{f}}_{>}$. Denote
$\mathbb{E}_{>}$, $\mathbb{E}_{k}$ and $\mathbb{E}_{kn}$ the expectation
operator over $\hat{\bm{f}}_{>}$, over $\hat{\bm{f}}_{k}$ and over
$\hat{\bm{f}}_{kn}$ resp. Under this notation,
$\mathbb{E}_{\bm{f}}=\mathbb{E}_{r}\mathbb{E}_{>}$. Using Eq. (125) and Eq.
(126),

$\displaystyle\mathbb{E}_{>}I_{2}=-2\hat{\bm{f}}_{\leq}^{T}\overline{\bm{D}}^{-1}$
$\displaystyle\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\hat{\bm{f}}_{\leq}-2\gamma\mathbb{E}_{>}\sum_{k>r}\hat{\bm{f}}_{k}^{T}$
$\displaystyle\left(\overline{\bm{\Lambda}}_{k}\overline{\bm{Z}}_{k}^{T}K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}\overline{\bm{Z}}_{k}\right)\hat{\bm{f}}_{k}+\Delta_{d}$
(147)

for some $\Delta_{d}\to 0$ in probability. The second term goes to zero since,
for each $k>r$ and $n\in E_{k}$

$\displaystyle\mathbb{E}_{kn}\hat{\bm{f}}_{kn}^{T}\left(\overline{\bm{\Lambda}}_{kn}\overline{\bm{Z}}_{kn}^{T}K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}\overline{\bm{Z}}_{kn}\right)\hat{\bm{f}}_{kn}$
$\displaystyle=$
$\displaystyle\hat{F}_{kn}^{2}/N_{kn}^{(d)}\operatorname{Tr}\overline{\bm{\Lambda}}_{kn}\left(\overline{\bm{Z}}_{kn}^{T}K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}\overline{\bm{Z}}_{kn}\right)$
$\displaystyle=$
$\displaystyle\hat{F}_{kn}^{2}/N_{kn}^{(d)}\gamma^{-1}m(\sigma_{kn}^{(d)})^{2}\operatorname{Tr}\left(\overline{\bm{Z}}_{kn}^{T}K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}\overline{\bm{Z}}_{kn}\right)$
$\displaystyle\leq$
$\displaystyle\gamma^{-1}\hat{F}_{kn}^{2}/N_{kn}^{(d)}m(\sigma_{kn}^{(d)})^{2}\|K^{(d)}_{\gamma}({\bm{X}},{\bm{X}})^{-1}\|_{\mathrm{op}}\operatorname{Tr}\left(\overline{\bm{Z}}_{kn}\overline{\bm{Z}}_{kn}^{T}\right)$
$\displaystyle\lesssim$
$\displaystyle\hat{F}_{kn}^{2}/N_{kn}^{(d)}m(\sigma_{kn}^{(d)})^{2}\operatorname{Tr}(\frac{1}{m}N_{kn}^{(d)}\frac{1}{N_{kn}^{(d)}}{\bm{Z}}_{kn}{\bm{Z}}_{kn}^{T})$
$\displaystyle=$
$\displaystyle\hat{F}_{kn}^{2}m(\sigma_{kn}^{(d)})^{2}\frac{1}{m}\operatorname{Tr}({\bm{I}}_{m}+\Delta_{kn}^{(d)})$
$\displaystyle\lesssim$ $\displaystyle
m(\sigma_{kn}^{(d)})^{2}(1+\|\Delta_{kn}^{(d)}\|_{\mathrm{op}})$
$\displaystyle=$ $\displaystyle
m/N_{kn}^{{(d)}}\hat{h}_{kn}^{2}(1+\|\Delta_{kn}^{(d)}\|_{\mathrm{op}})$
$\displaystyle\sim$ $\displaystyle
d^{s_{r}-s_{k}}\hat{h}_{kn}^{2}(1+\|\Delta_{kn}^{(d)}\|_{\mathrm{op}})$

Clearly, the sum over $k>r$ and $n\in[E_{k}]$ of the above is bounded by
$\lesssim d^{-\delta_{0}}$ in probability. Thus

$\displaystyle\mathbb{E}_{>}I_{2}=-2\left(\hat{\bm{f}}_{\leq}^{T}\overline{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\hat{\bm{f}}_{\leq}\right)+\Delta_{d}$
(148)

#### Estimate $I_{3}$.

Again, we use the fact that cross terms have mean zero

$\displaystyle\mathbb{E}_{>}I_{3}$
$\displaystyle=f^{T}({\bm{X}}){\bm{H}}f({\bm{X}})=m\left(\hat{\bm{f}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}^{T}{\bm{H}}\overline{\bm{Z}}_{\leq}\hat{\bm{f}}_{\leq}+\sum_{k>r}\mathbb{E}_{k}\hat{\bm{f}}_{k}^{T}\overline{\bm{Z}}_{k}^{T}{\bm{H}}\overline{\bm{Z}}_{k}\hat{\bm{f}}_{k}\right)$
(149)
$\displaystyle=\hat{\bm{f}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\hat{\bm{f}}_{\leq}+\sum_{k>r}\mathbb{E}_{k}\hat{\bm{f}}_{k}^{T}\overline{\bm{Z}}_{k}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{k}\hat{\bm{f}}_{k}+\Delta_{d}$
(150) $\displaystyle\equiv I_{3,1}+I_{3,2}+\Delta_{d}$ (151)

Note that

$\displaystyle I_{3,2}$
$\displaystyle=\sum_{k>r}\sum_{n}\mathbb{E}_{kn}\operatorname{Tr}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\hat{\bm{f}}_{kn}{\bm{f}}_{kn}^{T}\overline{\bm{Z}}_{kn}\overline{\bm{Z}}_{kn}^{T}$
(152)
$\displaystyle=\sum_{k>r}\sum_{n}\hat{F}_{kn}^{2}/N_{k}^{(d)}\operatorname{Tr}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{kn}\overline{\bm{Z}}_{kn}^{T}$
(153)
$\displaystyle=\sum_{k>r}\sum_{n}\hat{F}_{kn}^{2}/m\operatorname{Tr}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}({\bm{I}}_{m}+\Delta_{kn}^{(d)})$
(154)
$\displaystyle=m^{-1}\sum_{k>r}\sum_{n}\hat{F}_{kn}^{2}\left(\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)(1+\Delta_{d})$
(155)
$\displaystyle=m^{-1}\hat{F}_{>r}^{2}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}+\Delta_{d}^{\prime}$
(156)

where $\Delta_{d}^{\prime}\to 0$ in probability since
$\|{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\|_{\mathrm{op}}\lesssim
1$ in probability.

#### Estimate $I_{4}$.

The $I_{4}$ is similar to $I_{3,2}$ above and we have

$\displaystyle
I_{4}=m^{-1}\sigma_{\epsilon}^{2}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}+\Delta_{d}$
(157)

#### All Together.

Combining all terms we have

$\displaystyle\mathbb{E}_{>}I_{1}+I_{2}+I_{3}+I_{4}$ (158)
$\displaystyle=\left\|\left({\bm{I}}-\overline{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{\bm{f}}_{\leq}\right\|_{2}^{2}+\left(1+m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{F}_{>r}^{2}+m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\sigma_{\epsilon}^{2}+\Delta_{d}$
(159)
$\displaystyle=\left\|\overline{\bm{D}}^{-1}\overline{\bm{\Lambda}}_{\leq}^{-1}\hat{\bm{f}}_{\leq}\right\|_{2}^{2}+\left(1+m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{F}_{>r}^{2}+m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\sigma_{\epsilon}^{2}+\Delta_{d}$
(160) $\displaystyle\equiv T_{1}+T_{2}+T_{3}+\Delta_{d}$ (161)

where

$\displaystyle T_{1}$
$\displaystyle=\left\|\overline{\bm{D}}^{-1}\overline{\bm{\Lambda}}_{\leq}^{-1}\hat{\bm{f}}_{\leq}\right\|_{2}^{2}$
(162) $\displaystyle T_{2}$
$\displaystyle=\left(1+m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{F}_{>r}^{2}$
(163) $\displaystyle T_{3}$
$\displaystyle=m^{-1}\operatorname{Tr}\overline{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\sigma_{\epsilon}^{2}$
(164)

As such, it remains to handle

$\displaystyle\mathbb{E}_{r}(T_{1}+T_{2}+T_{3})=\mathbb{E}_{r}T_{1}+T_{2}+T_{3},.$
(165)

###  C.4 Reduction II: Reducing Traces to Integrals

Recall that $\overline{\bm{\Lambda}}_{\leq}$ is a diagonal matrix with
elements $\gamma^{-1}m(\sigma_{kn}^{(d)})^{2}$, whose multiplicity is
$N_{kn}^{(d)}$ for $k\leq r$ and $n\in[E_{k}]$. When $k=r$,
$\gamma^{-1}m(\sigma_{rn}^{(d)})^{2}\sim 1$, otherwise
$\gamma^{-1}m(\sigma_{kn}^{(d)})^{2}\sim d^{s_{r}-s_{k}}$. For convenience, we
let

$\displaystyle N_{<}\equiv N_{<}^{(d)},\quad N_{=}\equiv N_{r}^{(d)},\quad
N_{\leq}\equiv N_{\leq}^{(d)}$ (166)

Therefore,

$\displaystyle\overline{\bm{\Lambda}}_{\leq}^{-1}=\begin{bmatrix}\Delta&0\\\
0&{\bm{R}}\end{bmatrix}$ (167)

where $\Delta$ is an $N_{<}\times N_{<}$ diagonal matrix whose entries are
$\gamma m^{-1}(\sigma_{kn}^{(d)})^{-2}\sim d^{-(s_{r}-s_{k})}$, and ${\bm{R}}$
is a $N_{=}\times N_{=}$ diagonal matrix whose entries are $\gamma
m^{-1}(\sigma_{rn}^{(d)})^{-2}\sim 1$. As such, we claim that we can replace
$\hat{\bm{f}}_{\leq}$ by $[\bm{0},\hat{\bm{f}}_{=}]$ in estimating $T_{1}$ and
replace the $\Delta$ in $\overline{\bm{\Lambda}}_{\leq}^{-1}$ by any
$\rho{\bm{I}}_{N_{<}}$ for any finite non-negative constant $\rho$ in
estimating $T_{2}$. The first claim is obvious as, by Lemma 3,

$\displaystyle\|\overline{\bm{D}}^{-1}\overline{\bm{\Lambda}}^{-1}_{\leq}[\bm{{\bm{f}}}_{<},\bm{0}]^{T}\|_{2}\leq\|\overline{\bm{D}}^{-1}\|_{\mathrm{op}}\|\overline{\bm{\Lambda}}^{-1}_{\leq}[\bm{{\bm{f}}}_{<},\bm{0}]^{T}\|_{2}\lesssim\lambda_{1}^{-1}(m\sigma_{r-1}^{2}\gamma^{-1})\|\bm{{\bm{f}}}_{<}\|_{2}\lesssim
d^{-(s_{r}-s_{k})}\to 0$ (168)

To prove the second claim regarding estimating $T_{2}$, denote

$\displaystyle\tilde{\bm{D}}=\begin{bmatrix}\rho{\bm{I}}_{N_{<}}&0\\\
0&{\bm{R}}\end{bmatrix}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$
(169)

We claim that, in probability,

$\displaystyle
m^{-1}\operatorname{Tr}\left((\overline{\bm{D}}^{-2}-\tilde{\bm{D}}^{-2})\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)$
(170) $\displaystyle=$ $\displaystyle
m^{-1}\operatorname{Tr}\left((\overline{\bm{D}}^{-1}(\overline{\bm{D}}^{-1}-\tilde{\bm{D}}^{-1})+(\overline{\bm{D}}^{-1}-\tilde{\bm{D}}^{-1})\tilde{\bm{D}}^{-1})\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\to
0$ (171)

We only bound the first term as the second term can be handled similarly.

$\displaystyle
m^{-1}\operatorname{Tr}\left(\overline{\bm{D}}^{-1}(\overline{\bm{D}}^{-1}-\tilde{\bm{D}}^{-1})\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)$
(172) $\displaystyle=$ $\displaystyle
m^{-1}\operatorname{Tr}\left((\overline{\bm{D}}^{-1}-\tilde{\bm{D}}^{-1})\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-1}\right)$
(173) $\displaystyle=$ $\displaystyle
m^{-1}\operatorname{Tr}\left(\overline{\bm{D}}^{-1}(\tilde{\bm{D}}-\overline{\bm{D}})\tilde{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-1}\right)$
(174) $\displaystyle=$ $\displaystyle
m^{-1}\operatorname{Tr}\left((\tilde{\bm{D}}-\overline{\bm{D}})\tilde{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}\right)$
(175)

Then we use the facts that (1) the upper right $N_{<}\times N_{<}$ block
matrix of $(\tilde{\bm{D}}-\overline{\bm{D}})$ is a diagonal matrix whose
entries are in $[0,1]$ and the three remaining block matrices are zeros, and
(2) all entries in
$\tilde{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}$
are bounded above by a constant (each matrix in
$\tilde{\bm{D}}^{-1}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\overline{\bm{D}}^{-2}$
has bounded operator norm333Recall that
$\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$ follows the Marchenko-
Pastur distribution.) to conclude that

$\displaystyle\left|m^{-1}\operatorname{Tr}\left(\overline{\bm{D}}^{-1}(\overline{\bm{D}}^{-1}-\tilde{\bm{D}}^{-1})\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\right|\lesssim
m^{-1}N_{<}$ (176)

Thus

$\displaystyle
T_{2}=\left(1+m^{-1}\operatorname{Tr}\tilde{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{F}_{>r}^{2}+\Delta_{d}$
(177)

which will be handled later.

It remains to handle $T_{1}$. We make two steps of reductions in estimating
$\mathbb{E}_{r}T_{1}$. The first one is to replace
$\overline{\bm{\Lambda}}_{\leq}^{-1}$ by

$\displaystyle\tilde{\bm{\Lambda}}_{\leq}^{-1}=\begin{bmatrix}0&0\\\
0&{\bm{R}}\end{bmatrix}$ (178)

and the second one is to replace
$\overline{\bm{\Lambda}}_{\leq}^{-1}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$
by

$\displaystyle{\bm{W}}\equiv\begin{bmatrix}{\bm{I}}_{N_{<}}&{\bm{B}}\\\
{\bm{B}}^{T}&{\bm{C}}\end{bmatrix}\equiv\begin{bmatrix}{\bm{I}}_{N_{<}}&\overline{\bm{Z}}_{<}^{T}\overline{\bm{Z}}_{=}\\\
\overline{\bm{Z}}_{=}^{T}\overline{\bm{Z}}_{<}&\overline{\bm{Z}}_{=}^{T}\overline{\bm{Z}}_{=}+{\bm{R}}\end{bmatrix}$
(179)

Here we have applied

$\displaystyle\overline{\bm{Z}}_{<}^{T}\overline{\bm{Z}}_{<}={\bm{I}}_{N_{<}}+\Delta_{d}$
(180)

The reason we could do so is exactly the same as we replaced
$\overline{\bm{D}}$ by $\tilde{\bm{D}}$ above as we only perturb the entries
in the upper $N_{<}\times N_{<}$ block by $O(1)$.

Note that ${\bm{W}}$ is symmetric and is also strictly positive definite, i.e.
the minimal eigenvalue of ${\bm{W}}$ is $\gtrsim 1$; see the proof in Sec.C.6.
Thus by the Schur complement,

$\displaystyle\mathbb{E}_{r}T_{1}$
$\displaystyle=\left\|\begin{bmatrix}{\bm{I}}_{N_{<}}&{\bm{B}}\\\
{\bm{B}}^{T}&{\bm{C}}\end{bmatrix}^{-1}\begin{bmatrix}0&0\\\
0&{\bm{R}}\end{bmatrix}\begin{bmatrix}0\\\
\hat{\bm{f}}_{=}\end{bmatrix}\right\|_{2}^{2}+\Delta_{d}$ (181)
$\displaystyle=\mathbb{E}_{r}\left\|\begin{bmatrix}-{\bm{B}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}\hat{\bm{f}}_{=}\\\
({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}\hat{\bm{f}}_{=}\end{bmatrix}\right\|_{2}^{2}+\Delta_{d}$
(182)

By the fact that $\hat{\bm{f}}_{=}$ is mean zero and isotropic, we have the
above equal to

$\displaystyle\mathbb{E}_{r}T_{1}$
$\displaystyle=\operatorname{Tr}\left({\bm{R}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{B}}^{T}{\bm{B}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}+{\bm{R}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}{\bm{R}}\right)\hat{F}_{r}^{2}/N_{=}+\Delta_{d}$
(183)
$\displaystyle=\operatorname{Tr}({\bm{R}}{\bm{C}}^{-2}{\bm{R}})\hat{F}_{r}^{2}/N_{=}+\Delta_{d}^{\prime}+\Delta_{d}$
(184)

where

$\displaystyle\Delta_{d}^{\prime}=\operatorname{Tr}$
$\displaystyle\left({\bm{R}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{B}}^{T}{\bm{B}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}\right)/N_{=}+$
(185)
$\displaystyle\operatorname{Tr}\left({\bm{R}}(({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}-{\bm{C}}^{-2}){\bm{R}}\right)/N_{=}\,.$
(186)

We claim that $\Delta_{d}^{\prime}\to 0$ in probability. For the first term,
we have

$\displaystyle\operatorname{Tr}\left({\bm{R}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{B}}^{T}{\bm{B}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}\right)/N_{=}$
(187) $\displaystyle=$
$\displaystyle\operatorname{Tr}\left(({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}^{2}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{B}}^{T}{\bm{B}}\right)/N_{=}$
(188) $\displaystyle=$
$\displaystyle\|({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{R}}^{2}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}\|_{\mathrm{op}}\operatorname{Tr}\left({\bm{B}}^{T}{\bm{B}}\right)/N_{=}$
(189) $\displaystyle\leq$
$\displaystyle\|({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}\|_{\mathrm{op}}\|{\bm{R}}^{2}\|_{\mathrm{op}}\|({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}\|_{\mathrm{op}}\operatorname{Tr}\left({\bm{B}}^{T}{\bm{B}}\right)/N_{=}$
(190) $\displaystyle\lesssim$
$\displaystyle\operatorname{Tr}\left({\bm{B}}^{T}{\bm{B}}\right)/N_{=}\sim
N_{<}/N_{=}\to 0\,$ (191)

in probability as $d\to\infty$. We have used

$\displaystyle\|{\bm{R}}\|_{\mathrm{op}}\lesssim 1$ (192)
$\displaystyle\|({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}\|_{\mathrm{op}}\leq\|{\bm{W}}^{-1}\|_{\mathrm{op}}\lesssim
1$ (193)
$\displaystyle\frac{1}{N_{=}}\operatorname{Tr}\left({\bm{B}}^{T}{\bm{B}}\right)\sim
N_{<}/N_{=}\,.$ (194)

The last one holds because ${\bm{B}}^{T}{\bm{B}}$ is a rank $N_{<}$ matrix
with operator norm $\lesssim 1$. Note that this also implies that
${\bm{B}}^{T}{\bm{B}}$ has at most $N_{<}$ many non-zero singular values,
which is upper bounded by $\lesssim 1$. Using Von Neumann’s trace
inequalities, for any matrix ${\bm{A}}$, we have

$\displaystyle|\operatorname{Tr}{\bm{A}}{\bm{B}}^{T}{\bm{B}}|\leq\sum_{j}\sigma_{j}({\bm{A}})\sigma_{j}({\bm{B}}^{T}{\bm{B}})\lesssim
N_{<}\|{\bm{A}}\|_{\mathrm{op}}$ (195)

where $\sigma_{j}({\bm{A}})$ is the $j$-th (in descending order) singular
value of a matrix ${\bm{A}}$. Now we proceed to control the second term. Note
that

$\displaystyle({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}-{\bm{C}}^{-2}$ (196)
$\displaystyle=$
$\displaystyle({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}-({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{C}}^{-1}+({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{C}}^{-1}-{\bm{C}}^{-2}$
(197) $\displaystyle=$
$\displaystyle({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}{\bm{B}}^{T}{\bm{B}}{\bm{C}}^{-1}+({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-1}{\bm{B}}^{T}{\bm{B}}{\bm{C}}^{-2}$
(198)

As such, by Eq. (195) we have

$\displaystyle|\operatorname{Tr}{\bm{R}}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}{\bm{B}}^{T}{\bm{B}}{\bm{C}}^{-1}{\bm{R}}|/N_{=}$
(199) $\displaystyle=$
$\displaystyle|\operatorname{Tr}{\bm{C}}^{-1}{\bm{R}}^{2}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}{\bm{B}}^{T}{\bm{B}}|/N_{=}$
(200) $\displaystyle\lesssim$ $\displaystyle
N_{<}/N_{=}\|{\bm{C}}^{-1}{\bm{R}}^{2}({\bm{C}}-{\bm{B}}^{T}{\bm{B}})^{-2}\|_{\mathrm{op}}\to
0\,.$ (201)

The other term can be bounded similarly. This finishes the proof of
$\Delta_{d}^{\prime}\to 0$ in probability. To sum up, we have the test error
to be

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},{\bm{h}})=\operatorname{Tr}\left({\bm{R}}^{2}{\bm{C}}^{-2}\right)/N_{=}\hat{F}_{r}^{2}+\left(1+m^{-1}\operatorname{Tr}\tilde{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}\right)\hat{F}_{>r}^{2}+$
(202) $\displaystyle
m^{-1}\sigma_{\epsilon}^{2}\operatorname{Tr}\tilde{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}+\Delta_{d}$
(203)

#### Generalization Error via Marchenko-Pastur

The next step is to reduce the traces to the integral form when $d\to\infty$.
That is evaluating the followings as $d\to\infty$,

$\displaystyle\operatorname{Tr}\left({\bm{R}}^{2}{\bm{C}}^{-2}\right)/N_{=}\quad\mathrm{and}\quad
m^{-1}\sigma_{\epsilon}^{2}\operatorname{Tr}\tilde{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$
(204)

We begin with the simpler case $E_{r}=1$ and then consider $E_{r}>1$.

####  The $E_{r}=1$ case.

I.e., there is only one eigenspace with eigenvalues $\sim d^{-s_{r}}$. This is
the case for one-hidden layer convolutional kernels and dot-product kernels.
In this case, $n=0$ and

$\displaystyle{\bm{R}}=(\overline{\xi}_{r}^{(d)})^{-1}{\bm{I}}_{N_{\leq}},\quad\mathrm{with}\quad\overline{\xi}_{r}^{(d)}=\gamma^{-1}m(\sigma_{rn}^{(d)})^{2}=\frac{m}{N_{=}}\hat{h}^{2}_{r}\gamma^{-1}\to\overline{\xi}_{r}=\alpha^{-1}\hat{h}^{2}_{r}\gamma^{-1}$
(205)

Choosing $\rho=(\overline{\xi}_{r}^{(d)})^{-1}$ and applying Theorem 1, we
have when444Note that $N_{\leq}=N_{=}(1+o(1))$ $N_{=}/m\to\alpha\in(0,\infty)$

$\displaystyle\overline{\xi}_{r}^{-2}\operatorname{Tr}\left({\bm{C}}^{-2}\right)/N_{=}$
$\displaystyle\longrightarrow\int(1+\overline{\xi}_{r}t)^{-2}\mu_{\alpha}(t)dt\,\,$
(206)
$\displaystyle\frac{N_{\leq}}{m}\frac{1}{N_{\leq}}\operatorname{Tr}\tilde{\bm{D}}^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$
$\displaystyle\longrightarrow\alpha\overline{\xi}_{r}^{2}\int{t}(1+\overline{\xi}_{r}t)^{-2}\mu_{\alpha}(t)dt$
(207)

Therefore,

$\displaystyle\mathrm{Err}({\bm{X}};\lambda,{\bm{F}},{\bm{h}})=\left(\hat{F}_{r}^{2}\cdot\int\frac{\mu_{\alpha}(t)}{(1+\overline{\xi}_{r}t)^{2}}dt+\hat{F}_{>r}^{2}\right)+\left(\hat{F}_{>r}^{2}+\sigma_{\epsilon}^{2}\right)\,\cdot\alpha\overline{\xi}_{r}^{2}\int\frac{t\mu_{\alpha}(t)}{(1+\overline{\xi}_{r}t)^{2}}dt+\Delta_{d}$
(208)

Both integrals have closed form formulas and they are computed in Sec.C.5.

#### The $E_{r}>1$ Case.

Recall that ${\bm{R}}$ is a diagonal matrix with entries
$\gamma(m(\sigma_{rn}^{(d)})^{2})^{-1}$ whose multiplicity is $N_{rn}^{(d)}$.
We assume the limiting density exist

$\displaystyle\gamma(m(\sigma_{rn}^{(d)})^{2})^{-1}\to\gamma\alpha\hat{h}_{rn}^{-2}\quad\mathrm{and}\quad
N_{rn}^{(d)}/\sum_{n\in[E_{k}]}N_{rn}^{(d)}\to\tau_{rn}$ (209)

and let $\nu_{\bm{h}}(r)$ denote this distribution. For convenience, we still
${\bm{R}}$ to represent a (sequence of) diagonal matrix with limiting spectral
$\nu_{\bm{h}}(r)$. By our assumptions on ${\bm{h}}$, the support of
$\nu_{\bm{h}}(r)$ is bounded away from $0$ and $\infty$. Thus, ignoring
vanishing correction term between
$\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$ and
$\overline{\bm{Z}}_{=}^{T}\overline{\bm{Z}}_{=}$, we need to compute the limit
of the following

$\displaystyle\frac{1}{N_{\leq}}\operatorname{Tr}{\bm{R}}^{2}({\bm{R}}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq})^{-2}={\bm{R}}^{1/2}(1+{\bm{R}}^{-1/2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}{\bm{R}}^{-1/2})^{-1}{\bm{R}}^{-1}(1+{\bm{R}}^{-1/2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}{\bm{R}}^{-1/2})^{-1}{\bm{R}}^{1/2}$
(210)
$\displaystyle\frac{1}{N_{\leq}}\operatorname{Tr}({\bm{R}}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq})^{-2}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}=\frac{1}{N_{\leq}}\operatorname{Tr}\left(({\bm{R}}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq})^{-1}-({\bm{R}}+\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq})^{-2}{\bm{R}}\right)$
(211)

.

To evaluate the limit, we may need extra assumptions on the eigenfunctions
$\phi_{knl}^{(d)}$ to ensure ${\bm{R}}$ and ${\bm{Z}}_{\leq}$ are
asymptotically free. Nevertheless, under the freeness assumption, computing
self-consistent equations that characterize the asymptotic values of the trace
objects in Eq. (210) and Eq. (211) is then straightforward using tools from
operator-valued free probability [Mingo and Speicher, 2017]. We do not
elaborate on the details here, but refer the reader so many related works for
examples of how to apply these tools [Far et al., 2006, Adlam and Pennington,
2020a, b, Tripuraneni et al., 2021a, b].

###  C.5 Computing the Integrals.

It remains to compute the above integrals. Note that

$\displaystyle\overline{\xi}_{r}^{2}\int\frac{t\mu_{\alpha}(t)}{(1+\overline{\xi}_{r}t)^{2}}dt=\overline{\xi}_{r}\left(\int\frac{\mu_{\alpha}(t)}{(1+\overline{\xi}_{r}t)}dt-\int\frac{\mu_{\alpha}(t)}{(1+\overline{\xi}_{r}t)^{2}}dt\right)$
(212)

As such we only need to compute, for $k=1$ and $k=2$,

$\displaystyle\zeta_{\alpha}(\xi;\alpha,k)=\int\frac{\mu_{\alpha}(t)}{(1+\xi
t)^{k}}dt$ (213)

Note that one only needs $\zeta_{\alpha}(\xi;\alpha,1)$ as
$\zeta_{\alpha}(\xi;\alpha,k)$ can be obtained from
$\zeta_{\alpha}(\xi;\alpha,k-1)$ by taking derivative w.r.t. $\overline{\xi}$.
Denote $b_{\pm}=(1\pm\sqrt{\alpha})^{2}$ and $\Delta=\alpha_{+}-\alpha_{-}$.
Then

$\displaystyle\mu_{\alpha}(t)=\left(1-\frac{1}{\alpha}\right)^{+}\delta_{0}(t)+\frac{\sqrt{(\alpha_{+}-t)(t-\alpha_{-})}}{2\pi\alpha
t}\bm{1}_{[\alpha_{-},\alpha_{+}]}(t)dt$ (214)

With
$b=(1+\overline{\xi}_{r}\alpha_{-})/(\overline{\xi}_{r}(\alpha_{+}-\alpha_{-}))$
and
$c=\overline{\xi}_{r}\alpha_{-}/(\overline{\xi}_{r}(\alpha_{+}-\alpha_{-}))=\alpha_{-}/(\alpha_{+}-\alpha_{-})$,

$\displaystyle\int(1+\overline{\xi}_{r}t)^{-k}\mu_{\alpha}(t)dt$ (215)
$\displaystyle=$
$\displaystyle\left(1-\frac{1}{\alpha}\right)^{+}+\int_{[\alpha_{-},\alpha_{+}]}(1+\overline{\xi}_{r}t)^{-k}\frac{\sqrt{(\alpha_{+}-t)(t-\alpha_{-})}}{2\pi\alpha
t}dt$ (216) $\displaystyle=$
$\displaystyle\left(1-\frac{1}{\alpha}\right)^{+}+\frac{1}{2\pi\alpha\overline{\xi}_{r}}(\overline{\gamma}(\alpha_{+}-\alpha_{-}))^{1-k}\int_{0}^{1}(b+t)^{-k}(c+t)^{-1}\sqrt{t(1-t)}dt$
(217)

Thanks to wolframalpha.com, we have, after doing some algebra,

$\displaystyle\int_{0}^{1}\frac{\sqrt{t(1-t)}}{(t+b)(t+c)}dt$
$\displaystyle=\pi\left(-1+\frac{1+b+c}{\sqrt{bc(1+b)(1+c)}}\right)$ (218)
$\displaystyle\int_{0}^{1}\frac{\sqrt{t(1-t)}}{(t+b)^{2}(t+c)}dt$
$\displaystyle=\frac{\pi}{2\sqrt{b^{2}+b}((b+c+2bc)+2\sqrt{(b+1)(c+1)bc})}$
(219)

###  C.6 Proof of Lemma 3.

Note that this lemma is trivial if $\lim_{d\to\infty}m/N_{\leq
r}^{(d)}=\alpha^{-1}>1$ as
$\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$ follows the Marchenko-
Pastur distribution, and the smallest eigenvalues is bounded from below by
$\alpha_{-}=(1-\sqrt{\alpha})^{2}$. When $\alpha^{-1}\leq 1$, we need to use
the regularization term $\overline{\bm{\Lambda}}_{\leq}^{-1}$. We provide the
details below.

Recall that
$\overline{\bm{Z}}_{\leq}^{T}=[\overline{\bm{Z}}_{<}^{T},\overline{\bm{Z}}_{r}^{T}]$,
where $\overline{\bm{Z}}_{<}^{T}$ is a $m\times N_{<r}^{(d)}$ matrix
consisting of low frequency modes and $\overline{\bm{Z}}_{r}^{T}$ is a
$m\times N_{r}^{(d)}$ is a matrix consisting of critical frequency modes. We
have $N_{<r}^{(d)}\sim d^{s_{r-1}}$, $N_{r}^{(d)}\sim N_{\leq r}^{(d)}\sim
d^{s_{r}}$ and

$\displaystyle\overline{\bm{Z}}_{<}^{T}\overline{\bm{Z}}_{<}={\bm{I}}_{N_{<r}^{(d)}}+\Delta_{d}$
(220)

where $\mathbb{E}\|\Delta_{d}\|_{op}\to 0$ as $d\to\infty$ in probability. Let
${\bm{u}}=[\beta_{<}{\bm{e}}_{<}^{T},\beta_{r}{\bm{e}}_{r}^{T}]^{T}$ be a unit
vector in ${\mathbb{R}}^{N_{\leq r}^{(d)}}$, where ${\bm{e}}_{<}$ and
${\bm{e}}_{r}$ are unit vectors in ${\mathbb{R}}^{N_{<r}^{(d)}}$ and
${\mathbb{R}}^{N_{r}^{(d)}}$ resp., and $\beta_{<}^{2}+\beta_{r}^{2}=1$. We
want to show that for some $\lambda_{1}>0$

$\displaystyle\lambda_{1}\leq{\bm{u}}^{T}D{\bm{u}}={\bm{u}}^{T}\overline{\bm{\Lambda}}_{\leq}^{-1}{\bm{u}}+{\bm{u}}^{T}\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}{\bm{u}}$
(221)

Note that the entries in $\overline{\bm{\Lambda}}_{\leq}^{-1}$ corresponding
to the critical-frequencies are $m(\sigma_{nr}^{(d)})^{2}\sim 1$. Thus there
is a constant $c>0$ such that

$\displaystyle{\bm{u}}^{T}\overline{\bm{\Lambda}}_{\leq}^{-1}{\bm{u}}\geq
c\beta_{r}^{2}$ (222)

In addition, if $C\equiv\|\overline{\bm{Z}}_{r}{\bm{e}}_{r}\|_{2}$ then $C\leq
2\alpha_{+}$ in probability. Thus by the triangle inequality,

$\displaystyle{\bm{u}}^{T}D{\bm{u}}$ $\displaystyle\geq
c\beta_{r}^{2}+\|\beta_{<}\overline{\bm{Z}}_{<}{\bm{e}}_{<}+\beta_{r}\overline{\bm{Z}}_{r}{\bm{e}}_{r}\|_{2}^{2}$
(223) $\displaystyle\geq
c\beta_{r}^{2}+(\|\beta_{<}\overline{\bm{Z}}_{<}{\bm{e}}_{<}\|_{2}-\|\beta_{r}\overline{\bm{Z}}_{r}{\bm{e}}_{r}\|_{2})^{2}$
(224) $\displaystyle\geq
c\beta_{r}^{2}+((1-\Delta_{d})|\beta_{<}|-C|\beta_{r}|)^{2}$ (225)

where $\Delta_{d}\to 0$ in probability. If
$C|\beta_{r}|<\frac{1}{2}|\beta_{<}|$, the above is greater than
$(1/2-\Delta_{d})\beta_{<}^{2}+c\beta_{r}^{2}\gtrsim 1$; otherwise
$C|\beta_{r}|\geq\frac{1}{2}|\beta_{<}|$ and ${\bm{u}}^{T}D{\bm{u}}\geq
c\beta_{r}^{2}\geq c(\frac{1}{2C}\beta_{<})^{2}$ and as a result

$\displaystyle{\bm{u}}^{T}D{\bm{u}}\geq
2c\beta_{r}^{2}/2\geq(c\beta_{r}^{2}+c(\frac{1}{2C}\beta_{<})^{2})/2\geq
c\min(1,\frac{1}{2C})^{2}/2\gtrsim 1$ (226)

The other direction is easier as both $\overline{\bm{\Lambda}}^{-1}_{\leq}$
and $\overline{\bm{Z}}_{\leq}^{T}\overline{\bm{Z}}_{\leq}$ have operator norms
bounded above.

###  C.7 Proof of Claim 1

The proof is split into two part: the ultra-high frequency parts $k\geq j_{0}$
and the median-high-frequency part, $r<k<j_{0}$. The first part is done by a
moment-based calculation and the second part is done by matrix concentration
[Vershynin, 2010].

#### Controlling the Ultra-high-frequency.

Recall that

$\displaystyle\Delta_{kn}^{(d)}\equiv\frac{1}{N_{kn}^{(d)}}Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}-{\bm{I}}_{m}\,.$
(227)

By Assumption (4.), the diagonals are zero and we have

$\displaystyle\Delta_{kn}^{(d)}$
$\displaystyle=[Z_{kn}({\bm{x}}_{i})^{T}Z_{kn}({\bm{x}}_{j})/N_{kn}^{(d)}]_{i,j\in[m],i\neq
j}$ (228)

Then

$\displaystyle\mathbb{E}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}^{2}\leq\mathbb{E}\|\Delta_{kn}^{(d)}\|_{\mathrm{F}}^{2}=$
$\displaystyle\mathbb{E}\sum_{i\neq
j}|Z_{kn}({\bm{x}}_{i})^{T}Z_{kn}({\bm{x}}_{j})/N_{kn}^{(d)}|^{2}$ (229)
$\displaystyle=$
$\displaystyle\frac{1}{(N_{kn}^{(d)})^{2}}\mathbb{E}\sum_{i\neq
j}\sum_{l,l^{\prime}}\phi^{(d)}_{knl}({\bm{x}}_{i})\phi^{(d)}_{knl}({\bm{x}}_{j})\phi^{(d)}_{knl^{\prime}}({\bm{x}}_{i})\phi^{(d)}_{knl^{\prime}}({\bm{x}}_{j})$
(230) $\displaystyle=$
$\displaystyle\frac{1}{(N_{kn}^{(d)})^{2}}\mathbb{E}\sum_{i\neq
j}\sum_{l}\phi^{(d)}_{knl}({\bm{x}}_{i})^{2}\phi^{(d)}_{knl}({\bm{x}}_{j})^{2}$
(231) $\displaystyle=$
$\displaystyle\frac{1}{N_{kn}^{(d)}}m(m-1)\leq\frac{1}{N_{kn}^{(d)}}m^{2}$
(232)

Recall that $E_{k}$ grows at most exponentially, i.e. $E_{k}\leq C^{k}$ for
some constant $C$. Thus, choosing $d$ large enough such that
$Cd^{-\delta_{0}/4}<1$ and summing over
$k>j_{0}\equiv[4s_{r}/\delta_{0}+4]+4$,

$\displaystyle\mathbb{E}\sum_{k>j_{0}}\sum_{n\in
E_{k}}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}\lesssim$
$\displaystyle\sum_{k>j_{0}}C^{k}(m^{2}/N_{kn}^{(d)})^{1/2}$ (233)
$\displaystyle\lesssim$ $\displaystyle\sum_{k>j_{0}}C^{k}d^{-s_{k}/2+s_{r}}$
(234) $\displaystyle\leq$
$\displaystyle\sum_{k>j_{0}}C^{k}d^{-k\delta_{0}/2+s_{r}}$ (235)
$\displaystyle\lesssim$
$\displaystyle\sum_{k>j_{0}}d^{-k\delta_{0}/4+s_{r}}\lesssim
d^{-j_{0}\delta_{0}/4+s_{r}}\lesssim d^{-\delta_{0}}$ (236)

#### Controlling the Median-high-frequency.

It remains to show, for some $\epsilon>0$

$\displaystyle\mathbb{E}\sum_{r<k\leq j_{0}}\sum_{n\in
E_{k}}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}\lesssim d^{-\epsilon}\,.$ (237)

As there are only finitely many terms in this sum, we only need to show that
for each $k$ and $n$,

$\mathbb{E}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}\lesssim d^{-\epsilon}\,.$

We use the following theorem from Vershynin [2010] regarding matrix
concentration to prove this claim.

######  Theorem 6 (Theorem 5.62 Vershynin [2010]).

Let ${\bm{A}}$ be a $N\times m$ ($N\geq m$) matrix whose columns $A_{j}$ are
independent isotropic random vectors in ${\mathbb{R}}^{N}$ with
$\|A_{j}\|_{2}=N$ a.s. Let $K$ be defined as

$\displaystyle K=\frac{1}{N}\mathbb{E}\max_{j\leq m}\sum_{i\in[m],i\neq
j}|A_{j}^{T}A_{i}|^{2}$ (238)

Then

$\displaystyle\mathbb{E}\|{\bm{A}}^{T}{\bm{A}}/N-{\bm{I}}_{m}\|_{\mathrm{op}}\lesssim\sqrt{K\log(m)/N}$
(239)

We apply this theorem to ${\bm{A}}=Z_{kn}({\bm{X}})^{T}$. The columns of
$Z_{kn}({\bm{X}})^{T}$ are $A_{j}=\phi^{(d)}_{kn}({\bm{x}}_{j})$, $j\in[m]$
which are independent. Let $N=N_{kn}^{(d)}$. By Assumption (4.),

$\displaystyle A_{j}^{T}A_{j}=\sum_{l}\phi^{(d)}_{knl}({\bm{x}}_{j})^{2}=N.$
(240)

We claim that $K\lesssim_{k,q}m^{1+\frac{1}{q}}$ for any $q\geq 1$. Indeed,
let

$\displaystyle B_{j}=\sum_{i\in[m],i\neq j}|A_{j}^{T}A_{i}|^{2}$ (241)

We then remove the maximal function by paying an $m^{1/q}$ factor

$\displaystyle K=\frac{1}{N}\mathbb{E}\max_{j\leq
m}B_{j}\leq\frac{1}{N}m^{1/q}|\mathbb{E}B_{j}^{q}|^{1/q}=\frac{1}{N}m^{1/q}\left|\mathbb{E}(\sum_{i\in[m],i\neq
j}|A_{j}^{T}A_{i}|^{2})^{2q/2}\right|^{1/q}$ (242)

Next we apply the Minkowski inequality to swap the $L^{2q}$-norm and the
$l^{2}$-norm,

$\displaystyle K\leq\frac{1}{N}m^{1/q}\sum_{i\in[m],i\neq
j}(\mathbb{E}|A_{j}^{T}A_{i}|^{2q})^{1/2q\times
2}\leq\frac{1}{N}m^{1/q+1}(\mathbb{E}|A_{j}^{T}A_{i}|^{2q})^{1/2q\times 2}\leq
C_{k,q}^{2}m^{1/q+1}$ (243)

if $(\mathbb{E}|A_{j}^{T}A_{i}|^{2q})^{1/2q\times 2}\leq C_{k,q}^{2}N$, which
is done by hypercontractivities below. Indeed, for ${\bm{x}}_{j}$ fixed,
$Z_{rn}({\bm{x}}_{j})^{T}Z_{rn}({\bm{x}}_{i})$ is a linear combination of
$\phi_{knl}^{(d)}$, Assumption (2.) gives

$\displaystyle\mathbb{E}_{{\bm{x}}_{i}}|A_{j}^{T}A_{i}|^{2q}=$
$\displaystyle\mathbb{E}_{{\bm{x}}_{i}}|Z_{rn}({\bm{x}}_{j})^{T}Z_{rn}({\bm{x}}_{i})|^{2q}$
(244) $\displaystyle\leq$
$\displaystyle\left(C_{k,q}(\mathbb{E}_{{\bm{x}}_{i}}|Z_{rn}({\bm{x}}_{j})^{T}Z_{rn}({\bm{x}}_{i})|^{2})^{1/2}\right)^{2q}$
(245) $\displaystyle=$ $\displaystyle C_{k,q}^{2q}N^{q}$ (246)

where we applied

$\displaystyle\mathbb{E}_{{\bm{x}}_{i}}|Z_{rn}({\bm{x}}_{j})^{T}Z_{rn}({\bm{x}}_{i})|^{2}=$
$\displaystyle\mathbb{E}_{{\bm{x}}_{i}}|\sum_{l}\phi_{knl}^{(d)}({\bm{x}}_{i})\phi_{knl}^{(d)}({\bm{x}}_{j})|^{2}$
(247) $\displaystyle=$
$\displaystyle\mathbb{E}_{{\bm{x}}_{i}}\sum_{ll^{\prime}}\phi_{knl}^{(d)}({\bm{x}}_{i})\phi_{knl}^{(d)}({\bm{x}}_{j})\phi_{knl^{\prime}}^{(d)}({\bm{x}}_{i})\phi_{knl^{\prime}}^{(d)}({\bm{x}}_{j})$
(248) $\displaystyle=$
$\displaystyle\sum_{l}\phi_{knl}^{(d)}({\bm{x}}_{j})^{2}$ (249)
$\displaystyle=$ $\displaystyle N$ (250)

Therefore with ${\bm{A}}=Z_{kn}({\bm{X}})^{T}$, we have

$\displaystyle\mathbb{E}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}=\mathbb{E}\|Z_{kn}({\bm{X}})Z_{kn}({\bm{X}})^{T}/N_{kn}^{(d)}-{\bm{I}}_{m}\|_{\mathrm{op}}\lesssim_{k,q}\sqrt{m^{1+1/q}\log
m/N_{kn}^{(d)}}\,.$ (251)

For each fixed $k>r$, $s_{k}-s_{r}\geq\delta_{0}$, by choosing $q$
sufficiently large (depending on $k$ and $\delta_{0}$), we have

$\displaystyle\mathbb{E}\|\Delta_{kn}^{(d)}\|_{\mathrm{op}}\lesssim_{k}d^{-\delta_{0}/2}\,.$
(252)

##  Appendix D Additional Plots

To simulate the learning curves for higher-order scalings, e.g. $r=4$, we must
chose $d$ small. As such, we are in a strong finite-size correction regime. In
this section, we vary $d$ to visualize the finite-size effect of the
predictions. Note that for larger $d$ ($=60$ here), we can only simulate up to
the quadratic scaling. For smaller $d$ ($d=10$), we observe noticeable finite-
size correction. However, the predictions match the simulations quite well.
When $d$ become larger $d=60$, the predicted learning curves match the
simulation almost perfectly.

Figure 5: Tiny $d=10$

Figure 6: Small $d=20$

Figure 7: Large $d=60$