Spectrum of inner-product kernel matrices in the polynomial regime and multiple descent phenomenon in kernel ridge regression

Introduction

Kernel methods are among the most popular tools in statistics and machine learning and have been extensively studied in the classical bias-variance trade-off setting [BTA11, Wai19]. Over the past few years, they have attracted a renewed interest because of their connection to neural networks in the ‘neural tangent kernel’ regime [JGH18, LL18, DZPS18, LXS+19, AZLS19, COB19]. Moreover, it was argued in [BMM18] that kernel methods share a number of surprising phenomena with deep learning, which are not explained by classical theory. This prompted a number of works to study kernel methods in the ‘overfitted regime’, which brought to light several interesting behavior: near optimality of interpolators and benign overfitting [LR20, GMMM21, BLLT20], self-induced regularization [GMMM21, LRZ20] and double descent of the prediction risk [MM22, HMRT22]. These phenomena appear in the high-dimensional regime, when both the number of samples $n$ and the dimensionality $d$ of the data are large [RZ19], and are not captured by previous approaches such as capacity/source conditions [CDV07]. This motivates the development of theory specific to kernel methods in high-dimension.

[EK10] showed that when $n,d\to\infty$ with $n/d=\Theta(1)$, the random matrix ${\bm{H}}$ can be approximated consistently in operator norm by its linearization $\overline{{\bm{H}}}_{\text{lin}}$ (i.e., $\|{\bm{H}}-\overline{{\bm{H}}}_{\text{lin}}\|_{\rm op}\rightarrow 0$ in probability):

1.2 Precise asymptotics of KRR prediction error in the polynomial regime

where $\|\cdot\|_{\mathcal{H}}$ is the RKHS norm associated to kernel $H_{d}$ in $L^{2}({\mathcal{A}}_{d})$. The test error (or prediction error) of KRR is given by

1.3 Equivalence with a Gaussian covariates model

A recent string of work started showing equivalence between non-linear regression models and simpler Gaussian covariates models in high-dimension [MM22, GLR+20, HL20]. These results hint at some general universality phenomena in high-dimensional models, where the test error only depends on the covariance of the features [MS22] and a few properties of the non-linearity.

Here, we will simply make the following observation: the kernel ridge regression model has the same asymptotic prediction error in the polynomial regime as a simpler linear regression model with Gaussian covariates. Consider a target function:

where $\{Y_{ks}\}$ is the polynomial basis that diagonalizes inner-product kernels on $L^{2}({\mathcal{A}}_{d})$, i.e., $Y_{ks}$ is an eigenfunction of the kernel operator with eigenvalue $\mu_{d,k}/B({\mathcal{A}}_{d},k)$. Let us now state the equivalent linear regression model: we are given $n$ i.i.d. pairs $({\bm{z}}_{i},y_{i})_{i\in[n]}$ with

We fit this model using ridge regression with ridge parameter $\lambda>0$:

Such models were studied in the overfitted regime in [BLLT20, TB20, RMR21]. Here, we show that kernel ridge regression has the same asymptotic test error as the Gaussian covariates model (7):

Under the same assumptions as Theorem 3, for any $\kappa,\psi>0$ and $n/d^{\kappa}\to\psi$ as $d,n\to\infty$, we have

2 Notations

Related work

The prediction error of kernel ridge regression in the linear high-dimensional regime was studied in [LR20, LLS21, BMR21] using the linearization of the kernel in this regime [EK10]. In particular, [LR20] points out that the minimum RKHS norm interpolating solution (KRR with ridge penalty $\lambda\to 0^{+}$) can still generalize well. Another line of work considers linear regression models with Gaussian or sub-Gaussian covariates as in Eq. (7), which are technically easier to study and allows to focus on the interaction between eigenvalue decay and target function in the prediction error [BLLT20, TB20, RMR21, CLKZ21].

Inner-product kernel matrices in the polynomial regime

We start by recalling some basic properties of functional spaces over ${\mathcal{A}}_{d}$ (see Appendix C for a complete exposition). Let $L^{2}({\mathcal{A}}_{d}):=L^{2}({\mathcal{A}}_{d},{\rm Unif})$ be the space of square-integrable functions on ${\mathcal{A}}_{d}$ with scalar product and norm denoted by $\langle\cdot,\cdot\rangle_{L^{2}}$ and $\|\cdot\|_{L^{2}}$ given by

where we denoted $\mu_{d,k}(h)=\xi_{d,k}(h)B({\mathcal{A}}_{d},k)$. For ease of notation, we will write $\xi_{d,k}:=\xi_{d,k}(h)$ and $\mu_{d,k}:=\mu_{d,k}(h)$. Note that $\xi_{d,k}\geq 0$ for any $k$, by assumption of $H_{d}$ being positive semi-definite.

We define the empirical Gegenbauer matrices ${\bm{Q}}_{k}=(Q_{k}^{(d)}(\langle{\bm{x}}_{i},{\bm{x}}_{j}\rangle))_{ij\in[n]}$. Note that with the above notations ${\bm{Q}}_{k}=B({\mathcal{A}}_{d},k)^{-1}{\bm{Y}}_{k}{\bm{Y}}_{k}^{\mathsf{T}}$. We can therefore decompose the empirical kernel matrix in terms of the Gegenbauer matrices:

2 Limiting spectral distribution of the empirical kernel matrix

We will make the following genericity assumption on the kernel functions $\{h_{d}\}_{d\geq 1}$:

as the sum of a low-rank spiked matrix with diverging eigenvalues plus a multiple of the identity matrix (the high-degree part of the kernel plays the role of a ‘self-induced ridge regularization’). In particular, [GMMM21, MMM21a] used this decomposition to show that kernel ridge regression with any target function $f_{*}\in L^{2}({\mathcal{A}}_{d})$ learns exactly the projection ${\mathsf{P}}_{\leq\lfloor\kappa\rfloor}f_{*}$ on degree-$\lfloor\kappa\rfloor$ polynomials and none of the high-degree part ${\mathsf{P}}_{>\lfloor\kappa\rfloor}f_{*}$.

Recall the definition of the Marchenko-Pastur distribution: given a parameter $\psi>0$,

where $\lambda_{\pm}=(1\pm\sqrt{\psi})^{2}$.

3 Proof of Theorem 2

Several sufficient conditions for the Marchenko-Pastur law have been proved in the literature for matrices without independent entries [BZ08, Ada11, PS11, O’R12, Yas16]. Here, we use a simple condition presented for example in [Yas16, Remark 2.2] (see also [BZ08, PS11]), which applies to random matrices with iid isotropic rows, and which will directly imply Theorem 2.

Note that ${\bm{B}}$ is a symmetric matrix and we have

Application: kernel ridge regression in the polynomial regime

where $\|\cdot\|_{\mathcal{H}}$ is the RKHS norm associated to kernel $H_{d}$ and we denoted ${\bm{y}}=(y_{1},\ldots,y_{n})$. The test error (or prediction error) of KRR is given by

where ${\bm{h}}({\bm{x}})=(h_{d}(\langle{\bm{x}},{\bm{x}}_{i}\rangle/d))_{i\in[n]}$. We will also consider the training error and the RKHS norm of the KRR solution, which are given by

In this section, we characterize the asymptotic test error with high probability over a class of random target functions. We will discuss in Section 4.2 how to extend these results to fixed target functions.

We assume the following distribution over the sequence of target functions $\{f_{*,d}\in L^{2}({\mathcal{A}}_{d})\}_{d\geq 1}$:

Let $\{f_{*,d}\in L^{2}({\mathcal{A}}_{d})\}_{d\geq 1}$ be a sequence of target functions with decomposition in the polynomial basis

where we recall that $\mu_{d,k}=\xi_{d,k}B({\mathcal{A}}_{d},k)$ with $\xi_{d,k}$ the $k$-th Gegenbauer coefficient of $h_{d}$ (see Eq. (10)).

We are now in position to state our main theorem:

where $r_{\psi}(-\zeta)=\int(x+\zeta)^{-1}\nu_{{\rm MP},\psi}({\rm d}x)$ is the Stieljes transform of the Marchenko-Pastur distribution (see Eq. (14)).

as $d\to\infty$, where the convergence in probability is over the randomness in ${\bm{X}},{\bm{\varepsilon}},f_{*}$. Furthermore,

A detailed proof of Theorem 3 can be found in Appendix B. Let us give some intuition on the asymptotic formula for the test error (19). It will be instructive to consider the contribution to the test error of the three subspaces

2 Pointwise asymptotic test error

In order to prove the asymptotic test error formula, we would need to show that

As shown in [HMRT22], this can be reduced to showing that

Acknowledgements

This work was supported by NSF through award DMS-2031883 and the Simons Foundation through Award 814639 for the Collaboration on the Theoretical Foundations of Deep Learning. We also acknowledge the NSF grant CCF-2006489 and the ONR grant N00014-18-1-2729.

References

Appendix A Proof of Proposition 1: the case of spherical harmonics

The proof of Proposition 1 will rely on an explicit representation of spherical harmonics in terms of the generalized spherical coordinate system in dimension $d$. See for example [Ave12, DX13].

where $0\leq\theta_{1}\leq 2\pi$ and $0\leq\theta_{i}\leq\pi$ for $i=2,\ldots,d-1$. The uniform probability measure on the unit sphere is given by

where $|\alpha^{j+1}|=\alpha_{j+1}+\ldots+\alpha_{d-1}$, $d_{j}=2|\alpha^{j+1}|+d-j+1$,

A proof of this proposition can be found for example in [DX13]. For completeness, we include here the proof with our notations and normalization choice.

In order to check that Eq. (27) is a homogeneous polynomial, recall that in the spherical coordinates (25), we have

is a polynomial of degree $\alpha_{j}$. We can further write $h_{\bm{\alpha}}(x_{1},x_{2}):=(x_{1}^{2}+x_{2}^{2})^{\alpha_{d-1}/2}g_{\bm{\alpha}}(\theta_{1})$ as the real part or the imaginary part of the polynomial $(x_{2}+ix_{1})^{\alpha_{d-1}}$ (or a constant if $\alpha_{d-1}=\alpha_{d}=0$), depending on $\alpha_{d-1},\alpha_{d}$. We deduce that

A.2 Proof of Proposition 1

We will show that both these terms are $o_{d}(1)$, which implies the concentration in probability of the quadratic form.

We proceed similarly than in the main text. Consider ${\bm{C}}$ the square matrix of size $B(B-1)$ such that for any $\bm{\alpha}\neq{\bm{\beta}}$ and $\bm{\gamma}\neq{\bm{\delta}}$,

By assumption $\|{\bm{A}}\|_{{\rm op}}\leq 1$. Hence, it is sufficient to show

A.3 Technical lemmas

First note that by Hölder inequality followed by hypercontractivity on the sphere (Lemma 14),

Consider the representation (27). If $d-1\in r(\bm{\alpha},{\bm{\beta}},\bm{\gamma},{\bm{\delta}})$, then the expectation (36) is simply . Assume that $d-1\not\in r(\bm{\alpha},{\bm{\beta}},\bm{\gamma},{\bm{\delta}})$, then from the bound (37), we can decompose

Let us decompose the expectation using the representation (27):

The different terms contribute as follows in the above product. First, we can’t have $\alpha_{d-1}>0$ and $\beta_{d-1}>0$ at the same time, hence

If $j\not\in S_{\bm{\alpha}}\cup S_{{\bm{\beta}}}$,

Combining these contributions in Eq. (44) yields

Hence to prove the lemma, it is sufficient to show that $|M_{\bm{\alpha},{\bm{\beta}}}-1|=O(d^{-1/2})$.

Expanding $M_{\bm{\alpha},{\bm{\beta}}}$ yields

where $c_{j}^{\bm{\alpha}}=1$ if $\alpha_{j}>0$ and $=0$ is $\alpha_{j}=0$ (similarly for $c_{j}^{{\bm{\beta}}}$). Note that on the first line, the product can be simplified by telescoping the terms and we obtain

Appendix B Proof of Theorem 3: asymptotic characterization of KRR

In this section, we focus on the test error (we will write $f_{*}=f_{*,d}$ for simplicity):

where we recall that the kernel ridge regression solution is given by

with ${\bm{H}}=(h(\langle{\bm{x}}_{i},{\bm{x}}_{j}\rangle/d))_{ij\in[n]}$, ${\bm{y}}=(y_{1},\ldots,y_{n})$ and ${\bm{h}}({\bm{x}})=(h(\langle{\bm{x}},{\bm{x}}_{i}\rangle/d))_{i\in[n]}$.

We will decompose ${\mathcal{A}}_{d}$ into three orthogonal subspaces and bound the risk along each of them:

Recall that we can decompose the inner-product kernel in terms of Gegenbauer polynomials associated to ${\mathcal{A}}_{d}$:

Recall that we denote ${\bm{Q}}_{k}=B({\mathcal{A}}_{d},k)^{-1}{\bm{Y}}_{k}{\bm{Y}}_{k}^{\mathsf{T}}$ the matrix of the $k$-th Gegenbauer polynomial evaluated on the inner-product of the inputs. We will further denote:

By Theorem 6 in [MMM21a] (see also Proposition 6 and Corollary 1 in Section B.3.2), the high-degree component of the kernel matrices satisfy

Under the assumptions of Theorem 3, we have:

The proofs of Propositions 3 and 4 can be found in Sections B.2 and B.3 respectively. The characterization of the test error in Theorem 3 follows directly from these two propositions.

B.2 Proof of Proposition 3

For the second term, notice that by Theorem 6 in [MMM21a] (see also Proposition 6 and Corollary 1 in Section B.3.2), we have

Similarly, by Eq. (58) of Lemma 5 with ${\bm{A}}={\mathbf{I}}_{n}$, we get

Combining Eqs. (49), (50) and (51) in Eq. (48) yields the first of the three contributions:

Let us now simplify $B_{22}$: applying Theorem 6 in [MMM21a], we have

we can use Lemma 6 and simplify the expression of the different terms:

Follow the assumptions and notations in the proof of Theorem 3. We have

Follow the assumptions and notations in the proof of Theorem 3. We have

The first bound (56) follows simply from Lemma 4. For the second bound 57, we have

For the third bound (58), we follow some of the notations introduced in the proof of Lemma 4. In particular, by Sherman-Morrison-Woodbury formula, we have

Follow the assumptions of Theorem 3. We have

B.3 Convergence to expectation: proof of Proposition 4

Proposition 4 is a direct implication of the following proposition:

Under the assumptions of Theorem 3, we have

This proposition is proved in Section B.3.1, while some more technical bounds in expectation (instead of in probability, as in Section A.3) are deferred to Section B.3.2.

This directly imply the claim in Proposition 4.

Let us bound each term separately. First,

where we used that ${\rm Tr}({\bm{Q}}_{k}\bm{\Xi}{\bm{Q}}_{l}\bm{\Xi})\leq n\|{\bm{Q}}_{k}\|_{{\rm op}}\|{\bm{Q}}_{l}\|_{{\rm op}}\|\bm{\Xi}\|_{{\rm op}}^{2}$, $\|\bm{\Xi}\|_{{\rm op}}\leq\lambda^{-1}$, and Corollary 1 and Lemma 9.

Combining the bounds Eqs. (68), (69) and (70) yields

where we used the bound (81) in Proposition 7. For the high degree part, notice that ${\bm{A}}^{-1}\preceq(\mu_{d,k}{\bm{Q}}_{k}+\lambda{\mathbf{I}})^{-1}$, and therefore, by Assumption 3, there exists $\delta>0$ such that

where we used the bound (80) in Proposition 7. The cross-terms can be bounded in a similar manner. We conclude that

The term $T_{5}$ can be bounded similarly as $T_{4}$ and $T_{6}$. The proposition follows by combining bounds (65), (71), (72), (75) and (76).

B.3.2 Technical results: bounds in expectation

In this section, we gather some $L^{q}$-bounds necessary for the proof of Proposition 5 (instead of bounds in probability, as in Section A.3). We first recall some concentration results on matrices of spherical harmonics (or Fourier basis) proved in [GMMM21] and [MMM21a].

Then, there exists a constant $C>0$ such that for any $t>0$,

The following is a reformulation of Proposition 3 proved in [GMMM21] (see also Proposition 4 in [MMM21a] for a more general proof).

In the case of the hypercube ${\mathcal{A}}_{d}={\mathscrsfs Q}^{d}$, the same result holds for ${\bm{Q}}_{d-k}$.

Note that for ${\mathcal{A}}_{d}={\mathscrsfs Q}^{d}$, using Eq. (77), we have ${\bm{Q}}_{d-k}={\bm{S}}{\bm{Q}}_{k}{\bm{S}}$ and therefore the equality $\|{\bm{Q}}_{d-k}-{\mathbf{I}}_{n}\|_{{\rm op}}=\|{\bm{Q}}_{k}-{\mathbf{I}}_{n}\|_{{\rm op}}$.

Follow the same setting as Proposition 6. We have for any fixed $q\in[d]$ and any constant $\delta>0$,

Denote $r=pq_{*}/2\geq 1$. By Jensen’s inequality,

and conclude by taking $r$ sufficiently large (i.e., $p$ sufficiently large). ∎

It will be useful to state the following bound, which is a direct consequence of the above results.

which combined with Eq. (78) concludes the proof. ∎

In the proof of Proposition 5, we will use the following bounds on matrix ${\bm{L}}$:

Assume the same setting as Proposition 5. For any fixed $\delta>0$, we have

Furthermore, for any fixed integer $q\geq 1$, we have

For the second bound, we use again the matrix identity of Lemma 11 on ${\mathcal{E}}^{c}$:

On the event ${\mathcal{E}}$, we use that

which concludes the proof of this proposition. ∎

B.4 Proof of the asymptotic formula of the training error and RKHS norm

The proof is very similar to the proof for the prediction error and we will simply outline the main steps. First recall that the training error is given by

The RKHS norm of the KRR solution $\hat{f}(\cdot;{\hat{\bm{a}}}_{\lambda})$ is given by

Similarly, for the RKHS norm, we have the decomposition:

For the first term, we have by Sherman-Morrison-Woodbury formula,

We see that it is sufficient to show that

Similarly, we decompose ${\bm{y}}^{\mathsf{T}}\bm{\Xi}{\bm{y}}$ as in the proof of Proposition 5 and bound each term separately:

Similarly to the proof in Section B.1, combining the convergence in probability of the expectation in step 1, and the bound on the variance in step 2 yields the results for the training error and RKHS norm of Theorem 3.

B.5 Auxiliary lemmas

From the assumption, there exists a constant $C>0$ such that

The following proposition is a simple modification of the proof of Theorem 5.48 in [Ver10]:

The proof follows from the same argument as in the proof of Theorem 5.45 in [Ver10]. We will denote $C>0$ a generic constant that only depends on $q$. By symmetrization, we have

where the $\varepsilon_{i}$ are $n$ independent Rademacher random variables.

By noncommutative Khinchine’s inequality, we have

Denoting $\delta=4C\sqrt{\frac{\Gamma\log(\min(n,B))}{n}}$, this implies Eq. (84). Equation (85) is a simple consequence of bound (84). ∎

Finally, the following lemma provides a useful matrix algebra identity:

where ${\bm{\Pi}}={\mathbf{I}}_{n}-[{\bm{V}}^{-1/2}{\bm{P}}][{\bm{V}}^{-1/2}{\bm{P}}]^{\dagger}$, ${\bm{P}}={\mathbf{I}}_{n}-({\bm{U}}^{\dagger})^{\mathsf{T}}{\bm{U}}^{\mathsf{T}}$ and ${\bm{U}}^{\dagger}$ is the Moore-Penrose inverse of ${\bm{U}}$ (here by assumption, ${\bm{U}}^{\dagger}=({\bm{U}}^{\mathsf{T}}{\bm{U}})^{-1}{\bm{U}}^{\mathsf{T}}$).

Denote for convenience ${\bm{S}}={\bm{V}}^{-1/2}$. This identity comes from the observation that

and by repeatedly applying Sherman-Morrison-Woodbury identity. First,

We can then use the identity again on its inverse:

By the definition of the pseudo-inverse, ${\bm{P}}={\mathbf{I}}-{\bm{U}}({\bm{U}}^{\mathsf{T}}{\bm{U}})^{-1}{\bm{U}}^{\mathsf{T}}={\bm{P}}^{2}$. We apply a third time the SMW formula:

Appendix C Technical background

The dimension of each subspace is given by

C.1.2 Gegenbauer polynomials

We will use the following properties of Gegenbauer polynomials

These properties imply that —up to a constant— $Q_{k}^{(d)}(\langle{\bm{x}},{\bm{y}}\rangle)$ is a representation of the projector onto the subspace of degree -$k$ spherical harmonics

then we have the following equation holds in $L^{2}([-\sqrt{d},\sqrt{d}],\tau^{1}_{d-1})$ sense

By rotational invariance, the space $V_{k}$ of homogeneous polynomials of degree $k$ is an eigenspace of $\mathscrsfs{H}_{d}$, and we will denote the corresponding eigenvalue by $\xi_{d,k}(h_{d})$. In other words $\mathscrsfs{H}_{d}f({\bm{x}})\equiv\sum_{k=0}^{\infty}\xi_{d,k}(h_{d}){\overline{\mathsf{P}}}_{k}f$. The eigenvalues can be computed via

C.1.3 Hermite polynomials

Here and below, for $P$ a polynomial, ${\rm Coeff}\{P(x)\}$ is the vector of the coefficients of $P$. As a consequence, for any fixed integer $k$, we have

where $\mu_{k}(\bar{\sigma})$ and $\xi_{d,k}(\bar{\sigma})$ are given in Eq. (98) and (94).

C.2 Functions on the hypercube

Fourier analysis on the hypercube is a well studied subject [O’D14]. The purpose of this section is to introduce some notations that make the correspondence with proofs on the sphere straightforward. For convenience, we will adopt the same notations as for their spherical case.

It is easy to verify that (notice that $x_{i}^{k}=x_{i}$ if $k$ is odd and $x_{i}^{k}=1$ if $k$ is even)

C.2.2 Hypercubic Gegenbauer

Notice that the right hand side only depends on $\langle{\bm{x}},{\bm{y}}\rangle$ and therefore these polynomials are uniquely defined. In particular,

Notice that by weak convergence of $\langle\bm{1},{\bm{x}}\rangle/\sqrt{d}$ to the normal distribution, we have also convergence of the (rescaled) hypercubic Gegenbauer polynomials to the Hermite polynomials, i.e., for any fixed $k$, we have

C.3 Hypercontractivity of the uniform distribution on the sphere and the hypercube

By Holder’s inequality, we have $\|f\|_{L^{p}}\leq\|f\|_{L^{q}}$ for any $f$ and any $p\leq q$. The reverse inequality does not hold in general, even up to a constant. However, for some measures, the reverse inequality will hold for some sufficiently nice functions. These measures satisfy the celebrated hypercontractivity properties [Gro75, Bon70, Bec75, Bec92].

Besides this classical result, we will also use the following simple observation:

Note that for any $S\subseteq[d]$, we have Y_{S}({\bm{x}})=\big{(}\prod_{i\in[d]}x_{i}\big{)}\cdot Y_{S^{c}}({\bm{x}}). Hence, we have

Finally, we have the following similar hypercontractivity property on the sphere: