Kernel Alignment Risk Estimator: Risk Prediction from Training Data

Introduction

Kernel Ridge Regression (KRR) is a widely used statistical method to learn a function from its values on a training set . It is a non-parametric generalization of linear regression to infinite-dimensional feature spaces. Given a positive-definite kernel function $K$ and (noisy) observations $y^{\epsilon}$ of a true function $f^{*}$ at a list of points $X=x_{1},\ldots,x_{N}$, the $\lambda$-KRR estimator $\hat{f}^{\epsilon}_{\lambda}$ of $f^{*}$ is defined by

Despite decades of intense mathematical progress, the rigorous analysis of the generalization of kernel methods remains a very active and challenging area of research. In recent years, many new kernels have been introduced for both regression and classification tasks; notably, a large number of kernels have been discovered in the context of deep learning, in particular through the so-called Scattering Transform , and in close connection with deep neural networks , yielding ever-improving performance for various practical tasks . Currently, theoretical tools to select the relevant kernel for a given task, i.e. to minimize the generalization error, are however lacking.

While a number of bounds for the risk of Linear Ridge Regression (LRR) or KRR exist, most focus on the rate of convergence of the risk: these estimates typically involve constant factors which are difficult to control in practice. Recently, a number of more precise estimates have been given ; however, these estimates typically require a priori knowledge of the data distribution. It remains a challenge to have estimates based on the training data alone, enabling one to make informed decisions on the choices of the ridge and of the kernel.

We introduce the Signal Capture Threshold (SCT) $\vartheta(\lambda,N,K,\pi)$, which is determined by the ridge $\lambda$, the size of the training set $N$, the kernel $K$, and the observations distribution $\pi$ (more precisely, the dependence on $\pi$ is only through its first two moments). We give approximations for the expectation and variance of the KRR predictor in terms of the SCT.

Decomposing $f^{*}$ along the kernel principal components of the data distribution, we observe that in expectation, the predictor $\hat{f}^{\epsilon}_{\lambda}$ captures only the signal along the principal components with eigenvalues larger than the SCT. If $N$ increases or $\lambda$ decreases, the SCT $\vartheta$ shrinks, allowing the predictor to capture more signal. At the same time, the variance of $\hat{f}^{\epsilon}_{\lambda}$ scales with the derivative $\partial_{\lambda}\vartheta$, which grows as $\lambda\to 0$, supporting the classical bias-variance tradeoff picture .

We give an explicit approximation for the expected MSE risk $R^{\epsilon}(\hat{f}^{\epsilon}_{\lambda})$ and empirical MSE risk $\hat{R}^{\epsilon}(\hat{f}^{\epsilon}_{\lambda})$ for an arbitrary continuous true function $f^{*}$. We find that, surprisingly, the expected risk and expected empirical risk are approximately related by

We introduce the Kernel Alignment Risk Estimator (KARE) as the ratio $\rho$ defined by

where $G$ is the Gram matrix of $K$ on the observations. We show that the KARE approximates the expected risk; unlike the SCT, it is agnostic of the true data distribution. This result follows from the fact that $\vartheta(\lambda)\approx 1/m_{G}(-\lambda)$, where $m_{G}$ is the Stieltjes Transform of the Gram matrix $\frac{1}{N}G$.

Empirically, we find that the KARE predicts the risk on the Higgs and MNIST datasets. We see empirically that our results extend extremely well beyond the Gaussian observation setting, thus supporting our universality assumption (see Figure 1).

Our proofs (see the Appendix) rely on a finite-size analysis of generalized Wishart matrices, in particular the complex Stieltjes transform $m_{G}(z)$, evaluated at $z=-\lambda$, and on fixed-point arguments.

2 Related Works

The theoretical analysis of the risk of KRR has seen tremendous developments in the recent years. In particular, a number of upper and lower bounds for kernel risk have been obtained in various settings: notably, convergence rates (i.e. without control of the constant factors) are obtained in general settings. This allows one to abstract away a number of details about the kernels (e.g. the lengthscale), which don’t influence the asymptotic rates. However, this does not give access to the risk at finite data size (crucial to pick e.g. the correct lengthscale or the NTK depth ).

A number of recent results have given precise descriptions of the risk for ridge regression , for random features , and in relation to neural networks . These results rely on the analysis of the asymptotic spectrum of general Wishart random matrices, in particular through the Stieltjes transform . The limiting Stieltjes transform can be recovered from the formula for the product of freely independent matrices . To extend these asymptotic results to finite-size settings, we generalize and adapt the results of .

While these techniques have given simple formulae for the KRR predictor expectation, approximating its variance has remained more challenging. For this reason the description of the expected risk in is stated as a conjecture. In only the bias component of the risk is approximated. In the expected risk is given only for random true functions (in a Bayesian setting) with a specific covariance. In , the expected risk follows from a heuristic spectral analysis combining a PDE approximation and replica tricks. In this paper, we approximate the variance of the predictor along the principal components, giving an approximation of the risk for any continuous true function.

The SCT is related to a number of objects from previous works, such as the effective dimension of , the companion Stieltjes transform of , and particularly the effective ridge of . The SCT can actually be viewed as a direct translation to the KRR risk setting of .

3 Outline

In Section 2, we first introduce the Kernel Ridge Regression (KRR) predictor in functional space (Section 2.1) and formulate its train error and risk for random observations (Section 2.2).

The rest of the paper is then devoted to obtaining approximations for the KRR risk. In Section 3,the Signal Capture Threshold (SCT) is introduced and used to study the mean and variance of the KRR predictor (Sections 3.1 and 3.2). An approximation of the SCT in terms of the observed data is then given (Section 3.4). In Section 4, the expected risk and the expected empirical risk are approximated in terms of the SCT and its derivative w.r.t. the ridge $\lambda$. The SCT approximation of Section 3.4, together with the estimates of Section 4.1, leads to an approximation of the KRR risk by the Kernel Alignment Risk Estimator (KARE).

Setup

The regression problem is now stated as follows: given noisy observations $y^{\epsilon}_{i}=o_{i}\left(f^{*}\right)+\epsilon e_{i}$ with i.i.d. centered noises $e_{1},\ldots,e_{N}$ of unit variance, how can one reconstruct $f^{*}$?

We call the $N\times N$ matrix $G=\mathcal{O}K\mathcal{O}^{T}$ the Gram matrix: in the classical setting, when the observations are $o_{i}=\delta_{x_{i}}$ (with $\delta_{x}(f)=f(x)$), $G$ is the usual Gram matrix, i.e. $G_{ij}=K(x_{i},x_{j})$.

2 Training Error and Risk

We consider the least-squares error (MSE loss) of the KRR predictor, taking into account randomness of: (1) the test point, random observation $o$ to which is added a noise $\epsilon e$ (2) the training data, made of $N$ observations $o_{i}$ plus noises $\epsilon e_{i}\sim\nu$, where $o,o_{1},\ldots,o_{n}\sim\pi$ and $e,e_{1},\ldots,e_{N}$ are i.i.d. The expected risk of the KRR predictor is thus taken w.r.t. the test and training observations and their noises. Unless otherwise specified, the expectations are taken w.r.t. all these sources of randomness.

For (fixed) observations $o_{1},\ldots,o_{N}$, the empirical risk or training error of the KRR predictor $\hat{f}^{\epsilon}_{\lambda}$ is

For a random observation $o$ sampled from $\pi$ and a noise $\epsilon e$ (where $e\sim\nu$ is centered of unit variance as before), the risk $R^{\epsilon}(\hat{f}^{\epsilon}_{\lambda})$ of the KRR predictor $\hat{f}^{\epsilon}_{\lambda}$ is defined by

From now on, we will assume that $\left\langle\cdot,\cdot\right\rangle_{S}$ is a scalar product; note that in the classical setting, when $o$ is the evaluation of $f^{*}$ at a point $x\in\Omega$ with $x\sim\sigma$, the $S$-norm is given by $\|f\|_{S}^{2}=\int_{\Omega}f(x)^{2}\sigma(dx)$.

The following three operators $\mathcal{C}\to\mathcal{C}$ are central to our analysis:

The KRR reconstruction operator $A_{\lambda}:\mathcal{C}\to\mathcal{C}$, the KRR Integral Operator $T_{K}:\mathcal{C}\to\mathcal{C}$, and its empirical version $T^{N}_{K}:\mathcal{C}\to\mathcal{C}$ are defined by

Note that in the noiseless regime (i.e. when $\epsilon=0$), we have $\hat{f}^{\epsilon}_{\lambda}\big|_{\epsilon=0}=A_{\lambda}f^{*}$. Also note that $A_{\lambda}$ and $T^{N}_{K}$ are random operators, as they depend on the random observations. The operator $T_{K}$ is the natural generalization to our framework of the integration operator $f\mapsto\int K(x,\cdot)f(x)\sigma(dx)$, which is defined with random observations $\delta_{x}$ with $x\sim\sigma$ in the classical setting.

The reconstruction and empirical integral operators are linked by $A_{\lambda}=T^{N}_{K}(T^{N}_{K}+\lambda I_{\mathcal{C}})^{-1}$, which follows from the identity $\left(\frac{1}{N}\mathcal{O}K\mathcal{O}^{T}+\lambda I_{N}\right)^{-1}\mathcal{O}=\mathcal{O}\left(\frac{1}{N}K\mathcal{O}^{T}\mathcal{O}+\lambda I_{\mathcal{C}}\right)^{-1}$. As $N\to\infty$, we have that $T^{N}_{K}\to T_{K}$, and it follows that

3 Eigendecomposition of the Kernel

4 Gaussianity Assumption

As far as one is concerned with the first two moments of the $A_{\lambda}$ operator, for large but finite $N$, one can assume that the observations $o_{1},\ldots,o_{N}$ are Gaussian, i.e. that for any tuple of functions $(f_{1},\ldots,f_{N})$, the vector $\left(o_{1}(f_{1}),\ldots,o_{N}(f_{N})\right)$ is a Gaussian vector.

Though our proofs use this assumption, the ideas in suggest a path to extend them beyond the Gaussian case, where our numerical experiments (see Figure 1) suggest that our results remain true.

Predictor Moments and Signal Capture Threshold

A central tool in our analysis of the KRR predictor $\hat{f}^{\epsilon}_{\lambda}$ is the Signal Capture Threshold (SCT):

For $\lambda>0$, the Signal Capture Threshold $\vartheta(\lambda)=\vartheta(\lambda,N,K,\pi)$ is the unique positive solution (see Section B.2 in the Appendix) to the equation:

In this section, we use $\vartheta(\lambda)$ and the derivative $\partial_{\lambda}\vartheta(\lambda)$ for the estimation of the mean and variance of the KRR predictor $\hat{f}^{\epsilon}_{\lambda}$ upon which the Kernel Alignment Risk Estimator of Section 4 is based.

The expected KRR predictor can be expressed in terms of the expected reconstruction operator $A_{\lambda}$

for a polynomial $\boldsymbol{P}_{0}$ with nonnegative coefficients and $\boldsymbol{P}_{0}(0)=0$.

More generally, if we decompose a true function $f^{*}$ along the principal components (i.e. eigenfunctions) of $T_{K}$, the signal along the $k$-th principal component $f^{(k)}$ is captured whenever the corresponding eigenvalue $d_{k}\gg\vartheta(\lambda)$ and lost when $d_{k}\ll\vartheta(\lambda)$.

2 Variance of the predictor

We now estimate the variance of $\hat{f}^{\epsilon}_{\lambda}$ along each principal component in terms of the SCT $\vartheta(\lambda)$ and its derivative $\partial_{\lambda}\vartheta(\lambda)$. Along the eigenfunction $f^{(k)}$, the variance is estimated by $V_{k}$, where

There is a constant $\boldsymbol{C}_{1}>0$ and a polynomial $\boldsymbol{P}_{1}$ with nonnegative coefficients and with $\boldsymbol{P}_{1}(0)=0$ such that

As shown in Section 4.1, understanding the variance along the principal components (rather than the covariances between the principal components) is enough to describe the risk.

3 Behavior of the SCT

The behavior of the SCT can be controlled by the following (agnostic of the exact spectrum of $T_{K}$)

moreover $\vartheta(\lambda,N)$ is decreasing as a function of $N$.

As $\lambda\to 0$, the above upper bound for $\partial_{\lambda}\vartheta$ becomes useless. Still, assuming that the spectrum of $K$ has a sufficiently fast power-law decay, we get:

If $d_{k}=\Theta(k^{-\beta})$ for some $\beta>1$, there exist $c_{0},c_{1},c_{2}>0$ such that for any $\lambda>0$

4 Approximation of the SCT from the training data

Likewise, we have $\partial_{\lambda}\vartheta\approx\left(\partial_{z}m_{G}(z)/m_{G}(z)^{2}\right)|_{z=-\lambda}$, as shown in the Appendix.

Risk Prediction with KARE

The Kernel Alignment Risk Estimator (KARE) $\rho$ is defined by

In the following, using Theorems 1 and 22, we give an approximation for the expected risk and expected empirical risk in terms of the SCT and the true function $f^{*}$. This yields the important relation (2) in Section 4.2, which shows that the KARE can be used to efficiently approximate the kernel risk.

The expected risk is approximated, in terms of the SCT and the true function $f^{*}$, by

There exists a constant $\boldsymbol{C}_{2}>0$ and a polynomial $\boldsymbol{P}_{2}$ with nonnegative coefficients and with $\boldsymbol{P}_{2}(0)=0$, such that we have

(Sketch; the full proof is given in the Appendix). From the bias-variance decomposition:

In a Bayesian setting, assuming that $f^{*}$ is random with zero mean and covariance kernel $\Sigma$, the optimal choices for the KRR predictor are $K=\Sigma$ and $\lambda=\nicefrac{{\epsilon^{2}}}{{N}}$ (see Section B.7 in the Appendix). When $K=\Sigma$ and $\lambda=\nicefrac{{\epsilon^{2}}}{{N}}$, the formula of Theorem 6 simplifies (see Corollary in the Appendix) to

The empirical risk (or train error) $\hat{R}^{\epsilon}(\hat{f}^{\epsilon}_{\lambda})=\lambda^{2}(y^{\epsilon})^{T}(\frac{1}{N}G+\lambda I_{N})^{-2}y^{\epsilon}$ can be analyzed with the same theoretical tools. Its approximation in terms of the SCT is given as follows:

There exists a constant $\boldsymbol{C}_{3}>0$ and a polynomial $\boldsymbol{P}_{3}$ with nonnegative coefficients and with $\boldsymbol{P}_{3}(0)=0$ such that we have

2 KARE: Kernel Alignment Risk Estimator

While the above approximations (Theorems 6 and 7) for the expected risk and empirical risk depend on $f^{*}$, their combination yields the following relation, which is surprisingly independent of $f^{*}$:

Since $\vartheta$ can be approximated from the training set (see Proposition 5), so can the expected risk. Assuming that the risk and empirical risk concentrate around their expectations, we get the KARE:

Note that both $\rho$ and $\varrho$ are invariant (as is the risk) under the simultaneous rescaling $K,\lambda\leadsto\alpha K,\alpha\lambda$.

The KARE can be used to optimize the risk over the space of kernels, for instance to choose the ridge and length-scale. The most popular kernel selection techniques are (see Figure 3):

Cross-validation: accurate estimator of the risk on a test set, but costly to optimize (the predictor must be recomputed for each kernel and differentiating it in the space of kernels is hard).

Kernel likelihood (Chapter 5 of ): efficient to optimize and takes into account the ridge, but not a risk estimator; unlike the risk, not invariant under the simultaneous rescaling $K,\lambda\leadsto\alpha K,\alpha\lambda$.

Classical kernel alignment : very efficient to optimize and scale invariant, but not a risk estimator, not sensitive to small eigenvalues and inadequate to select hyperparameters such as the ridge.

The KARE combines the best features of the three above techniques:

it can be computed efficiently on the training data, and optimized over the space of kernels;

like the risk, it is invariant under the simultaneous rescaling $K,\lambda\leadsto\alpha K,\alpha\lambda$;

it is sensitive to the small Gram matrix eigenvalues and to the ridge $\lambda$.

Conclusion

In this paper, we introduce new techniques to study the Kernel Ridge Regression (KRR) predictor and its risk. We obtain new precise estimates for the test and train error in terms of a new object, the Signal Capture Threshold (SCT), which identifies the components of a true function that are being learned by the KRR: our estimates reveal a remarkable relation, which leads one to the Kernel Alignment Risk Estimator (KARE). The KARE is a new efficient way to estimate the risk of a kernel predictor based on the training data only. Numerically, we observe that the KARE gives a very accurate prediction of the risk for Higgs and MNIST datasets for a variety of classical kernels.

Broader Impact

This work is fundamental and may be used in any research area using Kernel methods, possibly leading to indirect social impacts. However, we do not predict any direct social impact.

Acknowledgements

The authors wish to thank A. Montanari and M. Wyart for useful discussions. This work is partly supported by the ERC SG CONSTAMIS. C. Hongler acknowledges support from the Blavatnik Family Foundation, the Latsis Foundation, and the the NCCR Swissmap.

References

Appendix

In Section A, we present the details for the numerical results presented in the main text (and in the Appendix) and we present additional experiments and some discussions.

In Section B, we present the proofs of the mathematical results presented in the main text.

Appendix A Numerical Results

For the MNIST dataset. We sample $N$ images of digits $7$ and $9$ from the MNIST training dataset (image size $d=24\times 24$, edge pixels cropped, all pixels rescaled down to $$and recentered around the mean value) and label each of them with$+1$and$-1$labels. We perform KRR with various ridge$\lambda$on this dataset with the selected kernel$k$times and calculate the MSE training error, risk, and the KARE for every trial ($k=10$for small$N$and$k=5$for$N=2000$). The risk is approximated using other$N_{2}=1000$ random samples of the MNIST training data.

For the Higgs Dataset. We randomly choose $N$ samples among those that do not have any missing features marked with $-999$ from the Higgs training dataset. The samples have $d=31$ features, and we normalize each feature column down to $$by dividing by the maximum absolute value observed among the selected samples. We replace the categorical labels ‘s’ and ‘b’ with regression values$+1$and$-1$respectively and perform KRR with various ridge$\lambda$. We repeat this procedure$k$times, which corresponds to sampling$k$different training datasets of$N=1000$samples to perform kernel regression, and calculate the MSE training error, the risk, and the KARE for every trial ($k=10$for small$N$and$k=5$for$N=1000$). The risk is approximated using other$N_{2}=1000$ random samples of the Higgs training data.

A.2 KARE predicts risk for various Kernels

A.3 KRR predictor in function space

A.4 KARE predicts risk in average for small NN

A.5 SCT and its behavior

for $k\geq 1$. In particular, we have $n_{d}(0)=1,n_{d}(1)={d\choose 1},n_{d}(2)={d\choose 2}+d,\ldots$. In general, $n_{d}(k)$ is the number of ways to partition $k$ into $d$ non-negative integers.

The true SCT is therefore approximated solving the following equation numerically

Note that in the Figure 2 in the main text, we limit the approximation to $k=10$ for $d=20$ because the multiplicity $n_{d}(k)$ grows polynomially with $d^{k}$.

Appendix B Proofs

Throughout our proofs, we will frequently rely on a polynomial analogue of the big-O notation, which we call big-P:

For two functions $f$ and $g$ (of one or several variables, defined on an arbitrary common domain $\mathcal{D}$), we write $f=\mathcal{P}(g)$ if $g$ is nonnegative over $\mathcal{D}$ and there exists a polynomial $\mathbf{P}$ with nonnegative coefficients and $\mathbf{P}(0)=0$ such that $|f|\leq\mathbf{P}(g)$ over $\mathcal{D}$.

Note that the big-O notation corresponds to the case when the polynomial $\mathbf{P}$ is of degree at most one.

B.1 Objects of Interest and general strategy

The central object of our analysis is the $N\times N$ Gram matrix $\mathcal{O}K\mathcal{O}^{T}$, in particular the related Stieltjes transform:

From now on, we denote $\vartheta(-z)$ by $\vartheta$. Note that here, in the Appendix, we use the resolvent notation: in particular the KRR reconstruction operator $A_{\lambda}$ is equal to $A(-\lambda)$.

Using the spectral decomposition of $K$, the entries of $\mathcal{O}K\mathcal{O}^{T}$ are given by:

where the sum converges absolutely (thanks to the trace assumption on $K$) and the entries of $A$ are then given by:

where $\mathcal{O}_{\cdot k}=\left(o_{i}(f^{(k)})\right)_{i=1,\ldots,N}$.

B.1.2 Shermann-Morrison Formula

As a result of Equations (8) and (9), the diagonal entries of the operator $A(z)=\frac{1}{N}K\mathcal{O}^{T}B(z)^{-1}\mathcal{O}$ are equal to

Another important observation is that the Stieltjes transform $m(z)$ and the $g_{k}(z)$ are closely related.

Dividing both sides by $N$ and using Equation (10), we obtain

B.2 Concentration of the Stieltjes Transform

where the well-posedness of the two infinite sums of the r.h.s is granted by the fact that:

being the difference of two absolutely convergent series, the second sum is also absolutely convergent.

Regarding the concentration of the $g_{k}$’s around $m$, we have the following result:

where $\boldsymbol{c}_{s}$ only depends on $s$.

since the derivative $g^{\prime}_{k}(z)$ of $g_{k}(z)$ is equal to $\frac{1}{N}\mathcal{O}_{\cdot k}^{T}B(z)^{-2}\mathcal{O}_{\cdot k}$. As a result, we obtain $\left|m-m_{(k)}\right|^{2s}=\frac{1}{N^{2s}}\frac{d_{k}^{2s}\left|g^{\prime}_{k}\right|^{2s}}{\left|1+d_{k}g_{k}\right|^{2s}}$. Using the fact that $\left|1+d_{k}g_{k}\right|\geq\left|d_{k}g_{k}\right|$ since $\Re\left(g_{k}\right)\geq 0$,

Let $\boldsymbol{B}=\left(B_{(k)},\overline{B_{(k)}},\ldots,B_{(k)},\overline{B_{(k)}}\right)$ and let us denote by $\boldsymbol{B}(i)$ the $i^{th}$ element of $\boldsymbol{B}$. Using Wick’s formula (Lemma 26), we have

where we recall that $\mathfrak{S}_{2s}^{\dagger}$ is the set of permutations with no fixed points and the product over $i$ is taken according to the order given by the cycle $c$ and does not depend on the starting point. Using the fact that the eigenvalues of $B_{(k)}$ are of the form $\nicefrac{{1}}{{\left(\lambda_{i}-z\right)}}$ with $\lambda_{i}\geq 0$,

Note that, since $\sigma\in\mathfrak{S}_{2s}^{\dagger},$ it has no fixed point, hence $\boldsymbol{c}(\sigma)\leq s$ and thus $K_{s}:=\sup_{N}\sum_{\sigma\in\mathfrak{S}_{2s}^{\dagger}}\frac{2^{2s-\boldsymbol{c}(\sigma)}}{N^{s-\boldsymbol{c}(\sigma)}}$ is finite. This yields the inequality

where $\boldsymbol{c}_{s}=2^{2s-1}\left[1+K_{s}\right]$. ∎

where $\boldsymbol{c}_{s}$ is the same constant as in Lemma 9.

The second bound is a direct consequence of the first one, Lemma 9 and convexity. It remains to prove the first bound. Recall Equation (13)

and hence, using a generalization of Cauchy-Schwarz inequality (Lemma 28), by:

where $\boldsymbol{c}_{1}$ is the constant in Lemma 9.

First bound: Following similar ideas to the one which provided Equation (13), notice that

B.3 Properties of the effective dimension and SCT

We begin with general properties on the Signal Capture Threshold $\vartheta$ (which depends on $\lambda,N$ and on the eigenvalues $d_{k}$ of $T_{K}$), valid for any kernel $K$.

moreover $\vartheta(\lambda,N)$ is decreasing as a function of $N$ and $\partial_{\lambda}\vartheta(\lambda,N)$ is decreasing as a function of $\lambda$.

Recall that $\vartheta(\lambda)$ is the unique positive real number such that

Differentiating Equation (7), the derivative $\partial_{\lambda}\vartheta(\lambda)$ is given by:

Using the fact that $T_{K}\left(T_{K}+\vartheta(\lambda)I_{\mathcal{C}}\right)^{-1}\leq I_{\mathcal{C}}$, one has

Inverting this inequality yields the desired inequalities.

In order to study the variation of $\vartheta(\lambda,N)$ as a function of $N$, we take the derivatives of Equation (7) w.r.t $\lambda$ and $N$, and notice that

In particular, since $\vartheta>\lambda$ and $\partial_{\lambda}\vartheta\geq 1$, we get that $\partial_{N}\vartheta(\lambda,N)<0$ hence $\vartheta(\lambda,N)$ is decreasing as a function of $N$.

Finally, we conclude by noting that since $\partial_{\lambda}\vartheta(\lambda,N)>0$, $\vartheta(\lambda,N)$ is an increasing function of $\lambda$ and thus, from the Equation (15) we have that $\partial_{\lambda}\vartheta(\lambda,N)$ is decreasing as a function of $\vartheta$ and thus as a function of $\lambda$.

B.3.2 Bounds under polynomial decay hypothesis

For any $\lambda>0$, the SCT is the unique solution of $\vartheta\left(\lambda,N\right)=\lambda+\frac{\vartheta\left(\lambda,N\right)}{N}\mathcal{N}(\vartheta(\lambda,N)).$ In particular, $\vartheta(0,N)$ is the unique solution of $\mathcal{N}(\vartheta(0,N))=N$.

Since $\mathcal{N}(t)$ is decreasing from $\infty$ to $0$, in order to study the asymptotic behavior of $\vartheta(0,N)$ as $N$ goes to infinity, one has to understand the rate of explosion of $\mathcal{N}(t)$ as $t$ goes to zero.

If $d_{k}=\Theta(k^{-\beta})$ with $\beta>1$, then $\mathcal{N}(t)=\Theta(t^{-\frac{1}{\beta}})$ when $t\to 0$.

If $d_{k}=\Theta(k^{-\beta})$ with $\beta>1$, then $\vartheta(0,N)=\Theta\left(N^{-\beta}\right)$.

With no assumption on the spectrum of $T_{K}$, the upper bound for the derivative of the SCT $\partial_{\lambda}\vartheta$ obtained in Proposition 3, becomes useless in the ridgeless limit $\lambda\to 0$. Yet, with the assumption of power-law decay of the eigenvalues of $T_{K}$ we can refine the bound with a meaningful one. In order to obtain this we first prove a technical lemma.

If $d_{k}=\Theta(k^{-\beta})$ with $\beta>1$, then $\sup_{N}\partial_{\lambda}\vartheta(0,N)<\infty$.

The derivative of the SCT with respect to $\lambda$ at $\lambda=0$ is given by:

Set $\alpha>1$, then for all $d_{k}\in[\alpha^{-1}t,\alpha t]$, we have that $\frac{d_{k}}{(t+d_{k})^{2}}\geq\frac{\alpha t}{(t+\alpha t)^{2}}=\frac{\alpha}{(1+\alpha)^{2}}\frac{1}{t}$. Thus,

Now, using Lemma 14, we are going to find a value of $\alpha$ such that $\#\{k\mid\alpha^{-1}\vartheta(0,N)<d_{k}<\alpha\vartheta(0,N)\}\geq cN$ for some universal constant $c$: this will conclude the proof.

$1\leq\partial_{\lambda}\vartheta(\lambda,N)\leq c$,

We start by proving the inequalities for the derivative of the SCT $\partial_{\lambda}\vartheta(\lambda,N)$. The left side of the inequality has already been proven in Proposition 3. For the right side, from Proposition 3, the derivative $\partial_{\lambda}\vartheta(\lambda,N)$ is decreasing in $\lambda$. In particular, by Lemma 15, $\partial_{\lambda}\vartheta(\lambda,N)\leq\sup_{N}\partial_{\lambda}\vartheta(0,N)<\infty.$ Thus, the right side holds with $c:=\sup_{N}\partial_{\lambda}\vartheta(0,N)$.

B.4 The Operator A⁡(z)A(z)

We have now the tools to describe the moments of the operator $A(z)$ which allow us to describe the moments of the predictor $\hat{f}_{\lambda}$.

using the big-P notation of Definition 5.

Note that in particular since the polynomial implicitly embedded in $\mathcal{P}$ vanishes at $0$, the right hand side tends to $0$ as $N\to\infty$.

Off-Diagonal terms: By a symmetry argument, we show that the off-diagonal terms are null. Consider the map $s_{k}:\mathcal{C}\to\mathcal{C}$ defined by $s_{k}:f\mapsto f-2\left\langle f,f^{(k)}\right\rangle_{S}f^{(k)}$, and note that $s_{k}(f^{(m)})=f^{(m)}$ if $m\neq k$ and $s_{k}(f^{(k)})=-f^{(k)}$. The map $s_{k}$ is a symmetry for the observations, i.e. for any observations $o_{1},\ldots,o_{N}$, and any functions $f_{1},\ldots,f_{N}$, the vector $(o_{i}(s_{k}(f_{i}))_{i=1,\ldots,N}$ and $(o_{i}(f_{i}))_{i=1,\ldots,N}$ have the same law. Thus, the sampling operator $\mathcal{O}$ and the operator $\mathcal{O}s_{k}$ have the same law, hence so do $A(z)$ and $A^{s_{k}}(z)$, where

Diagonal terms: Using Equation 10, we have

From this, using the fact that $\Re(g_{k})>0$, we obtain

Using Proposition 11, we can bound the first fraction by

Finally, putting everything together, we get:

B.4.2 Variance

thus, recalling that $g_{(k)}=\frac{1}{N}\mathcal{O}_{.,k}^{T}B_{(k)}(z)^{-1}\mathcal{O}_{.,k}$, we have

Thus, we obtain the following formula for the off-diagonal entry:

where $\boldsymbol{a}_{s},\boldsymbol{b}_{s}$ only depend on $s$.

where $\boldsymbol{c}_{1}$ is as in Proposition 11.

We use Proposition 11 and Lemma 30 twice to obtain

Using Formula (10) for the diagonal entries of $A$, we have:

Using Proposition 10, the absolute value of the r.h.s. can now be bounded by

Using the same notation as in the proof of Proposition 20,

We can bound the terms in the r.h.s. of the above by applying Proposition 10 and Lemma 18:

Let $\boldsymbol{r}_{s}=\max\{2^{2s-1}\boldsymbol{c}_{s},\boldsymbol{a}_{s}\}$ and $\boldsymbol{t}_{s}=\max\{2^{2s-1}\boldsymbol{c}_{s},\boldsymbol{b}_{s}\}$; then putting the pieces together we have

And finally, putting all the pieces together, we have

We can now describe the variance of the predictor. The variance of the predictor along the eigenfunction $f^{(k)}$ is estimated by $V_{k}$, where

There is a constant $\boldsymbol{C}_{1}>0$ such that, with the notation of Definition 5, we have

Using the law of total variance, we decompose the variance with respect to the observations $\mathcal{O}$ and the vector of noise $E=(e_{1},\dots,e_{N})^{T}$

Since the randomness is now only on $A$ through $\mathcal{O}$, from now on, we will lighten the notation by sometimes omitting the $\mathcal{O}$ dependence in the expectations.

We first show how the approximation $V_{k}(f^{*},\lambda,N,\epsilon)$ appears, and then establish the bounds which allow one to study the quality of this approximation.

Approximations: Decomposing the true function along the principal components $f^{*}=\sum_{k=1}^{\infty}b_{k}f^{(k)}$ with $b_{k}=\left\langle f^{(k)},f^{*}\right\rangle_{S}$, we have

Combining Equations 20 and 21, we obtain the approximation

Second term: To approximate, we apply Cauchy’s inequality to Equation (17) of Theorem 17:

By using the fact that $1\leq\left|\partial_{\lambda}\vartheta(\lambda)\right|$ (see Proposition 3), we have that

Finally, by putting the bounds for the two terms together we have

B.5 Expected Risk

Similarly to the proof of Theorem 22, we explain how the approximation of the expected arises, then we establish the bounds which allow one to study the quality of this approximation.

Thus the variance term is approximately equal to:

Noting that from Equation 15, we have $(\partial_{\lambda}\vartheta(\lambda)-1)=\frac{\partial_{\lambda}\vartheta(\lambda)}{N}\sum_{k=1}^{\infty}\frac{d_{k}^{2}}{(\vartheta(\lambda)+d_{k})^{2}},$ we get:

Hence, we get the following approximation of the variance term:

Putting the approximations of the bias and variance terms together, we obtain:

Now, we explain how to quantify the quality of the approximations, and thus how to get the bound stated in the theorem. Recall that, using the bias-variance decomposition, we split the expected risk into two terms, the bias term and the variance term. We show now that:

Combining the two inequations, and using the fact that $1\leq\partial_{\lambda}\vartheta(\lambda)$, we then get the desired inequality.

We decompose the true function $f^{*}$ into $f^{*}=\sum_{k=1}^{\infty}b_{k}f^{(k)}$ for $b_{k}=\left\langle f^{*},f^{(k)}\right\rangle_{S}$, and obtain

By the triangular inequality, we get that

Variance term: For the second term, recall that $(\partial_{\lambda}\vartheta(\lambda)-1)=\frac{\partial_{\lambda}\vartheta(\lambda)}{N}\sum_{k=1}^{\infty}\frac{d_{k}^{2}}{(\vartheta(\lambda)+d_{k})^{2}},$ and that

Using Theorem 22, we can control the terms in the first series: there is a constant $\boldsymbol{C}_{1}>0$ such that

whereas for the second series, as explained already above, we have

Finally, putting the pieces together, we conclude. ∎

B.6 Expected Empirical Risk

The expected empirical risk can be approximated as follows:

A small computation allows one to show that:

Using the definition of $y^{\epsilon}$ and the fact that the noise on the labels is centered and independent from the observations, this yields:

Similarly to the proof of Theorem 22, we explain how the approximation of the expected empirical risk appears, then we establish the bounds which allow one to study the quality of this approximation.

The second term can be approximated using Proposition 11 and Lemma 27: this yields

Hence, putting the two approximations together, the expected empirical risk is approximated by:

Now, we explain how to quantify the quality of the approximations, and thus how to get the bound stated in the theorem. Recall that, we split the expected empirical risk into two terms.

First term: We have already seen in Theorem 22 that by applying Lemma 27 to Equation (17) of Theorem 17 we get

Second Term: Using Proposition 11 and Lemma 27:

B.7 Bayesian Setting

Differentiating w.r.t. $M_{x}$, we obtain that the above error is minimized when

In other terms, in this Bayesian setting, the KRR predictor with kernel $K=\Sigma$ and ridge $\lambda=\frac{\epsilon^{2}}{N}$ minimizes the expected squared error at all points $x$.

Using Theorem 6, we obtain the following approximation of the expected risk for a general kernel $K$ and ridge $\lambda$:

For a random true function of zero mean and covariance kernel $\Sigma$ the expected risk is approximated by

This formula can be further simplified. First note that differentiating both sides of Equation 7 w.r.t. to $\lambda$, we obtain that

Secondly, differentiating both sides of Equation 7, we obtain, writing $K(\tau)=K+\tau(\Sigma-K)$

Putting everything together, we obtain that

B.8 Technical Lemmas

For any family $\boldsymbol{A}=\left(A^{(1)},\ldots,A^{(k)}\right)$ of $k$ square matrices of same size, any permutation $\sigma\in\mathfrak{S}_{k}$, we define:

where the product inside the trace is taken following the order given by the cycle and, by the cyclic property, does not depend on the starting point (see ). For example if $k=4$ and $\sigma$ is the product of transpositions $(1,3)(2,4),$

If $\boldsymbol{A}=\left(A^{(1)},\ldots,A^{(k)}\right)$ is a family of $k$ square symmetric random matrices of size $P$ independent from a standard Gaussian vector $w$ of size $P$, we have

Furthermore, if $w$ and $v$ are independent Gaussian vectors of size $P$ and independent from $\boldsymbol{A}$, then

therefore, it is sufficient to prove that

B.8.2 Bound on derivatives

Given a bound on a holomorphic function, one can obtain a bound on its derivative.

The inequality follows by considering $r=-\frac{1}{2}\Re(z)$ and using the fact that $F$ is decreasing. ∎

B.8.3 Generalized Cauchy-Schwarz inequality

Another result that we will use is the following generalization of the Cauchy-Schwarz inequality, which is a consequence of Hölder’s inequality.

For complex random variables $a_{1},...,a_{s}$, we have

The proof is done using an induction argument. The initialization, i.e. when $s=1$, is trivial.

For the induction step, assume that the result is true for $s$ terms and let us prove it for $s+1$ terms. By Hölder’s inequality applied for $p=s+1$ and $q=\frac{s+1}{s}$, we obtain:

where the second inequality is obtained by the induction hypothesis. ∎

B.8.4 Control on fixed points