# On High-Dimensional Gaussian Comparisons For Cross-Validation

Nicholas Kissel  Carnegie Mellon University  Jing Lei  Carnegie Mellon
University

###### Abstract

We derive high-dimensional Gaussian comparison results for the standard
$V$-fold cross-validated risk estimates. Our result combines a recent
stability-based argument for the low-dimensional central limit theorem of
cross-validation with the high-dimensional Gaussian comparison framework for
sums of independent random variables. These results give new insights into the
joint sampling distribution of cross-validated risks in the context of model
comparison and tuning parameter selection, where the number of candidate
models and tuning parameters can be larger than the fitting sample size. As a
consequence, our results provide theoretical support for a recent
methodological development that constructs model confidence sets using cross-
validation.

##  1 Introduction

Cross-validation (Stone, 1974; Allen, 1974; Geisser, 1975) is among the most
popular procedures for estimating the out-of-sample predictive performance of
statistical models fitted on data sets randomly sampled from a population.
Generally speaking, cross-validation estimates the out-of-sample prediction
accuracy by fitting and assessing a fitted model on separate subsets of data.
One of the most common forms of cross-validation is $V$-fold cross-validation,
where data are partitioned into $V$ folds (sets) of identical size; then, each
fold is used to assess the error of the model fitted using the other $V-1$
folds. Finally, the average of all $V$ estimates is used to create the cross-
validation risk estimate.

Cross-validation is commonly used in statistical learning problems wherein
researchers either compare the cross-validated risk of multiple models or
compare a cross-validated risk against some baseline method with known risk.
See Picard and Cook (1984); Arlot and Celisse (2010) for examples. The
popularity and simplicity of cross-validation has inspired numerous research
articles seeking to better understand its theoretical properties. In
particular, positive results have been established for parameter estimation
following the model selected by cross-validation, including risk consistency
and parameter estimation consistency. See Stone (1977b); Homrighausen and
McDonald (2017); Chetverikov et al. (2016); Celisse (2014) for various
examples from linear regression to nonparametric density estimation problems.

Despite the consistency results established for parameter estimation,
understanding the model selection properties of cross-validation has been a
challenging task, with most existing results being negative. In the early work
of Stone (1977a), it is shown that cross-validation is similar to AIC, and
hence prone to choosing overfitted linear regression models. Such an
overfitting tendency of cross-validation has been further studied in Shao
(1993); Zhang (1993); Yang (2007), which show that in the classical regime,
cross-validation often produces inconsistent model selection unless a very
unrealistic train-validate ratio is used. Indeed, these ratios are so extreme
that they can never be satisfied by standard V-fold cross-validation.
Furthermore, the unsatisfactory model selection performance of cross-
validation has been widely observed in practice, and many heuristic or
context-specific adjustments have been proposed, such as Efron and Tibshirani
(1997); Tibshirani and Tibshirani (2009); Yu and Feng (2014).

The model selection inconsistency of cross-validation can be understood as an
instance of the “winner’s curse.” Since the cross-validated risk of each model
is still a random variable, a particular model may have the smallest cross-
validated risk because its realized random fluctuation happens to be small
while the true optimal model has a much larger fluctuation. Such an intuition
calls for a more precise understanding of the sampling distribution of cross-
validated risks. A main challenge in studying the sampling distribution of
cross-validated risk is the global and heterogeneous dependence among each
individual empirical loss function. Bousquet and Elisseeff (2002) proved
convergence of cross-validated risk to the corresponding population quantity
under an expected leave-one-out loss stability condition. The population
target of cross-validated risk and its variability is further studied in Bates
et al. (2021).

In this work, we study the simultaneous fluctuations of the cross-validated
risks of many models around their mean values. In particular, we establish
high-dimensional Gaussian comparison results for the cross-validated risk
vector indexed by a collection of models, whose cardinality can potentially be
very large. Our main contributions are two fold. First, we extend the low-
dimensional central limit theorem by Austern and Zhou (2020) to the high-
dimensional case, combining their cross-validation error analysis with the
high-dimensional Gaussian comparison framework by Chernozhukov et al. (2013).
Second, we provide theoretically justifiable model selection confidence sets
using cross-validation, answering an open question left in the methodological
work Lei (2020).

Our theoretical development extends and merges two lines of current research:
central limit theorems for cross-validation and high-dimensional Gaussian
comparisons. Low-dimensional central limit theorems for cross-validation have
been developed recently by Austern and Zhou (2020) and Bayle et al. (2020).
While these low-dimensional results provide very useful insights for the
estimation of prediction risks of individual models, they cannot be used to
construct simultaneous confidence sets when many candidate models are being
compared. This is of particular interest because in practice cross-validation
is often used to compare and select from a large collection of models or
tuning parameters. Therefore, in order to understand the behavior of cross-
validation in selecting from many models, it is necessary to consider the
joint sampling distribution of the cross-validated risks. In the case of sums
of independent random vectors, high-dimensional Gaussian comparison has been
developed in the milestone work of Chernozhukov et al. (2013) (see also
Bentkus, 2005). Since then, similar results have been developed for
$U$-statistics (Chen, 2018) and stochastic processes with weak dependence such
as mixing or spatial process (Kurisu et al., 2021; Chang et al., 2021).
However, these extensions do not cover the cross-validation case, where each
term in the summation is dependent of each other with a similar magnitude of
dependence, violating the sparsity of dependence ($U$-statistics) and fast
decaying dependence (mixing and spatial processes). In fact, a different
extension of the Gaussian comparison result is needed for cross-validation,
which borrows the martingale decomposition and stability conditions in Austern
and Zhou (2020).

##  2 Preliminaries

Consider iid data $\mathbf{X}=(X_{0},X_{1},...,X_{n})$ with
$X_{i}\in\mathcal{X}$. We would like to simultaneously study the performance
of $p$ learning algorithms through the framework of $V$-fold cross-validation.

For notation simplicity, we assume $V$ evenly divides $n$. For $v\in[V]$, let
$I_{v}=\\{n(v-1)/V+1,...,nv/V\\}$ be the index that corresponds to the $v$th
fold of data. Let $\tilde{n}=n(1-1/V)$ be the training sample size used in
$V$-fold cross-validation. For each $r\in[p]$, let
$\ell_{r}(\cdot;\cdot):\mathcal{X}\times\mathcal{X}^{\tilde{n}}\mapsto\mathbb{R}$
be a loss function. Intuitively, we should think
$\ell_{r}(x_{0};(x_{1},...,x_{\tilde{n}}))$ as the loss function evaluated at
$x_{0}$ of a fitted model using training data $(x_{1},...,x_{\tilde{n}})$.
Here the index $r$ denotes a particular model or tuning parameter value. This
notation covers both supervised and unsupervised learning.

  1. 1.

In supervised learning, each point can be thought of as
$x=(y,z)\in\mathcal{Y}\times\mathcal{Z}$ where $z$ is a vector of covariates
and $y$ is the response variable. The loss function can often be written more
concretely as

$\ell_{r}(x_{0};(x_{1},...,x_{\tilde{n}}))=\rho(y_{0},\hat{f}_{r}(z_{0}))\,,$

where $\hat{f}_{r}(\cdot):\mathcal{Z}\mapsto\mathcal{Y}$ is a regression
function that predicts $y$ from $z$, trained using the $r$th model/tuning
parameter with input data $(x_{1},...,x_{\tilde{n}})$, and
$\rho(\cdot,\cdot):\mathcal{Y}^{2}\mapsto\mathbb{R}$ is a loss function
measuring the quality of predicting $y$ using $\hat{f}_{r}(z)$, such as
squared loss, $0$-$1$ loss, and hinge loss.

  2. 2.

In unsupervised learning,

$\ell_{r}(x_{0};(x_{1},...,x_{\tilde{n}}))=\rho(x_{0};\hat{f}_{r})\,,$

where $\hat{f}_{r}$ is a function describing the distribution of $X$ trained
from the $r$th model/tuning parameter with the input data
$(x_{1},...,x_{\tilde{n}})$, and the function $\rho$ is a loss function
assessing the agreement of the sample point $x_{0}$ and the fitted probability
model $\hat{f}_{r}$. Examples of $\rho$ include the negative likelihood in
density estimation and the proportion of total variance explained in dimension
reduction.

In model selection and parameter tuning, a particularly interesting scenario
is when the number of models being compared is large.

For each $r$, the $V$-fold cross-validated risk is

$\hat{R}_{{\rm
cv},r}=n^{-1}\sum_{i=1}^{n}\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})\,.$ (1)

where $\mathbf{X}_{-v}$ denotes the sub-vector of $\mathbf{X}$ excluding the
$v$th fold, and $v_{i}\in[V]$ is such that $i\in I_{v_{i}}$.

It is natural to expect $\hat{R}_{{\rm cv},r}$ to approximate the true average
risk of the fitted model:

$\tilde{R}_{r}=\frac{1}{V}\sum_{v=1}^{V}R_{r}(\mathbf{X}_{-v})\,,$

where

$R_{r}(\mathbf{X}_{-v})=\mathbb{E}\left[\ell_{r}(X_{0};\mathbf{X}_{-v})|\mathbf{X}_{-v}\right]\,,$

is the true risk of the $r$th model fitted using input data $\mathbf{X}_{-v}$.

The quantity $R_{r}(\mathbf{X}_{-v})$ still depends on the input data
$\mathbf{X}_{-v}$ and hence is a random variable itself. It would be natural
to consider its expected value

$R_{r}^{*}=\mathbb{E}R_{r}(\mathbf{X}_{-v})=\mathbb{E}\ell_{r}(X_{0};\mathbf{X}_{-1})\,.$

Statistical inference tasks such as model comparison would also involve
uncertainty quantification of such point estimates of risks, and we would hope
to establish central limit theorems of the form

$\frac{\sqrt{n}(\hat{R}_{{\rm cv},r}-\mu_{r})}{\sigma_{r}}\rightsquigarrow
N(0,1)$

with $\mu_{r}$ being either $\tilde{R}_{r}$ or $R_{r}^{*}$ and some
appropriate scaling $\sigma_{r}$ (Austern and Zhou, 2020; Bayle et al., 2020).

In the context of model comparison or tuning parameter selection, such
individual normal approximation would have limited practical use. For example,
in our numerical example in Section 4.1, individual confidence intervals fail
to simultaneously cover the targets when when $p$ is moderately large. To
cover this gap between theory and practice, we seek to establish a high-
dimensional Gaussian approximation in a similar fashion as in Chernozhukov et
al. (2013):

$\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\max_{1\leq r\leq
p}\sqrt{n}(\hat{R}_{{\rm cv},r}-\mu_{r})\leq
z\right)-\mathbb{P}\left(\max_{1\leq r\leq p}Y_{r}\leq
z\right)\right|\rightarrow 0$ (2)

for some centered Gaussian random vector $\mathbf{Y}=(Y_{1},...,Y_{p})$ with
matching covariance.

##  3 Main results

In this section, we establish a high-dimensional Gaussian approximation result
with random centering. In particular, we prove (2) with
$\mu_{r}=\tilde{R}_{r}$. In the following subsections, we present and discuss
the assumptions required for this result and provide its full statement as a
theorem in Section 3.3.

###  3.1 Symmetry and moment conditions on the loss function

The idea of cross-validation relies on independence and symmetry among data
points. We consider the following symmetry and moment conditions on the loss
functions $\ell_{r}$.

######  Assumption 1 (Symmetry and moment condition on $\ell_{r}$).

For each $r\in[p]$, the loss function $\ell_{r}(\cdot;\cdot)$ satisfies

  1. (a)

$\ell_{r}(x_{0};x_{1},...,x_{\tilde{n}})$ is symmetric in
$(x_{1},...,x_{\tilde{n}})$.

  2. (b)

$\mathbb{E}[\ell_{r}(X_{1};\mathbf{X}_{-1})-R_{r}(\mathbf{X}_{-1})]^{2}\geq\underline{\sigma}^{2}$
for some constant $0<\underline{\sigma}$.

Part (b) essentially assumes that the randomly centered cross-validated loss
function has non-degenerate conditional variance. This makes intuitive sense,
as we would expect the resulting confidence interval to have length at the
scale of $1/\sqrt{n}$. For example, if $\ell_{r}$ is a regression residual,
then this lower bound is at least as large as the prediction risk of the ideal
regression function. In the additional assumptions below, we will also have
the upper bound on the variance term.

###  3.2 Stability and tail conditions

A key consideration from Austern and Zhou (2020) in their low dimensional
central limit theorems for cross-validation is the stability of the loss
function and the average risk when one input sample point is replaced by an
iid copy. In the high dimensional case, we need the loss function to be stable
in a uniform sense across all $p$ candidate models indexed by $r\in[p]$. Thus,
we will consider stability conditions in the form of stronger tail
inequalities instead of the moment conditions used for the low dimensional
case. Such stronger tail conditions are common in high dimensional central
limit theorem literature, such as in Chernozhukov et al. (2013).

We use sub-Weibull concentration to describe the required tail behaviors of
random variables.

######  Definition 1 (Sub-Weibull Random Variables).

Let $K$ be a positive number, we say a random variable $X$ is $K$-sub-Weibull
($K$-SW for short) if there are positive constants $(a,b,\alpha)$ such that

$\mathbb{P}\left(\frac{|X|}{K}\geq t\right)\leq
ae^{-bt^{\alpha}}\,,~{}~{}\forall~{}t>0\,.$

This definition generalizes the well-known sub-exponential and sub-Gaussian
distributions, and has been systematically introduced in Vladimirova et al.
(2020); Kuchibhotla and Chakrabortty (2018).

######  Remark 1.

Unlike common practices in the literature, our notation of the sub-Weibull
tail inequality only focuses on the scaling $K$. We do not keep track of the
constants $a,b,\alpha$, which can vary from one instance to another as long as
they stay bounded and bounded away from zero. It is easy to check that our
notion of sub-Weibull is invariant under constant scaling: if $X$ is $1$-SW,
then $X$ is $c$-SW for all positive constant $c$. Aside from the scaling $K$,
the second (and only) important parameter in sub-Weibull tail inequality is
the exponent $\alpha$. In the literature, it is more common to write
$(K,\alpha)$-SW. Our proof can be adapted to keep explicit track of the
constant $\alpha$ in each instance at the cost of more complicated
bookkeeping, but that does not qualitatively change the results.

In our theoretical developments, the dependence on logarithm terms may be
complicated, as it involves the sub-Weibull constant $\alpha$, which may vary
between lines. For brevity of presentation, we absorb such logarithm terms
into the $\tilde{O}(\cdot)$ notation. Where $A\leq\tilde{O}(B)$ means that
there are positive constants $c_{1},c_{2}$ independent of $(n,p)$ such that
$A\leq c_{1}\log^{c_{2}}(n+p)B$.

To introduce the stability conditions, let $X_{i}^{\prime}$ be iid copies of
$X_{i}$ for $1\leq i\leq n$ and $\mathbf{X}^{i}$ be the vector obtained by
replacing $X_{i}$ with $X_{i}^{\prime}$. For any function $f(\mathbf{X})$,
define $\nabla_{i}f(\mathbf{X})=f(\mathbf{X})-f(\mathbf{X}^{i})$.

The main stability conditions involved in our normal approximation bounds are
the following.

######  Assumption 2.

There exists $\epsilon_{\ell}\in(0,1]$ such that

  1. (a)

For all $i\in[\tilde{n}]$, $r\in[p]$,
$\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$ is $n^{-1/2}\epsilon_{\ell}$-sub-
Weibull.

  2. (b)

For all $1\leq i<j\leq\tilde{n}$, $r\in[p]$,
$\nabla_{j}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$ is
$n^{-3/2}\epsilon_{\ell}$-sub-Weibull.

  3. (c)

For all $r$, $\ell_{r}(X_{0},\mathbf{X}_{-1})$ is $1$-sub-Weibull.

Assumption 2 requires that the first order difference
$\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$ has a scaling no larger than
$\epsilon_{\ell}n^{-1/2}$, the second order difference
$\nabla_{j}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$ has a scaling no larger
than $\epsilon_{\ell}n^{-3/2}$, and the original loss function $\ell_{r}$ has
a constant scaling. The sub-Weibull tail ensures that with high probability
all such quantities will not exceed their scalings by more than a poly-
logarithm factor. We require $\epsilon_{\ell}\leq 1$ for notation simplicity
and without much loss of generality, because the approximation error bound in
the main theorems becomes meaningless if $\epsilon_{\ell}>1$.

We further remark that the scaling assumption on the second order difference
$\nabla_{j}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$ is stronger than that in
Austern and Zhou (2020) by a factor of $\sqrt{n}$. This is due to a
fundamental difference between low dimensional CLT and high dimensional
Gaussian comparison, where the former only requires controlling the second
moment of error terms, while the latter requires controlling the supremum of
many such error terms. More specifically, define the randomly centered loss at
$X_{i}$

$K_{r,i}=\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})-R_{r}(\mathbf{X}_{-v_{i}})\,,$
(3)

and

$D_{r,i}=\sum_{j\notin I_{v_{i}}}\nabla_{i}K_{r,i}\,.$ (4)

A key result in the low dimensional CLT is that
$\|D_{r,i}\|_{2}\lesssim\epsilon_{\ell}$ provided
$\|n^{1/2}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})\|_{2}\leq\epsilon_{\ell}$
and
$n\|\nabla_{j}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})\|_{2}\leq\epsilon_{\ell}$.
However, in the high-dimensional regime, we need to simultaneously control
$D_{r,i}$ for all $1\leq r\leq p$, which cannot be guaranteed by a vanishing
second moment on each individual term. Our condition can be relaxed to
requiring a similar $\nabla_{j}\nabla_{i}\ell_{r}(X_{0},\mathbf{X}_{-1})$
being $n^{-1}\epsilon_{\ell}$-sub-Weibull, provided we can further assume that
$\frac{D_{r,i}}{\|D_{r,i}\|_{2}}$ is $1$-sub-Weibull. While this additional
assumption certainly seems reasonable in many situations, we choose to work
with the stronger condition on the second order difference as presented in
Assumption 2, which allows for a more streamlined presentation. Nevertheless,
the stability conditions in Assumption 2 are still practically plausible since
we should typically expect each $\nabla$ operator to reduce the scale of the
loss function by a factor up to $n$.

Here we provide two examples with simple intuition about the plausibility of
the stability conditions required in Assumption 2. Rigorous derivations seem
possible by invoking convexity and smoothness in M-estimation (Example 1, see
also Proposition 4 of Austern and Zhou (2020)), and random matrix theory
(Example 2). Another example satisfying such stability conditions is bagged
estimates (Chen et al., 2022).

######  Example 1 (Classical M-Estimator).

Consider a parametric loss function
$\ell(x_{0};x_{1},...,x_{\tilde{n}})=\ell(x_{0};\hat{\theta})$ with
$\hat{\theta}$ a parameter estimated from the input data
$(x_{1},...,x_{\tilde{n}})$. Under classical parametric regularity conditions
such $\hat{\theta}$ can satisfy the stability required in Assumption 2. For
example, if $\hat{\theta}$ is asymptotically linear (Tsiatis, 2006, Chapter
3), then we have

$\hat{\theta}=\theta_{0}+\frac{1}{\tilde{n}}\sum_{i=1}^{\tilde{n}}\varphi(x_{i})+o_{P}(n^{-1/2})\,.$

Then, the stability bound for $\nabla_{i}\ell$ is directly satisfied. The
second order stability condition seems plausible if the $o_{P}(n^{-1/2})$ has
weak dependence on each single data point.

######  Example 2 (Penalized Least Squares).

Now consider a high dimensional ridge-regression where we have paired sample
points $x_{i}=(z_{i},y_{i})$:

$\hat{\beta}=\arg\min{\tilde{n}}^{-1}\sum_{i=1}^{\tilde{n}}(y_{i}-z_{i}^{T}\beta)^{2}+\lambda\|\beta\|_{2}^{2}\,.$

When the dimensionality of $x_{i}$ is comparable or smaller than the sample
size, it is possible to argue that the empirical covariance matrix will be
well-conditioned with high probability, and hence changing any one sample
point will incur an $O(n^{-1})$ change in $\hat{\beta}$. For a simpler
argument under stronger assumptions, if the sample points are bounded and
$\lambda\gg n^{-1/2}$, then the stability requirement on $\nabla_{i}\ell$
holds. If $\lambda\gg n^{-3/4}$ then the stability requirement on
$\nabla_{j}\nabla_{i}\ell$ also holds.

###  3.3 Main theorem with random centering

Our first main result is a Gaussian comparison with random centering. In order
to state the result, we need to specify the covariance of the Gaussian vector
$\mathbf{Y}$. Using the notation
$K_{r,0,v}=\ell_{r}(X_{0};\mathbf{X}_{-v})-R_{r}(\mathbf{X}_{-v})$, define

$\sigma_{rs,v}=\mathbb{E}(K_{r,0,v}K_{s,0,v}|\mathbf{X}_{-v})\,,~{}~{}1\leq
r,s\leq p\,,$

to be the conditional variance/covariance of the loss functions given the
fitted model using input data $\mathbf{X}_{-v}$, and let

$\sigma_{rs}=\mathbb{E}\sigma_{rs,1}$

be the expected value of $\sigma_{rs,v}$. Let $\Sigma=[\sigma_{rs}]_{1\leq
r,s\leq p}$ be the corresponding $p\times p$ expected conditional covariance
matrix, and $\Sigma_{v}$ the corresponding random covariance matrix with
entries $\sigma_{rs,v}$.

We will show that $\sigma_{rs,v}\approx\sigma_{rs}$ (Lemma B.3) and
$\sqrt{n}(\hat{\mathbf{R}}_{\rm cv}-\tilde{\mathbf{R}})$ behaves like a
centered Gaussian vector with covariance matrix $\Sigma$.

######  Theorem 3.1 (High-dimensional CLT for Cross-validation with random
centering).

Assume Assumptions 1 and 2 hold, then we have

$\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\max_{1\leq r\leq
p}\sqrt{n}(\hat{R}_{{\rm cv},r}-\tilde{R}_{r})\leq
z\right)-\mathbb{P}\left(\max_{1\leq r\leq p}Y_{i}\leq
z\right)\right|\leq\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,,$ (5)

for $\mathbf{Y}=(Y_{1},...,Y_{p})\sim N(0,\Sigma)$.

######  Remark 2.

Theorem 3.1 implies that the Gaussian approximation error is small if
$\epsilon_{\ell}\lesssim n^{-c}$ for some constant $c>0$. The result of
Theorem 3.1 can be easily extended to the quantity $\max_{r}|\hat{R}_{{\rm
cv},r}-\tilde{R}_{r}|$ by applying Theorem 3.1 to the augmented vector
$(\hat{\mathbf{R}}_{\rm
cv}-\tilde{\mathbf{R}},\tilde{\mathbf{R}}-\hat{\mathbf{R}}_{\rm cv})$ with the
corresponding Gaussian vector $(\mathbf{Y},-\mathbf{Y})$.

##  4 Simultaneous Confidence Bands for Cross-Validated Risk

In this section we consider various statistical inference tools following from
Theorem 3.1, including constructing simultaneous confidence bands of the
average fitted risks and possible ways to construct confidence sets of the
“optimal” model.

###  4.1 Confidence bands

Following Theorem 3.1, we consider the coordinate-wise standardized process

$\sqrt{n}\Lambda^{-1/2}(\hat{\mathbf{R}}_{\rm cv}-\tilde{\mathbf{R}})\,,$

where

$\displaystyle\Lambda=$ $\displaystyle{\rm diag}(\sigma_{11},...,\sigma_{pp})$

is the diagonal submatrix of $\Sigma$.

Let $\hat{\Lambda}$ and $\hat{\Sigma}$ be the natural plug-in estimates (i.e.,
the average of all the within-fold empirical covariance matrices) of $\Lambda$
and $\Sigma$. In particular, let $\bm{\ell}=(\ell_{1},...,\ell_{p})$, and
$\hat{\Sigma}_{v}$ be the empirical covariance of
$\left\\{\bm{\ell}(X_{i};\mathbf{X}_{-v}):i\in I_{v}\right\\}$. Then
$\hat{\Lambda}$ and $\hat{\Sigma}$ can be the aggregated estimate.

$\displaystyle\hat{\Sigma}=$
$\displaystyle\frac{1}{V}\sum_{v=1}^{V}\hat{\Sigma}_{v}\,,~{}~{}~{}\hat{\Lambda}={\rm
diag}(\hat{\sigma}_{11},\dots,\hat{\sigma}_{pp})\,.$ (6)

Given a nominal type I error level $\alpha\in(0,1)$, the following fully data-
driven procedure computes a simultaneous normalized confidence band for the
vector $\tilde{\mathbf{R}}$ with asymptotic coverage $1-\alpha$.

Let $\hat{z}_{\alpha}$ be the upper $1-\alpha$ quantile of $\|Z\|_{\infty}$,
with $Z\sim N(0,\hat{\Lambda}^{-1/2}\hat{\Sigma}\hat{\Lambda}^{-1/2})$. Given
estimates $(\hat{\Lambda},\hat{\Sigma})$, $\hat{z}_{\alpha}$ can be
approximated efficiently using Monte-Carlo methods. Such a Monte-Carlo
approximation error can be controlled by combining the standard Dvoretsky-
Kiefer-Wolfowitz inequality and anti-concentration of Gaussian maxima. So we
use the theoretical value $\hat{z}_{\alpha}$ for brevity, which corresponds to
the limiting case of infinite Monte-Carlo sample size.

The simultaneous confidence band for cross-validated risk is

$\displaystyle\widehat{\rm CI}_{r}=\left[\hat{R}_{{\rm
cv},r}-\frac{\hat{\sigma}_{rr}^{1/2}\hat{z}_{\alpha}}{\sqrt{n}}\,,~{}\hat{R}_{{\rm
cv},r}+\frac{\hat{\sigma}_{rr}^{1/2}\hat{z}_{\alpha}}{\sqrt{n}}\right]\,,~{}~{}\text{for
each }r\in[p]\,,$ (7)

where $\hat{\sigma}_{rr}$ is the $r$th diagonal entry of $\hat{\Sigma}$ in
(6).

######  Corollary 4.1.

Under Assumptions 1 and 2, the confidence intervals constructed in (7) satisfy

$\mathbb{P}\left(\tilde{R}_{r}\in\widehat{\rm
CI}_{r}\,,~{}\forall~{}r\in[p]\right)\geq
1-\alpha-\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,.$

###  4.2 Model confidence set

Now, we consider the model/tuning selection problem. Let

$r^{*}=\arg\min_{r}\tilde{R}_{r}$

be the index of the candidate model with the smallest average fitted risk. We
hope to use the Gaussian comparison to construct a confidence set of $r^{*}$.
A simple way to do so is directly using (4.1):

$\hat{\mathcal{A}}_{0}=\left\\{r:\hat{R}_{{\rm
cv},r}-\hat{\sigma}_{rr}^{1/2}\hat{z}_{\alpha}/\sqrt{n}\leq\min_{s}\hat{R}_{{\rm
cv},s}+\hat{\sigma}_{ss}^{1/2}\hat{z}_{\alpha}/\sqrt{n}\right\\}\,.$ (8)

It is a direct consequence of Corollary 4.1 that

$\mathbb{P}(r^{*}\in\hat{\mathcal{A}}_{0})\geq
1-\alpha-\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,.$

However, $\hat{\mathcal{A}}_{0}$ often unnecessarily contains too many models
as it ignores the correlations among the coordinates of $\hat{\mathbf{R}}_{\rm
cv}$.

Following the idea in Lei (2020), we instead consider the following difference
based method, which takes into account the correlations of the cross-validated
risks. For each $r$, consider the risk difference vector $(\hat{R}_{{\rm
cv},r}-\hat{R}_{{\rm cv},s}:s\neq r)$, and apply the above framework to
$\ell^{(r)}_{s}\coloneqq\ell_{r}-\ell_{s}$ to test whether $\hat{R}_{{\rm
cv},r}-\hat{R}_{{\rm cv},s}>0$ for some $s\neq r$. Here we are considering a
one-sided hypothesis, so instead of the two-sided confidence band in 7, we
consider the one-sided version.

For each $r\in[p]$, consider $p-1$ difference loss functions
$\ell^{(r)}_{s}=\ell_{r}-\ell_{s}$ for $s\in[p]\backslash\\{r\\}$. Now apply
the cross-validation normal approximation theory to the $p-1$ loss functions
$(\ell^{(r)}_{s}:1\leq s\leq p,~{}s\neq r)$. Let $\hat{z}^{(r)}_{\alpha}$ be
the $1-\alpha$ quantile of the maximum of the corresponding
$(p-1)$-dimensional Gaussian vector with estimated covariance
$[\hat{\Lambda}^{(r)}]^{-1/2}\hat{\Sigma}^{(r)}[\hat{\Lambda}^{(r)}]^{-1/2}$,
where $\hat{\Sigma}^{(r)}$ and $\hat{\Lambda}^{(r)}$ are counterparts of
$\hat{\Sigma}$ and $\hat{\Lambda}$ when applied to the $(p-1)$ difference loss
functions $(\ell_{r}-\ell_{s}:s\in[p]\backslash\\{r\\})$. Then our model
confidence set is

$\hat{\mathcal{A}}=\left\\{r:\sup_{s\neq r}\hat{R}_{{\rm cv},r}-\hat{R}_{{\rm
cv},s}-[\hat{\sigma}^{(r)}_{ss}]^{1/2}\hat{z}^{(r)}_{\alpha}/\sqrt{n}\leq
0\right\\}\,.$ (9)

######  Proposition 4.2.

If the difference loss functions $\ell^{(r)}_{s}=\ell_{r}-\ell_{s}$ satisfy
Assumption 2 for all $r,s$, then

$\mathbb{P}\left(r^{*}\in\hat{\mathcal{A}}\right)\geq
1-\alpha-\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,.$

######  Remark 3.

Proposition 4.2 resolves an outstanding question about the theoretical
justification of the V-fold cross-validation with confidence (CVC) method Lei
(2020). Such a theoretical guarantee requires the difference loss functions to
satisfy Assumption 2, which is non-trivial. Because if the loss functions
$\ell_{r}$ and $\ell_{s}$ are highly correlated, their difference can be
singular and hence violate the stability requirement. In particular even if
$\nabla_{i}\ell_{r}\asymp n^{-1}$, for the stability condition to hold for the
difference loss we will need $\|\ell_{r}-\ell_{s}\|_{2}\gg n^{-1/2}$. In fact,
if the variance of $\ell_{r}-\ell_{s}$ vanishes sufficiently slowly , then it
is straightforward to modify the proof of Theorem 3.1 to accommodate this
scenario. Such a slow vanishing requirement on $\ell_{r}-\ell_{s}$ precludes
the case that both model $r$ and model $s$ produces $\sqrt{n}$-consistent
estimates. This intuition agrees with Yang (2007), which suggests that cross-
validation may not be model selection consistent if both candidate models are
$\sqrt{n}$-consistent. A simple illustration of this issue is given in Section
2.3 of Lei (2020).

###  4.3 Numerical Experiments

In this subsection we numerically verify the claim of Theorem 3.1 as well as
the model confidence sets considered in Section 4.1. We use $V=5$ in all
simulations and all plotted values are averaged over 1000 generated data sets.

#### Simultaneous coverage vs marginal coverage.

We first investigate the simultaneous coverage of confidence bands for the
cross-validated risk. To do so, we generate a predictor matrix $X\sim N(0,I)$
and response vector $Y=X\beta+\epsilon$ where $\epsilon\sim N(0,\sigma^{2}I)$
with $\sigma^{2}=\|\beta\|_{2}^{2}/\nu$ and $\beta$ is a sparse
$d$-dimensional vector with the first $s$ entries being 1 and the remainder
being 0. We set $d=20$, $s=5$, and $\nu=1000$. We then fit lasso regressions
across a grid of 50 regularization parameters and generate confidence bands.

Figure 1:  Simultaneous coverage of point-wise confidence intervals and the
confidence band from 7. The horizontal dashed lines represent the nominal
level given by $1-\alpha$.

Figure 1 shows the simultaneous coverage of a confidence band generated by all
point-wise confidence intervals (left) and the confidence band as specified in
7 at various values of $n$ and $\alpha$. We see that the latter method has
coverage much closer to the nominal level than the former. Therefore, the
point-wise procedure is insufficient for providing the correct simultaneous
coverage, suggesting that the simultaneous adjustment is indeed necessary.
This simulation also suggests that the coverage of the simultaneous band is
not overly conservative.

Figure 2:  Coverage of CV risk confidence intervals for forward selection
terminated at different values (left) and lasso at different regularization
values (right). The dashed horizontal line marks the nominal level given by
$1-\alpha$.

#### The importance of stability.

The impact of stability on the quality of Gaussian approximation of cross-
validated risks has been experimented in Austern and Zhou (2020). Here we
provide additional empirical results on this front. Consider confidence
intervals for forward selection in the same setting but with $d=10$ and
$\alpha=0.05$. Specifically, we look at the point-wise interval coverage of
forward selection terminated at different model sizes–one less than $s$, one
equal to $s$, and one larger than $s$. The left plot in Figure 2 shows that at
3 steps, the one-dimensional $95\%$-confidence interval based on Theorem 3.1
under covers regardless of the sample size, while at 5 and 7 steps the
coverage converges to nominal level as the sample size increases. This
observation is consistent with the stability condition. Before $s$ is reached,
forward selection is quite unstable in this setting, as the non-zero entries
of $\beta$ have the same magnitude. Therefore, the algorithm is equally likely
to pick any subset of $s-1$ non-zero coordinates, and changing the value of
one data point can incur a change of selected variables with non-negligible
chance, resulting in instability of the loss function. When forward selection
reaches exactly $s$ steps, it selects the correct subset with overwhelming
probability, leading to a stable loss function. When forward selection selects
one more variable, the index of the additional selected variable is not stable
but the fitted coefficient is very close to zero, so that the fitted loss
function is still stable. In contrast, the lasso algorithm is continuous in
the input data, hence the stability is much easier to hold for all values of
penalty parameters, as suggested by the right plot in Figure 2.

#### The advantage of difference-based model confidence sets.

Now we look at the model confidence set performance as applied to the lasso on
data with increasing $n$. Specifically, we are studying the size
$|\mathcal{A}|$ and coverage $\mathbb{P}(r^{*}\in{\mathcal{A}})$ for

  1. (1)

the naïve method $\mathcal{A}=\hat{\mathcal{A}}_{0}$ as defined in (8), and

  2. (2)

the difference based method $\mathcal{A}=\hat{\mathcal{A}}$ as defined in (9).

This simulation setting is again similar with the value of $s$ fixed at 5, but
$d$ grows at rate $n/10$. This time, to make our grid, we first find
$\lambda_{\text{max}}=\frac{1}{n}\|X^{T}Y\|_{\infty}$, then the grid is
$\lambda_{\text{max}}2^{i}/\sqrt{1-1/V}$ for $i\in\\{0,\dots,-9\\}$. The re-
scaling of $\sqrt{1-1/V}$ is done since the choice of $\lambda$ in lasso is
inversely proportional to the square root of the training sample size, and
$\lambda_{\max}$ may be a bit too small for the reduced sample size in
$V$-fold cross-validation. The left plot of Figure 3 shows that the difference
based method produces considerably smaller sets while maintaining coverage of
at least 0.95, supporting the intuition that the difference based method is
able to take into account the joint randomness of the cross-validated risks.
On the right plot of Figure 3, the empirical coverage of the naïve method is
always overly conservative, while the coverage of the difference based method
does get close to nominal level for certain values of $n$. Intuitively, such a
fluctuation of coverage as $n$ varies can be explained by whether the problem
is close to the boundary of the null hypothesis. More specifically, let
$\bm{\delta}_{r}=(\tilde{R}_{r}-\tilde{R}_{s}:s\neq r)$. In the difference
based method, the null hypothesis for a candidate $r$ is that
$\bm{\delta}_{r}$ is non-positive in each coordinate, which is a composite
null hypothesis. The Gaussian approximation is derived precisely for the
extreme point of the null hypothesis, where all coordinates of
$\bm{\delta}_{r}$ are zero. In practice, we will never really be working in
this scenario. Therefore, the supremum-based confidence set will be too
conservative if $\bm{\delta}_{r}$ has many large negative coordinates, and
will be nearly exact if most of the coordinates are close to $0$. In our
simulation, when $n$ increases, the relative performance of different tuning
parameters also changes. Indeed we do observe that when $n\leq 2500$, the
small sample size cannot quite distinguish the two best $\lambda$ values that
perform nearly equally well.

Figure 3:  Size and coverage of model confidence set procedures with
$\alpha=0.05$.

##  5 On deterministic centering

So far, we have focused on the high dimensional Gaussian comparison of cross-
validated risk with random centering, where the mean vector
$\tilde{\mathbf{R}}$ is data-dependent. This leads to the following question:
can we establish Gaussian comparison results with fixed centering? It is
natural to expect the fixed centering to be
$\mathbf{R}^{*}=\mathbb{E}(\hat{\mathbf{R}}_{\rm cv})$. Also, the
corresponding scaling should be based on the total variance
$\mathbb{E}\left[\bm{\ell}(X_{0};\mathbf{X}_{-1})\bm{\ell}(X_{0};\mathbf{X}_{-1})^{T}\right]-\mathbf{R}^{*}{{\mathbf{R}}^{*}}^{T}$,
which in addition to the variance term $\Sigma$ considered in the random
centering case above, must also take into account the variability caused by
the randomness of $\mathbf{R}(\mathbf{X}_{-1})$. Such a fixed centering
central limit theorem for cross-validated risks has been studied in Austern
and Zhou (2020) in the low-dimensional case. Our development here extends
their result to the high-dimensional case with a more streamlined proof.

###  5.1 Risk stability

In order to study the randomness in the risk function
$R_{r}(\mathbf{X}_{-1})$, we need the following stability conditions on the
risk functions.

######  Assumption 3 (Risk stability).

There exists a constant $\epsilon_{R}\in\mathbb{R}^{+}$ such that
$n\nabla_{i}R_{r}(\mathbf{X}_{-1})$ is $\epsilon_{R}$-SW for all $i$.

A remarkable difference between the risk stability and loss stability is in
the scaling factors. For the first order differencing operators $\nabla_{i}$,
the scaling factor of $n$ instead of $\sqrt{n}$ makes the risk stability
apparently harder to control than the loss stability. It is also considered
and briefly discussed in Austern and Zhou (2020) in the low-dimensional case.
Here we give more explanation of this key condition. At first, it seems
unreasonable to assume that $\epsilon_{R}$ is very small, as
$\nabla_{i}R_{r}(\mathbf{X}_{-1})$ generally should not be smaller than $1/n$.
However, a closer inspection suggests that the risk function (taking
expectation of the loss function $\ell_{r}(X_{0},\mathbf{X}_{-1})$ over the
evaluating point $X_{0}$) is usually much more stable than the loss function
itself, as taking conditional expectation usually increases stability. In fact
such an increase of stability can be quite substantial. For example, assume
that the loss function takes a parametric form:
$\ell(x_{0};x_{1},...,x_{\tilde{n}})=\ell(x_{0};\hat{\theta})$ where
$\hat{\theta}=\hat{\theta}(x_{1},...,x_{\tilde{n}})$ is a fitted parameter
from the input data $(x_{1},...,x_{\tilde{n}})$. Then

$\nabla_{i}R(\hat{\theta})\approx\left(\frac{dR}{d\theta}\bigg{|}_{\hat{\theta}}\right)\nabla_{i}\hat{\theta}\,,$

which should be much smaller than $\nabla_{i}\hat{\theta}$ if
$R(\theta)\coloneqq\mathbb{E}\ell(X_{0};\theta)$ is flat at $\hat{\theta}$.
This is usually the case when $\hat{\theta}$ is in a small neighborhood of the
optimal parameter value with minimum risk $\theta^{*}$.

Furthermore, we remark that our theoretical development does not require
$\epsilon_{R}$ to vanish asymptotically. Instead, we only need $\epsilon_{R}$
to be dominated by other vanishing terms such as $\epsilon_{\ell}$ and
$1/\sqrt{n}$.

###  5.2 Gaussian comparison with deterministic centering

Finding the covariance matrix for deterministic centering starts by
identifying the contribution of randomness from each single sample point. We
start by writing

$n\hat{R}_{{\rm cv},r}=\sum_{i=1}^{n}\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})\,.$

The part in the above sum that involves $X_{i}$ is

$\displaystyle\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})+\sum_{j\notin
I_{v_{i}}}\ell_{r}(X_{j};\mathbf{X}_{-v_{j}})\,.$ (10)

It is clear that $X_{i}$ plays two different roles in $\hat{R}_{{\rm cv},r}$:
(i) as the evaluation point in $\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})$, (ii) as
one of the $\tilde{n}$ fitting sample points in each of
$\ell_{r}(X_{j};\mathbf{X}_{-v_{j}})$ for $j\neq I_{v_{i}}$. The randomness
contributed by $X_{i}$ as an evaluating point should be captured by the
variability of the average loss function

$\bar{\ell}_{r}(X_{1})=\mathbb{E}\left[\ell_{r}(X_{1};\mathbf{X}_{-1})|X_{1}\right]\,,$

and the randomness contributed by $X_{1}$ as a fitting sample point should be
captured in the function

$\bar{R}_{r}(X_{1})=\mathbb{E}\left[R(\mathbf{X}_{-V})|X_{1}\right]\,.$

Let $\bm{\Phi}=(\phi_{rs}:1\leq r,s\leq p)$ be the covariance matrix given by

$\displaystyle\phi_{rs}\coloneqq$ $\displaystyle{\rm
Cov}\left[\bar{\ell}_{r}(X_{1})+\tilde{n}\bar{R}_{r}(X_{1})\,,~{}\bar{\ell}_{s}(X_{1})+\tilde{n}\bar{R}_{s}(X_{1})\right]\,~{}~{}1\leq
r,s\leq p\,.$ (11)

We assume that the marginal variance terms are bounded and bounded away from
$0$. The following assumption is analogous to Assumption 1(c).

######  Assumption 4.

There exist positive constants $\underline{\phi}$ and $\bar{\phi}$ such that
$\underline{\phi}\leq\phi_{rr}\leq\bar{\phi}$ for each $r\in[p]$.

######  Theorem 5.1 (Deterministic Gaussian Comparison).

Assume Assumptions 1, 2, 3 and 4 hold, then we have

$\displaystyle\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\sqrt{n}\max(\hat{\mathbf{R}}_{\rm
cv}-\mathbf{R}^{*})\leq z\right)-\mathbb{P}\left(\max\mathbf{Y}\leq
z\right)\right|$
$\displaystyle\leq\tilde{O}\left([\epsilon_{\ell}(1+\epsilon_{R})]^{1/3}+n^{-1/8}(1+\epsilon_{R})^{3/4}\right)$

for $\mathbf{Y}\sim N(0,\bm{\Phi})$.

###  5.3 Deterministic variance estimation

We now consider the problem of estimating $\phi_{rs}$. This problem has been
considered in Austern and Zhou (2020). We believe the estiamte stated in their
text is off by a factor of 2, and also only covers the case of two-fold cross-
validation. Our result below corrects the scaling and covers the general
$V$-fold case in a multivariate setting.

As suggested in (11) and Theorem 5.1, the covariance $\phi_{rs}$ is
essentially the sum of the marginal variability of each $X_{i}$. Indeed, we
have the following result

######  Theorem 5.2 (Marginal variance approximation).

Under Assumptions 1, 2, 3 and 4, we have

$\frac{n^{2}}{2}\mathbb{E}\left[\nabla_{i}\hat{R}_{{\rm
cv},r}\nabla_{i}\hat{R}_{{\rm cv},s}|\mathbf{X}_{-v_{i}}\right]-\phi_{rs}$

is $\epsilon_{\ell}(1+\epsilon_{R})$-SW.

Theorem 5.2 implies that we can simply estimate
$\mathbb{E}[\nabla_{1}\hat{R}_{{\rm cv},r}\nabla_{1}\hat{R}_{{\rm
cv},s}|\mathbf{X}_{v_{1}}]$ to approximate $\phi_{rs}$. This leads to the
following procedure, which requires a hold-out set of iid sample points
$X_{1}^{\prime},...,X_{m}^{\prime}$ from the same distribution, that are not
involved in any cross-validation folds. In practice, one can choose a small
but diverging value of $m=n^{a}$ with $a\in(0,1)$, then use $n-m$ sample
points for the $V$-fold cross-validation and $m$ hold-out sample points for
variance estimation.

For $i\in[n]$ and $j\in[m]$, define

$\hat{\mathbf{R}}_{{\rm cv}}^{i,j}$

to be the cross-validation risk vector obtained by replacing $X_{i}$ with
$X_{j}^{\prime}$. Then Theorem 5.2 implies the following.

######  Corollary 5.3.

Define

$\displaystyle\hat{\phi}_{rs}=\frac{n^{2}}{m}\sum_{j=1}^{m/2}\left(\hat{R}_{{\rm
cv},r}^{1,2j-1}-\hat{R}_{{\rm cv},r}^{1,2j}\right)\left(\hat{R}_{{\rm
cv},s}^{1,2j-1}-\hat{R}_{{\rm cv},s}^{2j}\right)\,.$ (12)

Then with probability at least $1-O((n+p)^{-1})$ we have

$\sup_{r,s}|\hat{\phi}_{rs}-\phi_{rs}|\leq\tilde{O}\left(\epsilon_{\ell}(1+\epsilon_{R})+m^{-1/2}\right)\,.$

The estimator in (12) estimates $\mathbb{E}[\nabla_{1}\hat{R}_{{\rm
cv},r}\nabla_{1}\hat{R}_{{\rm cv},s}|\mathbf{X}_{-v_{1}}]$ by taking empirical
average over $m/2$ conditional iid samples given the fitting data
$\mathbf{X}_{-v_{1}}$, which is supported by Lemma C.1 and Lemma C.2. In
practice, we can possibly also use

$\hat{\phi}_{rs}=\frac{n^{2}}{2m}\sum_{j=1}^{m}\left(\hat{R}_{{\rm
cv},r}-\hat{R}_{{\rm cv},r}^{i_{j},j}\right)\left(\hat{R}_{{\rm
cv},s}-\hat{R}_{{\rm cv},s}^{i_{j},j}\right)\,.$

which perturbs different entries instead of just the first one.

Corollary 5.3 provides an entry-wise error bound of the covariance estimation,
which is good enough for Gaussian approximation of the supremum, as
demonstrated in Corollary 4.1. The same kind of inference procedures
considered in Section 4 can be carried over to the deterministic centering
case, which is omitted here as there is little additional insight.

##  6 Discussion

Since its first appearance, high-dimensional Gaussian comparisons have found
wide applications in statistical inference problems, and have been extended
and improved by many authors. In addition to the extensions to dependent data
mentioned above, sharper results on the Gaussian comparison of independent
sums have been obtained in recent literature. For example, see Deng and Zhang
(2020); Kuchibhotla et al. (2021); Lopes (2020); Kuchibhotla and Rinaldo
(2020). In our work, the goal is to develop an asymptotic Gaussian comparison
to serve the purpose of statistical inference. Thus we did not attempt to
obtain the optimal Berry-Esseen type of convergence rates. Our proof uses the
Slepian interpolation as in the original work Chernozhukov et al. (2013), and
it seems possible to obtain better rates of convergence if the more refined
techniques are used.

The main motivation and application is understanding the joint randomness of
many cross-validated risks, and to provide theoretical foundations for
uncertainty quantification of cross-validation based model selection. Our
theory is particularly relevant to the “cross-validation with confidence”
method (Lei, 2020), where one uses the asymptotic Gaussian comparison to
construct a confidence set that contains the best model with a prescribed
confidence level. This method is connected to the literature of model
selection confidence set (Hansen et al., 2011), which has been studied using
various sequential hypothesis testing based approaches (Gunes and Bondell,
2012; Ferrari and Yang, 2015; Jiang et al., 2008). We expect the theory
outlined in this paper to be useful in developing a new model confidence set
estimator using cross-validation with both provable validity guarantees and
good practical performance.

##  Appendix A More notation, definition, and basic properties

###  A.1 Definition and properties of sub-Weibull concentration

######  Definition 2 (sub-Weibull).

Let $K,\alpha$ be positive numbers. We say a random variable $X$ is
$(K,\alpha)$-sub-Weibull (or $(K,\alpha)$-SW) if any of the following holds:

  1. 1.

There exists constant $a$ such that $\mathbb{P}\left(\frac{|X|}{K}\geq
t\right)\leq ae^{-t^{\alpha}}$, for all $t>0$.

  2. 2.

There exists constant $c$ such that $\|X\|_{q}\leq cKq^{1/\alpha}$ for all
$q\geq 1$.

The equivalence of these two definitions can be found in, for example, Theorem
2.1 of Vladimirova et al. (2020). The constants $a,c$ in the definition above
are not important, and are used here so that the two definitions have the same
$(K,\alpha)$ pair.

######  Proposition A.1 (Basic properties of sub-Weibull random variables).

If $X_{i}$ is $(K_{i},\alpha_{i})$-SW, for $i=1,2,$, then

  1. 1.

$X_{1}X_{2}$ is
$(K_{1}K_{2},\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}})$-SW.

  2. 2.

$X_{1}+X_{2}$ is $(K_{1}\vee K_{2},\alpha_{1}\wedge\alpha_{2})$-SW.

The following theorem controls the tail integral of sub-Weibull random
variables.

######  Lemma A.2.

If $Y$ is $(K,\alpha)$-sub-Weibull, then there exists constant $c>0$
independent of $K$ such that $\mathbb{E}[|Y|\mathds{1}(|Y|\geq wK)]\leq
cKw\exp(-w^{\alpha})$ for any $w\geq 1$.

###### Proof of Lemma A.2.

Without loss of generality, assume $Y\geq 0$ and $K=1$. Let $f(y)$ be the
density function of $Y$.

$\displaystyle\mathbb{E}[Y\mathds{1}(Y\geq w)]=$
$\displaystyle\int_{y=w}^{\infty}yf(y)dy$ $\displaystyle=$
$\displaystyle\int_{y=w}^{\infty}\int_{u=0}^{y}duf(y)dy$ $\displaystyle=$
$\displaystyle\int_{u=0}^{w}\int_{y=w}^{\infty}f(y)dydu+\int_{u=w}^{\infty}\int_{y=u}^{\infty}f(y)dydu$
$\displaystyle=$ $\displaystyle w\mathbb{P}(Y\geq
w)+\int_{u=w}^{\infty}\mathbb{P}(Y\geq u)du$ $\displaystyle\leq$
$\displaystyle wa\exp(-w^{\alpha})+\int_{u=w}^{\infty}a\exp(-u^{\alpha})du$
$\displaystyle=$ $\displaystyle
wa\exp(-w^{\alpha})+\frac{a}{\alpha}\int_{v=w^{\alpha}}^{\infty}\exp(-v)v^{\frac{1}{\alpha}-1}dv\,.$

When $\alpha\in(0,1]$, since $w\geq 1$ we have, by (Gabcke, 1979, Prposition
4.4.3)

$\displaystyle\int_{v=w^{\alpha}}^{\infty}\exp(-v)v^{\frac{1}{\alpha}-1}dv\leq\frac{1}{\alpha}e^{-w^{\alpha}}w^{1-1/\alpha}\,.$

So that $\mathbb{E}[Y\mathds{1}(Y\geq w)]\leq
a(1+\alpha^{-2})w\exp(-w^{\alpha})$.

When $\alpha>1$, $v^{1/\alpha-1}\leq 1$ on $[w,\infty)$ since $w\geq 1$, so we
have

$\displaystyle\int_{v=w^{\alpha}}^{\infty}\exp(-v)v^{\frac{1}{\alpha}-1}dv\leq\exp(-w^{\alpha})\,.$

So $\mathbb{E}[Y\mathds{1}(Y\geq w)]\leq 2aw\exp(-w^{\alpha})$ . ∎

The following lemma is a sub-Weibull version of martingale concentration
inequality, showing that the scaling of a martingale with stationary sub-
Weibull increments scales at the speed of $\sqrt{n}$, where $n$ is the time
horizon.

######  Lemma A.3.

Let $M=\sum_{i=1}^{n}M_{i}$ where the sequence $(M_{i})_{i=1}^{n}$ satisfies

  1. 1.

martingale property: $\mathbb{E}(M_{i}|M_{j}:1\leq j<i)=0$ for all $2\leq
i\leq n$, and $\mathbb{E}M_{1}=0$.

  2. 2.

sub-Weibull tail: $\sup_{i}\|M_{i}\|_{q}\leq cKq^{1/\alpha}$ for some
$c,\alpha>0$ and all $q\geq 1$.

Then we have, for $\alpha^{\prime}=\frac{2\alpha}{2+\alpha}$ and a positive
constant $c^{\prime}$,

$\|M\|_{q}\leq c^{\prime}\sqrt{n}Kq^{1/\alpha^{\prime}}\,,~{}~{}\forall q\geq
1\,.$

###### Proof of Lemma A.3.

By Theorem 2.1 of Rio (2009), we have for any $q\geq 2$

$\displaystyle\left\|M\right\|_{q}\leq\left[(q-1)\sum_{i=1}^{n}\|M_{i}\|_{q}^{2}\right]^{1/2}\leq\left[C(q-1)q^{2/\alpha}nK^{2}\right]^{1/2}\leq
C^{1/2}q^{\frac{2+\alpha}{2\alpha}}\sqrt{n}K\,.$

where $C$ is a constant depending only on $c$, and the second inequality
follows from the assumption $\|M_{i}\|_{q}\leq cKq^{1/\alpha}$. ∎

##  Appendix B Proof for random centering

###  B.1 Notation

We first collect some notation for the proof. For the ease of presentation, we
use $\mathbf{W}=(W_{r}:1\leq r\leq p)$ to denote the centered and scaled
random vector for which we would like to establish normal approximation. Thus,
in the proof, the symbol $W_{r}$ may refer to different objects in the proofs
of different theorems. In particular, for the random centering/scaling case
(Theorem 3.1), $W_{r}=\sqrt{n}(\hat{R}_{{\rm cv},r}-\tilde{R}_{r})$, while in
the deterministic centering/scaling case $W_{r}=\sqrt{n}(\hat{R}_{{\rm
cv},r}-R_{r}^{*})$

Recall the following notation:

  * •

$\mathbf{X}_{-v}$: the $\tilde{n}$ ($=n(1-1/V)$) subvector of $\mathbf{X}$
excluding those in index $I_{v}$.

  * •

$\mathbf{X}^{-i}$: the $(n-1)$ subvector of $\mathbf{X}$ excluding the $i$th
entry.

  * •

$\mathbf{X}^{i}$: the iid vector of $\mathbf{X}$ with $i$th entry being
$X_{i}^{\prime}$, an iid copy of $X_{i}$.

For random objects $(U,V)$ and function $f$ acting on $(U,V)$, we will also
use the notation $\mathbb{E}_{U}f(U,V)=\mathbb{E}[f(U,V)|V]$. For example
$\mathbb{E}_{X_{i}}f(\mathbf{X})=\mathbb{E}[f(\mathbf{X})|\mathbf{X}^{-i}]$.

###  B.2 Proof of Theorem 3.1

###### Proof of Theorem 3.1.

Recall the notation:

$W_{r}=\sqrt{n}(\hat{R}_{{\rm
cv},r}-\tilde{R}_{r})=\sum_{i=1}^{n}\frac{1}{\sqrt{n}}K_{r,i}\,.$

Consider the leave-one-out version of $W_{r}$:

$W_{r}^{i}=\sum_{j\neq
i}K_{j}(\mathbf{X}^{i})/\sqrt{n}=W_{r}-\frac{1}{\sqrt{n}}K_{r,i}-\frac{1}{\sqrt{n}}D_{r,i}\,.$

The plan is to use Slepian’s interpolation which smoothly bridges between
$\mathbf{W}$ and the corresponding Gaussian vector $\mathbf{Y}$. In order to
do so, we consider an intermediate object

$\hat{\mathbf{Y}}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{\mathbf{Y}}_{i}\coloneqq\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\Sigma_{v_{i}}^{1/2}\varepsilon_{i}$

with $\varepsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,I_{p})$,
$\Sigma_{v}=[\sigma_{rs,v}]_{1\leq r,s\leq p}$ being the conditional
covariance matrix of $\mathbf{K}(X_{0};\mathbf{X}_{-v})$ given the fitting
data $\mathbf{X}_{-v}$, and $v_{i}$ the fold id of sample point $i$.

Define the interpolating vector, for $t\in(0,1)$

$Z_{r}(t)=\sqrt{t}W_{r}+\sqrt{1-t}\sum_{i}\hat{Y}_{r,i}/\sqrt{n}\,,$

and the corresponding leave-one-out version

$\displaystyle Z_{r}^{i}(t)=\sqrt{t}W_{r}^{i}+\sqrt{1-t}\sum_{j\neq
i}\hat{Y}_{r,i}/\sqrt{n}$

which satsifies

$Z_{r}(t)-Z_{r}^{i}(t)=Z_{r,i}(t)+D_{i}(t)$

with

$\displaystyle Z_{r,i}(t)=$
$\displaystyle\sqrt{t}K_{i}/\sqrt{n}+\sqrt{1-t}\hat{Y}_{r,i}/\sqrt{n}\,,$
$\displaystyle D_{r,i}(t)=$ $\displaystyle\sqrt{t}D_{r,i}/\sqrt{n}\,.$

Let $h:\mathbb{R}^{p}\mapsto\mathbb{R}$ be such that for all $z\in\mathbb{R}$.
Define

$\displaystyle\sum_{r,s=1}^{p}|\partial_{r}\partial_{s}h(z)|=$ $\displaystyle
M_{2}(h)\,,$ (13)
$\displaystyle\sum_{r,s,u=1}^{p}|\partial_{r}\partial_{s}\partial_{u}h(z)|=$
$\displaystyle M_{3}(h)\,.$ (14)

Because
$\mathbb{E}[h(\mathbf{W})-h(\hat{\mathbf{Y}})]=\mathbb{E}\int_{0}^{1}\frac{dh(\mathbf{Z}(t))}{dt}dt$,
the main step in the proof is to control
$\mathbb{E}\frac{dh(\mathbf{Z}(t))}{dt}$.

By Taylor expansion:

$\displaystyle\mathbb{E}\frac{dh(\mathbf{Z}(t))}{dt}$
$\displaystyle=\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}[\partial_{r}h(\mathbf{Z}^{i}(t))Z^{\prime}_{r,i}(t)]$
$\displaystyle~{}~{}+\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\left(Z_{s,i}(t)+D_{s,i}(t)\right)Z_{r,i}^{\prime}(t)]$
$\displaystyle~{}~{}+\sum_{u=1}^{p}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\bigg{\\{}[Z_{s,i}(t)+D_{s,i}(t)][Z_{u,i}(t)+D_{u,i}(t)]$
$\displaystyle~{}~{}~{}~{}\times\left[\int_{0}^{1}(1-v)\partial_{u}\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t)+v\mathbf{Z}_{i}(t))dv\right]Z^{\prime}_{r,i}(t)\bigg{\\}}.$
(15)

The first term in (15)

$\displaystyle\mathbb{E}[\partial_{r}h(\mathbf{Z}^{i}(t))Z^{\prime}_{r,i}(t)]=\mathbb{E}\mathbb{E}_{X_{i},\varepsilon_{i}}[\partial_{r}h(\mathbf{Z}^{i}(t))Z^{\prime}_{r,i}(t)]=\mathbb{E}\left\\{\partial_{r}h(\mathbf{Z}^{i}(t))\mathbb{E}_{X_{i},\varepsilon_{i}}[Z^{\prime}_{r,i}(t)]\right\\}=0\,.$

The second term in (15) consists of two parts. First,

$\displaystyle\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))Z_{s,i}(t)Z_{r,i}^{\prime}(t)]$
$\displaystyle=$
$\displaystyle(2n)^{-1}\mathbb{E}\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\left[(\sqrt{t}K_{s,i}+\sqrt{1-t}\hat{Y}_{s,i})\left(\frac{K_{r,i}}{\sqrt{t}}-\frac{\hat{Y}_{r,i}}{\sqrt{1-t}}\right)\right]$
$\displaystyle=$
$\displaystyle(2n)^{-1}\mathbb{E}\left\\{\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\mathbb{E}_{X_{i},\varepsilon_{i}}\left[(\sqrt{t}K_{s,i}+\sqrt{1-t}\hat{Y}_{s,i})\left(\frac{K_{r,i}}{\sqrt{t}}-\frac{\hat{Y}_{r,i}}{\sqrt{1-t}}\right)\right]\right\\}$
$\displaystyle=$ $\displaystyle 0\,,$

by construction of $\hat{\mathbf{Y}}_{i}$. Now the second term in (15) reduces
to

$\displaystyle\sum_{s,r}\sum_{i}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))D_{s,i}(t)Z_{r,i}^{\prime}(t)]\,.$

By Lemma B.2, $D_{s,i}$ is $\epsilon_{\ell}$-SW. In
$Z_{r,i}^{\prime}(t)=(K_{r,i}/\sqrt{t}-\hat{Y}_{r,i}/\sqrt{1-t})/(2\sqrt{n})$,
$K_{r,i}$ is $1$-SW by Assumption 2, and $\hat{Y}_{r,i}\stackrel{{\scriptstyle
d}}{{=}}\sigma_{rr,v_{i}}\varepsilon_{r,i}$ is also $1$ sub-Weibull as
$\sigma_{r,v_{i}}$ is $1$-sub-Weibull according to the proof of Lemma B.3.
Therefore, $D_{s,i}(t)Z_{r,i}^{\prime}(t)$ is
$n^{-1}\epsilon_{\ell}\eta_{t}$-SW, where

$\eta_{t}=t^{-1/2}\vee(1-t)^{-1/2}\,.$ (16)

Now for any $\tau>0$, by Lemma A.2

$\displaystyle\sum_{s,r}\sum_{i}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))D_{s,i}(t)Z_{r,i}^{\prime}(t)]$
(17) $\displaystyle=$
$\displaystyle\sum_{s,r}\sum_{i}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))D_{s,i}(t)Z_{r,i}^{\prime}(t)\mathds{1}(|D_{s,i}(t)Z_{r,i}^{\prime}(t)|\leq
n^{-1}\tau\epsilon_{\ell}\eta_{t})]$
$\displaystyle~{}~{}+\sum_{s,r}\sum_{i}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))D_{s,i}(t)Z_{r,i}^{\prime}(t)\mathds{1}(|D_{s,i}(t)Z_{r,i}^{\prime}(t)|>n^{-1}\tau\epsilon_{\ell}\eta_{t})]$
$\displaystyle\leq$ $\displaystyle
n^{-1}\tau\epsilon_{\ell}\eta_{t}\sum_{s,r}\sum_{i}\mathbb{E}|\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))|$
$\displaystyle~{}~{}+M_{2}\sum_{s,r}\sum_{i}\mathbb{E}[|D_{s,i}(t)Z_{r,i}^{\prime}(t)|\mathds{1}(|D_{s,i}(t)Z_{r,i}^{\prime}(t)|>n^{-1}\tau\epsilon_{\ell}\eta_{t})]$
$\displaystyle\lesssim$
$\displaystyle\tau\epsilon_{\ell}\eta_{t}M_{2}+np^{2}e^{-\tau^{c}}\,.$

By choosing $\tau=c_{1}\log^{c_{2}}(n+p)$ with appropriate choices of
constants $c_{1},c_{2}$ independent of $(n,p)$, (17) is bounded by
$\tilde{O}(M_{2}\epsilon_{\ell}\eta_{t})$.

The third term in (15) is similarly controlled: Let

$Q_{rsu,i}=\int_{0}^{1}(1-v)\partial_{u}\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t)+v\mathbf{Z}_{i}(t))dv$

and

$T_{rsu,i}=[Z_{s,i}(t)+D_{s,i}(t)][Z_{u,i}(t)+D_{u,i}(t)]Z^{\prime}_{r,i}(t)\,.$

By definition of $M_{3}$ we have

$\sum_{r,s,u}|Q_{rsu,i}|\leq M_{3}$

and $T_{rsu,i}$ is $n^{-3/2}\eta_{t}$-SW. Thus the third term is controlled by

$\tilde{O}(n^{-1/2}M_{3}\eta_{t})\,.$

Since $\eta_{t}$ is integrable on $(0,1)$, we have shown that

$|\mathbb{E}h(\mathbf{W})-\mathbb{E}h(\hat{\mathbf{Y}})|\leq\tilde{O}(\epsilon_{\ell}M_{2}+n^{-1/2}M_{3})\,.$

Combining Lemma B.4 and the anti-concentration result111The anti-concentration
result there is for Gaussian processes. However, our $\hat{\mathbf{Y}}$ is a
Gaussian mixture, because $\hat{\mathbf{Y}}$ is Gaussian only when
conditioning on $\mathbf{X}$. We can condition on $\mathbf{X}$, provided that
$\sup_{r,v}\sigma_{rr,v}\leq\tilde{O}(1)$ with high probability. This can be
established if $\sigma_{rr,v}$ is $1$-SW, which is implied by the proof of
Lemma B.3. (Chernozhukov et al., 2013, Lemma 2.1), we have for any $\beta>0$

$\displaystyle\sup_{z}\left|\mathbb{P}\left(\max_{r}W_{r}\leq
z\right)-\mathbb{P}(\max_{r}\hat{Y}_{r}\leq z)\right|\leq$
$\displaystyle\tilde{O}\left(\epsilon_{\ell}\beta^{2}+n^{-1/2}\beta^{3}+\beta^{-1}\right)$
$\displaystyle\leq$ $\displaystyle\tilde{O}\left(\epsilon_{\ell}^{1/3}\vee
n^{-1/8}\right)\,.$ (18)

where the last inequality follows by choosing
$\beta=\min(\epsilon_{\ell}^{-1/3},n^{1/8})$.

To get the final approximation, let

$\Delta=\max_{r,s,v}|\sigma_{rs,v}-\sigma_{rs}|\,,$ (19)

and event

$\mathcal{E}=\\{\Delta\leq c_{1}\log^{c_{2}}(n+p)\\}\,,$ (20)

with appropriately chosen constants $c_{1},c_{2}$ such that, according to
Lemma B.3,

$\mathbb{P}(\mathcal{E})\geq 1-n^{-1}\,.$

Then

$\displaystyle\mathbb{P}\left(\max_{r}\hat{Y}_{r}\leq z\right)\leq$
$\displaystyle\mathbb{P}\left(\max_{r}\hat{Y}_{r}\leq
z|\mathcal{E}\right)+\mathbb{P}\left(\mathcal{E}^{c}\right)$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left(\max_{r}Y_{r}\leq
z\right)+\tilde{O}(n^{-1/6}+\epsilon_{\ell}^{1/3})+n^{-1}\,,$ (21)

where the last inequality uses Theorem 2 of Chernozhukov et al. (2015) between
$\hat{\mathbf{Y}}$ and $\mathbf{Y}$. On the other hand we have

$\displaystyle\mathbb{P}\left(\max_{r}\hat{Y}_{r}\leq z\right)\geq$
$\displaystyle\mathbb{P}\left(\max_{r}\hat{Y}_{r}\leq
z|\mathcal{E}\right)\mathbb{P}(\mathcal{E})$ $\displaystyle\geq$
$\displaystyle\left[\mathbb{P}\left(\max_{r}Y_{r}\leq
Z\right)-\tilde{O}(n^{-1/6}+\epsilon_{\ell}^{1/3})\right](1-n^{-1})$ (22)
$\displaystyle\geq$ $\displaystyle\mathbb{P}\left(\max_{r}Y_{r}\leq
Z\right)-\tilde{O}(n^{-1/6}+\epsilon_{\ell}^{1/3})\,.$ (23)

The claimed result follows by combining (18), (21), and (23). ∎

###  B.3 Proof of Corollary 4.1

###### Proof of Corollary 4.1.

First, define event $\mathcal{E}_{1}$ on the space of $\mathcal{X}^{n}$ as the
subset consisting all samples of size $n$ such that

$\displaystyle\sup_{v}\|\hat{\Sigma}_{v}-\Sigma\|_{\infty}\lesssim\tilde{O}\left(\frac{1}{\sqrt{n}}+\epsilon_{\ell}\right)\,,$

where the constants $c_{1},c_{2}$ in the $\tilde{O}$ notation is omitted. Then
combining Lemma B.3 and standard sub-Weibull concentration of iid sums we have

$\mathbb{P}(\mathcal{E}_{1})\geq 1-n^{-1}$

with appropriate choice of universal constants in $\tilde{O}(\cdot)$. Here the
$n^{-1/2}$ term comes from $\|\hat{\Sigma}_{v}-\Sigma_{v}\|_{\infty}$ and the
$\epsilon_{\ell}$ term comes from $\|\Sigma_{v}-\Sigma\|_{\infty}$.

Let

$\displaystyle\delta_{0}=$
$\displaystyle\left\|\sqrt{n}\hat{\Lambda}^{-1/2}(\hat{\mathbf{R}}_{\rm
cv}-\tilde{\mathbf{R}})\right\|_{\infty}\,,$ $\displaystyle\delta_{1}=$
$\displaystyle\left\|\sqrt{n}\Lambda^{-1/2}(\hat{\mathbf{R}}_{\rm
cv}-\tilde{\mathbf{R}})\right\|_{\infty}\,,$ $\displaystyle\delta_{2}=$
$\displaystyle\|\mathbf{Y}\|_{\infty}\,,~{}~{}\mathbf{Y}\sim
N(0,\Lambda^{-1/2}\Sigma\Lambda^{-1/2})\,,$ $\displaystyle\delta_{3}=$
$\displaystyle\|\tilde{\mathbf{Y}}\|_{\infty}\,,~{}~{}\tilde{\mathbf{Y}}\sim
N(0,\hat{\Lambda}^{-1/2}\hat{\Sigma}\hat{\Lambda}^{-1/2})\,.$

On $\mathcal{E}_{1}$ we have

$|\delta_{0}-\delta_{1}|\leq\tilde{O}(n^{-1/2}+\epsilon_{\ell})\,.$

Define $\mathcal{E}_{2}$ be the event that
$\left\|\sqrt{n}(\hat{\mathbf{R}}_{\rm cv}-\tilde{\mathbf{R}})\right\|\leq
2\sqrt{\log(n+p)}$. Then Theorem 3.1 implies that

$\mathbb{P}(\mathcal{E}_{2})\geq
1-n^{-1}-\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,.$

Then we have the following approximation.

$\displaystyle\mathbb{P}(\delta_{0}\leq t)\leq$
$\displaystyle\mathbb{P}(\delta_{0}\leq
t\,,~{}\mathcal{E}_{1}\cap\mathcal{E}_{2})+\mathbb{P}(\mathcal{E}_{1}^{c})+\mathbb{P}(\mathcal{E}_{2}^{c})$
$\displaystyle\leq$ $\displaystyle\mathbb{P}(\delta_{1}\leq
t+|\delta_{1}-\delta_{0}|\,,~{}\mathcal{E}_{1}\cap\mathcal{E}_{2})+\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left[\delta_{1}\leq
t+\tilde{O}(n^{-1/2}+\epsilon_{\ell})\right]+\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left[\delta_{2}\leq
t+\tilde{O}(n^{-1/2}+\epsilon_{\ell})\right]+\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left[\delta_{3}\leq
t+\tilde{O}(n^{-1/2}+\epsilon_{\ell})\big{|}\mathcal{E}\right]+\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})$
$\displaystyle\leq$ $\displaystyle\mathbb{P}\left[\delta_{3}\leq
t\right]+\tilde{O}(n^{-1/8}+\epsilon_{\ell}^{1/3})\,,$

where the third inequality holds because on $\mathcal{E}_{1}$
$|\delta_{1}-\delta_{0}|\leq\tilde{O}(n^{-1/2}+\epsilon_{\ell})$; the fourth
inequality holds by applying Theorem 3.1 to the scaled loss functions
$\ell_{r}/\sigma_{rr}^{1/2}$; the fifth inequality holds because conditioning
on the event $\mathcal{E}$ the two Gaussian vectors have covariance matrices
differing at most $\tilde{O}(n^{-1/2}+\epsilon_{\ell})$ and applying Theorem 2
of Chernozhukov et al. (2015); the last inequality holds by anti-concentration
of Gaussian maxima (Chernozhukov et al., 2013, Lemma 2.1).

The corresponding lower probability bound of $\mathbb{P}(\delta_{0}\leq t)$
can be obtained similarly. ∎

###  B.4 Auxiliary lemmas

######  Lemma B.1 (Properties of the difference operator).

Let $f$, $g$ be two functions of the vector
$(X_{1},...,X_{n},X_{1}^{\prime},...,X_{n}^{\prime})$ such that for some
$j\in[n]$, $\mathbb{E}_{X_{j}}g=0$ and $f$ is independent of $X_{j}^{\prime}$,
then

$\mathbb{E}[fg]=\mathbb{E}[(\nabla_{j}f)g]\,.$

###### Proof of Lemma B.1.

Let $f^{j}$ be the iid version of $f$ with input $X_{j}$ replaced by
$X_{j}^{\prime}$. Now It suffices to show that $\mathbb{E}[gf^{j}]=0$, which
holds true since
$\mathbb{E}[gf^{j}]=\mathbb{E}[f^{j}(\mathbb{E}_{X_{j}}g)]=0$. ∎

######  Lemma B.2.

Under Assumption 2, for all $i\in[n]$, $r\in[p]$, $D_{r,i}$ is
$\epsilon_{\ell}$-SW.

###### Proof of Lemma B.2.

Let $F_{k,i}$ be the sigma field generated by $F_{k}$ and $X_{i}^{\prime}$ for
$k\in[n]$, and $F_{0,i}$ be the sigma field genreated by $X_{i}^{\prime}$.
Because $\mathbb{E}(\nabla_{i}K_{j}|X_{i}^{\prime})=0$ for all $j\neq
I_{v_{i}}$, we have the following martingale sum representation

$D_{r,i}=\sum_{k=1}^{n}\mathbb{E}(D_{r,i}|F_{k,i})-\mathbb{E}(D_{r,i}|F_{k-1,i})=\sum_{k=1}^{n}\mathbb{E}(\tilde{\nabla}_{k}D_{r,i}|F_{k,i})\,,$

where $\tilde{\nabla}_{k}$ is the same operator as $\nabla_{k}$ expect that it
replaces $X_{k}$ by $X_{k}^{\prime\prime}$, a further iid copy. This is to
make sure that the difference operator $\tilde{\nabla}_{i}$ does not interfere
with $X_{i}^{\prime}$, which is already involved in $D_{r,i}$.

For each $k\in[n]$, if $k\in\\{i,j\\}$ we have
$\tilde{\nabla}_{k}\nabla_{i}K_{j}$ is $(\epsilon_{\ell}n^{-1/2},\alpha)$-SW
by part 1 of Assumption 2 and closure of sub-Weibull tails under additions.

For each $k\in[n]\backslash\\{i,j\\}$ we have
$\tilde{\nabla}_{k}\nabla_{i}K_{j}$ is either $0$ (if $k\in I_{v_{j}}$) or
$(\epsilon_{\ell}n^{-3/2},\alpha)$-SW (if $k\neq i,j$, by part 2 of Assumption
2).

So overall we conclude that $\tilde{\nabla}_{k}D_{r,i}$ is
$(\epsilon_{\ell}n^{1/2},\alpha)$-SW. The claimed result follows from Lemma
A.3. ∎

######  Lemma B.3.

Under Assumption 2, $\sigma_{rs,1}-\sigma_{rs}$ is $\epsilon_{\ell}$-SW.

###### Proof of Lemma B.3.

Without loss of generality, we work with $v=1$. First write
$\sigma_{rs,1}-\sigma_{rs}$ as the sum of martingale increments

$\displaystyle\sigma_{rs,1}-\sigma_{rs}=\sum_{i=1}^{\tilde{n}}\mathbb{E}(\nabla_{i}\sigma_{rs,1}|F_{i})\,.$
(24)

Next we control each $\mathbb{E}(\nabla_{i}\sigma_{rs,1}|F_{i})$. Let
$\|\cdot\|$ be any $L_{q}$ norm with $q\geq 1$,

$\displaystyle\|\mathbb{E}(\nabla_{i}\sigma_{rs,1}|F_{i})\|\leq\|\nabla_{i}\sigma_{rs,1}\|$
$\displaystyle\leq$
$\displaystyle\left\|\mathbb{E}_{0}\left(K_{r}(X_{0},\mathbf{X}_{-1})K_{s}(X_{0},\mathbf{X}_{-1})-K_{r}(X_{0},\mathbf{X}_{-1}^{i})K_{s}(X_{0},\mathbf{X}_{-1}^{i})\right)\right\|$
$\displaystyle\leq$
$\displaystyle\left\|K_{r}(X_{0},\mathbf{X}_{-1})K_{s}(X_{0},\mathbf{X}_{-1})-K_{r}(X_{0},\mathbf{X}_{-1}^{i})K_{s}(X_{0},\mathbf{X}_{-1}^{i})\right\|$
$\displaystyle\leq$
$\displaystyle\left\|\left[\nabla_{i}K(X_{0},\mathbf{X}_{-1})\right]K_{s}(X_{0},\mathbf{X}_{-1})\right\|+\left\|K_{r}(X_{0},\mathbf{X}_{-1}^{i})\left[\nabla_{i}K_{s}(X_{0},\mathbf{X}_{-1})\right]\right\|\,.$

Then it follows from Assumption 2 and Proposition A.1 that
$\mathbb{E}(\nabla_{i}\sigma_{rs,1}|F_{i})$ is $\epsilon_{\ell}n^{-1/2}$-SW.
Furhter applying Lemma A.3 to the martingale sum (24) we conclude that
$\sigma_{rs,1}-\sigma_{rs}$ is $\epsilon_{\ell}$-SW. ∎

######  Lemma B.4 (Bridging between smooth function and CDF of maximum).

For any $\beta>0$, there exists a function
$h=h_{\beta}:\mathbb{R}^{p}\mapsto\mathbb{R}$, such that for any random vector
$\mathbf{Z}\in\mathbb{R}^{p}$

$\mathbb{P}\left(\max_{r}Z_{r}\leq
t\right)\leq\mathbb{E}h(\mathbf{Z})\leq\mathbb{P}\left(\max_{r}Z_{r}\leq
t+\frac{\log p+1}{\beta}\right)$

and, for some universal constant $C$,

$\displaystyle M_{2}(h)=$
$\displaystyle\sup_{z\in\mathbb{R}^{p}}\sum_{r,s=1}^{p}|\partial_{r}\partial_{s}h(z)|\leq
C\beta^{2}\,,$ $\displaystyle M_{3}(h)=$
$\displaystyle\sup_{z\in\mathbb{R}^{p}}\sum_{r,s,u=1}^{p}|\partial_{r}\partial_{s}\partial_{u}h(z)|\leq
C\beta^{3}\,.$

###### Proof of Lemma B.4.

See Lemma A.5 and Corollary I.1 of Chernozhukov et al. (2013). ∎

##  Appendix C Proof for deterministic centering

###  C.1 Preparation

In preparation for the proof, we first take a closer look at some intermediate
quantities involved in forming the asymptotic covariance term. Intuitively,
the variance contributed by $X_{i}$ in the term
$\ell_{r}(X_{j};\mathbf{X}_{-v_{j}})$ for $j\notin I_{v_{i}}$ is from the
variability of $R_{r}(\mathbf{X}_{-v_{j}})$. Based on this intuition, we can
reduce (10) to

$\displaystyle\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})+\sum_{j\notin
I_{v_{i}}}R_{r}(\mathbf{X}_{-v_{j}})\eqqcolon g_{r,i}\,.$ (25)

The quantities in (25) still involve many sample points $X_{j}$ ($j\neq i$).
In order to pinpoint the variance contributed by $X_{i}$ alone, we consider
the following difference versions of $g_{r,i}$:

$\mathbb{E}(g_{r,i}|F_{i})-\mathbb{E}(g_{r,i}|F_{i-1})\,,\text{ and
}g_{r,i}-\mathbb{E}(g_{r,i}|\mathbf{X}^{-i})\,,$

where $F_{i}$ is the $\sigma$-field generated by $(X_{1},...,X_{i})$ and
$\mathbf{X}^{-i}=(X_{1},...,X_{i-1},X_{i+1},...,X_{n})$. The reason to
consider these two differences is rather technical, where the former allows us
to express $g_{r,i}$ as the sum of a sequence of martingale increments
$\\{\mathbb{E}(g_{r,i}|F_{i})-\mathbb{E}(g_{r,i}|F_{i-1}):~{}i=1,...,n\\}$,
and the latter provides $\mathbb{E}(g_{r,i}|\mathbf{X}^{-i})$ as a leave-one-
out approximation to $g_{r,i}$ with a manageable difference.

Let

$C_{rs,i}=\left[\mathbb{E}(g_{r,i}|F_{i})-\mathbb{E}(g_{r,i}|F_{i-1})\right]\left[g_{s,i}-\mathbb{E}(g_{s,i}|\mathbf{X}^{-i})\right]\,.$
(26)

In Lemma C.4 we will show that

$\mathbb{E}C_{rs,i}\approx\phi_{rs}\,.$ (27)

A key step in the proof is to ensure that this covariance is indeed
contributed mostly by $X_{i}$, which amounts to controlling

$\mathbb{E}(C_{rs,i}|\mathbf{X}^{-i})-\mathbb{E}(C_{rs,i})\,.$

It can be shown that
$\|\mathbb{E}(C_{rs,i}|\mathbf{X}^{-i})-\mathbb{E}(C_{rs,i})\|_{2}$ is small
using the standard Efron-Stein inequality. However, the high-dimensionality
requires some uniform bound of the realized values
$\mathbb{E}(C_{rs,i}|\mathbf{X}^{-i})-\mathbb{E}(C_{rs,i})$ over the triplet
$(r,s,i)$. This is established using our sub-Weibull conditions in Lemma C.3.

###  C.2 Proof of Theorem 5.1

###### Proof of Theorem 5.1.

Throughout this proof, the notation may be different from that in the proof of
Theorem 3.1, especially for $W$ and $Z$, due to the different centering and
scaling.

Let $\mathbf{Y}=(Y_{1},...,Y_{p})=\sum_{i=1}^{n}n^{-1/2}\varepsilon_{i}$,
where $\varepsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\bm{\Phi})$.
Define $\bm{\varepsilon}^{i}$ as the vector $(\varepsilon_{i}:1\leq i\leq n)$
with the $i$th element $\varepsilon_{i}$ replaced by its iid copy
$\varepsilon_{i}^{\prime}$. Note that each
$\varepsilon_{i}=(\varepsilon_{1,i},...,\varepsilon_{p,i})$ is itself a
$p$-dimensional vector.

Let $F_{i}$ be the $\sigma$-field generated by $\\{(X_{j}):j\leq i\\}$, and
for any function $f$ acting on $\mathbf{X}$, define

$\nabla_{i}f=f(\mathbf{X})-f(\mathbf{X}^{i})$

and

$\Delta_{i}f=\mathbb{E}(f|F_{i})-\mathbb{E}(f|F_{i-1})=\mathbb{E}(\nabla_{i}f|F_{i})\,.$
(28)

Define quantities

$\displaystyle W_{r}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left\\{\ell_{r}(X_{i},\mathbf{X}_{-v_{i}})-R_{r}\right\\}\,,$
$\displaystyle W_{r}^{i}$ $\displaystyle=\mathbb{E}(W_{r}|\mathbf{X}_{-i})\,.$

where $W_{r}$ is the deterministically centered quantity for which we would
like to establish Gaussian comparison, and $W_{r}^{i}$ is the corresponding
leave-one-out version.

We then have the following useful facts

$\displaystyle W_{r}-W_{r}^{i}=$ $\displaystyle
W_{r}-\mathbb{E}(W_{r}(\mathbf{X}^{i})|\mathbf{X})=\mathbb{E}(\nabla_{i}W_{r}|\mathbf{X})\,,$
(29)

and

$\displaystyle\nabla_{i}W_{r}=$ $\displaystyle
n^{-1/2}(\nabla_{i}g_{r,i}+D_{r,i})$ (30)

with $D_{r,i}$, $g_{r,i}$ defined in (4) and (25), respectively.

Similarly, for $t\in(0,1)$, consider interpolating variable

$Z_{r}(t)=\sqrt{t}W_{r}+\sqrt{1-t}Y_{r}\,,$

the leave-one-out version

$Z_{r}^{i}(t)=\sqrt{t}W_{r}^{i}+\sqrt{1-t}\sum_{j\neq
i}\varepsilon_{r,j}/\sqrt{n}\,.$

and the martingale increment with respect to the filtration $\\{F_{i}:1\leq
i\leq n\\}$,

$\Delta_{i}Z_{r}(t)=\sqrt{t}(\Delta_{i}W_{r})+\sqrt{1-t}\varepsilon_{r,i}/\sqrt{n}\,.$

Consider $h:\mathbb{R}^{p}\mapsto\mathbb{R}$ with quantities $M_{2}$ and
$M_{3}$ defined as in (13) and (14). In the following we will focus on
controlling $\mathbb{E}\left[h(\mathbf{W})-h(\mathbf{Y})\right]$, where the
bold font symbols represent the corresponding $p$-dimensional vectors:
$\mathbf{W}=(W_{1},...,W_{p})$, $\mathbf{Y}=(Y_{1},...,Y_{p})$, and
$\mathbf{Z}^{i}(t)=(Z_{1}^{i}(t),...,Z_{p}^{i}(t))$, etc. The only exception
is $\mathbf{X}$, which corresponds to the collection of $n$ iid samples
$(X_{1},...,X_{n})\in\mathcal{X}^{n}$.

Write $\frac{d\Delta_{i}Z_{r}(t)}{dt}=Z^{\prime}_{r,i}(t)$. By Taylor
expansion:

$\displaystyle\mathbb{E}\frac{dh(\mathbf{Z}(t))}{dt}$
$\displaystyle=\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}[\partial_{r}h(\mathbf{Z}^{i}(t))Z^{\prime}_{r,i}(t)]$
$\displaystyle~{}~{}+\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\left(Z_{s}(t)-Z_{s}^{i}(t)\right)Z_{r,i}^{\prime}(t)]$
$\displaystyle~{}~{}+\sum_{u=1}^{p}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\bigg{\\{}[Z_{s}(t)-Z_{s}^{i}(t)][Z_{u}(t)-Z_{u}^{i}(t)]$
$\displaystyle~{}~{}~{}~{}\times\left[\int_{0}^{1}(1-v)\partial_{u}\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t)+v\mathbf{Z}_{i}(t))dv\right]Z^{\prime}_{r,i}(t)\bigg{\\}}.$
(31)

The first term equals $0$ because $\mathbf{Z}^{i}(t)$ does not involve
$(X_{i},\varepsilon_{i})$, and
$\mathbb{E}_{X_{i},\varepsilon_{i}}Z_{r,i}^{\prime}(t)=0$.

The second term can be written as

$\displaystyle\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}[\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\left(Z_{s}(t)-Z_{s}^{i}(t)\right)Z_{r,i}^{\prime}(t)]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left(\frac{\Delta_{i}W_{r}}{\sqrt{t}}-\frac{\varepsilon_{r,i}}{\sqrt{n}\sqrt{1-t}}\right)\left[\sqrt{t}\left(W_{s}-W_{s}^{i}\right)+\sqrt{1-t}\frac{\varepsilon_{s,i}}{\sqrt{n}}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[\Delta_{i}W_{r}\left(W_{s}-W_{s}^{i}\right)-\phi_{rs}/n\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[\mathbb{E}(\nabla_{i}W_{r}|F_{i})\mathbb{E}(\nabla_{i}W_{s}|\mathbf{X})-\phi_{rs}/n\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[C_{rs,i}-\phi_{rs}+B_{rs,i}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}\,,$
(32)

where
$C_{rs,i}=\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})\mathbb{E}(\nabla_{i}g_{s,i}|\mathbf{X})$
is the same as defined in (26) and

$\displaystyle B_{rs,i}=$
$\displaystyle\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})\mathbb{E}(D_{s,i}|\mathbf{X})+\mathbb{E}(D_{r,i}|F_{i})\mathbb{E}(\nabla_{i}g_{s,i}|\mathbf{X})+\mathbb{E}(D_{r,i}|F_{i})\mathbb{E}(D_{s,i}|\mathbf{X})\,.$

In (32), the first equation follows by construction of $Z_{r,i}^{\prime}(t)$
and $Z_{s}^{i}$; the second equation follows by taking conditional expectation
over $\varepsilon_{i}$ and the definition of $\phi_{rs}$; the third and fourth
equations follow from (29) and (30), respectively.

By Assumptions 2, 3 and B.2 we have
$B_{rs,i}/[(1+\epsilon_{R})\epsilon_{\ell}]$ is sub-Weibull, using the same
argument as in (17), we have

$\displaystyle\frac{1}{2n}\mathbb{E}\sum_{rs,i}|B_{rs,i}\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))|\leq\tilde{O}\left[(1+\epsilon_{R})\epsilon_{\ell}M_{2}\right]\,.$
(33)

Define $\bar{C}_{rs,i}=\mathbb{E}(C_{rs,i}|\mathbf{X}^{-i})$. Now we are left
with the term

$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[C_{rs,i}-\phi_{rs}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[\bar{C}_{rs,i}-\phi_{rs}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[\bar{C}_{rs,i}-\mathbb{E}C_{rs,i}+\mathbb{E}C_{rs,i}-\phi_{rs}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$

where the equality holds by taking conditional expectation given
$\mathbf{X}^{-i}$ and realizing $\mathbf{Z}^{i}(t)$ is independent of $X_{i}$.
By Lemma C.3, $\bar{C}_{rs,i}-\mathbb{E}C_{rs,i}$ is
$\epsilon_{\ell}(1+\epsilon_{R})$-sub-Weibull. Repeating the truncation
argument used in (17) we get

$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{\left[\bar{C}_{rs,i}-\mathbb{E}C_{rs,i}\right]\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))\right\\}$
$\displaystyle\leq$
$\displaystyle\tilde{O}\left[\epsilon_{\ell}(1+\epsilon_{R})M_{2}\right]\,.$
(34)

We still need to control $\mathbb{E}C_{rs,i}-\phi_{rs}$, which is provided by
Lemma C.4. Thus we obtain

$\displaystyle\frac{1}{2n}\sum_{s=1}^{p}\sum_{r=1}^{p}\sum_{i=1}^{n}\mathbb{E}\left\\{|\mathbb{E}C_{rs,i}-\phi_{rs}|\times|\partial_{s}\partial_{r}h(\mathbf{Z}^{i}(t))|\right\\}\leq
c\epsilon_{\ell}M_{2}\,,$ (35)

for some universal constant $c$.

Combining (33), (34), and (35) into (32) we conclude that the second term in
(31) is upper bounded in absolute value by (using the simplifying assumption
$\epsilon_{\ell}<1$.)

$\tilde{O}\left[\epsilon_{\ell}(\epsilon_{R}+1)M_{2}\right]\,.$ (36)

The third term in (31) can be controlled using the following equation

$\displaystyle Z_{r}(t)-Z_{r}^{i}(t)=$
$\displaystyle\frac{1}{\sqrt{n}}\mathbb{E}\left[g_{r,i}\sqrt{t}+D_{r,i}\sqrt{t}+\varepsilon_{r,i}\sqrt{1-t}\Big{|}\mathbf{X}\right]\,,$
(37)

which holds by combining (29) and (30), and, by Assumptions 2, 3 and B.2, is
$n^{-1/2}(1+\epsilon_{R})$-sub-Weibull.

Similarly, by (28) and (30) we have

$Z_{r,i}^{\prime}(t)=\frac{1}{2\sqrt{n}}\left(\frac{\mathbb{E}(\nabla_{i}g_{r,i}+D_{r,i}|F_{i})}{\sqrt{t}}-\frac{\varepsilon_{r,i}}{\sqrt{1-t}}\right)\,,$

and hence

$\frac{Z_{r,i}^{\prime}(t)}{n^{-1/2}\eta_{t}(1+\epsilon_{R})}$

is sub-Weibull, where $\eta_{t}$ is defined in (16).

Putting together the sub-Weibull properties of $Z_{r}(t)-Z_{r}^{i}(t)$ and
$Z_{r,i}^{\prime}(t)$, we have

$\frac{[Z_{s}(t)-Z_{s}^{i}(t)][Z_{u}(t)-Z_{u}^{i}(t)]Z_{r,i}^{\prime}(t)}{n^{-3/2}\eta_{t}(1+\epsilon_{R})^{3}}$

is sub-Weibull. Therefore, applying the truncation argument in (17) again in
the third term of (31), we obtain an upper bound of

$\displaystyle\tilde{O}\left\\{n^{-1/2}\eta_{t}(1+\epsilon_{R})^{3}M_{3}\right\\}\,.$
(38)

Combining (36) and (38) with (31) and integrate the latter over $t\in(0,1)$ we
obtain

$\displaystyle|\mathbb{E}[h(\mathbf{W})-h(\mathbf{Y})]|\leq\tilde{O}\left[\epsilon_{\ell}(1+\epsilon_{R})M_{2}+n^{-1/2}(1+\epsilon_{R})^{3}M_{3}\right]\,.$
(39)

Again, using Lemma B.4 and the anti-concentration result (Chernozhukov et al.,
2013, Lemma 2.1), we have for any $\beta>0$

$\displaystyle\sup_{z}\left|\mathbb{P}\left(\max_{r}W_{r}\leq
z\right)-\mathbb{P}(\max_{r}Y_{r}\leq z)\right|$ (40)
$\displaystyle\leq\tilde{O}\left[\epsilon_{\ell}(1+\epsilon_{R})\beta^{2}+n^{-1/2}(1+\epsilon_{R})^{3}\beta^{3}+\beta^{-1}\right]$
$\displaystyle\leq\tilde{O}\left([\epsilon_{\ell}(1+\epsilon_{R})]^{1/3}\vee
n^{-1/8}(1+\epsilon_{R})^{3/4}\right)\,$ (41)

where the last inequality follows by choosing
$\beta=\min([\epsilon_{\ell}(1+\epsilon_{R})]^{-1/3},n^{1/8}(1+\epsilon_{R})^{-3/4})$.

∎

###  C.3 Proof of variance estimation with deterministic centering (Theorem
5.2)

The claimed result in Theorem 5.2 follows directly from the two following
lemmas.

######  Lemma C.1.

Under Assumptions 1, 2, 3 and 4,

$\sup_{1\leq r,s\leq
p}\left|\frac{n^{2}}{2}\mathbb{E}\left(\nabla_{i}\hat{R}_{{\rm
cv},r}\nabla_{i}\hat{R}_{{\rm
cv},s}\right)-\phi_{rs}\right|=O\left(\epsilon_{\ell}(1+\epsilon_{R})\right)\,.$

Lemma C.1 shows that the variance of $\nabla_{i}\hat{\mathbf{R}}_{{\rm cv}}$
is entry-wise close to the true covariance matrix $\Phi$. The next result
further reduces this variance term to the proportion contributed solely by
$X_{i}$.

######  Lemma C.2.

Under Assumptions 1, 2, 3 and 4,

$n^{2}\sup_{1\leq r,s\leq p}\left|\mathbb{E}\left[\nabla_{i}\hat{R}_{{\rm
cv},r}\nabla_{i}\hat{R}_{{\rm
cv},s}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}\left[\nabla_{i}\hat{R}_{{\rm
cv},r}\nabla_{i}\hat{R}_{{\rm
cv},s}\right]\right|\lesssim\tilde{O}\left[\epsilon_{\ell}(1+\epsilon_{R})\right]\,.$

###### Proof of Lemma C.1.

Let $f_{r}=n\hat{R}_{{\rm
cv},r}=\sum_{i=1}^{n}\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})$. For $-1\leq j\leq
n$, $j\neq i$, define

$\displaystyle E_{j,i}=\left\\{\begin{array}[]{ll}F_{0}\,,&j=-1\\\
\sigma(X_{i},X_{i}^{\prime})\,,&j=0\\\
\sigma(X_{1},...,X_{j},X_{i},X_{i}^{\prime})\,,&1\leq j\leq
n\,.\end{array}\right.$

and

$\displaystyle E_{j,i}^{-}=\left\\{\begin{array}[]{ll}E_{j-1,i}\,,&j\neq
i+1\,,\\\ E_{i-1,i}\,,&j=i+1\,.\end{array}\right.$

Then $\\{E_{j,i}:-1\leq j\leq n,~{}j\neq i\\}$ is a filtration.

Use notation $\ell_{r,i}=\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})$,
$\bar{\ell}_{r,i}=\bar{\ell}_{r}(X_{i})$,
$\ell_{r,i}^{\prime}=\ell_{r}(X_{i}^{\prime};\mathbf{X}_{-v_{i}})$,
$\bar{\ell}_{r,i}^{\prime}=\bar{\ell}_{r}(X_{i}^{\prime})$, and
$\bar{R}_{r,i}=\bar{R}_{r}(X_{i})$.

Using the decomposition

$\nabla_{i}f_{r}=\nabla_{i}K_{r,i}+\nabla_{i}\mathcal{R}_{r,i}+D_{r,i}$

we get

$\displaystyle\nabla_{i}f_{r}\nabla_{i}f_{s}=$
$\displaystyle\nabla_{i}K_{r,i}\nabla_{i}K_{s,i}+\nabla_{i}\mathcal{R}_{r,i}\nabla_{i}\mathcal{R}_{s,i}+\nabla_{i}K_{r,i}\nabla_{i}\mathcal{R}_{s,i}+\nabla_{i}\mathcal{R}_{r,i}\nabla_{i}K_{s,i}$
$\displaystyle+D_{r,i}\nabla_{i}K_{s,i}+D_{s,i}\nabla_{i}K_{r,i}+D_{r,i}\nabla_{i}\mathcal{R}_{s,i}+D_{s,i}\nabla_{i}\mathcal{R}_{r,i}+D_{r,i}D_{s,i}\,.$
(42)

For the first term

$\displaystyle\mathbb{E}\nabla_{i}K_{r,i}\nabla_{i}K_{s,i}$ $\displaystyle=$
$\displaystyle\mathbb{E}\nabla_{i}\ell_{r,i}\nabla_{i}\ell_{s,i}$
$\displaystyle=$ $\displaystyle{\rm
Cov}\left\\{\mathbb{E}\left[\nabla_{i}\ell_{r,i}|X_{i},X_{i}^{\prime}\right],~{}\mathbb{E}\left[\nabla_{i}\ell_{s,i}|X_{i},X_{i}^{\prime}\right]\right\\}$
$\displaystyle+\mathbb{E}\left\\{{\rm
Cov}\left[\nabla_{i}\ell_{r,i},~{}\nabla_{i}\ell_{s,i}|X_{i},X_{i}^{\prime}\right]\right\\}$
$\displaystyle=$ $\displaystyle 2{\rm
Cov}(\bar{\ell}_{r,i},\bar{\ell}_{s,i})+\mathbb{E}\left\\{{\rm
Cov}\left[\nabla_{i}\ell_{r,i},~{}\nabla_{i}\ell_{s,i}|X_{i},X_{i}^{\prime}\right]\right\\}$

where the second term is upper bounded by

$\displaystyle\left|\mathbb{E}\left\\{{\rm
Cov}\left[\nabla_{i}\ell_{r,i},~{}\nabla_{i}\ell_{s,i}|X_{i},X_{i}^{\prime}\right]\right\\}\right|$
$\displaystyle\leq$ $\displaystyle\frac{1}{2}\mathbb{E}\left\\{{\rm
Var}\left[\nabla_{i}\ell_{r,i}|X_{i},X_{i}^{\prime}\right]+{\rm
Var}\left[\nabla_{i}\ell_{s,i}|X_{i},X_{i}^{\prime}\right]\right\\}$
$\displaystyle\leq$ $\displaystyle 2\mathbb{E}\left\\{{\rm
Var}(\ell_{r,i}|X_{i})+{\rm Var}(\ell_{s,i}|x_{i})\right\\}$
$\displaystyle\leq$ $\displaystyle\sum_{j\neq
I_{v_{i}}}\mathbb{E}(\nabla_{j}\ell_{r,i})^{2}+\mathbb{E}(\nabla_{j}\ell_{s,i})^{2}$
$\displaystyle\leq$ $\displaystyle 2\epsilon_{\ell}^{2}$

where the second last step used Efron-Stein inequality.

So we conclude

$\left|\mathbb{E}\nabla_{i}K_{r,i}\nabla_{i}K_{s,i}-2{\rm
Cov}(\bar{\ell}_{r,1},\bar{\ell}_{s,1})\right|\leq 2\epsilon_{\ell}^{2}$ (43)

For the second term in (42), using the martingale decomposition

$\displaystyle\nabla_{i}\mathcal{R}_{r,i}=\sum_{0\leq j\leq n,j\neq
i}\mathbb{E}(\nabla_{i}\mathcal{R}_{r,i}|E_{j,i})-\mathbb{E}(\nabla_{i}\mathcal{R}_{r,i}|E_{j,i}^{-})$

we have, by orthogonality between martingale increments,

$\displaystyle\mathbb{E}\nabla_{i}\mathcal{R}_{r,i}\nabla_{i}\mathcal{R}_{s,i}$
$\displaystyle=$ $\displaystyle\mathbb{E}\sum_{0\leq j\leq n,j\neq
i}\left[\mathbb{E}(\nabla_{i}\mathcal{R}_{r,i}|E_{j,i})-\mathbb{E}(\nabla_{i}\mathcal{R}_{r,i}|E_{j,i}^{-})\right]\left[\mathbb{E}(\nabla_{i}\mathcal{R}_{s,i}|E_{j,i})-\mathbb{E}(\nabla_{i}\mathcal{R}_{s,i}|E_{j,i}^{-})\right]$
$\displaystyle=$ $\displaystyle 2\tilde{n}^{2}{\rm
Cov}(\bar{R}_{r}(X_{1}),\bar{R}_{s}(X_{1}))+\mathbb{E}\sum_{1\leq j\leq
n,j\neq
i}\left[\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{r,i}|E_{j,i})\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{s,i}|E_{j,i})\right]$

and the remainder term satisfies

$\displaystyle\left|\mathbb{E}\sum_{1\leq j\leq n,j\neq
i}\left[\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{r,i}|E_{j,i})\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{s,i}|E_{j,i})\right]\right|$
$\displaystyle\leq$ $\displaystyle\frac{1}{2}\sum_{1\leq j\leq n,j\neq
i}\left[\|\nabla_{j}\nabla_{i}\mathcal{R}_{r,i}\|_{2}^{2}+\|\nabla_{j}\nabla_{i}\mathcal{R}_{s,i}\|_{2}^{2}\right]$
$\displaystyle\leq$ $\displaystyle\epsilon_{\ell}^{2}$

Hence we have

$\left|\mathbb{E}\nabla_{i}\mathcal{R}_{r,i}\nabla_{i}\mathcal{R}_{s,i}-2\tilde{n}^{2}{\rm
Cov}(\bar{R}_{r,1},\bar{R}_{s,1})\right|\leq\epsilon_{\ell}^{2}\,.$ (44)

For the third term in (42), using the martingale decomposition of
$\nabla_{i}K_{r,i}$ and $\nabla_{i}\mathcal{R}_{s,i}$, we have

$\displaystyle\mathbb{E}\nabla_{i}K_{r,i}\nabla_{i}\mathcal{R}_{s,i}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathbb{E}(\nabla_{i}K_{r,i}|X_{i},X_{i}^{\prime})\mathbb{E}(\nabla_{i}\mathcal{R}_{s,i}|X_{i},X_{i}^{\prime})\right]+\sum_{1\leq
j\leq n,j\neq
i}\mathbb{E}\left[\mathbb{E}(\nabla_{j}\nabla_{i}K_{r,i})\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{s,i})\right]$
$\displaystyle=$ $\displaystyle 2\tilde{n}{\rm
Cov}(\bar{\ell}_{r,1},\bar{R}_{s,1})+\sum_{1\leq j\leq n,j\neq
i}\mathbb{E}\left[\mathbb{E}(\nabla_{j}\nabla_{i}K_{r,i})\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{s,i})\right]$

where the remainder term satisfies

$\displaystyle\left|\sum_{1\leq j\leq n,j\neq
i}\mathbb{E}\left[\mathbb{E}(\nabla_{j}\nabla_{i}K_{r,i})\mathbb{E}(\nabla_{j}\nabla_{i}\mathcal{R}_{s,i})\right]\right|$
$\displaystyle\leq$ $\displaystyle\sum_{1\leq j\leq n,j\neq
i}\|\nabla_{j}\nabla_{i}K_{r,i}\|_{2}\|\nabla_{j}\nabla_{i}\mathcal{R}_{s,i}\|_{2}$
$\displaystyle\leq$ $\displaystyle\epsilon_{\ell}^{2}\,.$

The fourth term can be bounded similarly. So we have

$\displaystyle\left|\nabla_{i}K_{r,i}\nabla_{i}\mathcal{R}_{s,i}+\nabla_{i}\mathcal{R}_{r,i}\nabla_{i}K_{s,i}-2\tilde{n}{\rm
Cov}(\bar{\ell}_{r,1}\bar{R}_{s,1})-2\tilde{n}{\rm
Cov}(\bar{R}_{r,1},\bar{\ell}_{s,1})\right|\leq 2\epsilon_{\ell}^{2}\,.$ (45)

For the other terms in (42), according to Lemma B.2, Assumption 2, and
Assumption 3, we have $\|D_{r,i}\|_{2}\lesssim\epsilon_{\ell}$,
$\|\nabla_{i}K_{r,i}\|_{2}\lesssim 1$, and
$\nabla_{i}\mathcal{R}_{r,i}\leq\epsilon_{R}$, so that the last five terms in
(42) are bounded by, up to constant factor, $\epsilon_{\ell}(1+\epsilon_{R})$.
The claimed result is proved. ∎

###### Proof of Lemma C.2.

We follow the notation in the proof of Lemma C.1. By the symmetry assumption
of $\ell$, we can assume $i>\tilde{n}$ without loss of generality.

Let $M=\nabla_{i}f_{r}\nabla_{i}f_{s}-\phi_{rs}$, and
$M_{j}=\mathbb{E}(M|F_{j})-\mathbb{E}(M|F_{j-1})=\mathbb{E}(\nabla_{j}M|F_{j})$
for $j=1,...,\tilde{n}$. The main task in the proof is to control
$\|M_{j}\|_{\psi_{\alpha}}$, which further reduces to controlling the norm of
$\nabla_{j}(\nabla_{i}f_{r}\nabla_{i}f_{s})$.

To begin with, we first write

$\nabla_{j}(\nabla_{i}f_{r}\nabla_{i}f_{s})=\nabla_{j}\nabla_{i}f_{r}\nabla_{i}f_{s}+\nabla_{i}f_{s}(\mathbf{X}^{j})\nabla_{j}\nabla_{i}f_{r}\,.$

Use the decomposition $\nabla_{i}f_{r}=\nabla_{i}g_{r,i}+D_{r,i}$ we have

$\displaystyle\left\|\mathbb{E}\left(\nabla_{i}f_{r}\nabla_{i}f_{s}|\mathbf{X}_{-v_{i}}\right)-\mathbb{E}(\nabla_{i}f_{r}\nabla_{i}f_{s})\right\|$
$\displaystyle\leq$
$\displaystyle\left\|\mathbb{E}\left[\nabla_{i}g_{r,i}\nabla_{i}g_{s,i}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}(\nabla_{i}g_{r,i}\nabla_{i}g_{s,i})\right\|$
$\displaystyle+\left\|\mathbb{E}\left[\nabla_{i}g_{r,i}D_{s,i}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}(\nabla_{i}g_{r,i}D_{s,i})\right\|$
$\displaystyle+\left\|\mathbb{E}\left[D_{r,i}\nabla_{i}g_{s,i}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}(D_{r,i}\nabla_{i}g_{s,i})\right\|$
$\displaystyle+\left\|\mathbb{E}\left[D_{r,i}D_{s,i}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}(D_{r,i}D_{s,i})\right\|\,.$

Using the fact that for any random variable $W$

$\|W-\mathbb{E}W\|_{q}\leq 2\|W\|_{q}\,,~{}~{}\forall~{}~{}q\geq 1\,,$

by assumption $D_{r,i}$ is $\epsilon_{\ell}$-SW, $\nabla_{i}g_{r,i}$ is
$(1+\epsilon_{R})$-SW, so the sum of the last three terms in the above
expression is $\epsilon_{\ell}(1+\epsilon_{R})$-SW.

Now for the first term, Let
$M=\mathbb{E}\left[\nabla_{i}g_{r,i}\nabla_{i}g_{s,i}|\mathbf{X}_{-v_{i}}\right]-\mathbb{E}(\nabla_{i}g_{r,i}\nabla_{i}g_{s,i})$
and
$M_{j}=\mathbb{E}(M|F_{j})-\mathbb{E}(M|F_{j-1})=\mathbb{E}(\nabla_{j}M|F_{j})$
for $1\leq j\leq\tilde{n}$. The main remaining task in the proof is to control
the tail of $M_{j}$.

For any $\ell_{q}$ norm $\|\cdot\|$ with $q\geq 1$,

$\displaystyle\|M_{j}\|=\|\mathbb{E}(\nabla_{j}M|F_{j})\|\leq\|\nabla_{j}M\|$
$\displaystyle=$
$\displaystyle\left\|\nabla_{j}\left[\mathbb{E}(\nabla_{i}g_{r,i}\nabla_{i}g_{s,i})|\mathbf{X}_{-v_{i}}\right]\right\|$
$\displaystyle\leq$
$\displaystyle\|\nabla_{j}(\nabla_{i}g_{s,i}\nabla_{i}g_{r,i})\|$
$\displaystyle=$
$\displaystyle\left\|\nabla_{j}\nabla_{j}g_{s,i}\nabla_{i}g_{s,i}+\nabla_{i}g_{s,i}(\mathbf{X}^{j})\nabla_{j}\nabla_{i}g_{s,i}\right\|$
$\displaystyle\leq$ $\displaystyle
2n^{-1/2}\epsilon_{\ell}(1+\epsilon_{R})\,.$

Then using Lemma A.3 we conclude $M$ is $\epsilon_{\ell}(1+\epsilon_{R})$-SW .
∎

###  C.4 Auxiliary lemmas

######  Lemma C.3 (Bounding $\bar{C}_{rs,i}-\mathbb{E}C_{rs,i}$).

Under the conditions in Theorem 5.1, $\bar{C}_{rs,i}-\mathbb{E}C_{rs,i}$ is
$\epsilon_{\ell}(1+\epsilon_{R})$-SW.

###### Proof of Lemma C.3.

For $j\neq i$, and function $f$ acting on $\mathbf{X}$, let
$\mathbf{X}^{j,-i}$ be the vector obtained by replacing $X_{j}$ in
$\mathbf{X}^{-i}$ with its iid copy $X_{j}^{\prime}$. Then by Jensen’s
inequality, we have, for $q\geq 1$

$\displaystyle\|\nabla_{j}\mathbb{E}(f|\mathbf{X}^{-i})\|_{q}^{q}=$
$\displaystyle\mathbb{E}\left\\{\mathbb{E}(f|\mathbf{X}^{-i})-\mathbb{E}\left[f(\mathbf{X}^{j,-i})|\mathbf{X}^{j,-i}\right]\right\\}^{q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\mathbb{E}\left[f(\mathbf{X})-f(\mathbf{X}^{j})|\mathbf{X}^{-i},X_{j}^{\prime}\right]\right\\}^{q}$
$\displaystyle\leq$ $\displaystyle\|\nabla_{j}f\|_{q}^{q}\,.$ (46)

Take $f$ to be $C_{rs,i}$, we have for $j\neq i$.

$\|\nabla_{j}\bar{C}_{rs,i}\|_{q}\leq\|\nabla_{j}C_{rs,i}\|_{q}\,.$

Next we control $\|\nabla_{j}C_{rs,i}\|_{q}$. By definition,

$\displaystyle\nabla_{j}C_{rs,i}=$
$\displaystyle\nabla_{j}\left[\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})\mathbb{E}(\nabla_{i}g_{s,i}|\mathbf{X})\right]$
$\displaystyle=$
$\displaystyle\left[\nabla_{j}\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})\right]\mathbb{E}(\nabla_{i}g_{s,i}|\mathbf{X})+\left[\nabla_{j}\mathbb{E}(\nabla_{i}g_{s,i}|\mathbf{X})\right]\left\\{\mathbb{E}\left[\nabla_{i}g_{r,i}(\mathbf{X}^{j})|F_{i}(\mathbf{X}^{j})\right]\right\\}\,.$
(47)

By Proposition A.1 and definition of $\epsilon_{\ell}$ and $\epsilon_{R}$,

$\displaystyle\nabla_{j}\nabla_{i}g_{r,i}=\nabla_{j}\nabla_{i}\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})+\sum_{k\notin
I_{v_{i}}}\nabla_{j}\nabla_{i}R(\mathbf{X}_{-v_{k}})~{}~{}\text{is}~{}~{}n^{-1/2}\epsilon_{\ell}\text{-SW}\,,$
(48)

and

$\nabla_{i}g_{r,i}=\nabla_{i}\ell_{r}(X_{i};\mathbf{X}_{-v_{i}})+\sum_{k\notin
I_{v_{i}}}\nabla_{i}R(\mathbf{X}_{-v_{k}})~{}~{}\text{is}~{}~{}(1+\epsilon_{R})\text{-SW}\,.$
(49)

Combining (48) and (49) with (47), and apply Proposition A.1 we conclude that
$\nabla_{j}C_{rs,i}$ is $n^{-1/2}\epsilon_{\ell}(1+\epsilon_{R})$-SW. For the
same reason as (46), we have $\nabla_{j}\bar{C}_{rs,i}$ is
$n^{-1/2}\epsilon_{\ell}(1+\epsilon_{R})$-SW. Then the desired result follows
from applying Lemma A.3 to the martingale sequence obtained by taking
conditional expectation of $\bar{C}_{rs,i}$ with respect to the filtration
$(F_{j,i}:1\leq 0\leq n,~{}j\neq i)$, where $F_{j,-i}$ is the sigma field
generated by $(X_{k}:1\leq k\leq j,~{}k\neq i)$. ∎

######  Lemma C.4 (Bounding $\mathbb{E}C_{rs,i}-\phi_{rs}$).

There exists a universal constant $c>0$ such that for all
$(r,s,i)\in[p]^{2}\times[n]$, $|\mathbb{E}C_{rs,i}-\phi_{rs}|\leq
c\epsilon_{\ell}$.

###### Proof of Lemma C.4.

Since $F_{i}$ is a sub $\sigma$-field of $\mathbf{X}$, and $\nabla_{j}g_{r,i}$
is centered, we have

$\displaystyle\mathbb{E}C_{rs,i}=\mathbb{E}\left[\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})\mathbb{E}(\nabla_{i}g_{s,i}|F_{i})\right]$
(50)

For $j\in\\{1,...,i-1\\}$, let $H_{j,i}$ be the $\sigma$-field generated by
$(X_{1},...,X_{j},X_{i})$. Define $H_{0,i}$ as the $\sigma$-field generated by
$X_{i}$, and $H_{-1,i}$ be the trivial $\sigma$-field. Then we can write the
martingale decomposition of $\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})$ as follows

$\displaystyle\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})=\sum_{j=0}^{i-1}\mathbb{E}(\nabla_{i}g_{r,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{r,i}|H_{j-1})\,.$

Apply this decomposition to both $\mathbb{E}(\nabla_{i}g_{r,i}|F_{i})$ and
$\mathbb{E}(\nabla_{i}g_{s,i}|F_{i})$ in (50), we get

$\displaystyle\mathbb{E}C_{rs,i}=$
$\displaystyle\mathbb{E}\sum_{j,k=0}^{i-1}\left[\mathbb{E}(\nabla_{i}g_{r,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{r,i}|H_{j-1})\right]\left[\mathbb{E}(\nabla_{i}g_{s,i}|H_{k})-\mathbb{E}(\nabla_{i}g_{s,i}|H_{k-1})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\sum_{j=0}^{i-1}\left[\mathbb{E}(\nabla_{i}g_{r,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{r,i}|H_{j-1})\right]\left[\mathbb{E}(\nabla_{i}g_{s,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{s,i}|H_{j-1})\right]\,,$
(51)

where the second equality holds because $H_{0}\subset H_{1}\subset
H_{2}\subset\cdots\subset H_{i-1}$ is a filtration.

When $j=0$,

$\displaystyle\mathbb{E}\left\\{\left[\mathbb{E}(\nabla_{i}g_{r,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{r,i}|H_{j-1})\right]\left[\mathbb{E}(\nabla_{i}g_{s,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{s,i}|H_{j-1})\right]\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\left[\mathbb{E}(g_{r,i}|X_{i})-\mathbb{E}g_{r,i}\right]\left[\mathbb{E}(g_{s,i}|X_{i})-\mathbb{E}g_{s,i}\right]\right\\}$
$\displaystyle=$ $\displaystyle\phi_{rs}\,.$ (52)

When $j>0$, we have, using Jensen’s inequality and (48),

$\displaystyle\|\mathbb{E}(\nabla_{i}g_{r,i}|H_{j})-\mathbb{E}(\nabla_{i}g_{r,i}|H_{j-1})\|_{2}=\|\mathbb{E}(\nabla_{j}\nabla_{i}g_{r,i}|H_{j})\|_{2}$
$\displaystyle\leq$ $\displaystyle\|\nabla_{j}\nabla_{i}g_{r,i}\|_{2}\lesssim
n^{-1/2}\epsilon_{\ell}\,.$ (53)

The claim follows by combining (52), (53) with (51). ∎

## References

  * Allen (1974) David M Allen.  The relationship between variable selection and data agumentation and a method for prediction.  _Technometrics_ , 16(1):125–127, 1974.
  * Arlot and Celisse (2010) Sylvain Arlot and Alain Celisse.  A survey of cross-validation procedures for model selection.  _Statistics surveys_ , 4:40–79, 2010.
  * Austern and Zhou (2020) Morgane Austern and Wenda Zhou.  Asymptotics of cross-validation.  _arXiv preprint arXiv:2001.11111_ , 2020.
  * Bates et al. (2021) Stephen Bates, Trevor Hastie, and Robert Tibshirani.  Cross-validation: what does it estimate and how well does it do it?  _arXiv preprint arXiv:2104.00673_ , 2021.
  * Bayle et al. (2020) Pierre Bayle, Alexandre Bayle, Lucas Janson, and Lester Mackey.  Cross-validation confidence intervals for test error.  _Advances in Neural Information Processing Systems_ , 33:16339–16350, 2020.
  * Bentkus (2005) Vidmantas Bentkus.  A lyapunov-type bound in rd.  _Theory of Probability & Its Applications_, 49(2):311–323, 2005.
  * Bousquet and Elisseeff (2002) Olivier Bousquet and André Elisseeff.  Stability and generalization.  _The Journal of Machine Learning Research_ , 2:499–526, 2002.
  * Celisse (2014) Alain Celisse.  Optimal cross-validation in density estimation with the ${L}^{2}$-loss.  _The Annals of Statistics_ , 42(5):1879–1910, 2014.
  * Chang et al. (2021) Jinyuan Chang, Xiaohui Chen, and Mingcong Wu.  Central limit theorems for high dimensional dependent data.  _arXiv preprint arXiv:2104.12929_ , 2021.
  * Chen et al. (2022) Qizhao Chen, Vasilis Syrgkanis, and Morgane Austern.  Debiased machine learning without sample-splitting for stable estimators.  _arXiv preprint arXiv:2206.01825_ , 2022.
  * Chen (2018) Xiaohui Chen.  Gaussian and bootstrap approximations for high-dimensional u-statistics and their applications.  _The Annals of Statistics_ , 46(2):642–678, 2018\.
  * Chernozhukov et al. (2013) Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  _The Annals of Statistics_ , 41(6):2786–2819, 2013.
  * Chernozhukov et al. (2015) Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Comparison and anti-concentration bounds for maxima of gaussian random vectors.  _Probability Theory and Related Fields_ , 162(1):47–70, 2015.
  * Chetverikov et al. (2016) Denis Chetverikov, Zhipeng Liao, and Victor Chernozhukov.  On cross-validated lasso in high dimensions.  _arXiv preprint arXiv:1605.02214_ , 2016.
  * Deng and Zhang (2020) Hang Deng and Cun-Hui Zhang.  Beyond gaussian approximation: Bootstrap for maxima of sums of independent random vectors.  _The Annals of Statistics_ , 48(6):3643–3671, 2020.
  * Efron and Tibshirani (1997) Bradley Efron and Robert Tibshirani.  Improvements on cross-validation: the $.632+$ bootstrap method.  _Journal of the American Statistical Association_ , 92(438):548–560, 1997.
  * Ferrari and Yang (2015) Davide Ferrari and Yuhong Yang.  Confidence sets for model selection by f-testing.  _Statistica Sinica_ , pages 1637–1658, 2015.
  * Gabcke (1979) Wolfgang Gabcke.  _Neue Herleitung und explizite Restabschätzung der Riemann-Siegel-Formel_.  PhD thesis, Georg-August-Universität Göttingen, 1979.
  * Geisser (1975) Seymour Geisser.  The predictive sample reuse method with applications.  _Journal of the American Statistical Association_ , 70(350):320–328, 1975.
  * Gunes and Bondell (2012) Funda Gunes and Howard D Bondell.  A confidence region approach to tuning for variable selection.  _Journal of Computational and Graphical Statistics_ , 21(2):295–314, 2012.
  * Hansen et al. (2011) Peter R Hansen, Asger Lunde, and James M Nason.  The model confidence set.  _Econometrica_ , 79(2):453–497, 2011.
  * Homrighausen and McDonald (2017) Darren Homrighausen and Daniel J McDonald.  Risk consistency of cross-validation with lasso-type procedures.  _Statistica Sinica_ , pages 1017–1036, 2017.
  * Jiang et al. (2008) Jiming Jiang, J Sunil Rao, Zhonghua Gu, and Thuan Nguyen.  Fence methods for mixed model selection.  _The Annals of Statistics_ , 36(4):1669–1692, 2008.
  * Kuchibhotla and Chakrabortty (2018) Arun Kumar Kuchibhotla and Abhishek Chakrabortty.  Moving beyond sub-gaussianity in high-dimensional statistics: Applications in covariance estimation and linear regression.  _arXiv preprint arXiv:1804.02605_ , 2018.
  * Kuchibhotla and Rinaldo (2020) Arun Kumar Kuchibhotla and Alessandro Rinaldo.  High-dimensional clt for sums of non-degenerate random vectors: $n^{-1/2}$-rate.  _arXiv preprint arXiv:2009.13673_ , 2020.
  * Kuchibhotla et al. (2021) Arun Kumar Kuchibhotla, Somabha Mukherjee, and Debapratim Banerjee.  High-dimensional clt: Improvements, non-uniform extensions and large deviations.  _Bernoulli_ , 27(1):192–217, 2021.
  * Kurisu et al. (2021) Daisuke Kurisu, Kengo Kato, and Xiaofeng Shao.  Gaussian approximation and spatially dependent wild bootstrap for high-dimensional spatial data.  _arXiv preprint arXiv:2103.10720_ , 2021.
  * Lei (2020) Jing Lei.  Cross-validation with confidence.  _Journal of the American Statistical Association_ , 115(532):1978–1997, 2020.
  * Lopes (2020) Miles E Lopes.  Central limit theorem and bootstrap approximation in high dimensions with near $1/\sqrt{n}$ rates.  _arXiv preprint arXiv:2009.06004_ , 2020.
  * Picard and Cook (1984) Richard R Picard and R Dennis Cook.  Cross-validation of regression models.  _Journal of the American Statistical Association_ , 79(387):575–583, 1984.
  * Rio (2009) Emmanuel Rio.  Moment inequalities for sums of dependent random variables under projective conditions.  _Journal of Theoretical Probability_ , 22(1):146–163, 2009.
  * Shao (1993) Jun Shao.  Linear model selection by cross-validation.  _Journal of the American statistical Association_ , 88(422):486–494, 1993.
  * Stone (1974) Mervyn Stone.  Cross-validatory choice and assessment of statistical predictions.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 111–147, 1974.
  * Stone (1977a) Mervyn Stone.  An asymptotic equivalence of choice of model by cross-validation and akaike’s criterion.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):44–47, 1977a.
  * Stone (1977b) Mervyn Stone.  Asymptotics for and against cross-validation.  _Biometrika_ , pages 29–35, 1977b.
  * Tibshirani and Tibshirani (2009) Ryan J Tibshirani and Robert Tibshirani.  A bias correction for the minimum error rate in cross-validation.  _The Annals of Applied Statistics_ , 3(2):822–829, 2009.
  * Tsiatis (2006) Anastasios A Tsiatis.  Semiparametric theory and missing data.  2006\.
  * Vladimirova et al. (2020) Mariia Vladimirova, Stéphane Girard, Hien Nguyen, and Julyan Arbel.  Sub-weibull distributions: Generalizing sub-gaussian and sub-exponential properties to heavier tailed distributions.  _Stat_ , 9(1):e318, 2020.
  * Yang (2007) Yuhong Yang.  Consistency of cross validation for comparing regression procedures.  _The Annals of Statistics_ , 35(6):2450–2473, 2007.
  * Yu and Feng (2014) Yi Yu and Yang Feng.  Modified cross-validation for penalized high-dimensional linear regression models.  _Journal of Computational and Graphical Statistics_ , 23(4):1009–1027, 2014.
  * Zhang (1993) Ping Zhang.  Model selection via multifold cross validation.  _The Annals of Statistics_ , pages 299–313, 1993.