An Equivalence Principle for the Spectrum of Random Inner-Product Kernel Matrices with Polynomial Scalings

Introduction

Our study of this model is motivated by recent problems in machine learning, statistics, and signal processing, where random matrices like (1.1) and their spectral properties play crucial roles. Examples include kernel methods (such as kernel-PCA and kernel-SVM ), covariance thresholding procedures , nonlinear dimension reduction , and probabilistic matrix factorization . Moreover, the closely-related non-Hermitian version of (1.1), where $A_{ij}=\frac{1}{\sqrt{n}}f_{d}(\sqrt{d}\boldsymbol{x}_{i}\boldsymbol{y}_{j})$ for two sets of vectors $\{\boldsymbol{x}_{i}\}_{i\leq n}$ and $\{\boldsymbol{y}_{j}\}_{j\leq p}$, appears in the random feature model , an interesting theoretical model for large random neural networks.

The remainder of the paper is organized as follows. In Section 2, we begin by examining the special case of polynomial kernel functions, with the equivalence principle formalized in Theorem 1 and a heuristic explanation provided in Section 2.2. We address the case of more general nonlinear functions in Section 2.3, stating the corresponding asymptotic characterizations in Theorem 2. Section 2.4 presents several numerical experiments to illustrate our theory, while Section 3 discusses related work in the literature. We dedicate Section 4 to the proof of Theorem 1 and gather auxiliary results in the appendix. In Section 5, we extend our findings to cases where data vectors are sampled from the isotropic Gaussian distribution rather than the spherical distribution. Finally, we conclude the paper in Section 6, discussing potential extensions of our results and open problems.

Before delving into the technical details, we first establish some notations employed throughout this paper.

Additionally, we assume that $\mathinner{\!\left\lvert\mathfrak{Re}(z)\right\rvert}\leq\tau^{-1}$ and $\tau\leq\eta\leq\tau^{-1}$, for some global constant $\tau>0$.

Probability distributions: $\mathsf{Unif}(\mathcal{S}^{d-1})$ denotes the uniform probability measure on $\mathcal{S}^{d-1}$. For two independent vectors $\boldsymbol{x},\boldsymbol{y}\sim\mathsf{Unif}(\mathcal{S}^{d-1})$, we define

Stochastic order notation: In our proof, we utilize a concept of high-probability bounds known as stochastic domination. This notion, first introduced in , provides a convenient way to account for low-probability exceptional events where some bounds may not hold. Consider two families of nonnegative random variables:

where $U^{(d)}$ is a possibly $d$-dependent parameter set. We say that $X$ is stochastically dominated by $Y$, uniformly in $u$, if for every (small) $\varepsilon>0$ and (large) $D>0$ we have

for sufficiently large $d\geq d_{0}(\varepsilon,D)$. If $X$ is stochastically dominated by $Y$, uniformly in $u$, we use the notation $X\prec Y$. Moreover, if for some complex family $X$ we have $\mathinner{\!\left\lvert X\right\rvert}\prec Y$, we also write $X=\mathcal{O}_{\prec}(Y)$. This stochastic order notation should not be mistaken for the conventional big $\mathcal{O}$ notation, which we will also use in this paper: for two deterministic sequences $X^{(d)}$ and $Y^{(d)}$, we write $X=\mathcal{O}_{\alpha}(Y)$, with some parameter $\alpha$, if $\mathinner{\!\left\lvert X\right\rvert}\leq C(\alpha)Y$ for all sufficiently large $d$. Here, the constant $C(\alpha)$ may depend on $\alpha$.

An Asymptotic Equivalence Principle

where $\left\{\mu_{k}\right\}_{k}$ is a set of expansion coefficients, and $q_{k}^{(d)}(x)$ denotes the $k$th Gegenbauer polynomial . The Gegenbauer polynomials, also known as the ultraspherical polynomials in the literature , form a set of orthogonal polynomial basis with respect to the probability measure $\mathsf{s}_{d}^{1}$ defined in (1.3). Specifically, we have $\mathsf{deg}\,q_{k}^{(d)}(x)=k$ and

where $\mathds{1}_{ij}$ is the Kronecker delta.

The coefficients of these polynomials can be determined by performing the Gram-Schmidt procedure on the monomial basis $\left\{0,1,x,x^{2},\ldots\right\}$ and by using the explicit formula for the moments of $\xi$ [see (D.5)]. The first few polynomials in the sequence are

Note that the coefficients of the Gegenbauer polynomials depend on the dimension $d$. To simplify the notation, we will suppress this dependence by writing throughout the paper

In Appendix A, we compile a list of properties of Gegenbauer polynomials that will be used in our analysis.

where $\left\{\boldsymbol{x}_{i}\right\}_{i\in[n]}$ is a collection of independent vectors drawn from $\mathsf{Unif}(\mathcal{S}^{d-1})$. Considering (2.1), we can express the inner-product kernel matrix in (1.1) as

The first component in the sum in (2.8) is a (shifted) Wishart matrix, constructed as follows:

For $z\in C_{+}$, let $s_{A}(z)$ and $s_{B}(z)$ denote the Stieltjes transforms of the ESDs of $A$ and $B$, respectively. The following theorem formalizes the asymptotic equivalence of the matrix models $A$ and $B$.

As $B$ is the linear combination of a (shifted) Wishart matrix and independent GOE matrices, its limiting eigenvalue density is given by an additive free convolution of a Marchenko-Pastur law with a semicircle law. Explicit formulas for the limiting density function

2 A Heuristic Explanation of the Equivalence Principle

The equivalence principle stated above has a simple heuristic explanation. To understand this, let us first revisit a crucial property of Gegenbauer polynomials. Given any $\boldsymbol{x}_{i},\boldsymbol{x}_{j}\in\mathcal{S}^{d-1}$,

where $\left\{Y_{k,a}(\boldsymbol{x}_{i})\right\}_{a\in[N_{k}]}$ represents a collection of orthonormal degree-$k$ spherical harmonics associated with $\boldsymbol{x}_{i}$, and

denotes the cardinality of the set. For a comprehensive explanation and the exact expression for $N_{k}$, we refer the reader to Appendix A. The above identity is valuable because it allows us to “linearize” the term $q_{k}(\sqrt{d}\,\boldsymbol{x}^{\mkern-1.5mu\mathsf{T}}\boldsymbol{y})$, transforming it into an inner product of two $N_{k}$-dimensional vectors comprised of spherical harmonics.

Using (2.13) and another identity, $q_{k}(\sqrt{d})=\sqrt{N_{k}}$, we can express the matrices $\left\{A_{k}\right\}$ in a factorized form:

where $S_{k}$ is an $N_{k}\times n$ matrix with entries consisting of the spherical harmonics, that is,

Heuristically, if we replace the entries of $S_{k}$ with i.i.d. standard Gaussians, we can expect the limiting ESD of $A_{k}$ to be characterized by the Marchenko-Pastur (MP) law with an aspect ratio parameter

Notice that this expression has the form of the matrix presented in (C.3) (in Appendix C). By further assuming that the family ${Y_{k,a}(\boldsymbol{x})}$ consists of not just uncorrelated but indeed independent random variables, $A$ adheres to the MP law in (C.4), meaning its limiting Stieltjes transform, denoted by $m(z)$, satisfies the equation:

3 General Nonlinear Kernels

In this subsection, we extend the equivalence principle established in Theorem 1 to encompass cases where the kernel $f_{d}$ in (1.1) is a general function, going beyond merely polynomials.

where $g\sim\mathcal{N}(0,1)$ and $\mathds{1}_{ij}$ is the Kronecker delta. The first four (normalized) Hermite polynomials are

where $\max_{i}\mathinner{\!\left\lvert c_{i,k}(d)\right\rvert}=\mathcal{O}_{k}(1/d)$.

In the subsequent discussion, we will establish an equivalence principle for the matrix in (1.1), given the following assumption on the function $f_{d}$.

Let $g\sim\mathcal{N}(0,1)$, and let $\left\{h_{k}(x)\right\}_{k}$ represent the Hermite polynomials defined above. The function $f_{d}$ in (1.1) satisfies the following conditions:

for some finite numbers $\left\{\mu_{k}\right\}$.

The sequence $\left\{\mu_{k}\right\}$ is square-summable, i.e.,

Let $w(x)$ denote the probability density functions of $\mathcal{N}(0,1)$, and let $w_{d}(x)$ denote the density function of the probability measure $\mathsf{s}_{d}^{1}$, as defined in (1.3). We have

When $f_{d}(x)=f(x)$ is a function that is independent of $d$, the following lemma offers simple sufficient conditions that can be used to verify that Assumption 1 holds.

Recall from Section 2.2 that the self-consistent equation (2.29) characterizes the free additive convolution of a (shifted) MP law and the semicircle law. Consequently, the statement of Theorem 2 can still be interpreted in terms of an equivalence principle: the limiting ESD of $A$ is equivalent to that of

Additionally, by using [9, Lemma C.1], we can conclude from the conditions (2.25) and (2.28) that

Given the $L$ found above, we define two functions

where $\hat{\mu}_{L}\coloneqq(\sigma^{2}-\sum_{0\leq k\leq L-1}\mu_{k}^{2})^{1/2}$. By using the property that the Gengenbauer polynomials are orthonormal [see (2.2)], we can write

To reach the inequality in (2.35), we have employed the characterizations presented in (2.30) and (2.31), which guarantee that the inequality holds for all $d\geq d(L,c)$, where $d(L,c)$ is some integer that may depend on $L$ and $c$. By (2.32), $\hat{\mu}_{L}^{2}\leq\eta^{4}c^{2}/24$ and $\mu_{L}^{2}\leq\eta^{4}c^{2}/24$. We can then further bound the right-hand side of (2.35) as

Let $\hat{A}$ be a matrix constructed according to (1.1) but with the function $f_{d}$ replaced by $\hat{f}_{d}$. Let $s_{\hat{A}}$ denote its Stieltjes transform. We can characterize $s_{\hat{A}}$ in two ways. On the one hand, since $\hat{f}_{d}$ is a linear combination of $L+1$ Gegenbauer polynomials, we can apply Theorem 1 to get

where the second inequality follows from (2.36). By the triangular inequality and estimates in (2.37) and (2.38), we have

for all sufficiently large $d$. Applying the Borel-Cantelli lemma (for a fixed $D>1$), we can then conclude that $s_{A}(z)$ converges to $m(z)$ almost surely. ∎

4 Numerical Experiments

In this subsection, we present numerical experiments that demonstrate the equivalence principle as stated in Theorems 1 and 2. We first examine the empirical spectral distribution (ESD) of individual polynomial component matrices $A_{k}$, as defined in (2.6). In the quadratic scaling regime, where $n$ is asymptotically proportional to $d/^{2}$ for a fixed /, $A_{2}$ is asymptotically equivalent to $B_{2}$ in (2.9). Consequently, the ESD of $A_{2}$ converges to the MP law, characterized by a density function given in (C.6). Both $A_{3}$ and $A_{4}$ are asymptotically equivalent to a GOE matrix, and their ESDs therefore converge to the standard semicircle law, as outlined in (2.19). As illustrated in Figure 1, the ESDs of $A_{2}$, $A_{3}$, and $A_{4}$ closely align with their respective limiting spectral densities.

In the second example, we consider linear combinations of the polynomial matrices $\left\{A_{k}\right\}$. Theorem 1 states that the ESD of any fixed linear combination of $A_{k}$ can be obtained by a free additive convolution between the MP law and the semicircle law. In Figure 2(a), we plot the ESD of

against the limiting spectral density. The theoretical curve (red solid line in the figure) is computed by solving the self-consistent equation (2.10) and then evaluating (2.12) numerically. In Figure 2(b), we show the ESD of another matrix in the form of

Here, $A_{5}^{M}$ and $A_{6}^{M}$ are the truncated version of $A_{5}$ and $A_{6}$, respectively. Specifically, for $M>0$, we have

and $A_{6}^{M}$ is defined similarly. Note that we apply this truncation to remove a small number of outliers in the entries of $A_{5}$ and $A_{6}$ that have very large magnitudes. The presence of these outliers intensifies the finite-size effect, causing the ESD to deviate from the limiting spectral density when the dimension $d$ is not very large. Due to the truncation step, the matrix in (2.42) is no longer a finite linear combination of $\left\{A_{k}\right\}$. Consequently, we employ Theorem 2 to compute the theoretical curve.

In the last example, we consider the ESD of the matrix in (1.1), where the nonlinear function is a soft-thresholding operator, i.e.,

Related Work

In this section, we discuss several related lines of work in the literature.

The linear asymptotic regime of the model in (1.1) was first investigated by Cheng and Singer , who established the weak limit of the ESD of the kernel matrix. The idea of “decorrelation” by expanding the nonlinear kernel function using an orthogonal polynomial basis was first introduced in that work and plays a significant role in our paper. The results of were obtained for data vectors sampled from isotropic Gaussian and spherical distributions. Do and Vu extended these results to other distributions (such as the Bernoulli case). The spectral norm of the kernel matrix was explored by Fan and Montanari , who also made the observation that the limiting distribution obtained in is the free additive convolution of an MP law and a semicircle law.

Asymptotics with polynomial scaling: Studies of high-dimensional statistical problems often focus on the linear asymptotic regime. In random matrix theory, the more general polynomial asymptotic regime was investigated by Bloemendal et al. , who demonstrated the local MP law for sample covariance matrices $X^{\mkern-1.5mu\mathsf{T}}X$, where $X$ is an $N\times n$ matrix with $\log N\asymp\log n$. In that work, the matrix $X$ is assumed to have independent entries. Although the kernel matrix in our work can also be written in the form of a (generalized) sample covariance matrix [see (2.15) and (2.20)], the matrix entries in our problem consist of spherical harmonics, which are uncorrelated but dependent random variables. As a result, the technical approach of cannot be directly applied here. In the context of kernel methods, exact asymptotics in the polynomial scaling regime were examined in Opper and Urbanczik and Bordelon et al. using nonrigorous statistical physics methods. The heuristics employed by these authors were similar to the scheme outlined in Section 2.2, specifically, treating the uncorrelated spherical harmonics as if they were independent standard normal random variables.

Non-Hermitian ensembles and random feature models: Lastly, we note that it is possible to extend the current study to a non-Hermitian version of (1.1). In this case, one would investigate an $n\times p$ matrix, with entries given by

Proof of the Main Result

This section is devoted to the proof of Theorem 1. We start by presenting a high-level outline of the proof in Section 4.1. The technical details are given in Sections 4.2 and 4.3, and in the appendix.

Note that, rather than mapping the indices of $A^{[i]}$ to $[n-1]$, we use $a,b\in[n]\setminus\left\{i\right\}$ to index the entries of $A^{[i]}$. Let $G(z)=(A-zI)^{-1}$ and $G^{[i]}(z)=(A^{[i]}-zI)^{-1}$ be the resolvents of $A$ and $A^{[i]}$, respectively.

We start our proof by applying Schur’s complement formula, which gives us

Recall that the Stieltjes transform of $A$ can be obtained as $s(z)=\frac{1}{n}\sum_{i}G_{ii}(z)$. The main technical step in our proof is to show that

where $\gamma_{a},\gamma_{b},\gamma_{c}$ are the three constants defined in Theorem 1. We will prove this estimate in Section 4.3 (see Proposition 2).

Using (4.3), we can now split the right-hand side of (4.2) into a leading term and an error term. Multiplying both sides of (4.2) by $G_{ii}(z)$ and averaging over $i$, we get

By (F.2) in the first step, and by applying (4.3) and the union bound in the second step,

Observe that the left-hand side of (4.4) has exactly the same functional form as the left-hand side of (2.10). The error estimate in (4.5) then implies that the Stieltjes transform $s(z)$ approximately satisfies the nonlinear equation in (2.10). By analyzing the uniqueness and stability of the solution $m(z)$ to this nonlinear equation (see Proposition 9 in Appendix G), we can then conclude that $s(z)\approx m(z)$.

2 Schur Complement: Reparameterization

We devote this and the next subsections to establishing the estimate in (4.3). Observe that, by construction, different coordinates of $A$ are statistically exchangeable. Thus, we only need to show (4.3) for a single index $i\in[n]$. We choose to do so for $i=n$.

and the entries of the minor $A^{[n]}$ are of the form

The random variables $\{A_{na}\}$ and $\{A^{[n]}_{ab}\}$ are weakly correlated. To make this correlation explicit, we use the following reparameterization. For every $a\in[n-1]$, let

By definition, $(R_{n}R_{n}^{\mkern-1.5mu\mathsf{T}}+\boldsymbol{x}_{n}\boldsymbol{x}_{n}^{\mkern-1.5mu\mathsf{T}})=I$, and hence $\boldsymbol{x}_{a}^{\mkern-1.5mu\mathsf{T}}\boldsymbol{x}_{b}=\boldsymbol{x}_{a}^{\mkern-1.5mu\mathsf{T}}(R_{n}R_{n}^{\mkern-1.5mu\mathsf{T}}+\boldsymbol{x}_{n}\boldsymbol{x}_{n}^{\mkern-1.5mu\mathsf{T}})\boldsymbol{x}_{b}$ for $a,b\in[n-1]$. We can then verify that $\sqrt{d}\,\boldsymbol{x}_{a}^{\mkern-1.5mu\mathsf{T}}\boldsymbol{x}_{b}=r(\xi_{a})r(\xi_{b})\sqrt{d-1}\widetilde{\boldsymbol{x}}_{a}^{\mkern-1.5mu\mathsf{T}}\widetilde{\boldsymbol{x}}_{b}+\xi_{a}\xi_{b}/\sqrt{d}$, where

The advantage of the new parameterizations in (4.9) and (4.11) is that the families $\left\{\xi_{a}\right\}$ and $\left\{\widetilde{\boldsymbol{x}}_{a}\right\}$ are independent, as we show below.

$\left\{\xi_{a}\right\}_{a<n}$ is an i.i.d. family with $\xi_{a}\sim\mathsf{s}_{d}^{1}$, where $\mathsf{s}_{d}^{1}$ is the probability measure defined in (1.3). $\left\{\widetilde{\boldsymbol{x}}_{a}\right\}_{a<n}$ is an i.i.d. family of random vectors with $\widetilde{\boldsymbol{x}}_{a}\sim\mathsf{Unif}(\mathcal{S}^{d-2})$. Moreover, $\left\{\xi_{a}\right\}$ and $\left\{\widetilde{\boldsymbol{x}}_{a}\right\}$ are mutually independent.

Fix $\boldsymbol{x}_{n}$, and let $E=[\boldsymbol{x}_{n}\,R_{n}]$. By construction, $E$ is an orthogonal matrix. Define

As $\boldsymbol{x}_{a}\sim_{\text{i.i.d.}}\mathsf{Unif}(\mathcal{S}^{d-1})$ and $E$ is independent of $\left\{\boldsymbol{x}_{a}\right\}_{a<n}$, we can conclude from the orthogonal invariant property of the spherical distribution that $\boldsymbol{v}_{a}\sim_{\text{i.i.d.}}\mathsf{Unif}(\mathcal{S}^{d-1})$. The statement of the lemma then follows from the observation that

where $\boldsymbol{v}_{a,\setminus 1}$ is the $(d-1)$-dimensional vector obtained from $\boldsymbol{v}_{a}$ after removing its first element. Then $\xi_{a}\sim\mathsf{s}_{d}^{1}$. Meanwhile, $\widetilde{\boldsymbol{x}}_{a}$ is uniformly distributed over $\mathcal{S}^{d-2}$, and is independent of $\xi_{a}$. Note that the above derivations are done for fixed $\boldsymbol{x}_{n}$ (i.e. when conditioned on $\boldsymbol{x}_{n}$). However, since the distributions of $\left\{\xi_{a}\right\}$ and $\left\{\widetilde{\boldsymbol{x}}_{a}\right\}$ are invariant to the choice of $\boldsymbol{x}_{n}$, we can conclude that $\left\{\xi_{a}\right\}$ and $\left\{\widetilde{\boldsymbol{x}}_{a}\right\}$ are indeed independent of $\boldsymbol{x}_{n}$. ∎

In what follows, we will use $\mathcal{E}_{t}(\xi_{1},\xi_{2})$ to denote any function in $\mathcal{E}_{t}$. The exact form of $\mathcal{E}_{t}(\xi_{1},\xi_{2})$ can change from one expression to another.

where $(k)_{t}=k!/(k-t)!$ denotes the falling factorial, and $\widetilde{q}_{k}(x)$ is the $k$th Gegenbauer polynomial in dimension $d-1$ as defined in (2.5).

For the first term on the right side of (4.13) with $t=0$, we define an $(n-1)\times(n-1)$ matrix $\widetilde{A}_{k}$, whose entries are given by:

For the first term on the right side of (4.13) with $t=k$, we define the vector

Using the expansion in (4.13), we can now write the matrix in (4.11) as

where $\Delta_{1},\ldots,\Delta_{5}$ are matrices defined as follows: for $a,b\in[n-1]$,

Later, we will demonstrate that the matrices $\left\{\Delta_{i}\right\}$ can all be considered as small noise terms that become negligible as $d\to\infty$. Moreover, by the construction of $\widetilde{A}_{k}$ in (4.14) and by Lemma 3, both $\widetilde{A}$ and its resolvent $\widetilde{G}(z)\coloneqq(\widetilde{A}-zI)^{-1}$ are independent of $\left\{\xi_{a}\right\}$. These observations lead us to consider the following expansion.

Let $G^{[n]}(z)$ and $\widetilde{G}(z)$ denote the resolvents of $A^{[n]}$ and $\widetilde{A}$, respectively. We have

where the second equality follows from the Woodbury matrix identity. By (4.9) and (4.16),

Substituting this into (4.22) then leads to the desired result. ∎

3 Schur Complement: the Leading Term

Next, we demonstrate that the right-hand side of (4.19) can indeed be separated into a leading term, whose form is provided in (4.3), and a small, random error term.

where $s(z)$ is the Stieltjes transform of the spectrum of $A$.

To prepare for the proof of Proposition 2, we first establish several intermediate results.

Let $\widetilde{s}(z)\coloneqq\frac{1}{n}\operatorname{tr}\widetilde{G}$ and $s^{[n]}(z)\coloneqq\frac{1}{n}\operatorname{tr}G^{[n]}$ be the Stieltjes transforms of $\widetilde{G}$ and $G^{[n]}$, respectively. We will establish (4.24) in three steps, by showing (a) $\chi_{\alpha\beta}(z)\approx\widetilde{s}(z)\delta_{\alpha\beta}$; (b) $\widetilde{s}(z)\approx s^{[n]}(z)$; and (c) $s^{[n]}(z)\approx s(z)$.

For step (a), we first recall the Ward identity (F.3), which gives us

where the inequality is due to (F.4). By the definitions presented in (4.16) and (4.20),

Applying (D.64b) in Proposition 7 with $\Psi_{c}=1/{\eta}$, we have

For step (b), we use the expansion formula in (4.17). By (F.8), (F.7), and the triangular inequality,

To bound $\Delta_{5}$, we first check from the definition in (4.12) that $\max_{a,b\in[n-1]}\mathinner{\!\left\lvert\mathcal{E}_{t}(\xi_{a},\xi_{b})\right\rvert}\prec 1/d$. Moreover, by its construction in (4.10), the function $r(\xi)\leq 1+\mathcal{O}(1/d)$ for all $\xi$. Combining these estimates then gives us $\mathinner{\!\left\lVert\Delta_{5}\right\rVert}_{\infty}\prec\frac{1}{d\sqrt{n}}$.

The matrix $\Delta_{2}$ in (4.18b) requires some additional care, as it contains terms with index $t=0$. First write

Applying the high-probability bound (D.45) and the property in (4.27), we get

which then implies that $\mathinner{\!\left\lVert\Delta_{2}\right\rVert}_{\infty}\prec\frac{1}{d\sqrt{n}}$. Substituting these bounds on $\mathinner{\!\left\lVert\Delta_{c}\right\rVert}_{\infty}$ for $1\leq c\leq 5$ into (4.26) then leads to

Now we move to step (c), where we compare $s^{[n]}(z)$ with $s(z)$. By the interlacing properties of the eigenvalues of $A^{[n]}$ and $A$, we have the following standard result (see [16, Lemma 7.5] for a proof):

where $C$ is an absolute constant. Finally, the statement of the lemma can be obtained by applying the triangular inequality and by using the estimates in (4.25), (4.29), and (4.30). ∎

Next, we show that the terms involving $\left\{\Delta_{c}\right\}_{1\leq c\leq 5}$ on the right-hand side of (4.19) are small. This is done in the following two lemmas.

Let $\alpha,\beta,k$ be integers such that $0\leq\alpha,\beta,k\leq L$ and $\alpha\neq\beta$. We have

where $\Lambda_{\beta}(\xi)$ denotes a diagonal matrix defined as

Write $u\coloneqq\boldsymbol{v}_{\alpha}^{\mkern-1.5mu\mathsf{T}}(\xi)\widetilde{G}(z)\Lambda_{\beta}(\xi)\widetilde{A}_{k}$. This is a vector with $n-1$ entries. In our proof, we will show that

of which (4.31) is then an immediate consequence. For each $a\in[n-1]$, define a matrix

We start by proving (4.33) for the special case when $\alpha>\beta=0$. Recall that $\boldsymbol{v}_{0}(\xi)=\boldsymbol{1}$. We then have $\mathinner{\!\left\lVert C_{a}\boldsymbol{v}_{0}(\xi)\right\rVert}\leq\|\widetilde{G}\|_{\mathsf{op}}\|\widetilde{A}_{k}\|_{\infty}\sqrt{n}\prec 1/\eta$, where the second step uses (F.1) and (D.44). Since $C_{a}$ is independent of $\boldsymbol{v}_{\alpha}(\xi)$, we can apply (D.64a) in Proposition 7 with $\Psi_{b}=1/\eta$ to get

By the same arguments, $\boldsymbol{v}_{0}^{\mkern-1.5mu\mathsf{T}}(\xi)C_{a}\boldsymbol{v}_{\beta}(\xi)\prec 1/\eta$ for $\beta>0$. In what follows, we assume $1\leq\alpha,\beta\leq L$ with $\alpha\neq\beta$. Using the Ward identity (F.3) in the second step, and applying (F.4), (D.44) in the third step, we get

By appealing to (D.64b) in Proposition 7 with $\Psi_{c}=1/\eta$, we have

One can check that the high-probability bounds in (4.34) and (4.35) are both uniform in $a\in[n-1]$. (To see this, note that we always bound $\max_{i}|(\widetilde{A}_{k})_{ia}|$ via $\|\widetilde{A}_{k}\|_{\infty}$ in the above derivations.) Applying the union bound over $a$ then gives us (4.33). ∎

By using (F.1) in the first step and (D.43) in the second step, we have $\|{G^{[n]}\boldsymbol{v}_{\beta}(\xi)}\|\leq\frac{1}{\eta}\mathinner{\!\left\lVert\boldsymbol{v}_{\beta}(\xi)\right\rVert}\prec\frac{\sqrt{n}}{\eta}$. The same arguments will also give us $\|{\boldsymbol{v}_{\alpha}^{\mkern-1.5mu\mathsf{T}}(\xi)\widetilde{G}(z)}\|\prec\sqrt{n}/\eta$. It follows that

Thus, one way to bound the left-hand side of (4.36) is to obtain estimates on the operator norm of $\Delta_{c}$. We start from $\Delta_{1}$ in (4.18a). Since it is a diagonal matrix,

where the second step is due to (D.43). For $\Delta_{2}$ in (4.18b), write $u_{k}\coloneqq\max_{a}\mathinner{\!\left\lvert r^{k}(\xi_{a})-1\right\rvert}$. Using the expansion in (4.28) and the triangular inequality, we have

The cases of $\Delta_{3}$ and $\Delta_{5}$ require a different approach. We start with $\Delta_{3}$. Recall the diagonal matrix notation introduced in (4.32). Write

From the construction of $\Delta_{3}$ in (4.18c), one can check that

We complete the proof by showing (4.36) for $\Delta_{5}$. By construction, $\Delta_{5}$ is the linear combination of a finite number of matrices \big{\{}\Delta_{5}^{(k,t)}\big{\}}, indexed by $k,t$. The entries of $\Delta_{5}^{(k,t)}$ are

Now consider the left-hand side of (4.36), with $\Delta_{c}$ replaced by $\Delta_{5}^{(k,t)}$. By using (4.50) and the notation introduced in (4.32) and (4.41), we have

We may now complete the proof of Proposition 2.

Observe that the different coordinates of $A$ are statistically exchangeable. Thus, it is sufficient to show (4.23) for any particular choice of the index $i$. We choose to do so for $i=n$. With the decomposition given in (4.19), our task boils down to verifying that the right-hand side of (4.19) concentrates around the leading term given in (4.23). To that end, we first show that the term $\theta(z)$ in (4.21) is close to

where $E=\mathfrak{Re}(z)$. This then gives us

We already have controls on the operator norm of $A$. Indeed, the estimates in (E.39) give us

Similarly, $\|\widetilde{A}\|_{\mathsf{op}}\prec 1$, as $\widetilde{A}$ is just the $(d-1)$-dimensional version of $A$. Substituting (4.57) and (4.58) into (4.55), and using the large deviation estimates in (4.24), we get

4 Proof of Theorem 1

where $\gamma_{a},\gamma_{b}$ are the constants defined in Theorem 1.

By (4.4), (4.5), (4.3), and the high-probability estimate given in Proposition 2, we have

where the error term $\Delta\coloneqq-\frac{1}{n}\sum_{i}Q_{i}(z)G_{ii}(z)$ satisfies the estimate

For any $0<\varepsilon<1/2$ and $D>0$, we have

In step (a), the second term on the right-hand side of the inequality follows from Proposition 9, where the error term $\omega$ in (G.2) is $\Delta/s(z)$; Step (b) holds for all sufficiently large $d$; Step (c) follows from the high-probability estimates given in (4.64) and (4.66). As $\varepsilon$ and $D$ can be chosen arbitrarily, we have verified (2.11). The almost sure convergence of $s(z)$ to $m(z)$ then follows from (2.11) and the Borel-Cantelli lemma.

The Isotropic Gaussian Model

Our approach is based on a simple coupling method, which allows us to define $J$ in the same probability space as the matrix $A$ in (1.1) and to write $J$ as a perturbation of $A$. To that end, we rewrite $\boldsymbol{g}_{i}$, for each $i\in[n]$, as

where $\nu_{i}\coloneqq\mathinner{\!\left\lVert\boldsymbol{g}_{i}\right\rVert}$ and $\boldsymbol{x}_{i}\coloneqq\boldsymbol{g}_{i}/\mathinner{\!\left\lVert\boldsymbol{g}_{i}\right\rVert}$. By the rotational invariance of the isotropic Gaussian distribution, $\boldsymbol{x}_{i}\sim\mathsf{Unif}(\mathcal{S}^{d-1})$, and $\boldsymbol{x}_{i}$ is independent of $\nu_{i}$. In our subsequent discussions, we will often use the following standard concentration result (see, e.g., [39, Theorem 3.1.1]): for any $t>0$,

where $c>0$ is some absolute constant. In other words, $\sqrt{d}(\nu_{i}-1)$ is a sub-Gaussian random variable whose sub-Gaussian norm is upper-bounded by an absolute constant. This then immediately implies that

Let $A$ and $J$ be the kernel random matrices defined in (1.1) and (5.1), respectively, where the nonlinear function

where $C_{L}$ is some constant that only depends on $L$.

For any $0\leq b\leq L$, we can expand the monomial $x^{b}$ as a linear combination of Gegenbauer polynomials, i.e.,

where $\left\{\beta_{b,k}\right\}$ are some expansion coefficients. Note that $\left\{\beta_{b,k}\right\}$ depend on the dimension $d$ [see, e.g., (2.3)], but we suppress this explicit dependence to streamline the notation. Recall the definition of the matrices $\left\{A_{k}\right\}$ in (2.6). Using (5.10) and after rearranging the terms, we can write (5.9) as

Here, the difference between $A$ and $J$ is represented by the sum of three $n\times n$ matrices, defined as

where $\left\{N_{k}\right\}$ are the constants given in (A.3). Next, we show that the perturbations by $\Delta_{1},\Delta_{2},\Delta_{3}$ only cause negligible changes in the Stieltjes transform of $A$. First, applying the triangular inequality gives us

By the estimates of Proposition 8, we have

To bound the expansion coefficients $\left\{\beta_{b,k}\right\}$, we let $\xi$ be a random variable sampled from the distribution $\mathsf{s}_{d}^{1}$ in (1.3). Using (5.10) and the orthogonality of the Gegenbauer polynomials $\left\{q_{k}(x)\right\}$ with respect to the law of $\xi$, we have

Substituting (5.16), (5.17), and (5.18) into (5.15) then yields

The operator norm of $\Delta_{2}$ can be bounded similarly:

Now we consider $\Delta_{3}$. Its operator norm is not small, but it has a low-rank structure. Recall from the decomposition in (2.15) that $\operatorname{rank}(A_{k}+\sqrt{\frac{N_{k}}{n}}I)\leq N_{k}$. It follows that

Given the estimates in (5.19), (5.21), and (5.22), we can reach the statement of the proposition by applying the perturbation inequalities (F.7) and (F.8) in Lemma 18 and the triangular inequality. ∎

As an immediate consequence of Proposition 3, we observe that Theorem 1 remains valid even when data vectors are sampled from the isotropic Gaussian distribution, as opposed to the spherical distribution. More specifically, the estimate provided in (2.11) continues to hold when we substitute $s_{A}(z)$ with $s_{J}(z)$, where $J$ represents the kernel matrix defined in (5.1) with $f_{d}$ corresponding to a linear combination of the first $L$ Gegenbauer polynomials, as expressed in (2.1).

Next, we establish a counterpart of Theorem 2 for the isotropic Gaussian case, under the following assumption on the function $f_{d}$ in (5.1).

Let $g\sim\mathcal{N}(0,1)$, and let $\left\{h_{k}(x)\right\}_{k}$ represent the Hermite polynomials defined in (2.22). The function $f_{d}$ in (5.1) satisfies the following conditions:

for some finite numbers $\left\{\mu_{k}\right\}$.

The sequence $\left\{\mu_{k}\right\}$ is square-summable, i.e.,

The conditions in Assumption 2 closely resemble those in Assumption 1. The only difference lies in the expression given in (5.26), where we replace the probability density $w_{d}(x)$ from (2.28) with a new density function, $\breve{w}_{d}(x)$. Furthermore, when $f_{d}=f$ is a function that remains independent of the dimension $d$, a simpler set of sufficient conditions can be employed, as stated in the following lemma.

The statement of Theorem 2 remains valid if we replace the matrix $A$ with the matrix $J$ as given in (5.1), and replace Assumption 1 with Assumption 2.

Our proof is a modification of the proof for Theorem 2, and we provide the details in Appendix H. ∎

In concluding this section, we highlight an important caveat concerning the results presented above. While changing the data distribution from spherical to isotropic Gaussian does not alter the weak limit of the empirical spectral distribution of the kernel matrix, the Gaussian model may indeed introduce additional spike eigenvalues that are absent in the original spherical case. To illustrate this, let us consider $J$ and $A$ as the kernel matrices defined in (5.1) and (1.1), respectively, with the nonlinear function chosen as $f_{d}(x)=x^{2}-1$. Furthermore, we assume that $n/d^{2}=+/o(1)$ for a constant $>/0$.

Notice that the nonlinear function in this case is equivalent (up to a constant) to the second-degree Gegenbauer polynomial $q_{2}(x)$, as seen in (2.3). Consequently, we can apply Proposition 8, resulting in:

This implies that, for any $\epsilon>0$, the spectrum of $A$ will be confined within the interval $[-d^{\epsilon},d^{\epsilon}]$ with high probability as $d\to\infty$. However, this is not the case for $J$. Although its bulk eigenvalues share the same weak limit as those of $A$, the matrix $J$ exhibits two spike eigenvalues of order $\mathcal{O}(\sqrt{d})$.

To see that, we first use the decomposition (5.2) to express

By defining $\delta_{i}\coloneqq\nu_{i}^{2}-1$ for $i\in[n]$, we can further expand (5.29):

where $E$ and $F$ are two $n\times n$ matrices representing the difference between $J$ and $A$. Observe that the contribution from $E$ is negligible. Indeed, from the estimate in (5.4), we obtain $\delta_{i}=\mathcal{O}_{\prec}(d^{-1/2})$, and by applying the union bound, $\max_{i\in[n]}\mathinner{\!\left\lvert\delta_{i}\right\rvert}=\mathcal{O}_{\prec}(d^{-1/2})$. Employing these estimates and the one in (5.27), we have

Next, we examine the second perturbation matrix $F$ in (5.30). $F$ has rank $2$, and its two nonnegative eigenvalues, denoted by $\lambda^{+}$ and $\lambda^{-}$, can be directly computed as

After substituting (5.34) into (5.32), the formula can be simplified as

Since $n/d^{2}\to$/, the two nonzero eigenvalues of $F$ are thus close to $\pm\sqrt{2\kappa d}$. Finally, given that $G=A+E+F$ with $\mathinner{\!\left\lVert A+E\right\rVert}_{\mathsf{op}}\prec 1$, standard eigenvalue perturbation arguments allow us to conclude that the matrix $J$ also has two large spike eigenvalues close to $\pm\sqrt{2\kappa d}$.

Summary and Discussions

We have established our results only for data vectors sampled from spherical or isotropic Gaussian distributions. In fact, Theorem 1 is expected to hold (with certain modifications in error bounds) for more general data distributions exhibiting reasonably fast decay at infinity. It is also worth noting that the error bounds in Theorem 1 are not optimal. Moreover, we currently constrain the imaginary part of the spectral parameter $z=E+i\eta$ to $\eta>\tau$ for some constant $\tau$. Although our proof approach can accommodate $\eta\geq d^{-c}$ for some small constant $c$ with more careful book-keeping, this is still far from the regime $\eta\sim n^{-1+\epsilon}$ required to reach individual eigenvalue locations.

Extending Theorem 1 to encompass the entire range $n^{-1+\epsilon}\leq\eta\leq 1$ may prove challenging due to the presence of multiscale structures in the eigenvalues. These structures emerge because the matrix $A$ in (2.7) is a sum of component matrices $\left\{A_{k}\right\}$ across different scales. As a result, accurately characterizing individual eigenvalues, including extremal ones, poses an interesting open problem. Related issues, such as eigenvector statistics and eigenvalue universality, also hinge on addressing this matter. The potential for multiscale structures sets the nonlinear model apart from standard mean field models like Wigner matrices or sample covariance matrices . We aim to tackle some of these questions in future papers.

Appendix A Gegenbauer Polynomials and Spherical Harmonics

In this appendix, we collect a few useful properties of the (normalized) Gegenbauer polynomials $\{q_{k}(x)\}_{k\geq 0}$ defined in (2.2). All of these results are standard, and their proofs can be found in e.g., .

We will repeatedly apply this recurrence relation in our proof of Proposition 1.

Using the Gram-Schmidt procedure, we can construct an orthonormal set of spherical harmonics. Let $Y_{k,i}(\boldsymbol{x})$, for $i\in[N_{k}]$, denote the $i$th degree-$k$ spherical harmonic in this orthonormal set. We then have

The Gengenbauer polynomials and spherical harmonics are deeply connected. In particular, we have the following identity, which can be viewed as a high-dimensional generalization of the classical addition theorem : for any $\boldsymbol{x},\boldsymbol{y}\in\mathcal{S}^{d-1}$,

In addition, for the special case of $\boldsymbol{x}=\boldsymbol{y}$, we have

The identity in (A.5) is useful because it allows us to “linearize” the term $q_{k}(\sqrt{d}\,\boldsymbol{x}^{\mkern-1.5mu\mathsf{T}}\boldsymbol{y})$ as an inner product of two vectors made of the spherical harmonics. Moreover, (A.4) implies that, if $\boldsymbol{x}\sim\mathsf{Unif}(\mathcal{S}^{d-1})$, the vector of spherical harmonics $(Y_{k,1}(\boldsymbol{x}),\ldots,Y_{k,N_{k}}(\boldsymbol{x}))$ is an isotropic random vector.

Using (A.5) and (A.6), we can rewrite the matrix $A_{k}$ defined in (2.6) as

where $S_{k}$ is a $N_{k}\times n$ matrix whose entries are the spherical harmonics, i.e.

and $\left\{\boldsymbol{x}_{i}\right\}_{i\in[n]}$ are the data vectors in the definition in (2.6).

Appendix B Proof of Proposition 1

We first state several simple properties of the set $\mathcal{E}_{m}$ introduced in (4.12).

where $r(\xi)$ is the function defined in (4.10).

The property (B.1a) follows immediately from the definition of $\mathcal{E}_{m}$. To show (B.1b), (B.1c), and (B.1d), we can assume without loss of generality that

for some $a,b$ such that $0\leq a,b\leq m$. The more general case, where $\mathcal{E}_{m}(\xi_{1},\xi_{2})$ is a linear combination of terms like (B.2), can be handled by using the linearity of $\mathcal{E}_{m+1}$ and $\mathcal{E}_{m+2}$.

The recurrent relation (A.1) implies that $\xi_{1}q_{a}(\xi_{1})=\sum_{i\in[0,a+1]}\mu_{i}q_{i}(\xi_{1})$ and $\xi_{2}q_{b}(\xi_{2})=\sum_{j\in[0,b+1]}\nu_{j}q_{i}(\xi_{1})$, where $\mu_{i}=\mathcal{O}(1),\nu_{i}=\mathcal{O}(1)$ are some fixed expansion coefficients. It follows that

where $\pi_{ij}=\max\left\{a+1,b+1\right\}/2-\max\left\{i,j\right\}/2\geq 0$. Thus, we $\frac{\xi_{1}\xi_{2}}{\sqrt{d}}\mathcal{E}_{m}(\xi_{1},\xi_{2})\in\mathcal{E}_{m+1}$.

The properties stated in (B.1c) can be proved similarly. Consider the function in (B.2). We have $\xi_{1}^{2}q_{a}(\xi_{1})=\sum_{i=0}^{a+2}\mu_{i}q_{i}(\xi_{1})$ for some expansion coefficients $\mu_{i}=\mathcal{O}(1)$. It follows that

where $\pi_{i}=\max\left\{a+2,b+2\right\}/2-\max\left\{i,b\right\}/2\geq 0$. Since $\mu_{i}/d^{\pi_{i}}=\mathcal{O}(1)$, we have $(\xi_{1}^{2}/d)\mathcal{E}_{m}(\xi_{1},\xi_{2})\in\mathcal{E}_{m+2}$. That $(\xi_{2}^{2}/d)\mathcal{E}_{m}(\xi_{1},\xi_{2})\in\mathcal{E}_{m+2}$ follows from analogous arguments.

Recall the definition of $r(\xi)$ in (4.10). We have

Note that $\mathcal{E}_{m}(\xi_{1},\xi_{2})\in\mathcal{E}_{m+2}$ by (B.1a). Applying the property in (B.1b) twice, we have $(\xi_{1}\xi_{2}/\sqrt{d})^{2}\mathcal{E}_{m}(\xi_{1},\xi_{2})\in\mathcal{E}_{m+2}$. By (B.1c), the last two terms of (B.3) also belong to $\mathcal{E}_{m+2}$. By the linearity of $\mathcal{E}_{m+2}$, we can conclude that the left-hand side of (B.1d) indeed belongs to $\mathcal{E}_{m+2}$. ∎

To show (B.4), we use the recurrent relation in (A.1), which gives us $\xi q_{k}(\xi)=a_{k}q_{k+1}(\xi)+a_{k-1}q_{k-1}(\xi)$. It follows that

Note that the last term on the right-hand side of above expression is in $\mathcal{E}_{k-1}$ and thus in $\mathcal{E}_{k+1}$ [by (B.1a)]. From (A.2), we have $a_{k}^{2}=(k+1)(1+\mathcal{O}(1/d))$ and $a_{k}a_{k-1}=\sqrt{k(k+1)}(1+\mathcal{O}(1/d))$. Thus, $a_{k}^{2}\frac{q_{k+1,k+1}}{d^{(k+1)/2}}=(k+1)\frac{q_{k+1,k+1}}{d^{(k+1)/2}}+\mathcal{E}_{k+1}(\xi_{1},\xi_{2})$, and similarly, $a_{k}a_{k-1}\frac{q_{k-1,k+1}+q_{k+1,k-1}}{d^{(k+1)/2}}=\sqrt{k(k+1)}\frac{q_{k-1,k+1}+q_{k+1,k-1}}{d^{(k+1)/2}}+\mathcal{E}_{k+1}(\xi_{1},\xi_{2})$.

Next, we show (B.5). Recall from the definition in (4.10) that $\sqrt{1-d^{-1}}r^{2}(\xi)=1-\xi^{2}/d$. Thus,

which then leads to (B.5), as $d/(d-1)=1+\mathcal{O}(1/d)$ and $[d/(d-1)]a_{k-1}a_{k}=\sqrt{k(k+1)}+\mathcal{O}(1/d)$. ∎

Next, we prove Proposition 1 by induction on the polynomial degree $k$. Throughout the proof, we use the shorthand notation introduced in (B.6). Recall that $q_{0}(x)=\widetilde{q}_{0}(x)=1$ and $q_{1}(x)=\widetilde{q}_{1}(x)=x$. It is then straightforward to verify the formula (4.13) for $k=0,1$. Specifically,

Now we carry out the induction. Assume that (4.13) holds for $k$ and $k-1$, with some $k\geq 1$. To prove it for $k+1$, we apply the recurrence relation in (A.1), which gives us

where in reaching the last step we have used (4.13) to expand $q_{k}(r_{1,2}x+\frac{\xi_{1}\xi_{2}}{\sqrt{d}})$ and $q_{k-1}(r_{1,2}x+\frac{\xi_{1}\xi_{2}}{\sqrt{d}})$.

On the right-hand side of (B.13) there are factors related to $x$ in the form of $x\widetilde{q}_{k-m}(x)$. They are polynomials of $x$, and can thus be rewritten as a linear combination of the orthogonal polynomials. To that end, we first recall from (2.5) that $\left\{\widetilde{q}_{k}(x)\right\}$ denote the orthogonal polynomials defined for dimension $d-1$. Similar to (A.1), they also satisfy a recurrence relation

Replacing all the factors of $x\widetilde{q}_{k-m}(x)$ in (B.13) by the right-hand side of (B.16), we can rewrite (B.13) as a linear combination of the orthogonal polynomials $\left\{\widetilde{q}_{k+1-m}(x)\right\}_{m=0}^{k+1}$, i.e.,

where $\left\{C_{m}(\xi_{1},\xi_{2})\right\}_{m=0}^{k+1}$ are some functions that only depend on $\xi_{1},\xi_{2}$ but not on $x$. Next, we identify the exact expressions for $C_{m}(\xi_{1},\xi_{2})$.

where the last equality follows from (A.2) and (B.15).

For each $m$ in the range $2\leq m\leq k+1$, we have

Using the formulas for $a_{k}$ and $\widetilde{a}_{k}$ in (A.2) and (B.15), we have

Substituting (B.23), (B.24), (B.25), (B.26) into (B.22), we have

Note that the expressions in (B.18), (B.21) and (B.27) exactly match the right-hand side of (4.13) for $k+1$. Thus, by induction, we can conclude that the formula (4.13) holds for all $k$.

Appendix C Review of the Marchenko-Pastur Law

In this appendix, we review some basic properties of the Marchenko-Pastur (MP) law that will be used in our work. Let $X$ be an $M\times n$ matrix whose entries $X_{\mu i}$ are independent complex-valued random variables satisfying

Additionally, $\mathinner{\!\left\lvert\sqrt{n}X_{\mu i}\right\rvert}$ has a sufficient number of bounded moments. We also assume that $M$ and $n$ satisfy the bounds

for some positive constant $C$, and define the aspect ratio

which may depend on $n$. Now consider an $n\times n$ matrix

where $\Sigma$ is a diagonal $M\times M$ matrix with matrix elements $\sigma_{i}$ on the diagonal. Note that the subtraction by $\frac{\operatorname{tr}\Sigma}{n}I$ in (C.3) makes sure that the diagonal elements of $R$ are approximately equal to 0. We do this to match the construction in (1.1).

Assuming that the law of the diagonal entries of $\Sigma$ is given by a probability law $\pi_{M}=M^{-1}\sum\delta_{\sigma_{i}}\to\pi$. The MP law asserts that the limiting Stieltjes transform of the eigenvalues of $R$, denoted by $m(z)$, satisfies the self-consistent equation

See, e.g., eq. (2.10) in Knowles and Yin . (Notice that our equation is slightly different due to a constant shift of the eigenvalues). In the special case of $\Sigma=tI$ and $R=tX^{\mkern-1.5mu\mathsf{T}}X-t\phi I$, the above formula reduces to

The limiting eigenvalue density associated with (C.5) is given by

Appendix D Moment Bounds and Concentration Inequalities

In this appendix, we derive several moment and concentration inequalities that will be used in our proof.

where $C(k,p)$ is some constant that only depends on $k$ and $p$.

Both the probability measure $\mathsf{s}_{d}^{1}$ and the Gegenbauer polynomials $q_{k}(x)$ depend on the dimension $d$. However, the upper bounds in (D.1) and (D.2) hold uniformly for all $d$.

The probability distribution of $\xi$ is given by

where $\Gamma(\cdot)$ is the gamma function. Thus, we have

where the second line uses the relationship between the beta function and the gamma function , and the last line is obtained by applying Stirling’s formula for the gamma function .

Next, we show (D.2). Two special cases are easy: for $k=0$, we have $\mathinner{\!\left\lVert q_{0}(\xi)\right\rVert}_{L^{p}}=1$; for $p=1$, $\mathinner{\!\left\lVert q_{k}(\xi)\right\rVert}_{L^{1}}\leq\mathinner{\!\left\lVert q_{k}(\xi)\right\rVert}_{L^{2}}=1$. Thus, in what follows we assume $k\geq 1$ and $p\geq 2$. Denote by $c_{k}$ the leading coefficient of the polynomial $q_{k}(\xi)$. Since $q_{k}(\xi)-c_{k}\xi^{k}$ is a polynomial of degree $k-1$, it can be written as linear combination of the lower order Gegenbauer polynomials, i.e.,

where $\left\{b_{k,i}\right\}_{i}$ are the expansion coefficients. By using the orthogonality of the Gegenbauer polynomials (see (2.2)), we have, for $0\leq i\leq k-1$,

Applying the Cauchy-Schwarz inequality, we have

for all $d\geq 2$. Using this estimate in (D.9) leads to

Starting from $\mathinner{\!\left\lVert q_{0}(\xi)\right\rVert}_{L^{p}}=1$, we can apply the above bound recursively to verify (D.2). ∎

In the discussion of the equivalence principle for general nonlinear kernel functions in Section 2.3 and Section 5, we examine three closely-related probability distributions: the standard normal distribution $\mathcal{N}(0,1)$, the probability measure $\mathsf{s}_{d}^{1}$ as defined in (1.3), and the distribution of the random variable

Here, $w(x)$ denotes the probability density function of $\mathcal{N}(0,1)$, $w_{d}(x)$ denotes the density function of the probability measure $\mathsf{s}_{d}^{1}$ as defined in (1.3), and $\breve{w}_{d}(x)$ represents the probability density function of the random variable $\breve{\xi}_{d}$ defined in (D.11).

We begin with (D.12). Let $M$ be a constant in the interval $[c_{3},\sqrt{d/2}]$. We split the integration in (D.12) into two parts:

In reaching the second step, we have used the condition that $\mathinner{\!\left\lvert f(x)\right\rvert}<c_{1}e^{c_{2}\mathinner{\!\left\lvert x\right\rvert}}$ for $x\geq M$, and we have also used the property that the densities $w(x)$ and $w_{d}(x)$ are both even functions.

A closed-form expression of $w_{d}(x)$ can be found in (D.3). By applying Stirling’s formula for the gamma function and Taylor’s expansion for $\log(1-x^{2}/d)$, it is straightforward to verify that $w_{d}(x)/w(x)=[1+\mathcal{O}(1/d)]e^{3x^{2}/(2d)}$ over the interval $\mathinner{\!\left\lvert x\right\rvert}\leq M$, where $M$ is some constant satisfying $M<\sqrt{d/2}$. Consequently, we have

where the final inequality follows from the Cauchy-Schwartz inequality and standard Gaussian tail bounds. By assumption, $\int f^{2}(x)w(x)\operatorname{d\!}x<\infty$. Upon substituting (D.16) and (D.17) into (D.15), we have

Since the constant $M$ can be chosen arbitrarily, we have then demonstrated (D.12).

In reaching the final step, we have used the inequality $\log(1+t^{2}/d)\leq t^{2}/d$, which then implies that $(1+t^{2}/d)^{-d/2}\geq e^{-t^{2}/2}$. It is straightforward to verify via Taylor’s expansion that

To bound the second integral on the right-hand side of (D.21), we use the inequality $\log(1+t^{2}/d)\geq(\log 2)t^{2}/d$ for $0\leq t\leq\sqrt{d}$. It follows that

where the final step uses standard Gaussian tail bounds. For the third integral on the right-hand side of (D.21), we obtain

By combining the bounds (D.22), (D.24), and (D.25) for the three integrals, we can establish that

Let $M_{d}$ be a positive number that depends on $d$. We divide the integral in (D.13) into two parts, following a similar approach as in (D.15). This results in:

where in obtaining the last inequality, we have utilized (D.26), the assumption that $\int f^{2}(x)w(x)\operatorname{d\!}x<\infty$, and the estimate provided in (D.17) to bound the first two terms on the right-hand side of (D.27). By applying the Cauchy-Schwarz inequality, we have:

In the following discussion, we derive several useful moment and high probability bounds for linear and quadratic functions of independent random variables. We first recall the following estimates, whose proof can be found in [16, Lemma 7.8, Lemma 7.9].

where $C(a,b,p)$ is a quantity that depends on $a$, $b$, and $p$, but not on $n$.

To bound the first term on the right-hand side, we use the standard decoupling technique. Observe that, for every $i\neq j$, the following identity holds:

where the sum ranges over all subsets of $[n]$, and $Z_{n}=2^{n-2}$. Applying this identity and the triangular inequality allows us to write

Note that the families of random variables $\left\{q_{a}(\xi_{i})\right\}_{i\in I}$ and $\left\{q_{b}(\xi_{j})\right\}_{j\in I^{c}}$ have mean-zero independent components with unit variances; they are also mutually independent. We can then apply (D.32b) to bound each term in the sum in (D.36). Moreover, the sum involves a total of $2^{n}-2$ terms. Thus,

Now we bound the second term on the right-hand side of (D.34). Write $X_{i}=q_{a}(\xi_{i})q_{b}(\xi_{i})-\delta_{ab}$ and $\sigma=\mathinner{\!\left\lVert X_{i}\right\rVert}_{L^{2}}$. Then $(X_{i}/\sigma)_{i}$ is a family of independent random variables that satisfy the condition of Lemma 13. Using (D.32a) gives us

By the triangular inequality (in the first line) and the Cauchy-Schwarz inequality (in the third line),

Substituting (D.37) and (D.42) into (D.34), and by applying the moment bound in (D.2), we complete the proof of the statement in (D.33). ∎

The moment estimates obtained above can be easily turned into high probability bounds, as follows.

where $r(\xi)$ is the function defined in (4.10).

Now we prove (D.45). The case of $k=0$ is trivial, so we assume $k\geq 1$. From the definition in (4.10), $r(\xi)=(1-\xi^{2}/d)^{1/2}+\mathcal{O}(1/d)$, and thus

Combining this deterministic bound with the high-probability bound in (D.43), we have that, for $\xi_{i}\sim\mathsf{s}_{d}^{1}$,

Applying (D.51) and the union bound then allows us to conclude (D.45). ∎

where the two functions $f_{1},f_{2}$ are such that