# Risk Convergence of Centered Kernel Ridge Regression with Large Dimensional
Data

Khalil Elkhalil, Abla Kammoun, Xiangliang Zhang, Mohamed-Slim Alouini and
Tareq Al-Naffouri The authors are with the Electrical Engineering Program,
King Abdullah University of Science and Technology, Thuwal, Saudi Arabia;
e-mails: {khalil.elkhalil, abla.kammoun, xiangliang.zhang, tareq.alnaffouri,
<EMAIL_ADDRESS>
###### Abstract

This paper carries out a large dimensional analysis of a variation of kernel
ridge regression that we call _centered kernel ridge regression_ (CKRR), also
known in the literature as kernel ridge regression with offset. This modified
technique is obtained by accounting for the bias in the regression problem
resulting in the old kernel ridge regression but with _centered_ kernels. The
analysis is carried out under the assumption that the data is drawn from a
Gaussian distribution and heavily relies on tools from random matrix theory
(RMT). Under the regime in which the data dimension and the training size grow
infinitely large with fixed ratio and under some mild assumptions controlling
the data statistics, we show that both the empirical and the prediction risks
converge to a deterministic quantities that describe in closed form fashion
the performance of CKRR in terms of the data statistics and dimensions.
Inspired by this theoretical result, we subsequently build a consistent
estimator of the prediction risk based on the training data which allows to
optimally tune the design parameters. A key insight of the proposed analysis
is the fact that asymptotically a large class of kernels achieve the same
minimum prediction risk. This insight is validated with both synthetic and
real data.

###### Index Terms:

Kernel regression, centered kernels, random matrix theory.

##  I Introduction

Kernel ridge regression (KRR) is part of kernel-based machine learning methods
that deploy a set of nonlinear functions to describe the real output of
interest [1, 2]. More precisely, the idea is to map the data into a high-
dimensional space $\mathcal{H}$, a.k.a. _feature space_ , which can even be of
infinite dimension resulting in a linear representation of the data with
respect to the output. Then, a linear regression problem is solved in
$\mathcal{H}$ by controlling over-fitting with a regularization term. In fact,
the most important advantage of kernel methods is the utilized _kernel trick_
or _kernel substitution_ [1], which allows to directly work with kernels and
avoid explicit use of feature vectors in $\mathcal{H}$.

Due to its popularity, a rich body of research has been conducted to analyze
the performance of KRR. In [3], a randomized version of KRR is studied with
performance guarantees in terms of concentration bounds. The work in [4]
analyzes the random features approximation in least squares kernel regression.
More relevant results can be found in [5] where upper bounds of the prediction
risk have been derived in terms of the empirical quadratic risk for general
regression models. Similarly for KRR models, an upper and lower bound on the
expected risk have been provided in [6] before being generalized to general
regularization operators in [7]. Therefore, most of the results related to the
performance analysis of KRR and related regression techniques are in the form
of upper or lower bounds of the prediction risk. In this work, we study the
problem from an asymptotic analysis perspective. As we will demonstrate in the
course of the paper, such an analysis brought about novel results that predict
in an accurate fashion prediction risks metrics. Our focus is on a variation
of KRR called centered kernel ridge regression (CKRR) that is built upon the
same principles of KRR with the additional requirement to minimize the _bias_
in the learning problem. This variation has been motivated by Cortes et al. in
[8] and [9, 10] where the benefits of centering kernels have been highlighted.
The obtained regression technique can be seen as KRR with centered kernels.
Moreover, in the high dimensional setting with certain normalizations, we show
that kernel matrices behave as a rank one matrix, thus centering allows to
neutralize this non-informative component and highlight higher order
components that retain useful information of the data.

To understand the behavior of CKRR, we conduct theoretical analysis in the
large dimensional regime where both the data dimension $p$ and the training
size $n$ tend to infinity with fixed ratio ($p/n\to\text{constant}$). As far
as _inner-product kernels_ are concerned, with mild assumptions on the data
statistics, we show using fundamental results from random matrix theory
elaborated in [11] and [12] that both the empirical and prediction risks
approach a deterministic quantity that relates in closed form fashion these
performance measures to the data statistics and dimensions. This important
finding allows to see how the model performance behaves as a function of the
problem’s parameters and as such tune the design parameters to minimize the
prediction risk. Moreover, as an outcome of this result, we show that it is
possible to jointly optimize the regularization parameter along with the
kernel function so that to achieve the possible minimum prediction risk. In
other words, the minimum prediction risk is always attainable for all kernels
with a proper choice of the regualrization parameter. This implies that all
kernels behave similarly to the linear kernel. We regard such a fact as a
consequence of the curse of dimensionality phenomenon which causes the CKRR to
be asymptotically equivalent to _centered linear ridge regression_. As an
additional contribution of the present work, we build a consistent estimator
of the prediction risk based on the training samples, thereby paving the way
towards optimal setting of the regularization parameter.

The rest of the paper is structured as follows. In section II, we give a brief
background on kernel ridge regression and introduce its centered variation. In
section III, we provide the main results of the paper related to the
asymptotic analysis of CKRR as well as the construction of a consistent
estimator of the prediction risk. Then, we provide some numerical examples in
section IV. We finally make some concluding remarks in section V.

### Notations: $\mathop{\mathbb{E}}\left[.\right]$ and
$\text{var}\left[.\right]$ stand for the expectation and the variance of a
random variable while $\to_{a.s.}$ and $\to_{prob.}$ respectively stand for
the almost sure convergence and the convergence in probability.
$\left\|.\right\|$ denotes the operator norm of a matrix and the $L_{2}$ norm
for vectors, $\mathop{\rm tr}\left[.\right]$ stands for the trace operator.
The notation $f=O\left(g\right)$ means that $\exists M$ bounded such that
$f\leq Mg$. We say that $f$ is $\mathcal{C}^{p}$ if the $p$th derivative of
$f$ exists and is continous.

##  II Background on kernel ridge regression

Let $\\{\left(\bm{x}_{i},y_{i}\right)\\}_{i=1}^{n}$ be a set of $n$
observations in $\mathcal{X}\times\mathcal{Y}$, where $\mathcal{X}$ denotes
the input space and $\mathcal{Y}$ the output space. Our aim is to predict the
output of new input points $\bm{x}\in\mathcal{X}$ with a reasonable accuracy.
Assume that the output is generated using a function
$f:\mathcal{X}\to\mathcal{Y}$, then the problem can be cast as a function
approximation problem where the goal is to find an estimate of $f$ denoted by
$\widehat{f}$ such that $\widehat{f}\left(\bm{x}\right)$ is close to the real
output $f\left(\bm{x}\right)$. In this context, the kernel learning problem is
formulated as follows

$\min_{f\in\mathcal{H}}\frac{1}{2}\sum_{i=1}^{n}l\left(y_{i},f\left(\bm{x}_{i}\right)\right)+\frac{\lambda}{2}\left\|f\right\|_{\mathcal{H}}^{2},$
(1)

where $\mathcal{H}$ is a reproducing kernel Hilbert space (RKHS),
$l:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}$ is a loss function and
$\lambda>0$ is a regularization parameter that permits to control overfitting.
Denoting by $\phi:\mathcal{X}\to\mathcal{H}$ a feature map that maps the data
points to the feature space $\mathcal{H}$, then we define
$k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ such that
${k}\left(\bm{x},\bm{x}^{\prime}\right)=\left\langle\phi\left(\bm{x}\right),\phi\left(\bm{x}^{\prime}\right)\right\rangle_{\mathcal{H}}$
for all $\bm{x},\bm{x}^{\prime}\in\mathcal{X}$ where ${k}$ is known as the
positive definite kernel corresponding to the feature map $\phi$. With these
definitions, the representer theorem [13, 14] shows that the minimizer of the
problem in (1) writes as
$f^{*}\left(\bm{x}\right)=\bm{\alpha}^{T}\phi\left(\bm{x}\right)$. Thus, we
can reformulate (1) as follows

$\min_{\bm{\alpha}\in\mathcal{H}}\frac{1}{2}\sum_{i=1}^{n}l\left(y_{i},\bm{\alpha}^{T}\phi\left(\bm{x}_{i}\right)\right)+\frac{\lambda}{2}\left\|\bm{\alpha}\right\|^{2}.$
(2)

When $l$ is the squared loss, the optimization problem in (2) can be
reformulated as

$\displaystyle\min_{\bm{\alpha}\in\mathcal{H}}\frac{1}{2}\left\|\bm{y}-\mathbf{\Phi}\bm{\alpha}\right\|^{2}+\frac{\lambda}{2}\left\|\bm{\alpha}\right\|^{2},$
(3)

where
$\mathbf{\Phi}=\left[\phi\left(\bm{x}_{1}\right),\cdots,\phi\left(\bm{x}_{n}\right)\right]^{T}$.
This yields the following solution
$\bm{\alpha}^{*}=\left(\mathbf{\Phi}^{T}\mathbf{\Phi}+\lambda\mathbf{I}\right)^{-1}\mathbf{\Phi}^{T}\bm{y}$,
where $\bm{y}=\\{y_{i}\\}_{i=1}^{n}$. Then, the output estimate of any data
point $\bm{s}$ is given by [1]

$\begin{split}\widehat{f}\left(\bm{s}\right)=\bm{\kappa}\left(\bm{s}\right)^{T}\left(\mathbf{K}+\lambda\mathbf{I}\right)^{-1}\bm{y},\end{split}$
(4)

where
$\bm{\kappa}\left(\bm{s}\right)=\\{{k}\left(\bm{s},\bm{x}_{i}\right)\\}_{i=1}^{n}$
is the information vector and $\mathbf{K}=\mathbf{\Phi}\mathbf{\Phi}^{T}$ with
entries $\mathbf{K}_{i,j}={\bm{\kappa}}\left(\bm{x}_{i},\bm{x}_{j}\right)$,
$1\leq i,j\leq n$. This is commonly known as the kernel trick which allows to
highly simplify the problem which boils down to solving a $n$-dimensional
problem. Throughout this paper, we consider the following data model

$\displaystyle y_{i}=f\left(\bm{x}_{i}\right)+\sigma\epsilon_{i},\quad
i=1,\cdots,n.$ (5)

where $f$ generates the actual output of the data and $\epsilon_{i}$ are
i.i.d. standard normal random variables with $\sigma^{2}$ assumed to be known.
We consider both the empirical (training) and the prediction (testing) risks
respectively defined as [15]

$\displaystyle\mathcal{R}_{\rm
train}\left(\widehat{f}\right)=\frac{1}{n}\mathbb{E}_{\bm{\epsilon}}\|\widehat{f}\left(\mathbf{X}\right)-f\left(\mathbf{X}\right)\|_{2}^{2},$
(6) $\displaystyle\mathcal{R}_{\rm
test}\left(\widehat{f}\right)=\mathop{\mathbb{E}}_{\bm{s}\sim\mathcal{D},\bm{\epsilon}}\left[\left(\widehat{f}\left(\bm{s}\right)-f\left(\bm{s}\right)\right)^{2}\right],$
(7)

where $\mathcal{D}$ is the data input distribution, $\bm{s}$ is taken
independent of the training data $\mathbf{X}$ and
$\bm{\epsilon}=\\{\epsilon_{i}\\}_{i=1}^{n}$. The above two equations
respectively measure the goodness of fit relative to the training data and to
new unseen data all in terms of the mean squared error (MSE).

###  II-A Centered kernel ridge regression

The concept of centered kernels dates back to the work of Cortes [8] on
learning kernels based on the notion of centered alignment. As we will show
later, this notion of centering comes naturally to the picture when we account
for the bias in the learning problem (also see the lecture notes by Jakkola
[16]). More specifically, we modify the optimization problem in (2) to account
for the bias as follows

$\min_{\alpha_{0},\bm{\alpha}\in\mathcal{H}}\frac{1}{2}\sum_{i=1}^{n}l\left(y_{i},\alpha_{0}+\bm{\alpha}^{T}\phi\left(\bm{x}_{i}\right)\right)+\frac{\lambda}{2}\left\|\bm{\alpha}\right\|^{2},$
(8)

where clearly we do not penalize the _offset_ (or the _bias_) $\alpha_{0}$ in
the regularization term . With $l$ being the squared loss, we immediately get
$\alpha_{0}^{*}=\bar{y}-\frac{1}{n}\bm{\alpha}^{T}\bm{\Phi}^{T}\mathbf{1},\quad\overline{y}=\frac{1}{n}\bm{1}^{T}\bm{y}.$
Substituting $\alpha_{0}^{*}$ in (8), we solve the centered optimization
problem given by

$\min_{\bm{\alpha}\in\mathcal{H}}\frac{1}{2}\left\|\mathbf{P}\left(\bm{y}-\mathbf{\Phi}\bm{\alpha}\right)\right\|_{2}^{2}+\frac{\lambda}{2}\bm{\alpha}^{T}\bm{\alpha},$
(9)

where $\mathbf{P}=\mathbf{I}_{n}-\frac{1}{n}\bm{1}_{n}\bm{1}_{n}^{T}$ is
referred as a _projection_ matrix or a _centering_ matrix [8, 16]. Finally, we
get

$\displaystyle\bm{\alpha}^{*}$
$\displaystyle=\left(\mathbf{\Phi}^{T}\mathbf{P}\mathbf{\Phi}+\lambda\mathbf{I}\right)^{-1}\mathbf{\Phi}^{T}\mathbf{P}\left(\bm{y}-\bar{y}\mathbf{1}\right)$
$\displaystyle\overset{(b)}{=}\mathbf{\Phi}^{T}\mathbf{P}\left(\mathbf{P}\mathbf{\Phi}\mathbf{\Phi}^{T}\mathbf{P}+\lambda\mathbf{I}_{n}\right)^{-1}\left(\bm{y}-\bar{y}\mathbf{1}_{n}\right)$
$\displaystyle=\mathbf{\Phi}^{T}\mathbf{P}\left(\mathbf{P}\mathbf{K}\mathbf{P}+\lambda\mathbf{I}_{n}\right)^{-1}\left(\bm{y}-\bar{y}\mathbf{1}_{n}\right)$
$\displaystyle=\mathbf{\Phi}^{T}\mathbf{P}\left(\mathbf{K}_{c}+\lambda\mathbf{I}_{n}\right)^{-1}\left(\bm{y}-\bar{y}\mathbf{1}_{n}\right),$

where $\mathbf{K}_{c}=\mathbf{P}\mathbf{K}\mathbf{P}$ is the _centered_ kernel
matrix as defined in [8, Lemma1] and $(b)$ is obtained using the Woodbury
identity. With some basic manipulations, the centered kernel ridge regression
estimate of the output of data point $\bm{s}$ is given by

$\widehat{f}_{c}\left(\bm{s}\right)=\bm{\kappa}_{c}\left(\bm{s}\right)^{T}\left(\mathbf{K}_{c}+\lambda\mathbf{I}\right)^{-1}\mathbf{P}\bm{y}+\bar{y}.$
(10)

Therefore, the feature map corresponding to $\mathbf{K}_{c}$ as well as the
information vector can be respectively obtained as follows

$\displaystyle\phi_{c}\left(\bm{s}\right)=\phi\left(\bm{s}\right)-\frac{1}{n}\sum_{i=1}^{n}\phi\left(\bm{x}_{i}\right),\>\bm{\kappa}_{c}\left(\bm{s}\right)=\mathbf{P}\bm{\kappa}\left(\bm{s}\right)-\frac{1}{n}\mathbf{P}\mathbf{K}\mathbf{1}_{n}.$
(11)

Throughout this paper, we consider inner-product kernels [1, 11] defined as
follows

$\displaystyle{k}\left(\bm{x},\bm{x}^{\prime}\right)=g\left(\bm{x}^{T}\bm{x}^{\prime}/p\right),\quad\forall\bm{x},\bm{x}^{\prime}\in\mathbb{R}^{p},$
(12)

and subsequently,
$\mathbf{K}=\left\\{g\left(\bm{x}_{i}^{T}\bm{x}_{j}/p\right)\right\\}_{i,j=1}^{n},$
where the normalization111This is equivalent to normalize all data points by
$\sqrt{p}$. This type of normalization has also been conisdered in [17]
following the heuristic of Jakkola. by $p$ in (12) is convenient in the large
$n,p$ regime as we will show later (also see [18] for similar normalization in
the analysis of LS-SVMs). In the following, we conduct a large dimensional
analysis of the performance of CKRR with the aim to get useful insights on the
design of CKRR. Particularly, we will focus on studying the empirical and the
prediction risks of CKRR which we define as

$\displaystyle\mathcal{R}_{\rm
train}=\frac{1}{n}\mathbb{E}_{\bm{\epsilon}}\sum_{i=1}^{n}\left|\widehat{f}_{c}({\bm{x}}_{i})-f({\bm{x}}_{i})\right|^{2}\hskip
14.22636pt,$ $\displaystyle\mathcal{R}_{\rm
test}=\mathbb{E}_{\bm{s}\sim\mathcal{D},{\bm{\epsilon}}}\left[\left|\widehat{f}_{c}({\bm{s}})-f({\bm{s}})\right|^{2}\right].$

The novelty of our analysis with respect to previous studies lies in that

  1. 1.

It provides a mathematical connection between the performance and the
problem’s dimensions and statistics resulting in a deeper understanding of
centered kernel ridge regression in the large $n,p$ regime.

  2. 2.

It brings insights on how to choose the kernel function $g$ and the
regularization parameter $\lambda$ in order to guarantee a good generalization
performance for unknown data.

As far as the second point is considered, we show later that both the kernel
function and the regularization parameter can be optimized jointly as a
consequence of the mathematical result connecting the prediction risk with
these design parameters. Our analysis does not assume a specific choice of the
inner-product kernels, and is valid for the following popular ones.

  * •

Linear kernels:
$k\left(\bm{x},\bm{x}^{\prime}\right)=\alpha\bm{x}^{T}\bm{x}^{\prime}/p+\beta$.

  * •

Polynomial kernels:
$k\left(\bm{x},\bm{x}^{\prime}\right)=\left(\alpha\bm{x}^{T}\bm{x}^{\prime}/p+\beta\right)^{d}$.

  * •

Sigmoid kernels:
$k\left(\bm{x},\bm{x}^{\prime}\right)=\tanh\left(\alpha\bm{x}^{T}\bm{x}^{\prime}/p+\beta\right)$.

  * •

Exponential kernels:
$k\left(\bm{x},\bm{x}^{\prime}\right)=\exp\left(\alpha\bm{x}^{T}\bm{x}^{\prime}/p+\beta\right)$.

##  III Main results

###  III-A Technical assumptions

In this section, we will present our theoretical results on the prediction
risk of CKRR by first introducing the assumptions of data growth rate, kernel
function $g$ and true function $f$. Without loss of generality, we assume that
the data samples $\bm{x}_{1},\cdots,\bm{x}_{n}\in\mathbb{R}^{p}$ are
independent such that
$\bm{x}_{i}\sim\mathcal{N}\left(\bm{0}_{p},\mathbf{\Sigma}\right),i=1,\cdots,n$,
with positive definite covariance matrix
$\mathbf{\Sigma}\in\mathbb{R}^{p\times p}$. Throughout the analysis, we
consider the large dimensional regime in which both $p$ and $n$ grow
simultaneously large with the following growth rate assumptions.

######  Assumption 1 (Growth rate).

As $p,n\to\infty$ we assume the following

  * •

Data scaling: $p/n\to c\in\left(0,\infty\right)$.

  * •

Covariance scaling: $\lim\sup_{p}\left\|\mathbf{\Sigma}\right\|<\infty$.

The above assumptions are standard to consider and allow to exploit the large
heritage of random matrix theory. Moreover, allowing $p$ and $n$ to grow large
at the same rate is of practical interest when dealing with modern large and
numerous data. The assumption treating the covariance scaling is technically
convenient since it allows to use important theoretical results on the
behavior of large kernel matrices [11, 12]. Under Assumption 1, we have the
following implications.

$\displaystyle\bm{x}_{i}^{T}\bm{x}_{i}/p$
$\displaystyle\to_{a.s.}\frac{1}{p}\mathop{\rm
tr}\mathbf{\Sigma}\triangleq\tau,\>i=1,\cdots,n.$ (13)
$\displaystyle\bm{x}_{i}^{T}\bm{x}_{j}/p$ $\displaystyle\to_{a.s.}0,\>i\neq
j,$ (14)

where $0<\tau<\infty$ due to the covariance scaling in Assumption 1. This
means that in the limit when $p\to\infty$, the kernel matrix $\mathbf{K}$ as
defined earlier has all its entries converging to a deterministic limit.
Applying a Taylor expansion on the entries of ${\bf K}$, and under some
assumption on the kernel function $g$, it has been shown in [11, Theorem 2.1]
that

$\displaystyle\left\|\mathbf{K}-\mathbf{K}^{\infty}\right\|\to_{prob.}0,$ (15)

where the convergence is in operator norm and $\mathbf{K}^{\infty}$ exhibits
nice properties and can be expressed using standard random matrix models. The
explicit expression of $\mathbf{K}^{\infty}$ as well as its properties will be
thoroughly investigated in Appendix A. We subsequently make additional
assumptions to control the kernel function $g$ and the data generating
function $f$.

######  Assumption 2 (Kernel function).

As in [11, Theorem 2.1], we shall assume that $g$ is $\mathcal{C}^{1}$ in a
neighborhood of $\tau$ and $\mathcal{C}^{3}$ in a neighborhood of 0. Moreover,
we assume that for any independent observations $\bm{x}_{i}$ and $\bm{x}_{j}$
drawn from $\mathcal{N}({\bf 0}_{p},{\bm{\Sigma}})$ and $k\in\mathbb{N}^{*}$,

$\mathbb{E}\left|g^{(3)}\left(\frac{1}{p}\bm{x}_{i}^{T}\bm{x}_{j}\right)\right|^{k}<\infty$

where $g^{(3)}$ is the third derivative of $g$.

######  Assumption 3 (Data generating function).

We assume that $f$ is $\mathcal{C}^{1}$ and polynomially bounded together with
its derivatives. We shall further assume that the moments of $f(\bm{x})$ and
its gradient are finite. More explicitly we need to have:

$\mathbb{E}_{\bm{x}\sim\mathcal{N}({\bf
0},{\bm{\Sigma}})}\left|f(\bm{x})\right|^{k}<\infty,$

and

$\displaystyle\mathop{\mathbb{E}}_{\bm{x}\sim\mathcal{N}({\bf
0},{\bm{\Sigma}})}\left\|\bm{\nabla}_{f}\left(\bm{x}\right)\right\|_{2}^{k}<\infty,\text{where}\>\bm{\nabla}_{f}\left(\bm{x}\right)=\left\\{\frac{\partial
f\left(\bm{x}\right)}{\partial x_{l}}\right\\}_{l=1}^{p}.$ (16)

As we will show later, the above assumptions are needed to guarantee a bounded
asymptotic risk and to carry out the analysis. Under the setting of
Assumptions 1, 2 and 3, we aim to study the performance of CKRR by
asymptotically evaluating the performance metrics defined in (6). Inspired by
the fundamental results from [11] and [12] in the context of spectral
clustering, then following the observations made in (13) and(14), it is always
possible to linearize the kernel matrix $\mathbf{K}$ around the matrix
$g\left(0\right)\bm{11}^{T}$ which avoids dealing with the original
intractable expression of $\mathbf{K}$. Note that the first component of the
approximation given by $g\left(0\right)\bm{11}^{T}$ will be neutralized by the
projection matrix $\mathbf{P}$ in the context of CKRR, which means that the
behavior of CKRR will be essentially governed by the higher order
approximations of $\mathbf{K}$. Consequently, one can resort to those
approximations to have an explicit expression of the asymptotic risk in the
large $p,n$ regime. This expression would hopefully reveal the mathematical
connection between the regression risk and the data’ statistics and dimensions
as $p,n\to\infty$.

###  III-B Limiting risk

With the above assumptions at hand, we are now in a position to state the main
results of the paper related to the derivation of the asymptotic risk of CKRR.
Before doing so, we shall introduce some useful quantities.

$\displaystyle\nu$ $\displaystyle\triangleq
g\left(\tau\right)-g\left(0\right)-\tau g^{\prime}\left(0\right).$

Also, for all $z\in\mathbb{C}$ at macroscopic distance from the eigenvalues
$\lambda_{1},\cdots,\lambda_{p}$ of
$\frac{1}{p}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{T}$, we define the _Stieltjes
transform_ of $\frac{1}{p}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{T}$ also known
as the _Stieltjes transform_ of the Marc̆enko-Pastur law as the unique
solution to the following fixed-point equation [19]

$\displaystyle{m}\left(z\right)$
$\displaystyle\triangleq-\left[cz-\frac{1}{n}\mathop{\rm
tr}\mathbf{\Sigma}\left(\mathbf{I}+{m}\left(z\right)\mathbf{\Sigma}\right)^{-1}\right]^{-1},$
(17)

where ${m}\left(z\right)$ in (17) is bounded as $p\to\infty$ provided that
Assumption 1 is satisfied. For ease of notation, we shall use ${m}_{z}$ to
denote ${m}\left(z\right)$ for all appropriate $z$. The first main result of
the paper is summarized in the following theorem, the proof of which is
postponed to the Appendix A.

######  Theorem 1 (Limiting risk).

Under Assumptions 1, 2 and 3 and by taking
$z=-\frac{\lambda+\nu}{g^{\prime}\left(0\right)}$ for kernel functions
satisfying $g^{\prime}\left(0\right)\neq 0$ 222The case of
$g^{\prime}\left(0\right)=0$ is asymptotically equivalent to take the sample
mean as an estimate of $f$ which is neither of practical nor theoretical
interest. and $z$ at macroscopic distance from the eigenvalues of
$\frac{1}{p}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{T}$, both the empirical and
the prediction risks converge in probability to a non trivial deterministic
limits respectively given by

$\displaystyle\mathcal{R}_{train}-\mathcal{R}_{train}^{\infty}\to_{prob.}0,$
(18)
$\displaystyle\mathcal{R}_{test}-\mathcal{R}_{test}^{\infty}\to_{prob.}0,$
(19)

where the expressions of $\mathcal{R}_{train}^{\infty}$ and
$\mathcal{R}_{test}^{\infty}$ are given in the top of the next page.

$\begin{split}\mathcal{R}_{train}^{\infty}&=\left(\frac{c\lambda{m}_{z}}{g^{\prime}\left(0\right)}\right)^{2}\frac{n\sigma^{2}+n\textbf{var}_{f}-n{m}_{z}\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]^{T}\mathbf{\Sigma}\left[\left(\mathbf{I}+{m}_{z}\mathbf{\Sigma}\right)^{-1}+\left(\mathbf{I}+{m}_{z}\mathbf{\Sigma}\right)^{-2}\right]\mathbf{\Sigma}\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]}{n-{m}^{2}_{z}\mathop{\rm
tr}\bm{\Sigma}^{2}\left(\mathbf{I}+{m}_{z}\bm{\Sigma}\right)^{-2}}+\sigma^{2}-2\sigma^{2}\frac{c\lambda{m}_{z}}{g^{\prime}\left(0\right)}.\end{split}$
(20)

$\displaystyle\mathcal{R}^{\infty}_{test}$
$\displaystyle=\frac{n\sigma^{2}+n\textbf{var}_{f}-n{m}_{z}\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]^{T}\mathbf{\Sigma}\left[\left(\mathbf{I}+{m}_{z}\mathbf{\Sigma}\right)^{-1}+\left(\mathbf{I}+{m}_{z}\mathbf{\Sigma}\right)^{-2}\right]\mathbf{\Sigma}\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]}{n-{m}^{2}_{z}\mathop{\rm
tr}\bm{\Sigma}^{2}\left(\mathbf{I}+{m}_{z}\bm{\Sigma}\right)^{-2}}-\sigma^{2}.$
(21)

Note that in the case where $\mathbf{\Sigma}=\mathbf{I}_{p}$, the limiting
risks in (20) and (21) can be further simplified as

$\displaystyle\mathcal{R}_{train}^{\infty}$
$\displaystyle=\left(\frac{c\lambda{m}_{z}}{g^{\prime}\left(0\right)}\right)^{2}$
$\displaystyle\times\frac{n\left(1+{m}_{z}\right)^{2}\left(\sigma^{2}+\textbf{var}_{f}\right)-n{m}_{z}\left(2+{m}_{z}\right)\left\|\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]\right\|^{2}}{n\left(1+{m}_{z}\right)^{2}-p{m}^{2}_{z}}$
$\displaystyle+\sigma^{2}-2\sigma^{2}\frac{c\lambda{m}_{z}}{g^{\prime}\left(0\right)}.$
$\displaystyle\mathcal{R}^{\infty}_{test}$
$\displaystyle=\frac{n\left(1+{m}_{z}\right)^{2}\left(\sigma^{2}+\textbf{var}_{f}\right)-n{m}_{z}\left(2+{m}_{z}\right)\left\|\mathop{\mathbb{E}}\left[\bm{\nabla}_{f}\right]\right\|^{2}}{n\left(1+{m}_{z}\right)^{2}-p{m}^{2}_{z}}$
$\displaystyle-\sigma^{2},$

where ${m}_{z}$ can be explicitly derived as in [20]

${m}_{z}=\frac{-\left(cz-c+1\right)-\sqrt{\left(cz-c-1\right)^{2}-4c}}{2cz}.$

######  Remark 1.

From Theorem 1 it entails that the limiting prediction risk can be expressed
using the limiting empirical risk in the following fashion.

$\displaystyle\mathcal{R}^{\infty}_{test}=\left(\frac{c\lambda{m}_{z}}{g^{\prime}\left(0\right)}\right)^{-2}\mathcal{R}^{\infty}_{train}-\sigma^{2}\left(\frac{g^{\prime}\left(0\right)}{c\lambda{m}_{z}}-1\right)^{2}.$
(22)

######  Lemma 1 (A consistent estimator of the prediction risk).

Inspired by the outcome of Theorem 1 summarized in Remark 1, we construct a
consistent estimator of the prediction risk given by

$\displaystyle\widehat{\mathcal{R}}_{test}=\left(\frac{c\lambda\widehat{{m}}_{z}}{g^{\prime}\left(0\right)}\right)^{-2}\mathcal{R}_{train}-\sigma^{2}\left(\frac{g^{\prime}\left(0\right)}{c\lambda\widehat{{m}}_{z}}-1\right)^{2},\>\text{with}\>\lambda>0,$
(23)

in the sense that
$\mathcal{R}_{test}-\widehat{\mathcal{R}}_{test}\to_{prob.}0,$ where ${m}_{z}$
can be consistently estimated as $\widehat{{m}}_{z}=\frac{1}{p}\mathop{\rm
tr}\left(\mathbf{XX}^{T}/p-z\mathbf{I}_{n}\right)^{-1},\>\mathbf{X}=\left[\bm{x}_{1},\cdots,\bm{x}_{n}\right]^{T}.$

###### Proof.

The proof is straightforward relying on the relation in (22) and the fact that

$\displaystyle\widehat{{m}}_{z}-{m}_{z}\to_{a.s.}0,$

as shown in [12, Lemma 1]. ∎

Since the aim of any learning system is to design a model that achieves
minimal prediction risk [15], the relation described in Lemma 1 by (23) has
enormous advantages as it permits to estimate the prediction risk in terms of
the empirical risk and hence optimize the prediction risk accordingly.

######  Remark 2.

One important observation from the expression of the limiting prediction risk
in (21) is that the information on the kernel (given by
$g^{\prime}\left(0\right)$ and $\nu$) as well as the information on $\lambda$
are both encapsulated in ${m}_{z}$ with
$z=-\frac{\lambda+\nu}{g^{\prime}\left(0\right)}$. This means that one should
optimize $z$ to have minimal prediction risk and thus jointly choose the
kernel $g$ and the regularization parameter $\lambda$. Moreover, it entails
that the choice of the kernel (as long as $g^{\prime}\left(0\right)\neq 0$) is
asymptotically irrelevant since a bad choice of the kernel can be compensated
by a good choice of $\lambda$ and vice-versa. This essentially implies that a
linear kernel asymptotically achieves the same optimal performance as any
other type of kernels 333This does not mean that all kernels will have the
same performance for a given regularization parameter but means that they will
achieve the same minimum prediction risk..

###  III-C A consistent estimator of the prediction risk

Although the estimator provided in Lemma 1 permits to estimate the prediction
risk by virtue of the empirical risk, it presents the drawback of being
sensitive to small values of $\lambda$. In the following theorem, we provide a
consistent estimator of the prediction risk constructed from the training data
$\left\\{\left(\bm{x}_{i},y_{i}\right)\right\\}_{i=1}^{n}$ and is less
sensitive to small values of $\lambda$.

######  Theorem 2 (A consistent estimator of the prediction risk).

Under Assumptions 1, 2 and 3 with $g^{\prime}\left(0\right)\neq 0$ and
$z=-\frac{\lambda+\nu}{g^{\prime}\left(0\right)}$, we construct a consistent
estimator of the prediction risk based on the training data such that

$\mathcal{R}_{test}-\widehat{\mathcal{R}}_{test}\to_{prob.}0,$
$\begin{split}\widehat{\mathcal{R}}_{test}&=\frac{1}{\left(cz\widehat{{m}}_{z}\right)^{2}}\left[\frac{1}{np}\bm{y}^{T}\mathbf{PX}\left(z\widetilde{\mathbf{Q}}_{z}^{2}-\widetilde{\mathbf{Q}}_{z}\right)\mathbf{X}^{T}\mathbf{P}\bm{y}+\textbf{var}\left(\bm{y}\right)\right]\\\
&-\sigma^{2},\end{split}$ (24)

with $\widetilde{\mathbf{Q}}_{z}$ is the resolvent matrix given by
$\widetilde{\mathbf{Q}}_{z}=\left(\frac{\mathbf{X}^{T}\mathbf{PX}}{p}-z\mathbf{I}_{p}\right)^{-1}$,
with ${\bf X}=\left[\bm{x}_{1},\cdots,\bm{x}_{n}\right]^{T}$. Moreover, in the
special case where $\mathbf{\Sigma}=\mathbf{I}_{p}$, the estimator reduces to

$\begin{split}\widehat{\mathcal{R}}_{test}&=\frac{n\left(1+\widehat{{m}}_{z}\right)^{2}\textbf{var}\left(\bm{y}\right)}{{n\left(1+\widehat{{m}}_{z}\right)^{2}-p\widehat{{m}}^{2}_{z}}}\\\
&-\frac{\widehat{{m}}_{z}\left(2+\widehat{{m}}_{z}\right)\left[\bm{y}^{T}\mathbf{P}\frac{\mathbf{X}\mathbf{X}^{T}}{n}\mathbf{P}\bm{y}-p\textbf{var}\left(\bm{y}\right)\right]}{n\left(1+\widehat{{m}}_{z}\right)^{2}-p\widehat{{m}}^{2}_{z}}-\sigma^{2}.\end{split}$
(25)

Theorem 2 provides a generic way to estimate the prediction risk from the
pairs of training examples
$\left\\{\left(\bm{x}_{i},y_{i}\right)\right\\}_{i=1}^{n}$. This allows using
the closed form expressions in (24) and (25)444The expression in (25) is
useful because it does not involve any matrix inversion unlike the one in
(24). with the same set of arguments in Remark 2 to jointly estimate the
optimal kernel and the optimal regularization parameter $\lambda$.

###  III-D Parameters optimization

We briefly discuss how to jointly optimize the kernel function and the
regularization parameter $\lambda$. As mentioned earlier, we exploit the
special structure in the expression of the consistent estimate
$\widehat{\mathcal{R}}_{test}$ where both parameters (the kernel function $g$
and $\lambda$) are summarized in $z$. We focus on the case where
$\mathbf{\Sigma}=\mathbf{I}_{p}$ due to the tractability of the expression of
$\widehat{\mathcal{R}}_{test}$ in (25). By simple calculations, we can show
that $\widehat{\mathcal{R}}_{test}$ is minimized when $\widehat{m}_{z}$
satisifies the equation

$\displaystyle
p\text{var}\left(\bm{y}\right)\left[p\widehat{m}_{z}^{2}+n\left(1+\widehat{m}_{z}\right)^{2}\right]$
$\displaystyle=A\left(n+n\widehat{m}_{z}+p\widehat{m}_{z}^{2}\right),$

where
$A=\bm{y}^{T}\mathbf{P}\frac{\mathbf{X}\mathbf{X}^{T}}{n}\mathbf{P}\bm{y}$,
which admits the following closed-form solution

$\begin{split}m^{\star}&=\frac{\sqrt{nA^{2}-4npA^{2}+8np^{2}A\text{var}\left(\bm{y}\right)-4np^{3}\text{var}^{2}\left(\bm{y}\right)}}{2\left(-pA+np\text{var}\left(\bm{y}\right)+p^{2}\text{var}\left(\bm{y}\right)\right)}\\\
&-\frac{nA-2pn\text{var}\left(\bm{y}\right)}{\left(-pA+np\text{var}\left(\bm{y}\right)+p^{2}\text{var}\left(\bm{y}\right)\right)}.\end{split}$
(26)

Then, look up $z^{\star}$ such that $\widehat{m}_{z^{\star}}=m^{\star}$.
Finally, choose $\lambda$ and $g\left(.\right)$ such that
$z^{\star}=-\frac{\lambda+\nu}{g^{\prime}\left(0\right)}$. In the general
case, it is difficult to get a closed from expression in terms of $z$ or
$\widehat{m}_{z}$, however it is possible to numerically optimize the
expression of $\widehat{\mathcal{R}}_{test}$ with respect to $z$. This can be
done using simple one dimensional optimization techniques implemented in most
softwares.

##  IV Experiments

###  IV-A Synthetic data

To validate our theoretical findings, we consider both Gaussian and Bernoulli
data. As shown in Figure 2, both data distributions exhibit the same behavior
for all the settings with different kernel functions. More importantly,
eventhough the derived formulas heavily rely on the Gaussian assumption, in
the case where the data is Bernoulli distributed, we have a good agreement
with the theoretical limits. This can be understood as part of the
universality property often encountered in many high dimensional settings.
Therefore, we conjecture that the obtained results are valid for any data
distribution following the model
$\bm{x}\sim\mathbf{\Sigma}^{\frac{1}{2}}\bm{z},$ where $\mathbf{\Sigma}$
satisfies Assumption 1 and $\\{\bm{z}_{i}\\}_{1\leq i\leq p}$ the entries of
$\bm{z}$ are i.i.d. with zero mean, unit variance and have bounded
moments555We couldn’t provide experiments for more data distributions due to
space limitations.. For more clarity, we refer the reader to Figure 1 as a
representative of Figure 2 when the data is Gaussian with $p=100$ and $n=200$.
As shown in Figure 1, the proposed consistent estimators are able to track the
real behavior of the prediction risk for all types of kernels into
consideration. It is worth mentioning however that the proposed estimator in
Lemma 1 exhibits some instability for small values of $\lambda$ due to the
inversion of $\lambda$ in (23). Therefore, it is advised to use the estimator
given by Theorem 2. It is also clear from Figure 1 that all the considered
kernels achieve the same minimum prediction risks but with different optimal
regularizations $\lambda$. This is not the case for the empirical risk as
shown in Figure 1 and (20) where the information on the kernel and the
regularization parameter $\lambda$ are decoupled. Hence, in contrast to the
prediction risk, the regularization parameter and the kernel can not be
jointly optimized to minimize the empirical risk.

###  IV-B Real data

As a further experiment, we validate the accuracy of our result over a real
data set. To this end, we use the real _Communities and Crime Data Set_ for
evaluation [21], which has 123 samples and 122 features. For the experiment in
Figure 3, we divide the data set to have $60\%$ training samples ($n=73$) and
the remaining for testing ($n_{test}=50$). The risks in Figure 3 are obtained
by averaging the prediction risk (computed using $n_{test}$) over $500$ random
permutation of the data. Although the data set is far from being Gaussian, we
notice that the proposed prediction risk estimators are able to track the real
behavior of the prediction risk for all types of considered kernels. We can
also validate the previous insight from Theorem 1 where all kernels almost
achieve the same minimum prediction risk.

Figure 1: CKRR risk with respect to the regularization parameter $\lambda$ on
Gaussian data
($\bm{x}\sim\mathcal{N}\left(0_{p},\\{0.4^{|i-j|}\\}_{i,j}\right)$, $n=200$
training samples with $p=100$ predictors, $\sigma^{2}=0.5$ for different types
of kernels. The data generating function is taken to be
$f\left(\bm{x}\right)=\sin\left(\bm{1}^{T}\bm{x}/\sqrt{p}\right)$. Figure 2:
CKRR risk with respect to the regularization parameter $\lambda$ on both
Gaussian and Bernoulli data (i.e.,
$\bm{x}\sim\mathcal{N}\left(0_{p},\\{0.4^{|i-j|}\\}_{i,j}\right)$ and
$\bm{x}=\\{0.4^{|i-j|}\\}_{i,j}^{\frac{1}{2}}\bm{z}$ with
$\bm{z}_{i}\sim_{i.i.d.}1-2\times\text{Bernoulli}\left(1/2\right)$
respectively). The noise variance is taken to be $\sigma^{2}=0.5$ and the data
generating function is
$f\left(\bm{x}\right)=\sin\left(\bm{1}^{T}\bm{x}/\sqrt{p}\right)$. Figure 3:
CKRR risk with respect to the regularization parameter $\lambda$ on the
_Communities and Crime Data Set_ where independent zero mean Gaussian noise
samples with variance $\sigma^{2}=0.05$ are added to the true response.

##  V Concluding Remarks

We conducted a large dimensional analysis of centered kernel ridge regression,
which is a modified version of kernel ridge regression that accounts for the
bias in the regression formulation. By allowing both the data dimension and
the training size to grow infinitely large in a fixed proportion and by
relying on fundamental tools from random matrix theory, we showed that both
the empirical and the prediction risks converge in probability to a
deterministic quantity that mathematically connects these performance metrics
to the data dimensions and statistics. A fundamental insight taken from the
analysis is that asymptotically the choice of the kernel is irrelevant to the
learning problem which asserts that a large class of kernels will achieve the
same best performance in terms of prediction risk as a linear kernel. Finally,
based on the asymptotic analysis, we built a consistent estimator of the
prediction risk making it possible to estimate the optimal regularization
parameter that achieves the minimum prediction risk.

## Appendix A

## Proof of Theorem 1

Here, we provide the details of the derivation for the prediction risk. The
analysis of the empirical risk follows in a very similar way and is thus
omitted. Before delving into the proof of Theorem 1, we shall introduce some
fundamental results on the asymptotic behavior of inner-product kernel
matrices established by El-Karoui [11, Theorem 2.1].

######  Theorem (Asymptotic behavior of inner product kernel random matrices).

Under the assumptions of Theorem 2.1 [11], the kernel matrix $\mathbf{K}$ can
be approximated by $\mathbf{K}^{\infty}$ in the sense that
$\|\mathbf{K}-\mathbf{K}^{\infty}\|\to 0$ almost surely in operator norm,
where

$\displaystyle\mathbf{K}^{\infty}=\left[g\left(0\right)+g^{\prime\prime}\left(0\right)\frac{\mathop{\rm
tr}\left(\mathbf{\Sigma}^{2}\right)}{2p^{2}}\right]\mathbf{1}\mathbf{1}^{T}+g^{\prime}\left(0\right)\frac{\mathbf{XX}^{T}}{p}+\nu\mathbf{I}_{n},$

where $\nu=g\left(\tau\right)-g\left(0\right)-\tau g^{\prime}\left(0\right)$.
A similar result can be found in [11] where the accuracy of ${\bf K}^{\infty}$
has been assessed as ${\bf K}={\bf K}^{\infty}+O_{\|.\|}(\frac{1}{\sqrt{n}}),$
where $O_{\|.\|}(\frac{1}{\sqrt{n}})$ denotes a matrix with spectral norm
converging in probability to zero with a rate $1/\sqrt{n}$.

Define

$\displaystyle\mathbf{Q}_{z}=\left(\mathbf{P}\frac{\mathbf{XX}^{T}}{p}\mathbf{P}-z\mathbf{I}_{n}\right)^{-1},\quad\widetilde{\mathbf{Q}}_{z}=\left(\frac{\mathbf{X}^{T}\mathbf{P}\mathbf{X}}{p}-z\mathbf{I}_{p}\right)^{-1}.$
(27)

Note that using the Woodbury identity, it is easy to show the following useful
relations

$\displaystyle\mathbf{Q}_{z}$
$\displaystyle=-\frac{1}{z}\mathbf{I}_{n}+\frac{1}{zp}{\bf
P}\mathbf{X}\widetilde{\mathbf{Q}}_{z}\mathbf{X}^{T}{\bf P},$ (28)
$\displaystyle\widetilde{\mathbf{Q}}_{z}$
$\displaystyle=-\frac{1}{z}\mathbf{I}_{p}+\frac{1}{zp}\mathbf{X}^{T}{\bf
P}{\bf Q}_{z}{\bf PX}.$ (29)

The above theorem has the following consequence

$\begin{split}&\left\|\left(\mathbf{K}_{c}+\lambda\mathbf{I}\right)^{-1}-\frac{1}{g^{\prime}\left(0\right)}\left[\mathbf{Q}_{z}+\frac{\frac{\nu}{g^{\prime}\left(0\right)}\mathbf{Q}_{z}\frac{1}{n}\bm{11}^{T}\mathbf{Q}_{z}}{1-\frac{\nu}{g^{\prime}\left(0\right)}\frac{1}{n}\bm{1}^{T}\mathbf{Q}_{z}\bm{1}}\right]\right\|\\\
&=O_{p}(\frac{1}{\sqrt{n}}),\end{split}$ (30)

where (LABEL:K_infty) is obtained by a simple application of the Sherman-
Morrison Lemma (inversion Lemma), along with the use of the resolvent identity
${\bf A}^{-1}-{\bf B}^{-1}={\bf A}^{-1}({\bf B}-{\bf A}){\bf B}^{-1}$, which
holds for any square invertible matrices ${\bf A}$ and ${\bf B}$. The proof of
the above Theorem follows from the application of a Taylor expansion of the
elements of $\frac{1}{p}{\bf XX}^{T}$ at the vicinity of their mean. Applying
the same approach for vector ${\bm{\kappa}}\left(\bm{s}\right)$, we get

${\bm{\kappa}}({\bf s})=g(0){\bf 1}+g^{\prime}(0)\frac{1}{p}{\bf X}{\bf
s}+\tilde{\bm{\kappa}}({\bf s}),$

where $\tilde{\bm{\kappa}}({\bf s})$ has elements

$\left[\tilde{\bm{\kappa}}({\bf
s})\right]_{i}=\frac{g^{\prime\prime}(0)}{2}\left(\frac{{\bf
x}_{i}^{T}s}{p}\right)^{2}+\frac{g^{3}(\xi_{i})}{6}\left(\frac{{\bf
x}_{i}^{T}{\bf s}}{p}\right)^{3},$ (31)

with $\xi_{i}\in\left[0,\frac{{\bf x}_{i}^{T}{\bf s}}{p}\right]$. Then, since
$\mathbb{E}\left|\frac{{\bf x}_{i}^{T}{\bf s}}{\sqrt{p}}\right|^{r}$ is
uniformly bounded in $p$ for all $r\in\mathbb{N}$, we have for all
$k\in\mathbb{N}$,

$\mathbb{E}\left\|\tilde{\bm{\kappa}}(\bf
s)\right\|_{2}^{k}=O(\frac{1}{n^{\frac{k}{2}}}).$ (32)

As shall be seen later, we need also to control
$\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}(\bm{s})\tilde{\bm{\kappa}}(\bm{s})^{T}\right]$.
This is performed in the following technical Lemma.

######  Lemma 2.

Let $\tilde{\bm{\kappa}}({\bm{s}})$ be as in (32). Then,

$\displaystyle\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}(\bm{s})\tilde{\bm{\kappa}}(\bm{s})^{T}\right]$
$\displaystyle=\frac{1}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\left(\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\right)^{2}{\bf
11}^{T}+O_{\|.\|}(\frac{1}{n^{\frac{3}{2}}}).$

Similarly, the following approximations hold true

$\displaystyle\mathbb{E}_{\bm{s}}\left[{\bf
X}{\bm{s}}\tilde{\bm{\kappa}}({\bm{s}})^{T}\right]=O_{\|.\|}(n^{-\frac{3}{2}}).$
$\displaystyle\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}({\bm{s}})\frac{{\bf
1}^{T}}{n}{\bf KP}\right]=\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}\frac{1}{n}\frac{{\bf 11}^{T}}{p}\frac{{\bf XX}^{T}}{p}{\bf
P}+O_{\|.\|}(n^{-\frac{3}{2}}).$

###### Proof.

To begin with, note that for ${\bf M}=\left\\{m_{ij}\right\\}_{i,j=1}^{n}$ a
random matrix whose elements satisfies, $m_{i,j}=O_{p}(n^{-\alpha})$ for some
$\alpha>0$, as $\|{\bf M}\|^{2}\leq\mathop{\rm tr}{\bf M}{\bf M}^{T}$, we have
$m_{i,j}=O_{p}(n^{-\alpha})\Rightarrow\|{\bf M}\|=O_{p}(n^{1-\alpha})$. We
first start by deriving
$\mathbb{E}_{\bm{s}}\left[\left[\tilde{\bm{\kappa}}({\bf
s})\right]_{i}\left[\tilde{\bm{\kappa}}({\bf s})\right]_{j}\right]$. We have

$\displaystyle\mathbb{E}_{\bm{s}}\left[\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{i}\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{j}\right]=\left(\frac{g^{\prime\prime}(0)}{2p^{2}}\right)^{2}\mathbb{E}_{\bm{s}}\left[{\bm{s}}^{T}{\bf
x}_{i}{\bf x}_{i}^{T}{\bm{s}}{\bm{s}}^{T}{\bf x}_{j}{\bf
x}_{j}^{T}{\bm{s}}\right]$
$\displaystyle+\frac{g^{\prime\prime}(0)}{2}\mathbb{E}_{\bm{s}}\left[\left(\frac{{\bf
x}_{i}^{T}{\bm{s}}}{p}\right)^{2}\left(\frac{{\bf
x}_{j}^{T}{\bm{s}}}{p}\right)^{3}\frac{g^{3}(\xi_{i})}{6}\right]$
$\displaystyle+\frac{g^{\prime\prime}(0)}{2}\mathbb{E}_{\bm{s}}\left[\left(\frac{{\bf
x}_{j}^{T}{\bm{s}}}{p}\right)^{2}\left(\frac{{\bf
x}_{i}^{T}{\bm{s}}}{p}\right)^{3}\frac{g^{3}(\xi_{j})}{6}\right]$
$\displaystyle+\mathbb{E}_{\bm{s}}\left[\frac{g^{3}(\xi_{i})}{6}\frac{g^{3}(\xi_{j})}{6}\left(\frac{{\bf
x}_{i}^{T}{\bm{s}}}{p}\right)^{3}\left(\frac{{\bf
x}_{j}^{T}{\bm{s}}}{p}\right)^{3}\right].$

Using Assumption 2, we can prove that $\mathbb{E}|g^{3}(\xi_{i})|^{r}$ is
bounded for all $r\in\mathbb{N}$. $\mathbb{E}\left[\left(\frac{{\bf
x}_{i}^{T}{\bm{s}}}{p}\right)^{2}\left|\frac{{\bf
x}_{j}^{T}{\bm{s}}}{p}\right|^{3}|\frac{g^{3}(\xi_{i})}{6}|\right]=O(n^{-\frac{5}{2}}).$
and
$\mathbb{E}\left[\left|\frac{g^{3}(\xi_{i})}{6}\right|\left|\frac{g^{3}(\xi_{j})}{6}\right|\left|\frac{{\bf
x}_{i}^{T}{\bm{s}}}{p}\right|^{3}\left|\frac{{\bf
x}_{j}^{T}{\bm{s}}}{p}\right|^{3}\right]=O(n^{-3}).$ Hence, by computing the
expectation over ${\bm{s}}$ of the first term, we obtain

$\displaystyle\mathbb{E}_{\bm{s}}\left[\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{i}\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{j}\right]=\frac{2}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\left(\frac{{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{j}}{p}\right)^{2}$
$\displaystyle+\frac{1}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\frac{1}{p^{2}}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{i}{\bf x}_{j}^{T}{\bm{\Sigma}}{\bf
x}_{j}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-\frac{5}{2}})}.$

When $i\neq j$, $\left(\frac{{\bf x}_{i}^{T}{\bm{\Sigma}}{\bf
x}_{j}}{p}\right)^{2}=O_{p}(n^{-1})$. Hence,

$\displaystyle\mathbb{E}_{\bm{s}}\left[\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{i}\left[\tilde{\bm{\kappa}}({\bm{s}})\right]_{j}\right]$
$\displaystyle=\frac{2}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\boldsymbol{\delta}_{i=j}\left(\frac{1}{p}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bm{x}}_{i}\right)^{2}$
$\displaystyle+\frac{1}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\frac{1}{p^{2}}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{i}{\bf x}_{j}^{T}{\bm{\Sigma}}{\bf x}_{j}$
$\displaystyle+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-\frac{5}{2}})}.$

Using (Proof of Theorem 1), we thus obtain

$\displaystyle\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}({\bm{s}})\tilde{\bm{\kappa}}({\bm{s}})^{T}\right]$
$\displaystyle=\frac{2}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\mathop{\rm
diag}\left\\{\left(\frac{1}{p}{\bf x}_{i}^{T}{\bm{\Sigma}}{\bf
x}_{i}\right)^{2}\right\\}_{i=1}^{n}$
$\displaystyle+\frac{1}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\left\\{\frac{1}{p}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf
x}_{i}\right\\}_{i=1}^{n}\left(\left\\{\frac{1}{p}{\bf
x}_{j}^{T}{\bm{\Sigma}}{\bf x}_{j}\right\\}_{j=1}^{n}\right)^{T}$
$\displaystyle+O_{\|.\|}(n^{-\frac{3}{2}}).$

It is easy to see that
$\left\|\frac{2}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\mathop{\rm
diag}\left\\{\left(\frac{1}{p}{\bf x}_{i}^{T}{\bm{\Sigma}}{\bf
x}_{i}\right)^{2}\right\\}_{i=1}^{n}\right\|={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-2})}$.
On the other hand, one can show that we can replace in the second term
$\frac{1}{p}{\bf x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{i}$ by
$\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}^{2}$. This is because

$\displaystyle\left\|\frac{1}{p^{2}}\left\\{\frac{1}{p}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{i}-\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\right\\}_{i=1}^{n}\left(\left\\{\frac{1}{p}{\bf
x}_{j}^{T}{\bm{\Sigma}}{\bf x}_{j}\right\\}_{j=1}^{n}\right)^{T}\right\|_{2}$
$\displaystyle\leq\frac{1}{p^{\frac{3}{2}}}\left\|\left\\{\frac{1}{p}{\bf
x}_{i}^{T}{\bm{\Sigma}}{\bf x}_{i}-\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\right\\}_{i=1}^{n}\right\|_{2}\left\|\left\\{\frac{1}{p^{\frac{3}{2}}}x_{j}^{T}{\bm{\Sigma}}{\bf
x}_{j}\right\\}_{j=1}^{n}\right\|_{2}$
$\displaystyle=O_{p}(n^{-\frac{3}{2}}).$

Putting all the above results together, we obtain

$\displaystyle\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}({\bm{s}})\tilde{\bm{\kappa}}({\bm{s}})^{T}\right]$
$\displaystyle=\frac{1}{p^{2}}\left(\frac{g^{\prime\prime}(0)}{2}\right)^{2}\left(\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\right)^{2}{\bf 11}^{T}$
$\displaystyle+O_{\|.\|}(n^{-\frac{3}{2}}).$

∎

Now using the approximation in (Theorem), we obtain

${\bm{\kappa}}_{c}({\bf s})=g^{\prime}(0)\frac{1}{p}{\bf PXs}+{\bf
P}\tilde{\bm{\kappa}}({\bf s})-\frac{1}{n}{\bf PK1}.$ (33)

######  Theorem (Asymptotic behavior of $\mathbf{Q}_{z}$ and
$\widetilde{\mathbf{Q}}_{z}$).

As in [12, Lemma 1], let Assumption 1 holds, then as $p\to\infty$ and all
$z\in\mathbb{C}\backslash\mathbb{R}_{+}$,

$\displaystyle\widetilde{\mathbf{Q}}_{z}$
$\displaystyle\leftrightarrow-\frac{1}{z}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1},\forall
z\in\mathbb{C}\backslash\textbf{supp}\left(\mathbf{\Sigma}\right),$ (34)

where $m_{z}$ is the unique stieltjes transform solution, for all such $z$, to
the implicit equation

$m_{z}=-\left(cz-\frac{1}{n}\mathop{\rm
tr}{\bm{\Sigma}}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}\right)^{-1},$

and the notation $\mathbf{A}\leftrightarrow\mathbf{B}$ means that as
$p\to\infty$, $\frac{1}{p}\mathop{\rm
tr}\mathbf{M\left(A-B\right)}\to_{a.s.}0$ and
$\bm{u}^{T}\left(\mathbf{A}-\mathbf{B}\right)\bm{v}\to_{a.s.}0$, for all
deterministic Hermitian matrices $\mathbf{M}$ and deterministic vectors
$\bm{u},\bm{v}$ of bounded norms. Moreover, from [22] and [23] and
$z\in\mathbb{C}\backslash\mathbb{R}_{+}$,

$\displaystyle\frac{1}{n}\mathop{\rm tr}{\bf MQ}_{z}+\frac{1}{z}\mathop{\rm
tr}{\bf
M}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}=\frac{1}{n}\psi_{n}(z),$
(35) $\displaystyle{\bm{u}}^{T}{\bf
Q}_{z}{\bm{v}}+\frac{1}{z}{\bm{u}}^{T}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}{\bm{v}}=\frac{1}{\sqrt{n}}h_{n}(z),$
(36)

where for all $k\in\mathbb{N}$, $\mathbb{E}\left|\psi_{n}(z)\right|^{k}$ and
$\mathbb{E}\left|h_{n}(z)\right|^{k}$ can be bounded uniformly in $n$ over any
compact at a macroscopic distance from the limiting support of
$\frac{1}{p}{\bf PXX}^{T}{\bf P}$.

######  Theorem 3 (An Integration by parts formula for Gaussian functionals).

[22] With $f$ satisfying Assumption 3 and for
$\bm{x}=\left[x_{1},\cdots,x_{p}\right]^{T}\sim\mathcal{N}\left(\mathbf{0}_{p},\mathbf{\Sigma}\right)$,
we have

$\displaystyle\mathbb{E}\left[x_{i}f\left(\bm{x}\right)\right]=\sum_{j=1}^{p}\left[\mathbf{\Sigma}\right]_{i,j}\mathop{\mathbb{E}}\left[\frac{\partial
f\left(\bm{x}\right)}{\partial x_{j}}\right],$ (37)

or equivalently,
$\mathbb{E}\left[\bm{x}f\left(\bm{x}\right)\right]=\mathbf{\Sigma}\mathop{\mathbb{E}}\left[\nabla
f\left(\bm{x}\right)\right].$

######  Theorem (Nash-Poincaré inequality).

[22] With $f$ satisfying Assumption 3 and for
$\bm{x}=\left[x_{1},\cdots,x_{p}\right]^{T}\sim\mathcal{N}\left(\mathbf{0}_{p},\mathbf{\Sigma}\right)$,
we have under the setting of the previous theorem, ${\rm var}({f}({\bf
x}))\leq\mathbb{E}\left[\nabla f({\bf x})^{T}{\bm{\Sigma}}\nabla f({\bf
x})\right].$

We shall also need the following differentiation formula. For
$i\in\left\\{1,\cdots,n\right\\}$ and $j\in\left\\{1,\cdots,p\right\\}$

$\frac{\partial\widetilde{\bf Q}_{z}}{\partial x_{ij}}=-\widetilde{\bf
Q}_{z}{\bf e}_{i}{\bf e}_{j}^{T}\frac{\bf PX}{p}\widetilde{\bf
Q}_{z}-\widetilde{\bf Q}_{z}\frac{{\bf X}^{T}{\bf P}}{p}{\bf e}_{j}{\bf
e}_{i}^{T}\widetilde{\bf Q}_{z},$ (38)

where ${\bf e}_{i}$ is the all zero vector with 1 at the $i$th entry. With
this background on the asymptotic behavior of kernel matrices, in the
following, we derive the limiting prediction and empirical risks. Recall that
the prediction risk writes as

$\displaystyle\mathcal{R}_{\rm test}$
$\displaystyle=\mathbb{E}_{{\bm{s},\epsilon}}\left[\left|{\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
Py}+\overline{y}-{f}({\bf s})\right|^{2}\right],$

where ${\bf y}=\left[y_{1},\cdots,y_{n}\right]^{T}={\bf f}({\bf
X})+{\bm{\epsilon}}$ and $\bar{y}=\frac{{\bf 1}^{T}}{n}{\bf f}({\bf
X})+\frac{1}{n}{\bm{\epsilon}}^{T}{\bf 1}$. Due to the independence of ${\bf
s}$ and ${\bm{\epsilon}}$, the prediction risk can be decomposed into a
variance and bias terms as $\mathcal{R}_{\rm test}=B+V,$ where

$\displaystyle V$
$\displaystyle=\mathbb{E}_{\bm{s},\epsilon}\left[\left|{\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
P}{\bm{\epsilon}}+\frac{1}{n}{\bm{\epsilon}}^{T}{\bf 1}\right|^{2}\right].$
$\displaystyle B$ $\displaystyle=\mathbb{E}_{\bf
s}\left[\left({\bm{\kappa}}_{c}({\bf s})^{T}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}{\bf Pf}({\bf X})+\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf
X})-f({\bf s})\right)^{2}\right].$

Now, computing the expectation over ${\bm{\epsilon}}$, we obtain

$\displaystyle V$
$\displaystyle=\sigma^{2}\mathbb{E}_{s}\left[{\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bm{\kappa}}_{c}({\bf s})\right]$
$\displaystyle+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-1})}$
$\displaystyle=\sigma^{2}\mathbb{E}_{s}\Biggl{[}{\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-2}{\bm{\kappa}}_{c}({\bf
s})$ $\displaystyle-\left|{\bm{\kappa}}_{c}({\bf s})^{T}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}\Biggr{]}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-1})}.$

Let us start by controlling the second term. Replacing ${\bm{\kappa}}_{c}({\bf
s})$ by (33),

$\displaystyle\mathbb{E}_{s}\left[\left|{\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}\right]$ $\displaystyle\leq
3g^{\prime}(0)^{2}\mathbb{E}_{\bm{s}}\left|\frac{1}{p}{\bf s}^{T}{\bf
X}^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}$
$\displaystyle+3\mathbb{E}_{\bm{s}}\left|\tilde{\bm{\kappa}}({\bf s})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-1})}.$

Computing the expectation over ${\bm{s}}$ of the term term of the last
inequality, we can show that

$\displaystyle\mathbb{E}_{\bm{s}}\left[\left|\frac{1}{p}{\bf s}^{T}{\bf
X}^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}\right]$ $\displaystyle=\frac{1}{p}\frac{{\bf
1}^{T}}{n}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf P}{\bf
X}{\bm{\Sigma}}{\bf X}^{T}{\bf P}}{p}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}{\bf 1}$
$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-1})}.$

On the other hand, from (32), we have

$\mathbb{E}\left|\tilde{\bm{\kappa}}({\bf s})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-1})}.$

The above approximations thus yield

$\mathbb{E}_{s}\left[\left|{\bm{\kappa}}_{c}({\bf s})^{T}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{{\bf
1}}{\sqrt{n}}\right|^{2}\right]=O_{p}(n^{-1}).$

It remains now to compute the first term. Using (32) along with (33), we have

$\displaystyle V$
$\displaystyle=\sigma^{2}\mathbb{E}_{\bm{s}}\left[\frac{1}{p}(g^{\prime}(0))^{2}{\bm{s}}^{T}{\bf
X}^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-2}{\bf P}{\bf
X}{\bm{s}}\right]+O_{p}(n^{-1})$
$\displaystyle=\sigma^{2}\frac{1}{p}(g^{\prime}(0))^{2}\mathop{\rm
tr}{\bm{\Sigma}}{\bf X}^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-2}{\bf P}{\bf X}+O_{p}(n^{-1}).$

From (LABEL:K_infty), $V$ can be approximated as

$\displaystyle V\stackrel{{\scriptstyle(a)}}{{=}}$
$\displaystyle\sigma^{2}\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}{\bf X}^{T}{\bf
P}{\bf Q}_{z}^{2}{\bf P}{\bf X}+O_{p}(n^{-1})$
$\displaystyle=\frac{\sigma^{2}}{p}\frac{\partial}{\partial
t}\left[\mathop{\rm tr}{\bm{\Sigma}}{\bf X}^{T}{\bf P}{\bf Q}_{t}{\bf P}{\bf
X}\right]_{t=z}+O_{p}(n^{-1})$
$\displaystyle{=}\sigma^{2}\frac{\partial}{\partial
t}\left[\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}{\bf X}^{T}{\bf P}{\bf
X}\widetilde{\bf Q}_{t}\right]_{t=z}+O_{p}(n^{-1})$
$\displaystyle=\sigma^{2}\frac{\partial}{\partial
t}\left[\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}+\frac{t}{p}\mathop{\rm
tr}{\bm{\Sigma}}\widetilde{\bf Q}_{t}\right]_{t=z}+O_{p}(n^{-1})$
$\displaystyle=\frac{\partial}{\partial
t}\left[\frac{t\sigma^{2}}{p}\mathop{\rm tr}{\bm{\Sigma}}{\widetilde{\bf
Q}_{t}}\right]_{t=z}+O_{p}(n^{-1}),$

where in $(a)$ the contribution of matrix
$\frac{\frac{\nu}{g^{\prime}\left(0\right)}\mathbf{Q}_{z}\frac{1}{n}\bm{11}^{T}\mathbf{Q}_{z}}{1-\frac{\nu}{g^{\prime}\left(0\right)}\frac{1}{n}\bm{1}^{T}\mathbf{Q}_{z}\bm{1}}$
has been discarded since it only induces terms of order $O_{p}(n^{-1})$. Using
(35), we thus have

$\displaystyle\frac{\partial}{\partial
t}\left[\frac{t\sigma^{2}}{p}\mathop{\rm tr}{\bm{\Sigma}}{\widetilde{\bf
Q}_{t}}\right]_{t=z}$ $\displaystyle=\frac{\partial}{\partial
t}\left[-\frac{\sigma^{2}}{p}\mathop{\rm tr}{\bm{\Sigma}}\mathop{\rm
tr}\left({\bf I}_{p}+m_{t}{\bm{\Sigma}}\right)^{-1}\right]_{t=z}$
$\displaystyle+O_{p}(n^{-1}).$

Taking the derivative over $z$, we find after simple calculations that

$\displaystyle\frac{\partial}{\partial
t}\left[\frac{t\sigma^{2}}{p}\mathop{\rm tr}{\bm{\Sigma}}{\widetilde{\bf
Q}_{t}}\right]_{t=z}$
$\displaystyle=\frac{\sigma^{2}m_{z}^{2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$
$\displaystyle\times\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}+O_{p}(n^{-1}).$

Putting all the above derivations together, we finally obtain

$V=\frac{\sigma^{2}m_{z}^{2}}{n-m_{z}^{2}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}+O_{p}(n^{-1}).$

Evaluation of the bias term.

To begin with, we first expand $B$ as

$\displaystyle B$ $\displaystyle=\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf
X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}{\bm{\kappa}}_{c}({\bf s}){\bm{\kappa}}_{c}({\bf
s})^{T}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf P}{\bf f}({\bf
X})\right]$ $\displaystyle+\mathbb{E}_{\bm{s}}\left|\frac{1}{n}{\bf 1}^{T}{\bf
f}({\bf X})-f({\bf s})\right|^{2}+2\mathbb{E}_{\bm{s}}\Biggl{[}{\bf f}({\bf
X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}{\bm{\kappa}}_{c}({\bf s})$
$\displaystyle\times\left(\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf X})-f({\bf
s})\right)\Biggr{]}.$

We will sequentially treat the above three terms. To begin with, we control
first
$\mathbb{E}_{\bm{s}}\left[{\bm{\kappa}}_{c}({\bm{s}}){\bm{\kappa}}_{c}({\bm{s}})^{T}\right]$
which we expand as

$\displaystyle\mathbb{E}_{\bm{s}}\left[{\bm{\kappa}}_{c}({\bm{s}}){\bm{\kappa}}_{c}({\bm{s}})^{T}\right]$
$\displaystyle=(g^{\prime}(0))^{2}\frac{1}{p^{2}}{\bf PX}{\bm{\Sigma}}{\bf
X}^{T}{\bf P}+g^{\prime}(0)\frac{1}{p}{\bf
PX}{\bm{s}}\tilde{\bm{\kappa}}(s)^{T}{\bf P}$
$\displaystyle+\frac{g^{\prime}(0)}{p}{\bf
P}\mathbb{E}_{\bm{s}}\left[\tilde{\bm{\kappa}}({\bm{s}}){\bf s}^{T}\right]{\bf
XP}+{\bf P}\mathbb{E}_{\bm{s}}\tilde{\bm{\kappa}}({\bf
s})\tilde{\bm{\kappa}}({\bf s})^{T}{\bf P}-\frac{1}{n}{\bf
P}\mathbb{E}_{\bm{s}}\tilde{\bm{\kappa}}({\bm{s}}){\bf 1}^{T}{\bf KP}$
$\displaystyle-\frac{1}{n}{\bf
PK1}\mathbb{E}_{\bm{s}}\tilde{\bm{\kappa}}({\bm{s}})^{T}{\bf
P}+\frac{1}{n^{2}}{\bf PK11}^{T}{\bf P}.$

Using (33) along with Lemma 2, we obtain

$\displaystyle\mathbb{E}_{\bm{s}}\left[{\bm{\kappa}}_{c}({\bm{s}}){\bm{\kappa}}_{c}({\bm{s}})^{T}\right]$
$\displaystyle=(g^{\prime}(0))^{2}\frac{1}{p^{2}}{\bf PX}{\bm{\Sigma}}{\bf
X}^{T}{\bf P}+\frac{1}{n^{2}}{\bf PK11}^{T}{\bf KP}$
$\displaystyle+O_{\|.\|}(n^{-\frac{3}{2}}).$

Replacing ${\bf K}$ by ${\bf K}_{\infty}$, we thus obtain

$\displaystyle\mathbb{E}_{\bm{s}}\left[{\bm{\kappa}}_{c}({\bm{s}}){\bm{\kappa}}_{c}({\bm{s}})^{T}\right]$
$\displaystyle=(g^{\prime}(0))^{2}\frac{1}{p^{2}}{\bf PX}{\bm{\Sigma}}{\bf
X}^{T}{\bf P}$ $\displaystyle+\frac{(g^{\prime}(0))^{2}}{n^{2}}{\bf
P}\frac{{\bf X}{\bf X}^{T}}{p}{\bf 11}^{T}\frac{{\bf XX}^{T}}{p}{\bf
P}+O_{\|.\|}(n^{-\frac{3}{2}}).$

From Assumption 3, we can prove that

$\left\|{\bf P}{\bf f}({\bf X})\right\|=O_{p}(\sqrt{n}).$

Hence, the first term in $B$ can be approximated as

$\displaystyle\mathbb{E}_{\bm{s}}\Biggl{[}{\bf f}({\bf X})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bm{\kappa}}_{c}({\bf
s}){\bm{\kappa}}_{c}({\bf s})^{T}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}$ $\displaystyle\times{\bf P}{\bf f}({\bf X})\Biggr{]}$
$\displaystyle=\left(g^{\prime}(0)\right)^{2}{\bf f}({\bf X})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf P}\frac{{\bf
X}{\bm{\Sigma}}{\bf X}^{T}}{p^{2}}{\bf P}$ $\displaystyle\times\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf P}{\bf f}({\bf
X})+\left(g^{\prime}(0)\right)^{2}{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}$ $\displaystyle\times\frac{1}{n^{2}}{\bf
P}\frac{{\bf XX}^{T}}{p}{\bf 11}^{T}\frac{{\bf XX}^{T}}{p}{\bf P}({\bf
K}_{c}+\lambda{\bf I}_{n})^{-1}{\bf Pf}({\bf X})$
$\displaystyle+O_{p}(n^{-\frac{1}{2}})$
$\displaystyle=Z_{1}+Z_{2}+O_{p}(n^{-\frac{1}{2}}).$

We will start by treating $Z_{1}$. From (LABEL:K_infty), we can show that

$Z_{1}=\frac{1}{p^{2}}{\bf f}({\bf X})^{T}{\bf P}{\bf X}\widetilde{\bf
Q}_{z}{\bm{\Sigma}}\widetilde{\bf Q}_{z}{\bf X}^{T}{\bf P}{\bf f}({\bf
X})+O_{p}(n^{-\frac{1}{2}}).$

To treat $Z_{1}$, we first decompose ${\bf P}{\bf f}({\bf X})$ as

${\bf P}{\bf f}({\bf X})={\bf f}({\bf X})-\mathbb{E}_{\bm{x}}(f(\bm{x})){\bf
1}+\mathbb{E}_{\bm{x}}(f(\bm{x})){\bf 1}-\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf
X}),$

where from the classical probability results, we have

$\left\|\mathbb{E}_{{\bf x}}(f({\bf x})){\bf 1}-\frac{1}{n}{\bf 1}^{T}{\bf
f}({\bf X}){\bf 1}\right\|_{2}=O_{p}(1).$

Hence, we can replace ${\bf P}{\bf f}({\bf X})$ by
$\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X})={\bf f}({\bf
X})-\mathbb{E}_{\bm{x}}(f(\bm{x})){\bf 1}$ up to an error
$O_{p}(n^{-\frac{1}{2}})$, thus yielding

$Z_{1}=\frac{1}{p^{2}}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X})^{T}{\bf
X}\widetilde{\bf Q}_{z}{\bm{\Sigma}}\widetilde{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})+O_{p}(n^{-\frac{1}{2}}).$

To treat $Z_{1}$, we shall resort to the following lemma.

######  Lemma 3.

Let ${\bf A}$ be a $p\times p$ symmetric matrix with a uniformly bounded
spectral norm. Let $f$ satisfy Assumption 3 and
$z\in\mathbb{C}\backslash\mathbb{R}_{+}$. Consider the following function.

$\displaystyle{h}({\bf X})$
$\displaystyle=\frac{1}{p^{2}}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})^{T}{\bf X}\widetilde{\bf Q}_{z}{\bf A}\widetilde{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X}).$

Then, for any $\delta>0$, ${\rm var}({h}({\bf X}))=O(n^{-1+\delta}).$

###### Proof.

The proof of Lemma 3 follows from the Nash-Poincaré inequality and the
differentiation formula in (38). Therefore, we omit technical details for
brievity. ∎

As per Lemma 3, $Z_{1}$ can be approximated as

$Z_{1}=\mathbb{E}\left[\frac{1}{p^{2}}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\widetilde{\bf Q}_{z}{\bm{\Sigma}}\widetilde{\bf
Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\right]+O_{p}(n^{-\frac{1}{2}+\delta}),$

for any $\delta>0$. Now let $\overline{\bf x}=\frac{1}{\sqrt{n}}{\bf
X}^{T}{\bf 1}$. The resolvent matrix $\widetilde{\bf Q}_{z}$ can be expressed
as

$\widetilde{\bf Q}_{z}=\left(\frac{1}{p}{\bf X}^{T}{\bf
X}-\frac{1}{p}\overline{\bf x}\overline{\bf x}^{T}-z{\bf I}_{p}\right)^{-1}.$

By the inversion Lemma, matrix $\widetilde{\bf Q}_{z}$ can be written as

$\widetilde{\bf Q}_{z}=\overline{\bf Q}_{z}+\frac{1}{p}\frac{\overline{\bf
Q}_{z}\overline{\bf x}\overline{\bf x}^{T}\overline{\bf
Q}_{z}}{1-\frac{1}{p}\overline{\bf x}^{T}\overline{\bf Q}_{z}\overline{\bf
x}},$

where $\overline{\bf Q}_{z}=\left(\frac{1}{p}{\bf X}^{T}{\bf X}-z{\bf
I}_{p}\right)^{-1}$. Hence $Z_{1}$ can be expanded as

$\displaystyle Z_{1}$
$\displaystyle=\frac{1}{p^{2}}\mathbb{E}\left[\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf
Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X})\right]$
$\displaystyle+\frac{1}{p^{3}}\mathbb{E}\left[\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}\overline{\bf x}\overline{\bf
x}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\left(1-\frac{1}{p}\overline{\bf x}^{T}\overline{\bf Q}_{z}\overline{\bf
x}\right)^{-1}\right]$
$\displaystyle+\frac{1}{p^{3}}\mathbb{E}\left[\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf
Q}_{z}\overline{\bf x}\overline{\bf x}^{T}{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\left(1-\frac{1}{p}\overline{\bf x}^{T}\overline{\bf Q}_{z}\overline{\bf
x}\right)^{-1}\right]$
$\displaystyle+\frac{1}{p^{4}}\mathbb{E}\Biggl{[}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}\overline{\bf x}\overline{\bf
x}^{T}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf Q}_{z}\overline{\bf
x}\overline{\bf x}^{T}\overline{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X})$
$\displaystyle\times\left(1-\frac{1}{p}\overline{\bf x}^{T}\overline{\bf
Q}_{z}\overline{\bf x}\right)^{-2}\Biggr{]}.$

We can show that the last three terms are $O(n^{-\frac{1}{2}+\delta})$. To
illustrate this, we will focus on the second term, as the derivations are
similar for the remaining quantities. By Cauchy-Schwartz inequality, we have

$\displaystyle\mathbb{E}\frac{1}{p^{3}}\left[\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}\overline{\bf x}\overline{\bf
x}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\left(1-\frac{1}{p}\overline{\bf x}^{T}\overline{\bf Q}_{z}\overline{\bf
x}\right)^{-1}\right]$
$\displaystyle\leq\sqrt{\mathbb{E}\left|\frac{1}{p^{\frac{3}{2}}}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}{\bf X}^{T}\frac{\bf
1}{\sqrt{n}}\right|^{2}}\sqrt{\mathbb{E}\left|\frac{1}{p^{\frac{3}{2}}}\frac{{\bf
1}^{T}}{\sqrt{n}}{\bf X}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf
Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\right|^{2}}.$

To treat the above bound, we first show that

$\mathbb{E}\left|\frac{1}{p^{\frac{3}{2}}}\frac{{\bf 1}^{T}}{\sqrt{n}}{\bf
X}\overline{\bf Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\right|^{2}=O(n^{-1+\delta}),$

for any $\delta>0$. Towards this end, we need the following control of the
variance of $\frac{1}{p^{\frac{3}{2}}}\frac{{\bf 1}^{T}}{\sqrt{n}}{\bf
X}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf Q}_{z}{\bf
X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf X})$ which we can be shown
to be $O(n^{-1+\delta})$.

######  Lemma 4.

Let $g$ satisfy Assumption 3 and $z\in\mathbb{C}\backslash\mathbb{R}_{+}$.
Consider the following function.

${g}({\bf X})=\frac{1}{p^{\frac{3}{2}}}\frac{{\bf 1}^{T}}{\sqrt{n}}{\bf
X}\overline{\bf Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X}).$

Then, for any $\delta>0$ ${\rm var}(g({\bf X}))=O(n^{-1+\delta}).$

The proof of Lemma 4 follows the same lines as Lemma 3 and is thus omitted.
With the control of the variance at hand, it suffices to show the following.

$\mathbb{E}\left[g({\bf X})\right]=O(n^{-\frac{1}{2}+\delta}).$ (39)

Let $\overline{\bf
Q}_{z,i}=\left(\sum_{k=1,k\notin\left\\{i\right\\}}^{n}\frac{1}{p}{\bf
x}_{k}{\bf x}_{k}^{T}-z{\bf I}_{p}\right)^{-1}$. We thus develop $g({\bf X})$
as

$\displaystyle\frac{1}{\sqrt{n}}\frac{1}{p^{\frac{3}{2}}}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}\left[{\bf
x}_{i}^{T}\overline{\bf Q}_{z}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf x}_{j})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\frac{1}{p^{\frac{3}{2}}}\sum_{i=1}^{n}\sum_{j=1,i\neq
j}^{n}\mathbb{E}\left[\frac{{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{j}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf x}_{j})}{1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{i}}\right]$
$\displaystyle+\frac{1}{\sqrt{n}}\frac{1}{p^{\frac{3}{2}}}\sum_{i=1}^{n}\mathbb{E}\left[\frac{{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{i}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf x}_{i})}{1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{i}}\right]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\frac{1}{\sqrt{n}}\frac{1}{p^{\frac{3}{2}}}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf x}_{j})$ $\displaystyle\times\left[\left(1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{i}\right)^{-1}-\left(1+\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}\overline{\bf Q}_{z,i}\right)^{-1}\right]$
$\displaystyle+O(n^{-\frac{1}{2}})$
$\displaystyle=\frac{1}{\sqrt{n}p^{\frac{3}{2}}}\mathbb{E}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}$ $\displaystyle\times\mathop{\rm
diag}\left\\{\frac{1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}{\bf
x}_{i}}{1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{i}}-\frac{1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}{\bf
x}_{i}}{1+\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}\overline{\bf
Q}_{z,i}}\right\\}_{i=1}^{n}$ $\displaystyle\times{\bf X}^{T}{\bf
1}+O(n^{-\frac{1}{2}}),$

where in (a), we use the fact that when $i\neq j$, $\mathbb{E}\left[{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf x}_{j})\left(1+\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}\overline{\bf
Q}_{z,i}\right)^{-1}\right]=0$. Using the fact that

$\displaystyle\left\|\mathop{\rm diag}\left\\{\frac{1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}{\bf x}_{i}}{1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{i}}-\frac{1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}{\bf x}_{i}}{1+\frac{1}{p}\mathop{\rm
tr}{\bm{\Sigma}}\overline{\bf Q}_{z,i}}\right\\}_{i=1}^{n}\right\|$
$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}O_{p}(n^{-\frac{1}{2}+\delta})}.$

we conclude that

$Z_{1}=\mathbb{E}\left[\frac{1}{p^{2}}\stackrel{{\scriptstyle\circ}}{{\bf
f}}({\bf X})^{T}{\bf X}\overline{\bf Q}_{z}{\bm{\Sigma}}\overline{\bf
Q}_{z}{\bf X}^{T}\stackrel{{\scriptstyle\circ}}{{\bf f}}({\bf
X})\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$

From the inversion lemma, we get

$\overline{\bf Q}_{z}{\bf x}_{i}=\frac{\overline{\bf Q}_{z,i}{\bf
x}_{i}}{1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf x}_{i}}.$

$Z_{1}$ writes as

$\displaystyle Z_{1}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{p^{2}}\mathbb{E}\left[\frac{\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bm{\Sigma}}\overline{\bf
Q}_{z,j}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{j})}{(1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{i})(1+\frac{1}{p}{\bf x}_{j}^{T}\overline{\bf Q}_{z,j}{\bf
x}_{j})}\right]$ $\displaystyle+O_{p}(n^{-\frac{1}{2}+\delta}).$

Using the fact that

$\left|\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bf
x}_{i}+\frac{1}{pz}\mathop{\rm tr}{\bm{\Sigma}}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-1}\right|^{k}=O(n^{-\frac{1}{2}}).$

Along with the relation

$1-\frac{1}{pz}\mathop{\rm tr}{\bm{\Sigma}}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-1}=-zcm_{z},$

we obtain

$\displaystyle Z_{1}$
$\displaystyle=c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{p^{2}}\mathbb{E}\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bm{\Sigma}}\overline{\bf
Q}_{z,j}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf x}_{j})$
$\displaystyle=c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\frac{1}{p^{2}}\mathbb{E}\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,i}{\bm{\Sigma}}\overline{\bf
Q}_{z,j}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf x}_{j})$
$\displaystyle+c^{2}z^{2}m_{z}^{2}\frac{1}{p^{2}}\sum_{i=1}^{n}\mathbb{E}\left[\left|\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i})\right|^{2}{\bf x}_{i}^{T}\overline{\bf
Q}_{z,i}{\bm{\Sigma}}\overline{\bf Q}_{z,i}{\bf
x}_{i}\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$ (40)

Now, we further proceed resorting to the inversion lemma

$\overline{\bf Q}_{z,i}=\overline{\bf Q}_{z,ij}-\frac{1}{p}\frac{\overline{\bf
Q}_{z,ij}{\bf x}_{j}{\bf x}_{j}^{T}\overline{\bf Q}_{z,ij}}{1+\frac{1}{p}{\bf
x}_{j}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{j}},$

which once plugged into (40) yields

$\displaystyle Z_{1}$
$\displaystyle=c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbb{E}\left[\frac{1}{p^{2}}\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}\overline{\bf
Q}_{z,ij}{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf x}_{j})\right]$
$\displaystyle-c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\mathbb{E}\Biggl{[}\frac{1}{p^{3}}\frac{\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{{\bf x}_{j}{\bf
x}_{j}^{T}}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}\overline{\bf
Q}_{z,ij}}{1+\frac{1}{p}{\bf x}_{j}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{j}}$
$\displaystyle\times{\bf x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{j})\Biggr{]}-c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\mathbb{E}\Biggl{[}\frac{1}{p^{3}}\frac{\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}}{1+\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{i}}$
$\displaystyle\times\overline{\bf Q}_{z,ij}{\bf x}_{i}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf
x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf x}_{j})\Biggr{]}$
$\displaystyle+c^{2}z^{2}m_{z}^{2}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\mathbb{E}\Biggl{[}\frac{1}{p^{4}}\frac{\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i}){\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{j}{\bf
x}_{j}^{T}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}\overline{\bf
Q}_{z,ij}}{(1+\frac{1}{p}{\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf
x}_{i})(1+\frac{1}{p}{\bf x}_{j}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{j})}$
$\displaystyle\times{\bf x}_{i}{\bf x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf
x}_{j}\stackrel{{\scriptstyle\circ}}{{f}}({\bf x}_{j})\Biggr{]}$
$\displaystyle\triangleq
Z_{11}+Z_{12}+Z_{13}+Z_{14}+O_{p}(n^{-\frac{1}{2}+\delta}).$

Let us first control $Z_{11}$. Taking the expectation over ${\bf x}_{i}$ and
${\bf x}_{j}$, we obtain

$\displaystyle Z_{11}$
$\displaystyle=c^{2}z^{2}m_{z}^{2}\frac{1}{p^{2}}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\mathbb{E}\left[\nabla
f\right]^{T}\mathbf{\Sigma}\mathbb{E}\left[\overline{\bf
Q}_{z,ij}{\bm{\Sigma}}\overline{\bf
Q}_{z,ij}\right]\mathbf{\Sigma}\mathbb{E}\left[\nabla f\right]$
$\displaystyle+c^{2}z^{2}m_{z}^{2}\frac{1}{p^{2}}\sum_{i=1}^{n}\mathbb{E}\left[\left|\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i})\right|^{2}{\bf x}_{i}^{T}\overline{\bf
Q}_{z,ij}\boldsymbol{\Sigma}\overline{\bf Q}_{z,ij}{\bf x}_{i}\right].$

It can be shown that for a vector ${\bf a}$ with uniformly bounded norm that

$\mathbb{E}\left[{\bf a}^{T}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}\overline{\bf
Q}_{z,ij}{\bf b}\right]=\mathbb{E}\left[{\bf a}^{T}\overline{\bf
Q}_{z}{\bm{\Sigma}}\overline{\bf Q}_{z}{\bf b}\right]+O(n^{-1}).$

Using the fact that $\mathbb{E}\left[\left|\frac{1}{p}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,ij}\boldsymbol{\Sigma}\overline{\bf
Q}_{z,ij}{\bf x}_{i}-\mathbb{E}\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}{\bf
Q}{\bm{\Sigma}}{\bf Q}\right|^{2}\right]=O(n^{-1})$, we thus obtain

$\displaystyle Z_{11}$
$\displaystyle=c^{2}z^{2}m_{z}^{2}\frac{1}{p^{2}}\sum_{i=1}^{n}\sum_{j=1,j\neq
i}^{n}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}\mathbb{E}\left[\overline{\bf
Q}_{z}{\bm{\Sigma}}\overline{\bf
Q}_{z}\right]{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]$
$\displaystyle+c^{2}z^{2}m_{z}^{2}\frac{1}{p}\sum_{i=1}^{n}{\bf
var}_{f}\mathbb{E}\left[\frac{1}{p}\mathop{\rm tr}{\bm{\Sigma}}{\bf
Q}{\bm{\Sigma}}{\bf Q}\right]+O(n^{-\frac{1}{2}}).$

By standard results from random matrix theory [23, 24], we get

$\displaystyle Z_{11}$ $\displaystyle=\frac{nm_{z}^{2}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-1}{\bm{\Sigma}}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-1}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\left({\bf I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}}$
$\displaystyle\times{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]+{\bf
var}_{f}\frac{m_{z}^{2}\mathop{\rm tr}{\bm{\Sigma}}^{2}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\left({\bf I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}}$
$\displaystyle+O(n^{-\frac{1}{2}}).$

Along the same arguments, we can show that

$\displaystyle Z_{12}=Z_{13}$ $\displaystyle=\frac{m_{z}^{3}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}\left({\bf I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}}$
$\displaystyle\times\mathbb{E}\left[\nabla f\right]^{T}{\bm{\Sigma}}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-1}{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]+O(n^{-\frac{1}{2}}).$

As for $Z_{14}$, using Hölder’s inequality, it can be bounded as

$\displaystyle\left|Z_{14}\right|$
$\displaystyle\leq\frac{c^{2}z^{2}m_{z}^{2}}{p^{2}\sqrt{p}}\sum_{i=1}\sum_{j=1,j\neq
i}\left(\mathbb{E}\left|\frac{1}{\sqrt{p}}{\bf x}_{i}^{T}\overline{\bf
Q}_{z,ij}{\bf x}_{j}\right|^{4}\right)^{\frac{1}{4}}$
$\displaystyle\times\left(\mathbb{E}\left|\frac{1}{\sqrt{p}}{\bf
x}_{j}^{T}\overline{\bf Q}_{z,ij}{\bm{\Sigma}}\overline{\bf Q}_{z,ij}{\bf
x}_{i}\right|^{4}\right)^{\frac{1}{4}}\left(\mathbb{E}\left|\frac{1}{\sqrt{p}}{\bf
x}_{i}^{T}\overline{\bf Q}_{z,ij}{\bf x}_{j}\right|^{4}\right)^{\frac{1}{4}}$
$\displaystyle\times\left(\mathbb{E}\left|\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{i})\stackrel{{\scriptstyle\circ}}{{f}}({\bf
x}_{j})\right|^{4}\right)^{\frac{1}{4}}$ $\displaystyle=O(n^{-\frac{1}{2}}).$

We thus conclude that

$\displaystyle Z_{1}$ $\displaystyle={\bf var}_{f}\frac{m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$
$\displaystyle+\frac{nm_{z}^{2}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}({\bf I}_{p}+m_{z}{\bf\Sigma})^{-1}{\bm{\Sigma}}({\bf
I}_{p}+m_{z}{\bf\Sigma})^{-1}{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]}{n-m_{z}^{2}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$ $\displaystyle-\frac{2m_{z}^{3}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}({\bf I}_{p}+m_{z}{\bf\Sigma})^{-1}$
$\displaystyle\times{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$

Now, we will treat the term $Z_{2}$. Similarly to $Z_{1}$, we can show that

$Z_{2}=\frac{1}{np^{2}}{\bf f}({\bf X})^{T}{\bf P}{\bf X}\widetilde{\bf
Q}_{z}{\bf X}^{T}\frac{{\bf 1}{\bf 1}^{T}}{n}{\bf X}\widetilde{\bf Q}_{z}{\bf
X}^{T}{\bf P}{\bf f}({\bf X})+O_{p}(n^{-\frac{1}{2}}).$

Using Lemma 4 along with (39), we thus obtain

$Z_{2}=O_{p}(n^{-\frac{1}{2}}).$

We have thus completed the treatment of the first term of the bias and shown
that

$\displaystyle\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bm{\kappa}}_{c}({\bf
s}){\bm{\kappa}}_{c}({\bf s})^{T}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}{\bf P}{\bf f}({\bf X})\right]$ $\displaystyle={\bf
var}_{f}\frac{m_{z}^{2}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$
$\displaystyle+\frac{nm_{z}^{2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}({\bf I}_{p}+m_{z}{\bf\Sigma})^{-1}$
$\displaystyle\times{\bm{\Sigma}}({\bf
I}_{p}+m_{z}{\bf\Sigma})^{-1}{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]-\frac{2m_{z}^{3}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}{n-m_{z}^{2}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$
$\displaystyle\times\mathbb{E}\left[\nabla f\right]^{T}{\bm{\Sigma}}({\bf
I}_{p}+m_{z}{\bf\Sigma})^{-1}{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$

The second term in the bias can be dealt with by noticing that

$\frac{1}{n}{\bf 1}^{T}f({\bf X})=\mathbb{E}f({\bf
x})+O_{p}(\frac{1}{\sqrt{n}}).$

Thus yielding

$\mathbb{E}_{\bm{s}}\left|\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf X})-f({\bf
s})\right|^{2}={\bf var}_{f}+O_{p}(n^{-\frac{1}{2}}).$

We now move to the last term in the bias. Using (33), we obtain

$\displaystyle 2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}\tilde{\bm{\kappa}}_{c}({\bm{s}})\left(\frac{1}{n}{\bf
1}^{T}{\bf f}({\bf X})-f({\bm{s}})\right)\right]$
$\displaystyle=2\mathbb{E}_{\bm{s}}\Biggl{[}{\bf f}({\bf X})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}\left\\{g^{\prime}(0)\frac{1}{p}{\bf PX}{\bm{s}}+{\bf
P}\tilde{\bm{\kappa}}_{c}({\bm{s}})-\frac{1}{n}{\bf PK}{\bf 1}\right\\}$
$\displaystyle\times\left(\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf
X})-f({\bm{s}})\right)\Biggr{]}$
$\displaystyle=-2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}g^{\prime}(0)\frac{1}{p}{\bf
PX}{\bm{s}}f({\bm{s}})\right]$ $\displaystyle+2\mathbb{E}_{\bm{s}}\left[{\bf
f}({\bf X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
P}\tilde{\bm{\kappa}}_{c}({\bm{s}}){\bm{s}}\left(\frac{1}{n}{\bf 1}^{T}{\bf
f}({\bf X})-f({\bm{s}})\right)\right]$
$\displaystyle-2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}\frac{1}{n}{\bf PK}{\bf
1}\left(\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf X})-f({\bm{s}})\right)\right].$

Since $\mathbb{E}_{\bm{s}}f({\bm{s}})-\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf
X})=O_{p}(n^{-\frac{1}{2}})$, we can replace in the two last terms
$\frac{1}{n}{\bf 1}^{T}{\bf f}({\bf X})$ by $\mathbb{E}_{\bm{s}}f({\bm{s}})$
with an error $O(n^{-\frac{1}{2}})$. In doing so, we obtain

$\displaystyle 2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}\tilde{\bm{\kappa}}_{c}({\bm{s}})\left(\frac{1}{n}{\bf
1}^{T}{\bf f}({\bf X})-f({\bm{s}})\right)\right]$ $\displaystyle=-2{\bf
f}({\bf X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}g^{\prime}(0)\frac{1}{p}{\bf
PX}{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]$
$\displaystyle+2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
P}\tilde{\bm{\kappa}}_{c}({\bm{s}})\left(\mathbb{E}_{\bm{s}}f({\bf
s})-f({\bm{s}})\right)\right]$ $\displaystyle+O_{p}(n^{-\frac{1}{2}}).$

Using (LABEL:K_infty) along with standard calculations as those used in Lemma
4, we obtain

$\displaystyle-2{\bf f}({\bf X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}g^{\prime}(0)\frac{1}{p}{\bf
PX}{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]$ $\displaystyle=-2{\bf f}({\bf
X})^{T}{\bf P}{\bf Q}_{z}\frac{1}{p}{\bf
PX}{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]+O_{p}(n^{-\frac{1}{2}})$
$\displaystyle=-\frac{2}{p}{\bf f}({\bf X})^{T}{\bf P}{\bf X}\widetilde{\bf
Q}_{z}{\bm{\Sigma}}\mathbb{E}\nabla f+O_{p}(n^{-\frac{1}{2}}).$

Replacing $\widetilde{\bf Q}_{z}$ by $\overline{\bf Q}_{z}$, and using similar
derivations as before, we obtain

$\displaystyle-2{\bf f}({\bf X})^{T}{\bf P}\left({\bf K}_{c}+\lambda{\bf
I}_{n}\right)^{-1}g^{\prime}(0)\frac{1}{p}{\bf
PX}{\bm{\Sigma}}\mathbb{E}\left[\nabla f\right]$
$\displaystyle=-2m_{z}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-1}{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$

Finally, using the same tools we can show that

$2\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf P}\left({\bf
K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
P}\tilde{\bm{\kappa}}_{c}({\bm{s}})\left(\mathbb{E}_{\bm{s}}f({\bf
s})-f({\bm{s}})\right)\right]$

will not contribute to the expression of the bias. In fact,

$\displaystyle\frac{2}{p^{2}}\mathbb{E}_{\bm{s}}\left[{\bf f}({\bf X})^{T}{\bf
P}\left({\bf K}_{c}+\lambda{\bf I}_{n}\right)^{-1}{\bf
P}\tilde{\bm{\kappa}}_{c}({\bm{s}})\left(\mathbb{E}_{\bm{s}}f({\bf
s})-f({\bm{s}})\right)\right]$ $\displaystyle=O_{p}(n^{-\frac{1}{2}+\delta}).$

We are now in position to estimate the bias term. Combining the results of all
derivations, for any $\delta>0$ we obtain

$\displaystyle B$ $\displaystyle=\frac{n{\bf var}_{f}}{n-m_{z}^{2}\mathop{\rm
tr}\mathop{\rm tr}{\bm{\Sigma}}^{2}({\bf I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}$
$\displaystyle-\frac{nm_{z}}{n-m_{z}^{2}\mathop{\rm tr}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}\left({\bf I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-1}$
$\displaystyle\times{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]-\frac{nm_{z}^{2}}{n-m_{z}^{2}\mathop{\rm tr}\mathop{\rm
tr}{\bm{\Sigma}}^{2}({\bf
I}_{p}+m_{z}{\bm{\Sigma}})^{-2}}\mathbb{E}\left[\nabla
f\right]^{T}{\bm{\Sigma}}$ $\displaystyle\times\left({\bf
I}_{p}+m_{z}{\bm{\Sigma}}\right)^{-2}\times{\bm{\Sigma}}\mathbb{E}\left[\nabla
f\right]+O_{p}(n^{-\frac{1}{2}+\delta}).$

This concludes the proof.

## Appendix B

## Proof of Theorem 2

The proof heavily relies on the same tools used earlier in the proof of
Theorem 1 and on the following observations

$\displaystyle\frac{1}{p^{2}}\bm{y}^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}\bm{y}$
$\displaystyle=\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)+\frac{1}{p^{2}}\bm{\epsilon}^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}\bm{\epsilon}$
$\displaystyle+O_{p}\left(n^{-\frac{1}{2}}\right)$
$\displaystyle=\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)+\frac{\sigma^{2}}{p}\mathop{\rm
tr}\frac{\mathbf{X}^{T}\mathbf{PX}}{p}\mathbf{\widetilde{Q}}_{z}$
$\displaystyle+O_{p}\left(n^{-\frac{1}{2}}\right)$
$\displaystyle=\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)+\sigma^{2}+z\frac{\sigma^{2}}{p}\mathop{\rm
tr}\mathbf{\widetilde{Q}}_{z}$
$\displaystyle+O_{p}\left(n^{-\frac{1}{2}}\right)$
$\displaystyle=\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)+\sigma^{2}-\frac{\sigma^{2}}{p}\mathop{\rm
tr}\mathbf{\Sigma}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}$
$\displaystyle+O_{p}\left(n^{-\frac{1}{2}}\right).$

It is also possible to show along the same lines as before that

$\displaystyle\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)$
$\displaystyle=-zm_{z}^{2}\mathop{\mathbb{E}}\left[\nabla_{f}\right]^{T}\mathbf{\Sigma}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}\mathbf{\Sigma}\mathop{\mathbb{E}}\left[\nabla_{f}\right]$
$\displaystyle+m_{z}\textbf{var}_{f}\frac{1}{p}\mathop{\rm
tr}\mathbf{\Sigma}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}+O_{p}\left(n^{-\frac{1}{2}+\delta}\right).$

Moreover, by computing

$\displaystyle\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}^{2}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)$
$\displaystyle=\frac{\partial}{\partial
z}\frac{1}{p^{2}}f\left(\mathbf{X}\right)^{T}\mathbf{PX}\mathbf{\widetilde{Q}}_{z}\mathbf{X}^{T}\mathbf{P}f\left(\mathbf{X}\right)$
$\displaystyle=\frac{\partial}{\partial
z}\Biggl{[}-zm_{z}^{2}\mathop{\mathbb{E}}\left[\nabla_{f}\right]^{T}\mathbf{\Sigma}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}\mathbf{\Sigma}\mathop{\mathbb{E}}\left[\nabla_{f}\right]$
$\displaystyle+m_{z}\textbf{var}_{f}\frac{1}{p}\mathop{\rm
tr}\mathbf{\Sigma}\left(\mathbf{I}+m_{z}\mathbf{\Sigma}\right)^{-1}\Biggr{]}+O_{p}\left(n^{-\frac{1}{2}+\delta}\right).$

With the above observation at hand, it is straightforward to show that

$\displaystyle\widehat{\mathcal{R}}_{test}=\mathcal{R}_{test}^{\infty}+O_{p}(n^{-\frac{1}{2}+\delta}),$

which is combined with Theorem 1 gives the claim of Theorem 2.

## References

  * [1] C. M. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006.
  * [2] T. Hofmann, B. Schölkopf, and A. Smola, “Kernel methods in machine learning,” _Annals of Statistics_ , vol. 36, no. 3, pp. 1171–1220, Jun. 2008\.
  * [3] A. E. Alaoui and M. W. Mahoney, “Fast randomized kernel ridge regression with statistical guarantees,” _NIPS_ , vol. 1, pp. 775–783, 2015.
  * [4] A. Rudi and L. Rosasco, “Generalization properties of learning with random features,” in _NIPS_ , 2017.
  * [5] S. Giguère, F. Laviolette, M. Marchand, and K. Sylla, “Risk bounds and learning algorithms for the regression approach to structured output prediction,” in _Proceedings of the 30th International Conference on International Conference on Machine Learning - Volume 28_ , ser. ICML’13. JMLR.org, 2013, pp. I–107–I–114. [Online]. Available: http://dl.acm.org/citation.cfm?id=3042817.3042831
  * [6] A. Caponnetto and E. De Vito, “Optimal rates for the regularized least-squares algorithm,” _Foundations of Computational Mathematics_ , vol. 7, no. 3, pp. 331–368, Jul 2007. [Online]. Available: https://doi.org/10.1007/s10208-006-0196-8
  * [7] L. H. Dicker, D. P. Foster, and D. Hsu, “Kernel ridge vs. principal component regression: Minimax bounds and the qualification of regularization operators,” _Electronic Journal of Statistics_ , vol. 11, no. 1, pp. 1022–1047, 2017.
  * [8] C. Cortes, M. Mohri, and A. Rostamizadeh, “Algorithms for Learning Kernels Based on Centered Alignment,” _Journal of Machine Learning Research_ , no. 13, pp. 795–828, 2012.
  * [9] Y. Lu, L. Wang, J. Lu, J. Yang, and C. Shen, “Multiple kernel clustering based on centered kernel alignment,” _Pattern Recognition_ , vol. 47, no. 11, pp. 3656 – 3664, 2014. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0031320314001903
  * [10] D. Cárdenas-Peña, D. Collazos-Huertas, and G. Castellanos-Dominguez, “Centered kernel alignment enhancing neural network pretraining for mri-based dementia diagnosis,” _Computational and Mathematical Methods in Medicine_ , 2016. [Online]. Available: https://doi.org/10.1155/2016/9523849
  * [11] N. El-Karoui, “The spectrum of kernel random matrices,” _The Annals of Statistics_ , vol. 38, no. 1, pp. 1–50, 2010.
  * [12] R. Couillet and F. Benaych-Georges, “Kernel Spectral Clustering of Large Dimensional Data,” _Electronic Journal of Statistics_ , vol. 10, pp. 1393–1454, 2016.
  * [13] G. Kimeldorf and G. Wahba, “Some results on Tchebycheffian spline functions,” _Journal of mathematical analysis and applications_ , vol. 33, no. 1, p. 82–95, 1971.
  * [14] B. Scholkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” _Springer_ , no. 1, p. 416–426, 2001.
  * [15] V. Cherkassky, X. Shao, F. M. Mulier, and V. N. Vapnik, “Model complexity control for regression using VC generalization bounds,” _IEEE Transactions on Neural Networks_ , vol. 10, no. 5, pp. 1075–1089, Sep 1999.
  * [16] T. Jaakkola, “Course materials for 6.867 machine learning, fall 2006. mit opencourseware,” _Massachusetts Institute of Technology._ , 2006.
  * [17] G. Valentini and T. G. Dietterich, “Bias-variance analysis of support vector machines for the development of svm-based ensemble methods,” _J. Mach. Learn. Res._ , vol. 5, pp. 725–775, Dec. 2004. [Online]. Available: http://dl.acm.org/citation.cfm?id=1005332.1016783
  * [18] Z. Liao and R. Couillet, “A Large Dimensional Analysis of Least Squares Support Vector Machines,” 2017. [Online]. Available: https://arxiv.org/abs/1701.02967
  * [19] R. Couillet and M. Debbah, _Random Matrix Methods for Wireless Communications_. Cambridge University Press, 2011.
  * [20] Z. Bai and J. W. Silverstein, _Spectral Analysis of Large Dimensional Random Matrices_. Springer, 2009.
  * [21] M. A. Redmond and A. Baveja, “Communities and crime data set,” _UCI Repository_ , 2002. [Online]. Available: http://archive.ics.uci.edu/ml/datasets/communities+and+crime
  * [22] W. Hachem, O. Khorunzhiy, P. Loubaton, J. Najim, and L. Pastur, “A New Approach for Mutual Information Analysis of Large Dimensional Multi-Antenna Channels,” _IEEE Transactions on Information Theory_ , vol. 54, no. 9, pp. 3987–4004, Sept 2008.
  * [23] W. Hachem, P. Loubaton, J. Najim, and P. Vallet, “On bilinear forms based on the resolvent of large random matrices,” _Annales de l’Institut Henri Poincaré- Probabilités et Statistiques_ , vol. 49, no. 1, pp. 36–63, 2013\.
  * [24] W. Hachem, M. Kharouf, J. Najim, and J. W. Silverstein, “A clt for information-theoretic statistics of non-centered gram ranodm matrices,” _Random Matrices, Theory and Applications_ , vol. 12, no. 2, 2012.