A Modern Theory of Cross-validation Through the Lens of Stability

Risk Estimation and Model Selection Properties of CV

In the theoretical study of cross-validation (CV), two of the most extensively examined properties are risk estimation and model selection.

Risk estimation refers to the ability to quantify the predictive risk of a fitted model on future, unseen data. Alternatively, one might consider the estimation risk relative to a true underlying parameter, which requires assuming the existence of a true model for the datagenerating process. In this article, we focus on predictive risk estimation, also known as generalization error estimation in the machine learning literature.

General results of CV risk consistency mostly rely on empirical process theory (Devroye, 1988;Vaart et al., 2006;Lecué and Mitchell, 2012). Results in specific regression settings include nonparametric regression and density estimation (Devroye and Lugosi, 2001;Györfi et al., 2006) and the Lasso (Homrighausen and McDonald, 2017).

Model selection concerns the ability to identify the best, or nearly the best, model from a set of candidates. This task is often divided into two variants: (i) identifying the true model that generated the data, and (ii) selecting the model with the best predictive performance, regardless of whether a true model exists within the candidate set. These are referred to as model selection for identification and model selection for estimation, respectively, in Arlot and Celisse (2010). Our focus in this article is primarily on the latter, as it requires fewer assumptions and subsumes the identification problem as a special case when the true model yields the smallest predictive risk.

The model selection properties of CV are generally more challenging to analyze. Early theoretical results indicate that conventional versions of CV-such as leave-one-out and K-fold CV-can be inconsistent even in simple linear regression settings (Zhang, 1993;Shao, 1993). These inconsistencies were later extended to more general frameworks in Yang (2007). The underlying reason is intuitive: standard CV procedures struggle to distinguish between the best parametric model and a slightly overfitting alternative, as both models tend to have estimation errors on the order of 1/ √ n. Evaluating their performance on a sample of size O(n) leads to a constant-order probability of selection error. To address this limitation, CV must be applied in a nonstandard way using a vanishing fraction of the data for model fitting and the remainder for risk evaluation.

Additionally, it has been widely observed that CV often favors models with slight under-smoothing. While such models may yield near-optimal predictive accuracy, they frequently suffer from reduced interpretability. To counteract this under-regularization tendency, several heuristic adjustments have been proposed, including the ".632 bootstrap" (Efron and Tibshirani, 1997) and the "1-standard-deviation rule" (Tibshirani and Tibshirani, 2009). Other examples and variants include Yu and Feng (2014) and Feng and Yu (2019).

Outline of This Article

In Section 2, we introduce the foundational setup of cross-validation along with the relevant stability conditions, and present the risk consistency and concentration results established by Bousquet and Elisseeff (2002). In addition, we provide a collection of general probabilistic tools concerning sub-Weibull random variables and their concentration properties. These tools can be used to extend the results of Bousquet and Elisseeff (2002) to settings involving more general, potentially unbounded loss functions. As such, they will serve as a key component in the analysis presented in subsequent sections.

In Section 3, we present results on the model selection consistency of cross-validation. We begin with batch mode, revisiting the consistency result of Yang (2007), and provide a simplified proof. As an extension, we introduce an online model selection framework and propose an online variant of cross-validation (Zhang and Lei, 2023). While this online method resembles classical leave-one-out cross-validation-which is known to be inconsistent in the batch setting-it can, somewhat unexpectedly, achieve consistent model selection under appropriate stability conditions.

Section 4 focuses on central limit theorems (CLTs) for cross-validation risk estimates. We begin by presenting the CLT results under two different centering schemes, originally developed by Bayle et al. (2020) and Austern and Zhou (2024), while streamlining and simplifying some of the proofs. In our treatment of the deterministic centering results, several technical conditions from the original work are either weakened or removed. We also unify the two centering approaches by identifying a condition under which they become equivalent. Finally, we discuss the extensions in Kissel and Lei (2022), which develop high-dimensional Gaussian approximation results in the spirit of Chernozhukov et al. (2013) and Chernozhukov et al. (2022).

Section 5 presents several important applications of the theoretical developments in Section 2 and Section 4 across various contexts. The first application introduces a model confidence set based on Gaussian comparisons of CV risk estimates, representing a methodological advancement that emphasizes constructing confidence intervals for the differences in model performance, rather than for each model individually. Building on this framework, we establish a model selection consistency result for selecting the most parsimonious model within the K-fold CV model confidence set, offering a principled correction for the well-documented overfitting tendency of CV. In Section 5.3, we leverage the CLT results to propose a novel method for high-dimensional mean testing, providing a simple alternative to selective inference techniques. Extending this line of inquiry, Section 5.4 presents a procedure for constructing confidence sets for the location of the minimum entry in a vector observed with noise. Finally, in Section 5.5, we establish an asymptotic coverage guarantee for the cross-conformal prediction method under appropriate stability conditions.

Section 6 presents several examples of stable estimators. Sections 6.1 and 6.2 focus on the stochastic gradient descent (SGD) algorithm under general regularity conditions, providing batch-mode stability results needed for theorems in Sections 2 and 4. Section 6.3 examines the stability of the difference between two competing estimators in a prototypical nonparametric regression setting. Finally, Section 6.4 explores the stability properties of various SGD variants in the online setting, as required for the online model selection consistency results in Section 3.3.

Other Related Topics

This article focuses on recent developments in cross-validation theory, particularly through the lens of stability conditions. While cross-validation is a widely used tool for model selection and risk estimation, it is by no means the only approach. Classical alternatives include regularized likelihood and data compression-based methods, such as used notation and symbols for the reader's convenience. The notation specifically related to the sample-splitting scheme of cross-validation is introduced in Section 2.1.

• We use (X i : i g 0) to denote an infinite sequence of i.i.d.

(independent and identically distributed) sample points, and D n = (X i : 1 f i f n) the dataset consisting of the first n sample points.

• X denotes the space of a single sample point, and X ¹n the n-fold product space, so that D n ∈ X ¹n .

• For a positive integer n, [n] := {1, 2, ..., n}.

• E X (•) means taking expectation over the randomness of X, conditioning on everything else.

• R d denotes the d-dimensional Euclidean space.

• ∥ • ∥ q denotes the L q -norm of a random variable. More concretely, ∥X∥ q = (E|X| q ) 1/q if 1 f q < ∞, and ∥ • ∥ ∞ is the essential supremum of a random variable. For a function h acting on a random variable X, ∥h∥ q,X = ∥h(X)∥ q .

• ∥ • ∥ denotes the Euclidean, Hilbert, or Banach-space norm, depending on the context, and ∥ • ∥ op denotes the operator norm of a matrix or linear operator.

• For two non-random sequences a n and b n such that b n g 0, , where c n is a finite-degree polynomial of log n and perhaps the logarithm of other model quantities such as the dimensionality of the data and/or the number of models being compared.

• P → denotes convergence in probability.

• ⇝ denotes convergence in distribution.

• For two random variables U , V defined on the same space and q g 1, d Wq (U, V ) denotes the Wasserstein-q distance between U and V .

• For two sequences of possibly random variables U n and V n such that V n g 0,

• For a positive integer m: S m-1 denotes the (m -1)-dimensional simplex in R m : {(w 1 , ..., w m ) ∈ R m : w r g 0, m r=1 w r = 1}.

• For a smooth function h defined on a Euclidean or Hilbert space, ḣ(•) and ḧ(•) denote the gradient and Hessian of h, respectively.

In applications of cross-validation, the estimation of f may be carried out on datasets with different sizes. To keep track of the potentially varying fitting sample size, we consider a collection of estimators indexed by different sample sizes. Let f : ∞ n=1 X ¹n → F be an estimation procedure that maps a sample of any size to the parameter space. We assume that f is permutation invariant (or symmetric).

Definition 2.1 (Symmetric Estimator). We say an estimator f :

Let ℓ : X × F → [-∞, ∞] be a loss function used to evaluate the quality of model fit. This notation covers both supervised learning and unsupervised learning.

Examples in supervised learning include regression and classification.

• Classification: X = (Y, Z) ∈ {0, 1} × R p , and f is a function mapping Z to [0, 1] that predicts the conditional probability of Y = 1 given Z. A common choice of ℓ in classification is the negative logistic likelihood ℓ ((y, z)

Examples in unsupervised learning include principal components analysis and mixture models.

• Principal components analysis (PCA):

where ∥ • ∥ denotes the Euclidean norm.

• Mixture models: X ∈ R p , f is a density function (e.g., in a reproducing kernel Hilbert space or a Gaussian mixture), and

The overall quality of a fitted model f under the loss function ℓ

Here the notation E X (•) means taking conditional expectation over the randomness of X only. We will often consider the risk of a fitted model f = f (D n ), which is obtained by applying the fitting procedure f to the input dataset D n . To simplify notation, we write

when there is no ambiguity on the fitting procedure.

Under such a setting, there are two natural quantities of interest. The first is R(D n ), the risk of the estimated parameter obtained by applying a particular fitting procedure f to the fitting dataset D n . This is a random quantity, because it depends on the random fitting dataset D n . Another quantity of interest is the population version of R(D n ):

which is a deterministic quantity that depends only on the fitting procedure f , the data distribution P X , and the sample size n.

Neither R(D n ) nor µ n can be straightforwardly estimated without additional assumptions, such as uniform convergence over the function class or, as we will see later this article, stability of f .

A common challenge in estimating the risk R(D n ) or the average risk µ n is to handle simultaneously the randomness in f (D n ) (estimation) and in approximating the risk of f (D n ) (evaluation). The cross-validation technique does this by splitting the dataset to create a holdout sample for empirical evaluation of the risk, and further reducing the estimation variability by repeated sample splitting. In this article, we focus on the commonly used K-fold cross-validation estimate, described as follows.

1. Fix K such that n te := n/K, the test subsample size, is an integer. 2Let n tr = nn te be the training subsample size.

2. For k = 1, ..., K, let I n,k = {n te (k -1) + 1, ..., n te k} be the indices of the kth fold of the data, D n,k = {X i : i ∈ I n,k } be the corresponding kth fold of the data, and

We may drop the dependence on n and K and simply write Rcv . We will consider the following two cases.

• K = O(1). This corresponds to the commonly used 5-fold or 10-fold CV.

The most typical case is K = n, which corresponds to the leave-one-out (LOO) CV.

While we might intuitively expect Rcv,n,K to approximate R(D n ) or µ n , it is important to note that, by construction, Rcv represents the risk of models trained on samples of size n tr rather than n, since each term in (2.2) is a loss evaluated at a parameter estimate fitted on a dataset of size n tr .

With a sample size n tr , the counterpart of µ n is simply µ ntr . But for the risk of the fitted model R(D), the CV procedure creates a total of K versions at the sample size n tr : {R(D n,-k ) : 1 f k f K}. Hence the counterpart of R(D n ) at sample size n tr targeted by the CV risk estimate Rcv,n,K in (2.2) should be their average:

It is straightforward to check that Rcv,n,K is unbiased for Rcv,n,K and µ ntr , in the sense that E( Rcv,n,K ) = E Rcv,n,K = µ ntr . The following picture describes the chain of approximation among the quantities of interest. Rcv,n,K

--→ Rcv,n,K

where the approximation in each arrow has its own interpretation.

(a) : the sampling variability in the evaluation of empirical risk;

(b,b ′ ) : the concentration of the risk of estimated models around the expected value;

As we will see in this article, all three approximations can be established using appropriate stability conditions.

Stability Conditions

Now we describe the specific notions of stability to be considered in this article. Intuitively, stability refers to a form of continuity of the fitting procedure f , such that f (D) ≈ f (D ′ ) whenever D ≈ D ′ . Because f is random, stability must be stated in a stochastic manner, such as via a norm or tail-probability bound on a random variable. We will consider two particular types of pairs (D, D ′ ) that are close to each other. In the first, D ′ is obtained by removing one element in D. Recall that D n = {X 1 , ..., X n } consists of a collection of i.i.d. sample points. We can use D n-1 to denote the dataset obtained by removing one sample point from D n . This leads to the leave-one-out stability.

a function that maps a dataset of any size and an auxiliary input in [0, 1] to a real number. For q ∈ [1, ∞], define the Leave-One-Out L q stability (LOO-L q ) of h at sample size n as

In this definition, we allow h to involve additional randomness besides the dataset D n , through U ∼ Unif(0, 1). There is no loss of generality by assuming U ∼ Unif(0, 1) since we can use U to generate random vectors of arbitrary distribution. In particular, this definition covers the loss function ℓ(X 0 , D n ), which involves the random evaluation point X 0 . For example, in Theorem 2.1 below, we will consider S loo q,n (ℓ(X 0 , •)), which equals ∥ℓ(X 0 , D n )ℓ(X 0 , D n-1 )∥ q , where the ℓ q norm is with respect to the joint randomness of the n + 1 independent sample points X 0 , ..., X n .

In another type of neighboring pairs (D, D ′ ), D ′ is obtained by replacing one entry of D by an i.i.d. copy. For each i ∈ {0, 1, ..., n} we let X ′ i be an i.i.d. copy of X i , independent of everything else. For any i ∈ [n] and function h : ∞ n=1 X ¹n → R 1 which may involve additional randomness, we define the perturb-one difference operator

where D i n is obtained by replacing the ith entry of D n by X ′ i . The perturb-one stability is defined as follows.

Definition 2.3 (Perturb-One Stability in L q -Norm). Let h : [0, 1] × ∞ n=1 X ¹n → R 1 be a function that maps a dataset of any size and an auxiliary input in [0, 1] to a real number. For q ∈ [1, ∞], define the Perturb-One L q (PO-L q ) stability of h at sample size n as

where U ∼ Unif(0, 1) is independent of D n .

Although from a practical point of view the difference is not substantial, mathematically the Perturb-One stability is weaker than the Leave-One-Out stability. Roughly speaking, Perturb-One Stability only requires that each individual sample point contributes a small proportion to the output, while Leave-One-Out stability also requires continuity with respect to sample size. In particular, a simple application of triangle inequality implies that the PO-L q stability is no larger than twice the LOO-L q stability . Proposition 2.1 (LOO Stability ⇒ PO Stability). S po q,n (h) f 2S loo q,n (h).

Conversely, PO-L q stability does not necessarily imply LOO-L q stability. For example, consider the following h which does not use the auxiliary input:

Clearly, h satisfies PO-L q stability at rate O(1/n) whenever ∥X 1 ∥ q < ∞ but h does not satisfy LOO-L q stability at the same rate.

Risk Consistency of LOOCV

Now we present a risk consistency result for the leave-one-out cross-validation under the LOO-L q stability condition of the loss function.

Theorem 2.1. Let 1 f q f r f ∞ be such that 1/q + 1/r = 1. Assume the estimator f is symmetric, the loss function has bounded L r -norm:

The proof of Theorem 2.1 is quite elementary and elegant. The key step is a "label-swapping trick" that makes use of the symmetry of data, which further leads to the application of the stability condition.

Proof. The proof is based on a direct expansion of the squared difference between Rcv,n,n and R(D n ):

Next, we expand each of the three terms.

For the first term, we have

For the second term, we have, by definition of Rcv,n,n , symmetry of f and exchangeability of sample points,

(2.8)

For the third term, note that the assumptions 1 f q f r f ∞ and 1/q + 1/r = 1 imply that r g 2, so

where the first inequality uses the fact that |ℓ

, and symmetry of f and the data points. Putting (2.7), (2.8), and (2.9) in (2.6), we obtain

(2.10)

Recall that D i n is the perturb-one dataset obtained by replacing the ith entry of D n by an i.i.d. copy X i . We have the following key equation (the "label-swapping trick"). (2.11) where the first equality holds due to exchangeability between X 1 and X ′ 1 , and the second equality holds due to exchangeability between X ′ 0 and X 1 . As a result,

By the LOO-L q stability of ℓ(X 0 , D n ) we have

and by triangle inequality

Thus Hölder's inequality implies that

A similar treatment of swapping X 2 , X ′ 2 , and X 0 in the term

leads to the same upper bound of 3Bϵ n . Therefore the RHS of (2.10) is bounded by 6Bϵ n , which leads to the claimed result.

The "label-swapping trick" (2.11) is a special case of a more general technique called the "free nabla trick" (Theorem 4.11). Besides its crucial role in proving CLTs for CV risk estimates, the free nabla trick also allows us to extend Theorem 2.1 by providing upper bounds on | Rcv,n,K -Rcv,n,K | for general values of K under a Perturb-One stability condition. An application of this upper bound is presented in Section 5.5, in the context of conformal prediction.

Concentration of CV Risk Estimate

Theorem 2.1 gives a general consistency result for the LOO-CV risk estimates. In this subsection, we study the concentration bounds for such estimates. The main challenge is that many mean-zero terms that disappear under expectation or are combined in the proof of Theorem 2.1, such as those in (2.7) and (2.8), have complicated correlations that contribute to the random fluctuation of Rcv,n,n around R(D n ). Therefore, we will need stronger stability conditions to establish exponential concentration inequalities for CV risk estimates.

The starting point is the following well-known bounded-difference concentration result.

Theorem 2.2 (McDiarmid's Inequality). For any function h that maps D n to a real number, if there exist constants c 1 , ..., c n such that |∇ i h| f c i almost surely for each i ∈ [n], then

The bounded difference condition in McDiarmid's inequality corresponds to the PO-L ∞ bound. Therefore, we consider the following perturb-one stability in the L ∞ -norm condition on the loss function of the fitted model:

(2.12)

Again, in this subsection, we focus on the LOO-CV, where K = n and D n,-i denotes the dataset obtained by removing the ith entry.

Under condition (2.12), for any

Therefore,

then the triangle inequality implies that Rcv,n,n satisfies the bounded difference assumption in McDiarmid's inequality with c i = 2Bn -1 + ϵ n . Therefore, we have the following result. Proposition 2.2. Assume ∥ℓ(X 0

1. If the loss function ℓ(•, •) and fitting procedure f satisfy (2.12), we have, for any t > 0,

2. If the loss function ℓ(•, •) and fitting procedure f satisfy S loo ∞,n (ℓ(X 0 , f (•))) f ϵ n , then we have, for any t > 0,

and

In other words, Rcv,n,n concentrates around µ n at an O(n -1/2 + n 1/2 ϵ n ) scale and a sub-Gaussian tail. This result can be extended to general K-fold cross-validation to bound | Rcv,n,Kµ ntr |.

Remark 2.1. The concentration bounds based on McDiarmid's inequality presented in Proposition 2.2 have a much stronger requirement on the stability bound ϵ n than in the consistency result Theorem 2.1. In the consistency result, Theorem 2.1 implies that Rcv,n,n -R(D n ) = o P (1) as long as ϵ n = o(1). In contrast, the concentration results in Proposition 2.2 require ϵ n = o(n -1/2 ) in order to have a vanishing bound. Such a difference reflects the subtle nature of the dependence between the cross-validated losses, which is a major challenge in uncertainty quantification for CV model evaluation and comparison.

Extensions to Sub-Weibull Stability

The concentration result Proposition 2.2 requires the perturb-one loss stability to be uniformly bounded, which may be too restrictive, as it precludes many common scenarios involving unbounded support, such as linear regression with Gaussian noise under the squared loss. Unlike the consistency result, which can avoid the uniform bound condition by leveraging different combinations of conjugate indices in Hölder's inequality, relaxing the L ∞ -norm stability assumption in the concentration result is less straightforward, as it typically requires extending the bounded-difference inequality to the case of infinite support. Here, we present an approach based on sub-Weibull tail conditions. Some results and properties of sub-Weibull random variables may be of more general interest.

Definition 2.4 (Sub-Weibull Random Variables). For a pair of fixed positive real numbers (», ³), we say a random variable Z is (», ³)-sub-Weibull (or (», ³)-SW) if

The definition of sub-Weibull random variables is concerned with the tail behavior. It involves two parameters: a scaling parameter » and an exponent parameter ³. A larger ³ corresponds to a heavier tail. In particular, a random variable Z is (», ³)-SW if and only if (Z/») 1/³ has a sub-exponential tail. The sub-Gaussian tail condition with standard deviation » corresponds to (», 1/2)-SW, and the subexponential condition corresponds to (», 1)-SW. It may also be helpful to interpret the extreme case of ³ = 0 as uniformly bounded.

Like sub-Gaussian and sub-exponential tail conditions, we can bound the L q -norms of a sub-Weibull random variable proportionally to the scale parameter.

Proposition 2.3 (Sub-Weibull Condition in L q -Norm). If a random variable Z satisfies (2.13) then

In most cases, we will only focus on the rates of the scale parameters of sub-Weibull residual terms and usually do not distinguish two scaling parameters differing by a constant factor. Therefore, we will use the two sub-Weibull conditions (2.13) and (2.14) interchangeably.

Proof of Proposition 2.3. For q g 1,

where Γ(z) = ∞ 0 t z-1 e -t dt is the gamma function. So for q g 1 we have

If 1 f k < ³, then by Jensen's inequality

So in either case we have

Now using k! g (k/e) k for all k g 1 we have

. We have by Markov's inequality applied to e (¼|Z|) 1/α for any t > 0

The sum and product of sub-Weibull random variables are also sub-Weibull, with the scale parameter preserved through the sum and product, and a slightly different exponent parameter.

Proof. For the first part, take r, s > 1 such that 1/r + 1/s = 1

where the last inequality follows from minimizing r ³ 1 s ³ 2 over all (r, s) such that 1/r + 1/s = 1, which is achieved at r

For the second part,

The third part is a direct consequence of Jensen's inequality.

To prove the fourth part, without loss of generality we assume Z 1 g 0 almost surely. Let q g 1 be a positive number. Consider two cases. First, if 1 f q < 1/´, then

´.

Second, if q g 1/´then

´.

Lemma 2.3 (Maximal Inequality).

Let Z 1 , ..., Z m be (», ³)-SW random variables. Then there exists a constant c ³ depending only on ³, such that max Vaart and Wellner (1996) implies that, for another c ′ ³ , max i Zi /(c ′ ³ log m) is sub-exponential, i.e., max i Zi is (c ′′ ³ log m, 1)-SW for another constant c ′′ ³ . Therefore, the claimed result follows from Part (4) of Proposition 2.4. Now we return to the bounded-difference concentration inequality. For a function h acting on a dataset D n consisting of i.i.d. sample points, we want to generalize McDiarmid's inequality by establishing concentration of h(D n ) around its mean value under sub-Weibull tail conditions of the perturb-one differences ∇ i h(D n ).

We begin with the standard martingale decomposition of h(D n ). With the convention D 0 = ∅, we have

(2.15)

Similar to the proof of McDiarmid's inequality, the summand {E(∇ i h(D n )|D i ) : 1 f i f n} is a martingale increment sequence adapted to the filtration of Ã-fields generated by {D i : 1 f i f n}. On the other hand, Part (3) of Proposition 2.4 implies that if all the differences

are also sub-Weibull. However, a direct application of Part (2) of Proposition 2.4 may not yield a useful upper bound, since the resulting scaling parameter may involve n in an undesirable way.

The following theorem offers a good upper bound for sums of martingale increments with sub-Weibull tails.

Theorem 2.4 (Rio's Inequality). Let q g 2 and (S n ) ng0 be a sequence of random variables with finite L q -norm. Let

We refer readers to Rio (2009) for the proof of Theorem 2.4. Using Rio's inequality we can prove a bounded difference inequality under sub-Weibull conditions. Theorem 2.5 (Sub-Weibull Difference Concentration). Let h be a realvalued function acting on X ¹n and D n a set of n independent sample points in

Using Rio's inequality and the sub-Weibull condition of

The claim holds since (2/q) f 2.

This result shows that h(D n ) deviates from the mean value on the correct order of » 2 1 + ... + » 2 n . The original McDiarmid's inequality can be stated as a special case of Theorem 2.5 with ³ i = 0 and c i = » i . Now we have a concentration inequality for cross-validation risk estimates under a relaxed version of perturb-one stability. Definition 2.5 (Perturb-One Sub-Weibull Stability). Let h : [0, 1] × ∞ n=1 X ¹n → R 1 be a function that maps a dataset of any size and an auxiliary input in [0, 1] to a real number. We say h satisfies the Perturb-One Sub-Weibull (PO-SW) stability with parameter (» n , ³)

Theorem 2.6 (Concentration of LOO-CV Risk Under Sub-Weibull Stability). Assuming i.i.d. data, if ℓ(X 0 , •) satisfies PO-SW stability with parameter (ϵ n , ³) at sample size n-1, and ℓ(X

Proof of Theorem 2.6. By definition

Using the triangle inequality and the sub-Weibull moment bound we have for q g 1

Then the claimed result follows from Theorem 2.5.

Theorem 2.6 shows that the LOO-CV risk estimate concentrates around the mean value at essentially the same rate as in Proposition 2.2, with the same dependence, up to constant factors, on (B, ϵ n ). We will revisit sub-Weibull stability in Section 4.3 when proving highdimensional Gaussian comparison for cross-validated risks.

Similar to Theorem 2.1 and Proposition 2.2, the concentration result under sub-Weibull stability in Theorem 2.6 can be extended to bound | Rcv,n,Kµ ntr | in general K-fold CV.

Bibliographic Notes

The point estimation target of cross-validation is studied in Bates et al. (2024). The random target Rcv,n,K also appeared as the "random centering" in central limit theorems of cross-validation (Austern and Zhou, 2024;Bayle et al., 2020), and was used as an intermediate inference target for the construction of model confidence sets using cross-validation (Lei, 2020;Kissel and Lei, 2022).

The main consistency and concentration results, Theorem 2.1 and Proposition 2.2, are based on Bousquet and Elisseeff (2002). Other important contributions to the study of stability and generalization error include Kutin and Niyogi (2002), who developed similar generalization error bounds in parallel to Bousquet andElisseeff (2002), andShalev-Shwartz et al. (2010) and Wang et al. (2016), who established the connection between stability and generalization error in a broader context of empirical risk minimization. More recently, Bousquet et al. (2020) obtained sharper generalization error bounds under uniform stability.

The stability conditions introduced in Definitions 2.2, 2.3 and 2.5 are most directly relevant to the theoretical study of cross-validation. Alternative stability conditions have been considered in the study of other predictive inference tools. For example, Barber et al. (2021) and Liang and Barber (2023) used a notion of "out-of-sample stability" and a high-probability version of leave-one-out stability to analyze several variants of conformal prediction methods.

More discussion and details about sub-Weibull random variables can be found in Vladimirova et al. (2020) and Kuchibhotla and Chakrabortty (2022). The combination of L q -norm upper bound for martingale increments (Theorem 2.4) and sub-Weibull difference (Theorem 2.5) originated from personal communication with Arun Kumar Kuchibhotla, and a similar form has appeared in Appendix B of Kuchibhotla et al. (2023). Alternative approaches to proving extensions of McDiarmid's inequality using high-probability bounded differences can be found in Kutin (2002) and Combes (2024).

In this article, we consider generic cross-validation methods applied to an arbitrary learning algorithm f under stability conditions. Parameter estimation and risk consistency results for cross-validation have been established for specific estimators under more specific conditions, such as the lasso (Homrighausen and McDonald, 2017;Chetverikov et al., 2021), ridge regression (Patil et al., 2021), and gradient-descent least squares regression (Patil et al., 2024).

In this section, we consider another key application of cross-validation: model and tuning parameter selection. The overarching goal of inference in this context is to choose, from a finite collection of candidate models, one that best fits the data.

We adopt the perspective of model selection for estimation (see Section 1.2 and references therein), in which we neither assume that the true model is included in the candidate set nor require the existence of a true model at all. Instead, each candidate model is viewed as a specific procedure for producing a function that predicts future data points, and a good model is one that yields accurate predictions. In other words, in our context of model selection, each "model" simply corresponds to a particular estimator of an underlying parameter of the data-generating distribution, which further leads to a procedure for predicting future, unseen data.

This framework imposes minimal assumptions on both the datagenerating process and the model-fitting procedures, capturing the flexibility and broad applicability of cross-validation. All that is required is an i.i.d. sample, a loss function, and a collection of fitting procedures. Notably, this setup also encompasses tuning parameter selection, and for simplicity, we will refer to both tasks under the umbrella of "model selection." 423

Using the notation and terminology in Section 2.1, assume we have a collection of symmetric estimators { f (r) : 1 f r f m}, where for each r ∈ M = {1, ..., m}, f (r) : ∞ n=1 X ¹n → F maps a dataset D n to a parameter f (r) (D n ). We would like to pick the r that gives the best predictive

The model selection consistency question we want to address is: When does rcv,n,K coincide with the "best candidate model"?

Pairwise Selection Consistency

Intuitively, for the cross-validation model selection procedure described above to successfully identify the "best model," two conditions must be met. First, the best model should exhibit superior predictive accuracy at the population level-a quality gap must exist. Second, the crossvalidation risk estimates must concentrate around their corresponding population risks, ensuring that the variability in the estimates does not obscure the true differences between models.

In this section, we follow the setting in Yang (2007), where two further simplifications are made. The first simplification is to consider only two models: M = {1, 2}. This pairwise setting facilitates the articulation of some key intuitions, and the results can be easily extended to any finite candidate set such that m = |M| = O(1). The second is a regression setting with the squared loss:

where the noise ϵ has a finite second moment and E(ϵ|Z) = 0.

Under the regression model (3.1) and the squared loss, it is often more convenient to study the excess risk of each fitted candidate model. For each r ∈ M define the excess risk

where ¶ (r) n := f (r) (D n )f * is the estimation error of the regression function.

Without loss of generality, assume that the first model (r = 1) is better than the other. Intuitively, this requires Ψ 2,n to be larger than Ψ 1,n . In order to give a precise mathematical notion of "model 1 is better than model 2", we need to take into account the following issues:

(i) The quantities Ψ 1,n and Ψ 2,n are random, so the comparison between the models must be stated in a stochastic way.

(ii) The separation between Ψ 2,n and Ψ 1,n needs to be large enough to ensure consistency.

With these issues in mind, the following "risk dominance" condition seems natural.

Assumption 3.1 (Risk Dominance). The two sequences of random variables (Ψ r,n : n g 1) satisfy Ψ 2,n > 0 almost surely and

where (x) + := max(x, 0), with the convention 0/0 = ∞.

In other words, the risk dominance condition requires that, with high probability, Ψ 2,n exceeds Ψ 1,n by at least a constant proportion of Ψ 1,n . This definition precludes the scenario of Ψ 2,n /Ψ 1,n → 1 in probability, where the two risk sequences become asymptotically equivalent.

We can derive some useful consequences from the risk dominance condition.

Proposition 3.1. (a) An equivalent statement of the risk dominance condition (3.3) is that for any ε > 0 there exists c ε > 0 such that

Proof of Proposition 3.1. Part (a) follows from a direct rearrangement of the tail probability of

The claim follows by picking c = c ε large enough in the above equation so that the probability is no larger than ε for all sufficiently large n. Part (b) follows from the fact that Ψ 2,n f Ψ 1,n implies

Part (c) follows from the following:

Now we state the following pairwise comparison consistency result for K-fold cross-validation, originally due to Yang (2007). Recall that in K-fold cross-validation, the training sample size is n tr = n(1 -1/K) and the test sample size is n te = nn tr = n/K. Theorem 3.1. Assume D n consists of i.i.d. data from the regression model (3.1) with ϵ independent of Z, Eϵ = 0, and

2 and all sample sizes n. If the excess risks Ψ r,n (r = 1, 2) satisfy Assumption 3.1 and

then, for K-fold CV risk estimates with K = O(1), we have

Before proving Theorem 3.1, we briefly discuss the additional conditions stated above. The boundedness condition ∥ f

is more or less standard, and is needed to establish rates of convergence for empirical approximations of the excess risks using the test sample. It can be relaxed to

and is not too heavy-tailed. The second condition n te Ψ 2,ntr P → ∞ is more nontrivial and specifies a requirement on the test and training split ratio determined by the convergence rate of the inferior model. Intuitively, the test sample evaluates the empirical risk with a standard error rate Ψ 2,ntr / √ n te . Then we need Ψ 2,ntr -Ψ 1,ntr , which is asymptotically the same order as Ψ 2,ntr , to dominate this risk evaluation error. We illustrate the effect of split ratio in a simple example below, which also reflects the tightness of the condition (3.4).

Further extensions, such as other sample splitting schemes and diverging K, and the implications of Theorem 3.1 in traditional parametric and nonparametric regression settings, are discussed in Remarks 3.1 and 3.2.

Example 3.1. Consider the following comparison of two models:

and f (2) are constant functions, where f (1) uses the constant 0 and f (2) uses the sample average in D n . In this case, f * = f (1) . By construction, Ψ 1,n = 0 for all n and Ψ 2,n = (n -1 n i=1 Y i ) 2 ≍ 1/n. Let Rss,r (r = 1, 2) be the contribution to the CV risk estimate from a single test sample fold, defined in (3.5) below. It is straightforward to verify that n te ( Rss,2 -Rss,1 ) = -2 Ȳte Ȳtr + Ȳ 2 tr , and P( Rss,2 -Rss,1 f 0) → 0 if and only if (3.4) holds.

Proof of Theorem 3.1. We drop the dependence on n whenever there is no confusion. We focus on such a single sample split with D tr and D te denoting the fitting and validation data, with corresponding index sets I tr , I te .

For r = 1, 2, use simplified notation

Observe that

For r = 1, 2, consider the sample-split risk estimate Rss,r defined as

(3.5) Further define,

Our argument is based on the following two facts about À i . Fact (i): Conditioning on D tr , (À i : i ∈ I te ) are i.i.d. with mean zero. Fact (ii), Rss,1 g Rss,2 if and only if

The plan is to use Tchebyshev's inequality to bound the probability of the event in (3.6). We first need to control the conditional variance of À i given D tr .

where the first inequality uses E(ϵ i |Z i ) = 0 and the boundedness of ¶ 1 + ¶ 2 , and the second inequality follows from Cauchy-Schwartz.

Let E be the event that Ψ 2,ntr -Ψ 1,ntr > 0. Then

Using Tchebyshev's inequality conditionally given D tr we have

where the second inequality follows from the conditional variance bound (3.7), the O P (•) term follows from part (c) of Proposition 3.1, and the o P (•) term follows from the theorem assumption (3.4).

As a result,

where the exchange of expectation and limit is guaranteed by the boundedness of P( Rss,1 > Rss,2 |D tr )1 E . Finally, part (b) of Proposition 3.1 guarantees that P(E c ) = o(1). Combining this with (3.9) we conclude that the RHS of (3.8) is o(1).

Remark 3.1 (Other variants of K-fold CV). Theorem 3.1 also applies to other CV variants than the standard K-fold CV, such as repeated random split, where I tr is repeatedly generated by uniformly sampling n tr indices from [n] without replacement, and I te = [n] \ I tr . Another example is "reversed K-fold CV", which uses the same K-fold split but

Online Model Selection

In this subsection, we introduce the online model selection problem, setting the stage for a variant of cross-validation tailored for the online setting.

So far, we have only considered the batch setting, where the fitting algorithm f takes in a static dataset D n of an arbitrary sample size n and outputs f (D n ) ∈ F. The online setting, which becomes increasingly popular in large-scale machine learning, differs from the batch setting with profound algorithmic and computational consequences. Definition 3.1 (Online Estimator). Let X be the sample space of a single data point, and F a parameter space. An online algorithm f = ( fn : n g 1) is an infinite sequence of mappings such that (a) f1 : X → F , (b) For i g 2, fi : F × X → F , (c) After obtaining fi , the sample points {X j : j f i} and previous estimates { fj : j < i} are not used again by the algorithm in the future.

An online estimator f = ( fi : i g 1) takes in a (possibly infinite) stream of sample points (X i : i g 1) in X , with an initial estimate f0 , and updates the estimate from fi to fi+1 = fi+1 ( fi , X i+1 ) after seeing the (i + 1)th sample point (i g 0). Part (c) ensures that the online algorithm satisfies the "single-pass" (or constant memory) property, meaning that the sample points and estimates are only processed at the time of sampling and will not be stored. This does not exclude the popular Polyak averaging estimate, which outputs the final estimate as the average of ( fi : i g 1), because we only need to store the current average estimate, which costs a constant amount of memory.

Online algorithms often come with great computational efficiency and play a crucial role in modern machine learning. Perhaps the most important example of such online estimators is stochastic gradient descent.

Example 3.2 (Online Stochastic Gradient Descent). Suppose the estimation target is f * = arg min E X φ(f, X), for an objective function φ(•, •). The stochastic gradient descent estimate fn estimates f * by minimizing the empirical average objective function n i=1 φ(f, X i ) using an online update scheme.

(i) Initialize f0 = 0

The sequence (³ i : i g 1) specifies the learning rate at each step.

To provide a more concrete example, and better connect to recent literature, we consider the following variant of "sieve-SGD" tailored for nonparametric regression under squared loss, where each sample point X = (Y, Z), with Y ∈ R 1 and F being a class of functions mapping from Z to R 1 . Let (ϕ j : j g 1) be a set of orthonormal basis functions in L 2,Z . The sieve-SGD initializes f0 = 0 and the update from fi-1 to fi is:

(3.10)

Here ³ i is the learning rate, J i ∈ N is the number of basis functions used at sample size i, and w is a shrinkage exponent. It has been shown in Zhang and Simon (2022) that with suitable choices of (³ i , J i ) and w, the sieve-SGD estimate nearly achieves the minimax optimal error rate over Sobolev classes of functions. When model selection is of interest, a fundamental difference emerges between the batch setting and the online setting. In the batch setting, the evaluation criterion for model selection is often tied to a particular sample size, which is typically n but could also be n tr such as in Kfold cross-validation. The batch cross-validation can approximate the quantities µ ntr and Rcv,n,K because the algorithm can access versions of f (D ntr ) through f (D n,-k ), and repeatedly evaluate their losses at many independent sample points. This is not possible in the online setting, where, as specified in Definition 3.1, for each sample size i, the loss function of the estimate fi at sample size i can only be evaluated once using ℓ(X i+1 , fi ). As a result, the model selection criterion in the online setting must have the following two distinct features.

1. The model selection criterion in the online setting needs to be asymptotic, as it must describe the behavior of an infinite sequence.

2. As a consequence, in the online setting, a candidate in the model selection problem must be a sequence itself, which specifies the update rule for each sample size.

Given these thoughts, suppose we have m online algorithms, ( f (r) i : i g 1) for r = 1, ..., m, we can define the online model selection target as the r * ∈ [m] such that, assuming its existence,

where Ψ r,i is the excess risk of the fitted model f (r) i , as defined in (3.2) in the squared loss regression setting.

Such an optimal r * naturally exists in many nonparametric regression settings. Consider the sieve-SGD example in (3.10). For simplicity, suppose we only need to choose the number of basis (J i : i g 1).

Existing theory Zhang and Simon (2022) suggests the optimal choice is J i = Bi Ä for a pair of positive constants (B, Ä ) that depends on the function class. Thus, the candidate set of m online algorithms correspond to m distinct combinations {(B (r) , Ä (r) ) : r = 1, ..., m}, where for each r ∈ [m], the corresponding online algorithm uses the update rule (3.10) with

In this example, an online model selection procedure should point us to a choice of sequence {J (r) i : i g 1} that will eventually become the best for each i when the sample size keeps increasing. In contrast, standard batch model selection methods only focus on a single sample size and do not suggest a good choice of model for a different sample size. For example, if batch model selection returns a choice of J = 10 at sample size 100, it is unclear whether this J comes from a (B, Ä ) combination of (1, 1/2) or ( √ 10, 1/4), which would lead to different choices of J for other values of sample size.

Assumption 3.3 (Estimator stability).

There exist positive constants c, b 1 , and b 2 such that, for j f i, r = 1, 2,

Recall that ∇ j is the perturb-jth-entry operator defined in (2.5). If we ignore the conditional expectation, which makes the L 2 -norm smaller, (3.13) roughly requires that the value of fi at a future sample point change no more than i -br when one of the fitting sample points is replaced by an i.i.d. copy. Intuitively, if each fitting sample point contributes equally to the estimate, then the influence of a single sample point on the outcome should be of order i -1 and such a stability condition can be plausible if b r f 1. The almost-sure bound condition can be relaxed to L q norm bounds with different choices of conjugate indices in Hölder's inequality, together with stronger moment conditions on the risk function. We provide examples satisfying (3.13) in Section 6.4.

We consider the following online rolling validation risk estimate.

The RV risk estimate has some promising features.

1. The RV risk estimate Ξ r,n resembles the leave-one-out cross-validation, as each fitted model f (r) i -fitted using all sample points up to time i-is evaluated at one independent sample point X i+1 .

2. The RV risk estimate Ξ r,n is fully compatible with the online learning protocol described in Definition 3.1, as one can obtain

3. When n becomes large, Ξ r,n is increasingly representative of the asymptotic behavior of Ψ r,n and È r,n , thereby serving the purpose of online model selection objective in (3.11).

The intuition is that, if the fitting procedure does not too heavily rely on any single sample point, as implied by the stability condition (3.13), then the contribution of each single sample point to Ψ r,i will vanish as i increases. Thus each single sample point still only contributes a small proportion of the sum Ξ r,n , making the sum concentrate around its expected value. This is formally stated in the following theorem.

Theorem 3.2 (Model Selection Consistency of Rolling Validation). Under the regression model (3.1) and the squared loss function, assume that the function space is uniformly bounded:

then the rolling validation risks Ξ r,n defined in (3.14) satisfy

where, recall the notation D i = (X j : 1 f j f i). The first three terms are mean zero noises, and the last term is the signal, which will be negative for large values of n. The plan is to derive upper bounds for each of the noise terms so that their sum is unlikely to override the signal. In this proof, since the claimed result is consistency, we will not keep track of all constant factors. To simplify notation, we will use "≲" to denote that the inequality holds up to a constant factor.

Weighted Rolling Validation

In the online model selection criterion (3.11) and the model quality assumption Assumption 3.2, it is clear that the better model is defined through its superior asymptotic prediction risk. When reflected in the rolling validation risk Ξ r,n , the expected model quality difference is

Under Assumption 3.2, this quantity diverges as n increases, and is driven by the contributions from large values of i. In other words, these assumptions suggest that the model quality is better reflected by the performance of f (r) i at larger values of i. This intuition leads to the following weighted version of rolling validation (wRV).

Ξ r,n,À :=

where À g 0 is a positive constant controlling the fraction of weight tilted towards the larger values of i.

The original rolling validation risk Ξ r,n defined in (3.14) corresponds to the special case of Ξ r,n,0 , and it is straightforward to check that the wRV risks are also fully compatible with the online learning protocol in Definition 3.1. More generally, we can write wRV as

where (w i : i g 1) is a non-decreasing sequence of positive weights such that lim

The specific choice of w i = i À is quite natural, as a constant fraction of all sample points will have a constant fraction of the maximum weight. For example, at time i, a (1 -2 -1/À ) proportion of sample points will have weights at least half as much as the maximum weight i À . These sample points are also the most recent and thus more representative of the asymptotic model quality than the earlier sample points. In contrast, the exponential weighting w i = e Ài will only assign non-trivial weights to a constant number of sample points, leading to an overly large variance in the wRV statistic.

Consistency of wRV. Interestingly, weighted rolling validation enjoys the same consistency properties, regardless of the value of À.

Corollary 3.3. Let Ξ r,n,À be the wRV statistic defined in (3.15). Then, for any À g 0, under the same conditions as in Theorem 3.2, we have

Intuitively, the weights i À affect the variance and model quality gap in the same manner and will cancel out in the final signal-to-noise ratio calculation. The proof follows exactly the same procedure as that of Theorem 3.2.

The effect and practical choice of ξ. Although Theorem 3.3 shows that wRV can be consistent for any nonnegative value of À, it does not reflect the benefit of using a positive value of À. In practice, the choice of À may still impact the finite sample bias-variance trade-off, and wRV with a positive À generally outperforms the original RV with À = 0. In general, a small À leads to smaller variance but also larger bias (one way to interpret the bias is the representativeness of the wRV statistic for the asymptotic model quality). A theoretically rigorous understanding of the effect of À requires a finite sample analysis of the wRV model selection error probability, which is an open problem. Here we provide a heuristic argument to demonstrate the "quick response" effect of a positive À in a noiseless scenario, and give a simple rule-of-thumb choice of À.

Assume that the expected risks have the form

We hope to have Ξ 1,i < Ξ 2,i soon after i exceeds i * . In the noiseless case, Ξ r,i,À reduces to its expected value, which equals (when r = 1)

provided that Àa ̸ = -1. A similar equation can be derived for r = 2. Thus, approximately the noiseless wRV would be able to correctly favor the better model (r = 1) when i is large enough such that

The quantity 1/(À + 1a) in the above equation can be interpreted as the asymptotic delay ratio caused by the rolling sum. It is apparent that although the delay ratio is a monotone decreasing function of À, the improvement vanishes as À increases. In practice, a rule-of-thumb choice of À is À RoT ≡ 1 , which controls the asymptotic delay ratio to be at most 1, while keeping the variance at a reasonable level. Empirically, the model selection accuracy of wRV often observes the most significant improvement when À increases from 0 to 1, and becomes relatively stable for À ∈ [1, 3].

Bibliographic Notes

The pairwise model selection consistency of K-fold CV in mean regression under the squared loss is largely based on Yang (2007). Our Theorem 3.1 is a simplified version of Theorem 1 in Yang (2007), which covers more general tail conditions. As shown in Example 3.1, the condition n te È 2,ntr → ∞ captures the correct sufficient condition in this simple example. In terms of difference, our result is stated without the "convergence at exact rate" condition, which is required in Yang's Theorem 1. This simplification comes from a more streamlined treatment. Earlier investigations of model selection consistency by cross-validation include Zhang (1993) and Shao (1993).

The online model selection problem is motivated by recent developments in nonparametric stochastic approximation, such as kernel SGD (Kivinen et al., 2004;Dieuleveut and Bach, 2016) and sieve SGD (Zhang and Simon, 2022;Quan and Lin, 2024). The rolling validation and weighted rolling validation results are based on Zhang and Lei (2023), which contains a detailed treatment of the diverging candidate set setting, as well as numerical examples demonstrating the practical advantage of wRV over traditional CV model selection.

This section focuses on mean regression with squared loss. A natural question is whether similar results can be obtained for more general learning problems and loss functions. In this direction, Yang (2006) developed CV model selection consistency in binary classification that parallels Section 3.1 and Yang (2007). Wang and Shen (2006) also attempted to obtain a similar extension, requiring an accurate estimate of the conditional probability of Y given X used for a resampling step. In contrast, the proofs of the consistency results in this section mainly rely on the concentration of loss functions around their mean values. The general tools developed in Section 2.3, Section 2.4, and Section 2.5 could be used to provide useful results to establish similar model selection guarantees for general loss functions and online algorithms.

The model selection consistency of cross-validation has also been studied outside of the i.i.d. case. For example, Racine (2000) developed a variant of cross-validation for consistent time series model selection. Some early explorations of methods similar to the weighted rolling validation in the time series forecasting setting include Wei (1992) and Ing (2007). In the network analysis literature, Chen and Lei (2018), Li et al. (2020), and Qin and Lei (2022) used different sample splitting and cross-validation techniques for community-structured network model selection.

The focus of this section shifts back to risk estimation, where the target of CV risk estimate Rcv,n,K , defined in (2.2), is either the average risk of the leave-one-fold-out estimates defined as Rcv,n,K in (2.3) or its population version µ ntr defined in (2.1). Consistency and exponential concentration of Rcv,n,K around these targets have been established under various stability conditions in Sections 2.3 to 2.5. In this section, we are interested in pushing these results one step further to establish corresponding central limit theorems. These results, which are of independent theoretical interest, also pave the way for new methodological developments in a variety of inference problems, including constructing model selection confidence sets. As shown in Section 5.2, these can be further developed into improved model selection methods that overcome the theoretical limitation of standard CV methods observed in Section 3.

Starting in this section, we will use a new notation. For i ∈ [n], let k i ∈ [K] be the fold index such that i ∈ I n,k i . So we can write the cross-validation risk estimate as a sum over i:

CLT with Random Centering

We first consider the "random centering" case, where the quantity Rcv,n,K is treated as the target of Rcv,n,K and hence the centering of the central limit theorem. Proposition 2.2 and Theorem 2.6 suggest that perturb-one stability of the loss function with rate o(1/ √ n) implies sub-Gaussian type concentration of Rcv,n,K at rate n -1/2 + o(1). In this section, we further refine the analysis to show that under essentially the same stability condition the CV risk estimate centers at Rcv,n,K with √ n-asymptotic normality.

Our starting point is an intuition that if the risk function ℓ(X 0 , D n ) is perturb-one stable at a rate that vanishes faster than 1/ √ n, then the variance contributed to ℓ(X 0 , D n ) by the fitting part D n will be upper bounded by o(1) using the Efron-Stein inequality.

Theorem 4.1 (Efron-Stein Inequality). Let h be a function that maps D n to a real number, where D n = (X i : 1 f i f n) consists of independent entries (not necessarily identically distributed or even defined in the same space). Then

The proof of the Efron-Stein inequality is elementary and is based on the telescoping decomposition (2.15). A full proof can be found in Boucheron et al. (2013, Chap. 3.1).

As a consequence, one can expect the summand term ℓ(X i , D n,-k i ) in Rcv,n,K to behave like its de-randomized version

The subindex n tr in the notation ℓ ntr (•) emphasizes that this function only depends on the training sample size. More generally, we can define, for any n g 1,

(4.1)

The following result quantifies the contribution to the total variance of the loss function ℓ(X i , D n,-k i ) from D n,-k i , i.e., the fitting randomness.

Proposition 4.1. Assume the fitting procedure is symmetric in the sense of Definition 2.1 and D n consists of i.i.d. sample points. Let

Proof of Proposition 4.1. Using the variance decomposition

Thus it suffices to show

Using the Efron-Stein inequality on the conditional distribution of ℓ(X 0 , D n ) given X 0 we have

The desired result follows by taking expected value on both sides.

Proposition 4.1 suggests that if the loss function is stable enough then the randomness in ℓ(X 0 , D n ) comes mostly from X 0 and it is reasonable to expect

which is the average of n i.i.d. random variables. However, (4.2) may not be accurate enough in general cases for the purpose of Gaussian approximation. We will need a slight modification on the LHS so that it focuses on the randomness contributed from the loss function evaluation, as detailed in the following theorem.

Theorem 4.2 (CV-CLT, Random Centering). Let (X i : i g 0) be an i.i.d. sequence of sample points in X and

where

for all ϵ > 0 and ℓ n (X 0 ), ln (X 0 ) defined in (4.1). Then for any value of K that evenly divides n and is allowed to change with n, the K-fold cross-validated risk Rcv,n,K satisfies

where the symbol "⇝" means "convergence in distribution", and Rcv,n,K is defined in (2.3).

Proof of Theorem 4.2. The proof proceeds by rigorously establishing a modified version of the "approximate linearity" of Rcv,n,K in (4.2).

The rest of the proof consists of two main steps. In the first step we show that

In fact, by construction, the LHS of (4.5) multiplied by

For the squared term (À i -Ài ) 2 we have, using the same argument as in Proposition 4.1,

where the second inequality follows from the fact

as implied by Theorem 4.12 in Section 4.4. As a result, the total contribution from the squared terms in (4.6) is o(nÃ 2 ntr ). For the cross term (À i -Ài )(À j -Àj ), if k i = k j , then by conditioning on D n,-k i and using conditional independence of X i and X j given D n,-k i we have E(À i -Ài )(À j -Àj ) = 0. Now we only need to control the cross-term contributed by (i, j) pairs such that j / ∈ I k i . The main trick here is the "free nabla" lemma (Theorem 4.11 in Section 4.4). By applying Theorem 4.11 twice we have, for j /

whose absolute value is o(Ã 2 ntr /n tr ). Thus the total contribution from the cross-terms in (4.6) is o(nÃ 2 ntr ). This establishes (4.5). It remains to show that

Ài ⇝ N (0, 1) , which, under (4.4), follows from the standard Lindeberg central limit theorem for triangular arrays (Durrett, 2019, Theorem 3.4.5).

Remark 4.1. The proof of Theorem 4.2 implies that

where d W,1 is the Wasserstein-1 distance. Moreover, if we assume a specific rate of convergence in (4.3), as well as a stronger integrability condition than (4.4), such as uniform boundedness of ∥[ℓ n (X 0 )µ n ]/Ã n ∥ 3 , then we can establish the Berry-Esseen-type finite sample error bounds of normal approximation in the Wasserstein-1 and Kolmogorov-Smirnov distances.

Remark 4.2. In most applications, it is reasonable to assume that there exists a pair of positive constants (c, c) such that

for all n. For example, in regression problems, the lower bound corresponds to the aleatoric uncertainty that comes from the observation noise, and cannot be reduced by the prediction function. The upper bound usually corresponds to the parameter space and the total variability of the random variables.

As another consequence of Proposition 4.1, one can expect to obtain a consistent estimate of Ã ntr using the sample variance of {ℓ(X i , D n,-k i ) : i ∈ [n]} as the variance is mainly contributed by the evaluation point X i .

CLT with Deterministic Centering

In this subsection, we develop a central limit theorem for the CV risk estimate with the deterministic centering: µ ntr = E Rcv,n,K .

The random centering CLT proved in the previous subsection provides a head start. A key intermediate step in that proof is, under the loss stability condition (4.3) and recalling that ln

which can be rewritten equivalently as

where

accounts for the difference between the random centering Rcv,n,K and the deterministic centering µ ntr .

The additional challenge in order to establish a deterministic centering CLT for Rcv,n,K is to characterize the randomness in G, a mean zero random variable that involves each sample point in a symmetric fashion. Unlike the loss function, where each sample point X i is evaluated exactly once in ℓ(X i , D n,-k i ), the randomness in G features the sum of the risks of parameters fitted from different leave-one-fold-out subsamples, and each X i contributes, in a symmetric fashion as other fitting sample points, to all but one such risks. Thus we will proceed by considering the projection of G onto each X i , which is involved in the parameter fitting part in R(D n,-k j ) for all j / ∈ I k i . Define

and ḡn (x) = g n (x)nµ n . (4.10)

The quantity g n (X i ) is the counterpart of ℓ n (X i ) when considering the contribution of X i to the fitting variability.

By construction, it is straightforward to check that

It is reasonable to expect the linear decomposition

If we can establish (4.11) with enough level of accuracy, then one can establish a CLT with deterministic centering from (4.8) using standard i.i.d. sums.

To further explore the plausibility of (4.11), we consider the following martingale decomposition of G:

with the convention D 0 = ∅.

The martingale increment E(∇ i G|D i ) represents the contribution to G from the ith sample point. We shall also see that, at least intuitively, when the loss function is stable, the choice of ordering of i in (4.12) does not matter. Now it seems natural to expect this increment to be close to ḡntr (X i ). Indeed we have, using another round of martingale decomposition with respect to the Ã-fields generated by (D j , X i ),

(4.13)

Now we arrive at the following approximate linear form for G:

It remains to control the residual term n i=1 Λ i . Two observations are useful.

For any

Proof of Theorem 4.4. According to the proof of Theorem 4.2, the first stability condition establishes (4.8). Now the second stability condition ensures that n -1/2 i Λ i = o P (Ä ntr ) according to (4.15). The result then follows from the standard Lindeberg CLT for the triangular array { lntr (X i ) + ḡntr (X i ) :

The asymptotic variance Ä 2 ntr is the total variance of the sum of two parts: ℓ ntr (X i ) and g ntr (X i ), which corresponds to the evaluation and estimation randomness contributed by the sample point X i , respectively. In contrast, the random centering CLT Theorem 4.2 only captures the evaluation part involving ℓ ntr (X i ), leaving the randomness of estimation to the random centering Rcv,n,K . While it is relatively straightforward to assume the boundedness of Var(ℓ ntr (X 1 )), the behavior of Var(ℓ ntr (X 1 )+ g ntr (X 1 )) is harder to track in general. The correlation between ℓ ntr (X 1 ) and g ntr (X 1 ) was analyzed in several special cases in Austern and Zhou (2024, Sec. 4), and it may be positive, negative, or zero. In particular, when f is obtained by an M-estimator whose objective function coincides with the loss function ℓ(x, f ) and satisfies regularity conditions such as convexity and smoothness, the correlation is shown to be zero. This makes a positive constant lower bound of Var(ℓ ntr (X 1 ) + g ntr (X 1 )) more or less reasonable. In order to assume an upper bound on Var(g ntr (X 1 )), it helps to consider the variability of n tr ∇ 1 R(D ntr ), because there are n tr terms in G that involves X 1 . An upper bound on Var[g ntr (X 1 )] can be obtained through the first-order risk stability ∇ 1 R(D ntr ):

(4.17)

In order for the upper bound in (4.17) to be O(1), it would require a risk stability rate of O P (1/n tr ), which seems substantially more stringent than the corresponding rate assumed for the loss stability in (4.3) as well as Part 1 of the stability assumption in Theorem 4.4. We argue that in many machine learning contexts, the risk can have much better stability than the loss function itself. As a first example, consider

More generally, suppose that F is a Hilbert space, and we have the first-order approximation

Usually ∇ i f corresponds to the stability of the loss function ℓ(X 0 ; f (D n )) if ℓ(X 0 , f ) is smooth in f , which could be as small as O(n -1 ). For risk minimization estimators, it is typically the case that dR(f ) df f = f (Dn) ≈ 0. Thus for such approximate risk minimization estimators it is reasonable to expect ∇ i R(D n ) to be much smaller than ∇ i ℓ(X 0 , D n ), making the requirement of |∇ i R(D n )| j n -1 plausible. This intuition also helps explain the condition

) required in Theorem 4.4. We will resume this discussion after presenting the variance estimation result below. Further discussion and examples of stable algorithms are presented in Section 6, especially in Section 6.3 for the risk stability.

Variance estimation with deterministic centering. It is unclear how to construct a consistent estimate of the asymptotic variance of Rcv,n,K without auxiliary sample points that are not in D n . Now suppose we do have access to some auxiliary sample points, independently drawn from the same distribution. The question we seek to answer here is how many such auxiliary sample points would be sufficient to guarantee a consistent estimate of the asymptotic variance of the CV risk estimate with deterministic centering.

If we have É(n) auxiliary sample points, then we can generate É(1) i.i.d. copies of Rcv,n,K and then obtain a consistent variance estimate from these i.i.d. copies. Now we show that one can also obtain a consistent variance estimate using only É(1) auxiliary sample points.

We begin with the proof of Theorem 4.4, which is based on the following more refined version of (4.8) and (4.11) with an explicit expression of the residual,

where

To study the variance of Rcv,n,K , it is natural to consider ∇ i Rcv,n,K . In particular we will consider i = n, and, assuming that

18) The convenience of choosing i = n in ∇ i is reflected in the last term of the equation above, where all (Λ i : i < n) disappear under the ∇ n operator. It is important to note that (4.18) also holds for arbitrary index i ∈ [n], with the last term becoming ∇ i Λ i,n , where Λ i,n is the counterpart of Λ n when one swaps index i and index n in the martingale decompositions in (4.12) and (4.13).

The first term in the RHS of (4.18) gives us the desired i.i.d. random variable with the desired asymptotic variance. Thus, if we can show that the other terms in the RHS of (4.18) are all sufficiently small with vanishing L 2 norms, then one can use the empirical variance of n∇ i Rcv,n,K with É(1) many distinct i's, which would only require É(1) auxiliary sample points.

Among the residual terms in the RHS of (4.18), the most nontrivial one to control is the second term ∇ n i / ∈I K (À i -Ài ). In the proof of Theorem 4.2 we have shown

thus we only need the ∇ n operator to cancel out the additional √ n factor in this upper bound. While this is intuitively plausible, we will need some additional assumptions for a rigorous result. We consider the following.

The two parts of Assumption 4.1 represent different possible strengthening of the stability condition (4.3). Part (a) is easy to hold when the number of folds K is a constant or grows with n at a slow rate. When K is a constant, Part (a) of Assumption 4.1 is equivalent to the original stability condition (4.3). On the other hand, part (b) of Assumption 4.1 requires a "second-order stability". Intuitively, if applying the perturbone difference operator ∇ i reduces the magnitude of ℓ(X 0 , D n ) by a factor of o(1/ √ n), then applying two such operators will reduce the magnitude by a factor of o(1/n). The advantage of Part (b) is that it can potentially cover the case of large values of K, such as K = n. But such an advantage comes at the cost of having to verify the second-order stability condition. We will discuss examples of second-order stability in Section 6.2.

Formally, assume we have auxiliary sample points X ′ 1 , ..., X ′ naux independent of D n with sample size n aux = É(1). Consider the following pseudo-estimate of the asymptotic variance of the CV risk with deterministic centering. If ∥ lntr (X 1 ) + ḡntr (X 1 )∥ 4 = O(Ä ntr ), then Theorem 4.5 implies consistent estimation with n aux = É(1). In practice, we can hold out n aux = √ n for variance estimation and conduct cross-validation on the remaining sample points.

Proof of Theorem 4.5. Using (4.18), we have

where ∆ i equals the sum of the last three terms in the RHS of (4.18).

The second and third terms on the RHS of the above equation has L 1 norm bounded by o(1) because ∥W i /Ä ntr ∥ 2 = √ 2 and ∥∆ i /Ä ntr ∥ 2 = o(1) uniformly over i according to Theorem 4.14. Thus it suffices to show the first term is o P (1). In fact, the first term has mean 0 and variance

Revisiting the risk stability conditions. The variance estimation argument above also provides further insights to the asymptotic variance Ä 2 ntr and the risk stability condition in Theorem 4.4. To begin with, we repeat the derivation of (4.12), (4.13), (4.14), (4.15), and (4.18) for the risk R(D n ). In particular, (4.14) becomes

Λi

where Λi is the counterpart of Λ i in (4.14). Also (4.18) becomes (suppose we put i = 1 in the final position in the first martingale decomposition (4.12))

where

∥ 2 by a similar argument used in (4.15). Therefore, we have the following result on the relationship between the variances of ∇ 1 R(D ntr ) and g ntr (X 1 ), which provides the reverse direction of (4.17).

Proof of Theorem 4.6. Equation (4.17) shows that the ratio is upper bounded by 1/2. On the other hand, recall that

, which provides the other direction.

Putting things together we can have the following more refined view of Theorem 4.4 as well as its relationship with Theorem 4.2. Proposition 4.2 (CV-CLT, Consolidated Version). Assuming an i.i.d. sequence of X i : i g 0 and a symmetric loss function ℓ(•, •).

that (4.22), (4.16) hold. Proof of Proposition 4.2. For the first part, it suffices to show that ∥ Rcv,n,Kµ ntr ∥ 2 = o(Ã ntr / √ n), which is in turn implied by ∥R(D ntr )µ ntr ∥ 2 = o(Ã ntr / √ n) using symmetry among the folds and triangle inequality. The condition n∥∇ 1

For the second part, it suffices to show that the conditions in Theorem 4.4 are satisfied. In particular, we only need to verify the second-order risk stability in Theorem 4.4. First, Theorem 4.6 ensures that ∥g ntr (X

Thus the second-order risk stability condition (part 2) of Theorem 4.4 is satisfied because by (4.22) we have

Proposition 4.2 provides a more transparent view of the role played by the risk stability quantities

CLT with deterministic centering. In the first part of Proposition 4.2, the condition n∥∇ 1 R(D n )∥ 2 = o(Ã n ) indicates that the fitting randomness in nR(D n ) contributed from a single sample point X 1 has much smaller variability than the evaluation ℓ n (X 1 ). In this case the second-order risk stability ∥∇ 1 ∇ 2 R(D n )∥ 2 is irrelevant and the deterministic centering essentially behaves like the random centering.

In the second part of Proposition 4.2, the condition n∥∇ 1 R(D n )∥ 2 = Ω(Ã n ) means that the variability in the fitting is at least as large as that in the evaluation. The condition Ä n = Ω(Ã n ) ensures that these two sources of variability do not cancel each other when we add up ℓ n (X 1 ) and g n (X 1 ). As mentioned in the discussion after Theorem 4.4, it is often plausible to assume a nonnegative correlation between ℓ n (X 1 ) and g n (X 1 ), and hence no cancellation of variance. In this case, the key condition on the second-order risk stability is that

So this is analogous to the original loss stability condition (4.3) in the sense that we obtain an o(n -1/2 ) factor from applying one more ∇-operator.

because the fitting variability dominates and we can have a higher tolerance of the instability.

High-dimensional Gaussian Comparison

The central limit theorems developed in Section 4.1 and Section 4.2 are concerned with a single estimator f . In the context of model selection and model comparison, it is necessary to study the simultaneous fluctuation of the CV risk estimates of multiple estimators. While the extension of the central limit theorem to multivariate cases with a finite and fixed number of estimators is straightforward, in this subsection we consider the situation with a potentially large number of estimators.

To begin with, we first introduce a high-dimensional Gaussian comparison result for sums of independent random vectors, initiated by Chernozhukov et al. (2013). A typical form of high-dimensional Gaussian comparison is to compare the distribution of the maximum of a random vector to that of a Gaussian vector with matching mean and variance. For a vector z = (z r : 1 f r f m) ∈ R m , let max z = max r∈[m] z r be the value of the maximum entry of z. The following result is a simplified version of Theorem 2.5 of Chernozhukov et al. (2022) adapted for our setting.

Theorem 4.7 (Gaussian Approximation, Simplified i.i.d. Case). Let (Z i :

for a constant C » depending on » only.

Gaussian comparison results like this provide probability-theoretic foundations for high-dimensional inference tasks, such as construction of simultaneous confidence intervals for high dimensional mean vectors.

When studying the joint randomness of the CV risk estimates of m estimators, it is more convenient to use vector notation. Let ℓ be the vector version of (ℓ (r) : 1 f r f m). Define ℓ n , ln , Rcv,n,K , and Rcv,n,K accordingly. Let S n be the m × m diagonal matrix with the rth diagonal entry being Ã (r)

n (X 0 )), and Γ n the correlation matrix of ℓ n (X 0 ):

Our goal is to characterize the simultaneous fluctuation of Rcv,n,K as an m-dimensional random vector. Theorem 4.2 implies that the mdimensional vector √ nS -1 ntr ( Rcv,n,K -Rcv,n,K ) is entry-wise asymptotically standard normal, provided that, among other regularity conditions, the loss stability condition (4.3) holds for each estimator f (r) . Therefore, it is natural to consider the corresponding multivariate Gaussian vector Y n ∼ N (0, Γ n ) as the target in the high-dimensional Gaussian comparison.

More specifically: Does (4.23) still hold when W n is replaced by

The following result provides a positive answer to this question under a strengthened version of the stability condition (4.3).

Theorem 4.8. Let (X i : i g 0) be a sequence of i.i.d. sample points, and ( f (r) : r ∈ [m]) a collection of symmetric estimators. Assume 1. there exists a positive constant ³ and a sequence ϵ n such that, for all r ∈ [m] and all n, √ nK

2. and there exists a finite positive constant » such that for all r and n l(r)

In Theorem 4.8, we allow the number of models m = m n to vary with n. Equations (4.25) and (4.26) are counterparts of (4.3) and (4.4) respectively. The sub-Weibull tail enables a simultaneous upper bound on many residual terms. The factor of √ K on the LHS of (4.25) is necessary here because the "free nabla trick" used in the proof of Theorem 4.2 no longer works for general L q norms and we can only use Rio's inequality within each test sample fold and resort to triangle inequality to combine the folds. Condition (4.25) is less stringent for smaller values of K, and becomes comparable to (4.3) if K = O(1). In practice, K-fold CV with a constant K, such as K = 5 and K = 10, is perhaps the most common version.

Remark 4.3. Theorem 4.8 is stated for the maximum entry of the marginally standardized error vector, which corresponds to one-sided simultaneous confidence bound. The same result can be easily extended to a two-sided version in that (4.28) holds true, under the same conditions, for the corresponding vectors of absolute values | Rcv,n,K -Rcv,n,K | and |Y n |. This can be achieved by considering m more loss functions ℓ (r+m) = -ℓ (r) for r ∈ [m] and applying Theorem 4.8 to these 2m loss functions.

Proof of Theorem 4.8. The proof is based on a more refined version of the asymptotic linear representation (4.5) using the sub-Weibull stability condition. Let

By Theorem 4.

15, each • n,r , the rth coordinate of • n , is (2ϵ ntr , ³+1)-SW. Consider i.i.d. random vectors Z i = S -1 ntr lntr (X i ). Condition (4.26) and Theorem 2.3 imply that max |Z 1 1) , (4.29) where the last o(1) statement comes from the first part of (4.27).

Under the second part of (4.27), one can find a ¶ n such that

where the first o(1) in the second inequality comes from (4.29), the second o(1) in the same line comes from the fact that max |• n | is (c ³ ϵ ntr (log m) ³+1 , ³ + 1)-SW as guaranteed by Theorem 2.3 and the choice of ¶ n , and the last inequality comes from a Gaussian anticoncentration result (Chernozhukov et al., 2022, Lemma J.3). Similarly we can obtain the other direction and conclude the proof.

High-dimensional Gaussian comparison with deterministic centering.

It is possible to extend Theorem 4.8 to deterministic centering. While it is desirable to have a result analogous to Proposition 4.2, the dichotomy of whether the validation uncertainty dominates the fitting uncertainty needs to be checked for each coordinate. Define µ n , g n (X 0 ), and ḡn (X 0 ) as the vector versions of µ n , g n (X 0 ), and ḡn (X 0 ), respectively. For simplicity, we consider the following two scenarios.

In the first scenario, where

n ϵ n , ³)-SW for all r ∈ [m] with ϵ n vanishing fast enough in the sense of (4.27), then Theorem 4.8 would hold for the deterministic centering without much modification. The following theorem is a deterministic centering counterpart of Part 1 of Proposition 4.2. Theorem 4.9. Assume the same conditions in Theorem 4.8.

n εn , 3)-SW for all r ∈ [m] with εn (log m) 3+1 = o(1), then

Proof of Theorem 4.9. For the first part, let •n = √ n max |S -1 ntr ( Rcv,n,K -µ ntr )|. The sub-Weibull condition assumed on ∇ 1 R (r) (D n ) implies that, by Theorem 2.5 and triangle inequality, each coordinate of

for some constant c 3 depending only on 3.

The condition εn (log m) 3+1 = o(1) ensures existence of a sequence ¶n such that ¶n = o(1/ √ log m) and (log m) 3+1/2 ϵ ntr = o( ¶n ). Therefore,

where the first inequality follows from union bound, the second inequality from Theorem 4.8 and the sub-Weibull tail of •n , the third from Gaussian anti-concentration (Lemma J.3 of Chernozhukov et al. (2022) since ¶n = o(1/ √ log m). The other direction can be established similarly.

In the second scenario, when ∇ 1 R (r) is large as in the second case in Proposition 4.2 for some values of r, the deterministic centering would lead to a different asymptotic correlation matrix, whose (r, s)-entry equals Corr( l(r)

The following theorem is a deterministic centering counterpart of Theorem 4.4.

n is (», ³)-SW with constants (», ³), and ϵ n such that (4.27) holds, then

where T n is an m × m diagonal matrix with the rth diagonal entry being

n , and Y n a centered Gaussian vector whose covariance matrix has (r, s)-entry Corr( l(r) ntr (X 0 ) + ḡ(r) ntr (X 0 ), l(s) ntr (X 0 ) + ḡ(s) ntr (X 0 )).

Proof of Theorem 4.10. Using the vector versions of (4.8), (4.9), and (4.14), we have

where •n is the analogous version of • n in the proof of Theorem 4.8 with S n replaced by T n . The first-order stability condition implies that •n is component-wise (c ³ ϵ n , ³ + 1/2)-SW, and the martingale decomposition of

The rest of the proof follows a similar argument as in Theorem 4.9.

A fully consolidated version of deterministic centering Gaussian comparison would involve studying the scenario where the risk stability condition ∇ 1 R (r) (D n ) j n -1 for some but not all coordinates r ∈ [m],

which is an open problem. Practically, the random centering Gaussian comparison is easier to use with easier-to-verify assumptions. In some applications, random centering Gaussian comparison provides sufficient theoretical support for non-trivial inference problems. We provide a couple of such examples in Section 5.

Variance and correlation estimation. In order to use the high-dimensional Gaussian comparison Theorem 4.8 to construct confidence sets for Rcv,n,K , one needs to estimate not only the marginal variances (Ã (r) ntr ) 2 for r ∈ [m], but also the pairwise correlation Γ ntr (r, s) in (4.24). The same argument in the proof of Theorem 4.3 can provide entry-wise L 2 -consistency of the sample variance and correlation estimates. However, such entry-wise consistency is not sufficient for high-dimensional inference, because when m = m n increases with n, the Gaussian comparison result is not a central limit theorem as there is no fixed limit of convergence in distribution. To make valid inference with the estimated variances and correlations, we need a more refined error analysis of the estimated variance and covariance, which will be detailed in Section 5.1 below.

Auxiliary Lemmas

In this subsection we collect the auxiliary lemmas used in earlier subsections.

The following lemma explains the "free nabla" trick.

Lemma 4.11 (Free Nabla Lemma). Let D = (X 1 , ..., X n , X ′ 1 , ..., X ′ n , W ) be independent random objects such that for each i ∈ [n], X i and X ′ i have the same distribution. Let h and g be two functions that act on subsets of D such that for some i ∈ [n] 1. g(D), h(D), and g(D)h(D) all have finite L 1 norms, 2. E X i g = 0, and 3. h does not involve X ′ i .

Then

The operator ∇ i in this lemma is defined in the usual sense since h does not involve X ′ i . In particular, ∇ i h is the difference between the original h and the one obtained by replacing X i by X ′ i in the calculation of h.

Proof of Theorem 4.11. Let h i be the version of h obtained by swapping X i and X ′ i in D. It suffices to show that E(gh i ) = 0. In fact, by construction h i does not involve X i since h does not involve X ′ i . Thus

The following lemma provides a convenient bound on ∇ i À j for i / ∈ I k j .

Lemma 4.12.

Proof of Theorem 4.12. Let D i n,-k j be the dataset obtained by replacing X i by X ′ i in D n,-k j while keeping everything else the same. Then

where the last step follows from

by first taking expectation over X j .

The following lemma is a special case of Rio's inequality Theorem 2.4 and generalizes Theorem 2.5. Lemma 4.13 (Rio's Inequality For Conditionally Independent Sequences). Let Z 1 , ..., Z n be random variables that are conditionally independent given a Ã-field C such that E(Z i |C) = 0 for all i. Then for any q g 2,

Proof of Theorem 4.13. The proof is identical to that of Theorem 2.5.

Theorem 4.13 can be used to prove various important consequences of the loss stability, collected in the following lemma.

The following holds for any fixed q g 2 and each i ∈ [n],

Proof of Theorem 4.14. For Part 1, the first inequality follows from the conditional version of Rio's inequality as applied to the martingale increment sequence

for 1 f i f n, as in the proof of Theorem 2.5. The second inequality follows from the first one and Jensen's inequality. The third one follows from triangle inequality.

In the proofs of Parts 2 and 3, without loss of generality we assume i = 1 and k 1 = 1.

For Part 2, observe that

where the second line follows from an application Theorem 4.11. For Part 3, consider k ̸ = 1 and observe that (∇ 1 À i : i ∈ I n,k ) is a martingale increment sequence with respect to the collection of nested Ã-fields:

where the first inequality follows from Theorem 4.13 because (∇ 1 À i : i ∈ I n,k ) are conditionally i.i.d. given (D n,-k , X ′ 1 ) and have conditional expectation 0. Then the desired result follows from the triangle inequality by summing over k ∈ {2, ..., K}.

For Part 4, fix a k ∈ [K], then using Theorem 4.13 and Part (1) we obtain

Thus the claimed result follows from triangle inequality.

The following lemma provides ℓ q norm bound for the random centering residuals.

Lemma 4.15. Let (X i : i g 0) be an i.i.d. sequence, f a symmetric estimator, and ℓ(•, •) a loss function.

Proof of Theorem 4.15. Use notation

and the result follows from Part 4 of Theorem 4.14.

Bibliographic Notes

The rigorous normal approximation of cross-validated statistics has been conjectured and regarded as an open problem in Yang (2006).

The univariate and low-dimensional central limit theorems for crossvalidated risk estimates are first established in Bayle et al. (2020) and Austern and Zhou (2024). The random centering CLT based on the first-order loss stability (4.3) presented in Theorem 4.2 largely follows the argument in Bayle et al. (2020). The same random centering CLT is also proved in Austern and Zhou (2024), which relied on an additional second-order loss stability as in Part (b) of Assumption 4.1. This secondorder stability condition turned out to be redundant for the random centering CLT. While Bayle et al. (2020) does not offer a deterministic centering CLT, Austern and Zhou (2024) provides one. Our deterministic centering CLT (Theorem 4.4) requires fewer assumptions and offers additional insights through the asymptotic linear representation of the G term as well as the consolidated version Proposition 4.2. In particular, unlike Austern and Zhou (2024), our deterministic centering CLT does not require the second-order loss stability or a specific bound on the first-order risk stability. Our consolidated version Proposition 4.2 also makes it clear that, when the risk variability is non-negligible, it is the relative size of the second-order risk stability to the first-order risk stability that matters. In addition to simplified assumptions, our proof of the deterministic centering CLT is also substantially simpler than that of Austern and Zhou (2024). Our proof clearly reflects the role played by stability, which is solely in the asymptotic linearity of the CV risk estimate. The deterministic asymptotic variance estimate presented in (4.19) is the same as the one proposed in Austern and Zhou (2024).

A deterministic-centering CLT for K-fold cross-validation is established by Li (2023) for regression under high-level assumptions, including (i) the third derivative of the loss with respect to the regression function's value is identically zero, and (ii) the estimated regression function admits an asymptotic linear representation. These assumptions are highly nontrivial: The third-derivative condition forces ℓ((z 0 , y 0 ), f ) to depend on f quadratically in f (z 0 ), and the asymptotic linear representation is an abstract property of the loss/estimator pair which requires verification. The developments in Section 4.2 show that the stability conditions in Theorem 4.4 and proposition 4.2 are sufficient to guarantee such asymptotic linear representations.

The high-dimensional Gaussian comparison for independent random vectors has been studied as early as in Bentkus (2003), with a major breakthrough in Chernozhukov et al. (2013) that allows the dimensionality to be exp(n a ) for small values of a. Many further developments have been made since Chernozhukov et al. (2013), with the most recent ones being Chernozhukov et al. (2022) and Chernozhukov et al. (2023). Readers are referred to those papers for a more comprehensive literature review on this topic. In our case we used results from Chernozhukov et al. (2022) as it does not require the asymptotic covariance matrix Γ n to be non-degenerate. It is possible to use more refined techniques, such as those in Chernozhukov et al. (2023), to obtain approximation error bounds with better dependence on n.

The high-dimensional Gaussian approximation for cross-validated risks presented in Section 4.3 is developed in Kissel and Lei (2022), which is motivated from the model selection literature (Lei, 2020). The proofs presented in this article are much simpler and more streamlined than those in Kissel and Lei (2022), and feature a similar reduction in required assumptions as in the univariate CLT. Our Gaussian comparison result for cross-validation can also be viewed as an instance of high-dimensional Gaussian comparison for dependent data. Existing results in the literature include vector-valued U -statistics (Chen, 2018) and vector-valued stochastic processes with mixing properties (Chang et al., 2024;Kurisu et al., 2024). These results are based on a short-range dependence structure, where the dependence between two random vectors decays quickly as the indices move apart from each other. However, such a structure does not hold for cross-validation, where the loss functions evaluated at any two sample points have a small but nonzero correlation due to the involvement of sample points in both evaluation and estimation. tuning parameter selection, the stochastic dominance condition becomes less plausible, and, in order to take into account the uncertainty in the CV risk estimates, it becomes more reasonable to aim at finding all "nearly optimal" models. The tools developed in Section 4.3 are suitable for this purpose.

From risk confidence set to model confidence set. Suppose we run K-fold CV for each of the m candidate estimators. Our goal is to learn

and in the case that there exist multiple good models with close to optimal values, find a confidence set Θ ¦ [m] such that

for some pre-specified Type I error level

´. It may not seem entirely natural to choose r * as the target parameter of interest, as it is the average risk of K fitted parameters obtained from K leave-one-fold-out estimates. This choice is aligned with Theorem 4.8, which provides simultaneous uncertainty quantification of Rcv,n,K -Rcv,n,K . Moreover, Theorem 4.9 suggests that if the risks are stable enough, then such a random centering result is equivalent to the deterministic centering with sample size n tr .

Let ´∈ (0, 1) be a nominal Type I error level and Σ a positive semi-definite matrix. Define

where | • | means entry-wise absolute value. Then according to Theorem 4.8, under the conditions therein, the following is an asymptotic (1 -´) simultaneous confidence setfoot_4 of Rcv,n,K :

where Γ ntr is the correlation matrix defined in (4.24). However, the confidence set in (5.4) is not practical as it involves unknown quantities: (Ã (r) ntr , Γ ntr , t ´,Γn tr ). We can estimate Ã ntr for each r as in Theorem 4.3. Similarly, we can estimate Γ ntr (r, s) by

(5.5)

After obtaining Γntr , we can numerically approximate t ´,Γ n tr using Monte Carlo methods. The quality of approximation depends on the Monte Carlo sample size and can be made arbitrarily small. In the following analysis we will ignore this approximation error and simply work under the assumption that t ´,Γ n tr is available. However, our analysis does explicitly account for the estimation error in Γntr and the difference between t ´,Γ n tr and t ´,Γn tr .

Theorem 5.1. If (4.25) and (4.26) hold and

then for any ´∈ (0, 1)

Recall that (4.25) and (4.26) are the stability and tail conditions of the loss ℓ (r) (X 0 , D n ), respectively. These conditions are stated uniformly over r ∈ [m]. Intuitively, a K-fold version of Theorem 2.6 suggests that each coordinate of Rcv,n,K concentrates around its expected value at scale O(n -1/2 +ϵ ntr ). The additional (log m) 3³+4 factor required in (5.6) is used to ensure (i) the concentration of Rcv,n,K holds for all entries at the same time, and (ii) the correlation matrix estimation is uniformly accurate over all entries.

Proof of Theorem 5.1. We use simplified notation by dropping the subindices whenever there is no confusion: R = Rcv,n,K , R = Rcv,n,K , Γ = Γntr , Γ = Γ ntr . Let Ŝ be the m × m diagonal matrix whose rth diagonal entry is Ã(r) ntr , and S its population version whose corresponding entry is Ã (r) ntr . Define

We first control |W 0 -W 1 |, observe that

By Theorem 4.15 and Proposition 2.4, each entry of

)-SW and hence, by Theorem 2.3, the maximum entry-wise absolute value of

). Combine this with Part 1 of Theorem 5.2 we have

Next we relate t ´,Γ to t ´,Γ . Under the assumption that (ϵ ntr + n -1/2 )(log m) 3³+4 = o(1), there exists a sequence ¶ n such that

. By the Gaussian-to-Gaussian comparison result (Proposition 2.1 of Chernozhukov et al. (2022)) we have, on E,

where c is a universal constant, and hence

Finally we arrive at the following inequalities

where the first inequality follows from union bound, the second from another union bound over (E, E c ) and the fact that t ´,Γ f t ˜,Γ on E, the third from Theorem 4.8, the fourth from the Gaussian anti-concentration (Lemma J.3 of Chernozhukov et al. (2022)), and the last line follows from the choice of ¶ n . A similar argument provides the other direction and concludes the proof.

Theorem 5.1 leads to a confidence set Θ satisfying (5.2):

(5.7) In other words, this confidence set includes all models whose confidence lower bound is not dominated by all the confidence upper bounds, because otherwise we are 1 -´confident that there exists an s such that R(s) cv,n,K < R(r) cv,n,K , as witnessed by the s that violates the inequality in (5.7).

Proposition 5.1. Under the same conditions in Theorem 5.1, we have

Remark 5.1. The model confidence set given in (5.7) and its theoretical guarantee, Proposition 5.1, do not require the minimizer of Rcv,n,K to be unique. The interpretation of the coverage guarantee in Proposition 5.1 is that it holds marginally for each individual r * that minimizes R(r) cv,n,K .

Model confidence set based on differences. The model confidence set Θ0 given in (5.7) is often too conservative as it compares the least optimistic prediction of R(s) cv,n,K with the most optimistic prediction of R(r) cv,n,K . Here the term "most optimistic" corresponds to the lower bound of the simultaneous confidence set and vice versa. In practice, the correlation between the loss functions calculated from different estimators tends to be positive, as they are targeting the same task of estimating the parameter. Therefore, it is highly unlikely that a large positive fluctuation of the CV risk estimate in one estimator is associated with a large negative fluctuation in another.

To take into account the likely positive correlation between the cross-validated risks, one approach is to consider the difference of CV risks and loss functions. For a given r, we want to test whether r ∈ arg min s R(s) cv,n,K , which is equivalent to test

To this end, consider the difference loss function

(5.8)

The difference based model confidence is obtained by constructing a confidence lower bound for

by applying Theorem 4.8 to the (m -1)-dimensional vectors (ℓ (r,s) (X i , D n,-k i ) : s ̸ = r). We include r in the model confidence set Θ if and only if this confidence lower bound is non-positive.

Similarly define ℓ (r,s) r,s) ntr be the empirical standard deviation estimated from {ℓ (r,s) (X i , D n,-k i ) : 1 f i f n}, and Γ(r) ntr the corresponding (m -1) × (m -1) empirical correlation matrix.

To construct a confidence lower bound for max s̸ =r R(r) cv,n,K -R(s) cv,n,K , we only need the quantile of the maximum of a one-sided Gaussian vector. For a positive semidefinite matrix Γ, define u ´,Γ = upper ´-quantile of max N (0, Γ) .

(5.9)

Then we have the following result. , respectively, and (5.6) holds, then, for each r

Moreover, the confidence set

The proof of Proposition 5.2 is almost identical to that of Theorem 5.1 and is omitted.

Remark 5.2. As we mentioned in the motivation of the difference based CV model confidence set, the loss functions ℓ (r) (•, •) and ℓ (s) (•, •) are often positively correlated. If the estimators f (r) and f (s) are very similar, such as the same algorithm with two slightly different tuning parameter values, then the correlation can even be close to 1. Such a high correlation may lead to a very small standard deviation Ã (r,s) n of the difference loss function, and hence pose additional challenge for the stability condition (4.25) to hold. Nevertheless, (4.25) is still conceptually similar to the original version in that it requires an o(1/ √ n) factor after taking a ∇ i operator in the parameter estimation part. We will provide an example of the difference loss stability in Section 6.3.

Stopping rule in sequential model building. For example, in order to determine the number of steps in forward stepwise regression, in each step one can construct a confidence set by comparing the current model and all models with one more variable, then if the current model is included in the confidence set, it means going one more step will not significantly increase the performance and the procedure should stop.

3. Improved single model selection. In many scenarios, it is straightforward to pick a single element from the model confidence set. For example, one can pick a model having a particular desired feature suggested by subject knowledge. In the common situation of nested models, it is often beneficial to select the most parsimonious model in the confidence set. The motivation is a simple reflection of Occam's razor: With similar performance, the simpler one wins.

In this section, we will focus on a concrete example related to the third item in the list above: improved single model selection. What we will see is that, under fairly standard conditions in a classical linear subset selection problem, the most parsimonious model in the CV model confidence set can yield consistent subset selection while the ordinary CV fails.

To proceed, we consider a fixed-dimensional subset selection problem, where the data D n consists of n i.i.d. copies of X = (Y, Z) with

where Z ∈ R p is the covariate vector, f ∈ R p the regression coefficient, and • a mean-0 noise, independent of Z. Let J = {0, 1} p be the collection of all possible subsets of [p]. The goal is to recover J * , the support of f . It is well known, as illustrated in Example 3.1, that K-fold cross-validation that chooses

is generally inconsistent, regardless of the value of K. Instead, consistency requires reversed K-fold with a diverging K. The following result shows that, under certain regularity conditions, the most parsimonious element in the CV model confidence set (5.7) is consistent for conventional K-fold CV. In particular, let Ĵcvc = arg min{|J| : J ∈ Θ0 } .

Here the subindex "cvc" refers to "Cross-Validation with Confidence" (Lei, 2020). By construction, it is straightforward to check that Ĵcv is always contained in Θ0 (and Θ) as long as t ´,Γ n tr g 0. Therefore, Ĵcvc is guaranteed to be at least as parsimonious as Ĵcv , and can be viewed as a correction of Ĵcv , which is known to be prone to overfitting in practice. In this regard, Ĵcvc is a principled approach to reflect a common intuition about CV: The ordinary CV tends to select a slightly over-fitting model, so the best model is usually achieved by one that is slightly more regularized than the ordinary CV choice.

Proposition 5.3. Assume the linear regression model with a fixed distribution of (Z, •) that does not change with n, such that

and all entries of Z and • are (1, ³)-SW for an ³ > 0. Let ℓ(x, f ) = (yz T f ) 2 be the squared loss. If, for each J, f (J) (D n ) is supported only on J and satisfies the stability condition

for a sequence ϵ n = o(1); and for each

(5.13)

for all n large enough, and ´= ´n satisfies ´n = o(1) and ´n > 1/n then

Here, the inner ∥ • ∥ denotes the Euclidean norm of a p-dimensional vector.

Remark 5.3. Proposition 5.3 is stated with a set of pretty strong but reasonable conditions to deliver the desired model selection consistency result. These conditions can be verified in the setting of Example 3.1 to show that the most parsimonious element in the confidence set can lead to consistent (or bounded error probability if ´is a fixed constant) model selection with a standard K-fold cross-validation scheme. Nearly all the conditions can be relaxed to some extent with lengthier bookkeeping and/or more refined arguments. We list some possible relaxations below.

1. The dimensionality p, the smallest non-zero entry of f , and the eigenvalues of Σ Z , can be allowed to change with n at slow enough rates. These may require faster rates of the stability term. In particular, let ¼ min and f min be the smallest eigenvalue of Σ Z and the smallest non-zero entry of f , respectively. Instead of (5.13), now assume that

to ensure (5.15), and p 1/2 ϵ n = o(1) to establish (5.16).

2. The stability condition can be relaxed to L 2 -norm instead of sub-Weibull, with the addition of the uniform integrability as in Theorem 4.3 (which implies the Lindeberg condition (4.4)).

3. The requirement of J * = arg min µ

n can be relaxed to

using the Gaussian anti-concentration inequality.

Moreover, it is also possible to obtain a similar consistency of Ĵcvc for the difference-based CV model selection confidence set Θ in (5.10). The conditions and arguments would be similar to those of Proposition 5.3, with everything stated for the difference of the estimated parameters.

Now we examine the practical plausibility of the two conditions (5.12) and (5.13). First, (5.13) corresponds to the standard root-n-consistency for correctly specified linear regression models. The maximum likelihood estimator and least squares estimator often satisfy this in the fixed-p setting. To examine the stability condition (5.12), let f be an estimate of f obtained from D n using the least squares method. For notational convenience we will consider ∇ n instead of ∇ 1 . Let Σn be the sample covariance matrix from D n , and Γn = n -1 n i=1 Z i Y i . Then we have, by the Sherman-Morrison formula,

where only the second terms in each factor involve (Z n , Y n ), both of which are of order O P (1/n). Thus ∇ n f = O P (n -1 ), making (5.12) plausible for any ϵ n k n -1/2 .

Proof of Proposition 5.3. Let J 0 = {J ¦ [p] : J * ª J} and J 1 = J \J 0 . For an arbitrary J ∈ J , using f and

which is (cn -1/2 ϵ n , ³ ′ )-SW for some universal constant c and ³ ′ depending only on ³.

Therefore, the K-fold version of Theorem 2.6 implies that R(J) cv,n,K = µ (J) ntr + O P (n -1/2 + ϵ n ) for each J. For J ∈ J 0 , the corresponding f (J) satisfies µ (J) ntr g Var(• 0 ) + C for a positive constant C depending on the minimum eigenvalue of Σ Z and the minimum nonzero entry of f , while µ

As a result, with probability tending to 1 we have

Γn tr , (5.14) because the fluctuation Ã(J) ntr t ´n, Γn tr / √ n in the confidence set is of order at most log(p/´n)/ √ n, which is not enough to override the gap between µ (J)

ntr and µ

ntr . As a result we have 1) .

(5.15)

Now for J ∈ J 1 , we have, letting

which is (cn -1 ϵ n , ³ ′ )-SW. Thus using the deterministic version of Theorem 5.1 and Theorem 4.9, we have

where Γntr,1 is the submatrix of Γntr restricted to J 1 . Because t ´n, Γn tr g t ´n, Γn tr ,1 , (5.16) also holds for t ´n, Γn tr . Now combine this with (5.15) we conclude that P( Ĵcvc = J * ) = 1o(1).

Testing Many Means

In this subsection, we demonstrate an application of the tools developed in Section 4 in a prototypical inference problem at the intersection of high-dimensional inference, multiple comparison, and selective inference. Let (X i : 1 f i f n) be an i.i.d. sample in R m with EX 1 = ϑ = (ϑ 1 , ..., ϑ m ) T . We are interested in testing H 0 : max ϑ f 0 , vs H a : max ϑ > 0 .

(5.17)

This problem has many connections to other well-studied problems. Here we briefly mention two examples. A more comprehensive literature review is given in Section 5.6. One example is marginal regression, where we have a response variable and an m-dimensional covariate vector. Let each X i,r be the product of the response and the (centered) rth covariate in the ith sample point; then the problem (5.17) corresponds to testing whether any covariate has positive correlation with the response, and it can be extended to testing any non-zero correlation. Such a test is useful, for example, in deciding the stopping rule in forward stepwise regression. Another example is evaluation of many fitted prediction models. Assume there are m regression or classification probability functions, fitted from external training data. We want to compare their performance on a target data distribution from which we have n i.i.d. samples points. Then each X i,r corresponds to the realized loss function of the rth prediction function on the ith sample point. The testing problem (5.17) corresponds to testing whether any of these m prediction models are better than a benchmark, assuming, after re-centering, that the benchmark has 0 risk.

Our proposed solution to the testing problem (5.17) is based on constructing a confidence lower bound for the parameter ¹ = max ϑ and shares some similarity with the problem of selective inference. To begin with, observe that (5.17) is equivalent to

ϑ s or an arbitrary element in arg max s∈[m] ϑ s if the maximum is achieved at multiple entries. The optimal index s * is usually unknown and can be regarded as an auxiliary parameter (or nuisance parameter), the knowledge of which would make the inference problem much easier. If s * were known, then one can apply any one-sample location test such as Student's t-test or its robust versions to solve (5.17). When s * is unknown, the ordinary one-sample location test applied to ŝ, defined as the argmax of the sample mean θ, will generally lead to biased inference due to the "double-dipping" or "selection bias" issue, meaning that the same data is used to select the index ŝ and to compute the test statistic. Thus any valid inference must take into account the randomness in ŝ. Existing approaches include the following.

• Selective Inference. This framework aims at making the p-values conditionally valid on the event ŝ * = s for each s ∈ [m]. These methods typically require non-trivial assumptions on the data distribution and the power can be less than ideal due to the conditioning.

• Simultaneous Inference. This framework relies on controlling the simultaneous fluctuation of the sample mean vector θ, such as in the high-dimensional Gaussian comparison result in a similar fashion as in Proposition 5.2. This approach does not require conditioning but the worst-case bound could be loose if there are many "obviously bad" coordinates.

Our approach aims at striking a balance between these two: We want to focus on the few entries that showed enough promise to be legitimate candidates for s * , which allows us to avoid having the typical √ log m additional width of the confidence bound as in the simultaneous inference. On the other hand, we also want to avoid conditioning on selection. The key idea is that conditioning is unnecessary if the selection is stable and the test statistic is cross-validated. In fact, we have more to offer: If the selection is stable and the test statistic is cross-validated, a confidence lower bound on max ϑ can be constructed from an asymptotically normal test statistic with a small conservative bias.

Cross-validated max-entry evaluation. The first idea for fixing the double-dipping issue with data-driven max-entry selection is to split the sample used for selecting ŝ and to evaluate ϑ ŝ. Consider a leave-one-out CV scheme:

Unfortunately this Tn is not generally asymptotically normal. To see why, observe that for each i, the quantity

The stability condition (4.3) does not hold for this ℓ(x, D). In particular, consider the case of m = 2, ϑ 1 = ϑ 2 = 0, X i ∼ N (0, I 2 ). The contribution of sample point j in the comparison between θ(-i)

if X j is replaced by an i.i.d. copy X ′ j . When this happens, the resulting change in

). Even such a narrow miss of o(•) to O(•) would lead to a failure of the CLT, as demonstrated in the simulation in Figure 5.1, where the histogram of Tn in (5.18) is plotted in green, which clearly deviates from the standard normal density curve (black curve). This also demonstrates the sharpness of the stability condition (4.3) in Theorem 4.2 -even relaxing o P (n -1/2 ) to O P (n -1/2 ) may lead to failure of the CLT! Stabilizing the maximum value. The main reason for the instability of X i,ŝ (-i) is that ŝ(-i) may change completely when a single data point is perturbed. Such an instability of discrete argmax can be treated using techniques developed in differential privacy (Dwork and Roth, 2014). In general, differential privacy requires a statistic to be "insensitive" to perturbations at any individual sample point, which is closely related to the first-order stability condition. Inspired by the "exponential mechanism" (McSherry and Talwar, 2007), we consider the following cross-validated softmax statistic:

with weights ŵ(-i)

) .

The exponential mechanism stabilizes the maximum value by "softening" the indicator of the maximum entry to a weight vector, ŵ(-i) , with larger entries receiving more weight. By construction, the weights are smooth functions of the sample means, and hence enjoy a similar level of stability as the sample means. The quantity

where s(-i) is sampled according to the weight ŵ(-i) , which is a smoothed version of ŝ(-i) .

The coefficient ¼ is a tuning parameter in the exponential mechanism. A larger ¼ puts more weight on large entries, which has less bias but is less stable. In the extreme case

Theorem 5.3 is a direct consequence of Theorem 4.2. The choice of notation R i has a direct correspondence with R(D n,-i ) in the context of Theorem 4.2.

The conditions in Theorem 5.3, although directly corresponding to those in Theorem 4.2, are still somewhat abstract in the context of the many means testing problem. In particular, Theorem 5.3 does not specify which ¼ would lead to the stability condition

We provide a set of more transparent sufficient conditions in the proposition below.

Proposition 5.4. Let Σ = E(X 1ϑ)(X 1ϑ) T be the covariance matrix. Assume

• There exist a pair of positive constants c 1 , c 2 such that w

Remark 5.4. In Theorem 5.3 and Proposition 5.4, we allow for the triangular-array setting where m may change with n. Moreover, all the theoretical results and methods developed below can be extended to K-fold CV for any value of K. It would be possible to relax the L ∞ -norm bound on X i,s to sub-Weibull type bound using more refined arguments.

The choice of ¼ n = o( √ n) is used to cover the worst-case scenario.

In most realistic situations, such as when there are a few coordinates whose ϑ s values are substantially larger than the rest, then the stability condition can be satisfied by larger values of ¼ n . We discuss data-driven choices of ¼ below after proving Proposition 5.4 and describe the test procedure.

Proof of Proposition 5.4. We first show that Ã n = Ω(1). In fact

To simplify notation, we write ŵ(-i)

as w s , and ŵ(-i,j)

obtained from replacing X j by X ′ j . We also write θ(-i) as θ, and θ′ the corresponding perturb-one version. By construction

Recall that n tr = n-1 since we are considering LOOCV. By construction we have for each s ∈ [m]

Similarly we can show w ′ s ws g exp(-4¼ n n -1 tr B), and the claimed result follows.

Remark 5.5. One may wonder whether the exponential weighting trick alone can fix the failure of CLT in the statistic Tn in (5.18) without cross-validation. That is, consider

Qi with Qi = s∈ [m] ws X i,s and ws ∝ exp(¼ θs ). The empirical distribution of Tn is plotted as the blue histogram in Figure 5.1, showing a clear shift from the desired standard normal curve. This demonstrates that both cross-validation and exponential weighting are necessary to achieve asymptotic normality for the (randomly centered) statistic T n in (5.19).

Inference. As seen in the proof of Proposition 5.4, we have Ã 2 n ≍ 1, and

Thus the conditions required in Theorem 4.3 are also satisfied, which leads to the following.

Corollary 5.4. Let Ãn be the empirical variance of (Q i : i ∈ [n]). Under the conditions in Proposition 5.4, Ãn /Ã n → 1 in probability, where

Moreover, for ´∈ (0, 1), let u ´be the (1-´)-quantile of the standard normal, then the statistic T n defined in (5.19) satisfies

(5.20)

Proof of Theorem 5.4. The consistency of variance estimation follows directly from Theorem 4.3. Now we prove the second claim.

where the inequality follows from the fact that R i f max ϑ, and the last equation follows from Theorem 5.3. Therefore, the following rejection rule is asymptotically valid for the testing problem (5.17) at nominal level

´,

Bias and power analysis. In the confidence lower bound in Theorem 5.4, the term due to the variance, which is Ãn u ´/√ n, does not involve the dimensionality m. This is in sharp contrast with the simultaneous inference approach, which, in the case of independent coordinates, would typically involve a √ log m factor associated with the variance term. Such a reduction in variance is gained at a price of a bias term max ϑ-n -1 i∈[n] R i . Empirically the cross-validated and exponentially weighted confidence lower bound in (5.20) with an adaptively chosen ¼ (see below) usually achieves a favorable bias-variance trade-off when compared to other methods. Here we provide a simple upper bound of this bias under a boundedness condition. Proposition 5.5. If |X i,s | f B < ∞ almost surely for all i, s, and

Then there exists a constant C such that

Proof of Proposition 5.5. Without loss of generality, assume

and s ∈ T if we pick C large enough in the definition of T . Therefore, on event E, we have

Proposition 5.5 implies that the test statistic T n is asymptotically powerful against alternatives with signal size É(log(m + n)/¼ n ). See Section 5.6 for further references on power analysis and numerical results.

Choosing the tuning parameter. As mentioned in the discussion after Remark 5.4, the requirement of

4 is for the worst-case scenario. For an example of a non-worst-case scenario, suppose the largest entry ϑ s * k ϑ s for s ̸ = s * , the selection of arg max θs would be stable even without the exponential mechanism. In general, the optimal choice of ¼ n should be as large as allowed by the stability condition, and is often a complex function of the mean vector ϑ and covariance matrix Σ.

Suppose we are given a finite collection L of candidate values of ¼, a plausible strategy is to select the largest one in this set such that the stability condition holds approximately.

By the insight in Proposition 4.1, the stability condition we want to check for each ¼ ∈ L is the following:

for j ̸ = i, where Ã 2 ntr is defined as in Theorem 5.3. The quantity Ã 2 ntr can be estimated using the sample variance of (Q i : 1 f i f n) as in Theorem 4.3. Here we only need to estimate

2 . We use a leave-two-out bootstrap estimate. For a distinct triplet

Thus we estimate

dently sampled index triplets. Here B is the bootstrap sample size. Denote this estimate as ¶2 (¼) for each ¼ ∈ L.

Finally, we choose ¼ as the largest value in a given candidate set such that, for some small constant ε (e.g., ε = 0.05),

where Ã2 (¼) is the corresponding estimate of Ã 2 ntr .

Let

In the above description of cross-conformal prediction, we omitted the rounding, ignored the 1/(n te + 1) term in the definition of pk (y, z), and changed the normalization factor from 1/(n te + 1) to 1/n te . These will only cause negligible differences in the resulting quantities and will disappear in the asymptotic analysis below. In return, we gain a perfect matching between the cross-conformal p-value pcc and a general cross-validation statistic considered in this article.

Due to the dependence caused by the fold rotation, the fold-specific conformal p-values pk (y, z) are often dependent and hence it is hard to claim the finite-sample distribution-free coverage guarantee (5.23) for the cross-conformal prediction set Ĉcc . However, in practice it is often observed that cross-conformal prediction sets offer accurate coverage with favorable efficiency (measured as the size of the prediction set) compared to the split conformal prediction set, thanks to a more efficient use of the data.

Here we use the tools developed in Section 2 and Section 4 to prove an asymptotic coverage guarantee for the cross-conformal prediction under two mild assumptions. In the rest of this subsection, we will use x to denote (y, z), and X i = (Y i , Z i ). Theorem 5.6. Under Assumption 5.1, if ∥∇ 1 s(X 0 ; D n )∥ q f ϵ n for q g 1 and a sequence ϵ n , then the cross-conformal p-value satisfies

As a consequence, if

In practice, it is often the case that ϵ n ≳ n -1 , thus the condition ϵ n = o(n -1/2-1/q ) is plausible for q > 2. The requirement of ϵ n = o(n -1/2-1/q ) can be further relaxed if we strengthen the stability condition to ∇ 1 s(X 0 ; D n ) being (ϵ n , ³)-sub-Weibull. In this case it suffices to require ϵ n = õ(n -1/2 ), where õ(•) means "o(•) up to a polylog factor".

Proof of Theorem 5.6. Define

Here we keep in mind that such an ℓ(•; •) depends on an external random quantity X n+1 .

The proof consists of two steps. The first step is to show that, in the same spirit of Theorem 4.2,

In fact, letting

≲n -1 + ϵ q/(1+q) n .

(5.24)

To bound EÀ i À j in the final step of the inequality above, observe that EÀ i À j ̸ = 0 only if k i ̸ = k j and in this case using the "free-nabla" lemma (4.11), we have, using the bounded density assumption,

where the last equation follows from Theorem 5.7.

The next step is to show that

By construction, for each individual k, s(X 0 ; D n,-k ) and s(X n+1 ; D n,-k ) have the same conditional distribution given D n,-k . As a result we have

We have, using the martingale decomposition and Efron-Stein (Theorem 4.1)

Therefore, for any pair of k, k ′ ∈ [K] we have

which further implies that

) .

(5.26) The claimed result follows from (5.24) and (5.26) because (5.25) implies that Q(X n+1 ; D n,-1 ) ∼ U (0, 1).

Lemma 5.7. Under Assumption 5.1, there exists a constant c 1 such that

Proof of Theorem 5.7. We use shorthand notation to denote s(

By definition

(5.27) because the binary difference is nonzero only if the gap s(X 0 ) -s(X n+1 ) is dominated by the perturbation, whose absolute value equals ∆. Now we prove Part (a) as follows. For any ϵ > 0 we have

where the first inequality follows from (5.27), the second from union bound, the third from the bounded conditional density of s(X 0 ) given everything except X 0 , and the last by choosing ϵ = ϵ q/(1+q) n . For Part (b), the proof is similar except that taking expectation over X 0 before taking the L 2 norm increases the stability. For any ϵ > 0 we have

Now using triangle inequality we have

The claimed result follows by choosing ϵ = ϵ q/(2+q) n

Testing Many Means

The inference of maximum values in a vector considered in Section 5.3 is connected to many threads in the literature. The sample-splitting and stabilization idea is partially inspired by the "adaptive bootstrap" approach in McKeague and Qian (2015) for detecting marginal correlations. The connection to selective inference is straightforward; a general introduction to selective inference can be found in Taylor and Tibshirani (2015). A more directly relevant selective inference work on inference for the max value is Reid et al. (2017).

It is natural to extend the maximum value inference to high-dimensional mean testing, for which a vast literature exists; see Huang et al. (2022) for a recent overview of this topic. A particularly relevant work is Chernozhukov et al. (2019), who used high-dimensional Gaussian comparison to achieve asymptotically valid inference.

The argmin confidence set problem considered in Section 5.4 also has many connections in the literature and is a subject of active development. The problem of identifying the optimal population in a finite collection from noisy observations has been studied as early as Gibbons et al. (1977) and Gupta and Panchapakesan (1979), with further developments such as Futschik and Pflug (1995), Hansen et al. (2011), Mogstad et al. (2024), Dey et al. (2024), and Kim and Ramdas (2025). Other related problems include discrete MLE (Hammersley, 1950;Choirat and Seri, 2012) and rank verification (Hall and Miller, 2009;Xie et al., 2009;Hung and Fithian, 2019).

The testing-many-means problem and the related argmin inference problem can also be viewed from a nuisance parameter perspective. Let S m-1 be the (m -1)-dimensional simplex. Then for any w ∈ S m-1 define ¹ w = ϑ T w f max ϑ. Then any valid confidence lower bound of ¹ w is also a valid confidence lower bound for max ϑ. The method developed in Section 5.3 uses a data-driven version of a particular w * , the indicator of the maximum entry of ϑ:

This idea can be taken one step further: the current target w * aims to maximize w T ϑ and does not account for the variance of w T X. A strictly more efficient target w is

where Σ = E(X -ϑ)(X -ϑ) T is the noise covariance matrix. The objective function in the definition of w above is known as the Sharpe ratio (Palomar, 2025). It can also be modified to other similar forms such as w = arg max

for some penalty parameter ¼. Then it would be possible to extend the procedure in Section 5.3 to such weight vectors w , with a potentially improved bias-variance trade-off.

The softmax statistic Q i in Section 5.3 is an example of making and using stable algorithms for statistical inference. The idea of using stable nuisance parameter estimates in semiparametric inference without de-biasing or cross-fitting has been explored in Chen et al. (2022). In particular, the use of soft-min for stabilized optimal policy choice has been studied in Chen et al. (2023). The presented method of data-driven choice of ¼ works by numerically verifying the desired stability level. An important theoretical and methodological research topic is to estimate and test the stability of a given algorithm. This is challenging; see Kim and Barber (2023) and Luo and Barber (2024) for some early results.

In the next section we present additional results on stability implied by regularization.

Conformal Prediction

For a general introduction to conformal prediction, we refer the reader to the books Vovk et al. (2005) and Angelopoulos et al. (2024). Some early references on conformal prediction for outlier detection and regression include Lei et al. (2013), Lei and Wasserman (2014), and Lei et al. (2018). The split conformal method is also known as inductive conformal prediction (Papadopoulos et al., 2002;Lei et al., 2014). The cross-conformal prediction method was introduced by Vovk (2015) and analyzed in Vovk et al. (2018), achieving a coverage level of approximately 1-2³. A variant of cross-conformal, called CV+, was introduced and analyzed in Barber et al. (2021), with asymptotic coverage 1 -³ under a similar but different set of stability and bounded density conditions. A key component in our proof for cross-conformal is the transfer of stability from the conformity score function to the ranking indicator function (Theorem 5.7), where the inequality (5.27) is adapted from a similar analysis in Hu and Lei (2024).

In this section we examine the stability of several general classes of estimators. To simplify discussion, we will focus on the case where the loss function ℓ(x, f ) is Lipschitz in f , so that the stability of loss function is implied by the stability of f itself. Although we do not specify the parameter space F, it would be helpful to think of it as either a finite-dimensional Euclidean space or a Hilbert space so that the gradients and Hessians can be defined in the usual sense.

So far, the least squares regression estimate in Section 5.2 is the only example in which stability has been checked. This example can be extended to the general class of asymptotically linear estimators. Recall that an estimator f is asymptotically linear with influence function h

where the o P (n -1/2 ) should be interpreted in terms of the norm of the residual vector. Then for such asymptotically linear estimators, we have

507 However, this is still heuristic, as o P (n -1/2 ) is still weaker than the L 2 and sub-Weibull stability conditions required in the previous sections. Typically, establishing such stability conditions is nontrivial and may require specific tools depending on the model and estimator. For example, random matrix theory may be useful in establishing L 2norm stability of the least squares estimate.

We will consider the more general and abstract approach of establishing stability of estimators via regularization. For example, penalized M -estimators such as ridge regression and support vector machines are known for their stability due to their strongly convex optimization formulation. In Sections 6.1 and 6.2, we make this more rigorous and general for the stochastic gradient descent algorithm. In Section 6.3, we illustrate the subtlety of the stability of the difference loss functions (see Section 5.1 for the context) in a prototypical nonparametric regression example, where the stability is achieved by another form of regularization: truncated basis expansion. In Section 6.4 we consider the stability conditions involved in the consistency analysis of online rolling validation in Section 3.3, focusing on a combination of stochastic gradient descent and nonparametric basis expansion, known as the sieve SGD (introduced in (3.10) in Section 3.2). Other related results in analyzing the stability of M-estimators, as well as alternative ways to construct stable estimators, such as resampling and bootstrap-aggregating, are discussed in Section 6.6.

Stability of Stochastic Gradient Descent

Consider the M-estimation problem

for some objective function φ : F × X → R, and the stochastic gradient descent algorithm (SGD):

1. Initialize f0 = 0 2. For t = 0, ..., n -1, update ft+1 = ft -³ t+1 φ( ft ; X t+1 ) , where (³ t : 1 f t f n) is a sequence of learning rates, and φ(•; x) is the gradient of φ(•; x) viewed as a function of f . 3. The final estimate f (D n ) = fn .

One may also consider the averaged final estimate f (D n ) = n -1 n t=1 ft . The stability results derived below can be carried over to the averaged case. The SGD estimate may appear to violate the symmetry condition in Definition 2.1. This issue can be resolved by randomizing the order of sample points in the SGD update sequence.

The motivation of analyzing (approximate) M-estimators obtained by SGD is quite natural. First, M-estimators cover a wide range of parametric and nonparametric estimators such as the maximum likelihood estimator, the least squares estimator, local polynomial regression, their robust counterparts, as well as machine learning tools such as kernel regression, support vector machine, and neural networks. Second, the SGD algorithm is arguably the most popular method for large-scale Mestimation. Under convexity and smoothness conditions on the objective function φ, the simple form of SGD allows precise characterization of the algorithmic behavior of the estimated parameter and hence facilitates a general and abstract stability analysis.

As seen in the ridge regression, the strong convexity of the objective function φ(•, x) plays an important role in the stability of the resulting estimator. Definition 6.1 (Strong convexity). Let µ be a positive number. A function h :

, where ∥ • ∥ denotes the Euclidean or Hilbert space norm. Assumption 6.1 (Smoothness and convexity of objective function). There exists positive constants (µ, C 0 , ´) such that for all x ∈ X

Part (c) of Assumption 6.1 is also known as "´-smoothness". These are among the standard assumptions in the analysis of the convergence and stability of SGD. In deep neural networks, it is unrealistic to assume convexity of φ(•; x), but we can assume that these conditions hold within a local convergence basin. Parts (a) and (c) require that µ f ´. When µ = ´, φ(•; x) is close to a quadratic function. The gap between µ/á nd 1 reflects how close φ(•; x) is to quadratic; the conditions are more stringent when a large value of µ/´is required.

Example 6.1 (Ridge Regression). Consider x = (y, z) and X ¦ B(r x ), F = B(r f ) (B(r) denotes a centered ball with radius r.)

Then φ(•; x) satisfies Assumption 6.1 with µ = ¼, C 0 = r 2 x (1 + r f ) + ¼r f , and ´= r 2

x + ¼.

In practice, the objective function φ(•, ; x) usually is not strongly convex without a regularizing penalty term such as the ridge penalty. It is still possible to establish similar stability results for SGD if E X φ(•; X) is strongly convex. We provide a detailed example in Section 6.4 for the sieve SGD estimator introduced in Section 3.2.

The main consequence of Assumption 6.1 is the following contractivity result of SGD. Proposition 6.1 (Contractivity of SGD). Assume φ satisfies parts (a), (c) of Assumption 6.1, and ³ ∈ (0, 2/(´+ µ)). The SGD update mapping

The proof of Proposition 6.1, though elementary, is non-trivial. We will not provide a proof here. See the bibliographic notes in Section 6.6 for further references. However, a one-line proof can be obtained, with a different contraction factor, if we assume that φ(•, x) is twice differentiable: Let φ(•, x) be the Hessian matrix, which is lower bounded by µI because of strong convexity. Then, by the mean value theorem we have

where [u, v] denotes the line segment joining u and v. Theorem 6.1 (First-order stability of SGD). Under Assumption 6.1, if the learning rate satisfies ³ t = ´-1 t -a for an a ∈ (0, 1) and

The intuition behind Theorem 6.1 is straightforward. If X i is perturbed to X ′ i , the SGD algorithm will remain the same up to time i, and future updates apply to the perturbed solution at time i. If i is small, X i appears early in the SGD iteration, and the perturbation on fi will shrink exponentially in later steps due to the regularity conditions in Assumption 6.1. If i is large, the learning rate is already very small, limiting the perturbation from changing X i to X ′ i . In classical theoretical convex stochastic approximation, the most popular learning rates are a ∈ (1/2, 1], ensuring t ³ t = ∞ and t ³ 2 t < ∞ (Nesterov, 2013;Bubeck, 2015). The range of a allowed by Theorem 6.1 covers ³ ∈ (1/2, 1) but not ³ = 1. This is due to the much weaker contraction upper bound, since t i=1 i -a grows polynomially in t when a < 1 while only logarithmically when a = 1. A potential sharper stability result for a = 1, or showing its absence, remains an open problem.

As mentioned in the intuition above, the proof of Theorem 6.1 is a direct consequence of the following lemma, which establishes the multiple-step contractivity of SGD updates under Assumption 6.1. Lemma 6.2. Under Assumption 6.1 and the conditions in Theorem 6.1, we have for all t g i

Proof of Theorem 6.2. By construction, ∥∇

. By contractivity of SGD (Proposition 6.1), for all t > i,

If i f n/2, according to Theorem 6.2 with t = n we have

Remark 6.1 (Extension to Multi-epoch SGD). In practice, SGD is more commonly implemented in multiple epochs. The results and proof technique of single-epoch SGD presented above may be useful in establishing the stability of multi-epoch SGD algorithms. For a simple example, consider a second epoch using the same learning rate sequence. Let f (e) t denote the estimate f at step t in epoch e with e ∈ {1, 2}. Let f (e),i t be the corresponding estimate with X i replaced by X ′ i in both epochs, and f (e),i t the estimate with X i in the first epoch and X ′ i in the second epoch. Then ∥∇ i

≲ n -a . The same argument extends to any finite number of epochs. A more refined stability analysis of multi-epoch SGD under general learning-rate schemes is an interesting direction for future work.

Second-order Stability

The second-order stability concerns the quantity

and is used in the deterministic centering CLT in Section 4.2. Although the discussion in Section 4.2 emphasizes that it is more important to compare

We first show how to control the second-order loss stability through the second-order parameter stability. Proposition 6.2. If the loss function ℓ(x, f ) : X ×F → R is c 0 -Lipschitz, and l(x, •) is c 1 -Lipschitz for all x ∈ X , then

where ϵ n,i,j = max{∥∇ j fn ∥ , ∥∇ i f j n ∥ , ∥∇ j f i n ∥}, and f i n ( f j n ) denotes the estimate obtained with X i (X j ) replaced by its i.i.d. copy.

One example of such ϵ n,i,j is given in Theorem 6.2, which implies that ϵ 2 n,i,j has order n -2a . The main result in this subsection is to show that ∥∇ i ∇ j fn ∥ can also be bounded by n -2a up to a logarithmic factor.

Proof of Proposition 6.2. Let f ij n be the estimate obtained with (X i , X j ) replaced by their i.i.d. copies. By the mean value theorem, there exists

The claimed result follows from the fact that

Our main additional condition for the second-order stability of SGD estimates is the smoothness of Hessian of φ(•; x). Assumption 6.2. The Hessian matrix of φ(•; x), denoted φ(•; x), is C 2 -Lipschitz for a constant C 2 and all x:

It may be possible to avoid such a smooth-Hessian condition using a similar co-coerciveness argument as in the proof of Proposition 6.1. Here we work under Assumption 6.2 for simplicity and transparency. Example 6.2. Ridge logistic regression: x = (y, z) ∈ B(r x ), and f ∈ B(r f ).

Then φ satisfies Assumption 6.1 and Assumption 6.2 with µ = ¼, C 0 = r x + ¼r f , ´= r 2 x /4 + ¼, and C 2 = r 3 x .

Theorem 6.3 (Second-Order Stability of SGD). Under the same assumptions as in Theorem 6.1, assume in addition that φ satisfies Assumption 6.2 such that ³ t µ f 1 for all t, then

for a constant c depending on (a, ´, µ, C 0 , C 2 ).

The proof of Theorem 6.3 is lengthy and deferred to Section 6.5. Similar to the proof of first-order stability, the key idea is to establish a contractivity result analogous to (6.1) for the second-order difference. This is enabled by the smoothness of the Hessian (Assumption 6.2). A technical challenge is handling the remainder term of the first-order Taylor expansion of φ, which must be tracked across all SGD iterations.

Difference Loss Stability

According to Remark 5.2, a technical condition needed for the validity of the difference-based CV model confidence set involves stability conditions on the standardized loss differences, which is used to establish high-dimensional Gaussian comparison for the cross-validated risk differences. This often becomes tricky because the loss differences usually have vanishing variances, and standardization means dividing the difference loss function by a vanishing standard deviation, which could make the stability condition less likely to hold. However, a more optimistic argument is that this condition retains its original intuition: taking a ∇ operator reduces the magnitude of the difference loss by more than a factor of 1/ √ n, even though the difference loss itself may be vanishing. In this subsection we illustrate this intuition through a simplified nonparametric regression problem using truncated basis expansion estimates.

Ignoring constant factors, we use the symbol "≲" to denote "upper bound up to a constant factor" and "a ≍ b" for "a ≲ b and b ≲ a".

Let (X i :

The assumptions on the data distribution are the following.

Assumption 6.3. 1. ϵ i is independent of Z i , has mean zero, finite variance Ã 2 ϵ , and is uniformly bounded by a constant;

2. Z ij is mean zero, has unit variance, is uniformly bounded by a constant, and the components are mutually independent;

3. The regression coefficients f j ≍ j -1+a 2 for some a > 0 for all j g 1.

The conditions in Assumption 6.3 are stronger than technically necessary and are made for a simpler presentation. The uniform bound on ϵ and Z can be relaxed to sub-Weibull tail conditions, and the independence between the entries of Z can be relaxed to E(Z ij |Z i1 , ..., Z i,j-1 ) = 0 for all j. The independence between ϵ and Z can be relaxed in a similar fashion. The assumed exact rate of f j in Part 3 is more stringent than the usual Sobolev ellipsoid condition which only assumes upper bound on |f j |. We need the lower bound on |f j | to track the bias and variance precisely. More discussion can be found in Section 6.6.

For each j, we can naturally estimate f j by

We consider a truncated estimate of f = (f j : j g 1). Let J be a positive integer,

where f (J)

The tuning parameter J controls the number of basis functions used in the truncated estimation.

The optimal choice of J depends on the sample size. For each candidate fitting procedure r, let J r be a positive integer varying with the sample size n such that lim n→∞ J r = ∞ at a speed specified by the estimator r. For example J r = n Är , where {Ä r : r ∈ [m]} is a collection of positive numbers corresponding to different growth rates of J r as the sample size n increases. We will use f (r) to denote f (Jr) for convenience and the estimators will be evaluated under the squared loss: ℓ((y, z), f ) = (y -f (z)) 2 . Now we are ready to quantify the stability conditions for the difference loss ℓ (r,s) (X 0 , D n ) and the difference risk R (r,s) (D n ). Proposition 6.3 (Difference Loss Stability). Assume the regression model (6.2) and the conditions for ϵ, Z, f in Assumption 6.3. If J r < J s and J r /J s either converges to 0 or stay bounded away from 0 as n increases, then there is a positive constant » depending only on the model such that (a)

The proof of Proposition 6.3 is deferred to Section 6.5. The idealized analytical form of the model and the estimator enables us to track bias and variance terms precisely, as well as their perturb-one stability.

The optimal choice of J with minimum risk is J * ≍ n 1/(1+a) . Thus, for the scaling parameter on the RHS of (6.4) and (6.5) to be o(1), one needs J s J a r j n = (J * ) 1+a . This is possible if both J s and J r are no larger than J * , or J s is slightly larger than J * but J r is sufficiently small. In other words, the stability condition is more likely to hold when comparing over-smoothing estimators. This condition is violated if J r g J * , suggesting that the stability condition required by the CLT may not hold for under-smoothing estimators. Such subtleties in the stability condition of the difference loss resemble the well-known behavior of the regular CV: it can eliminate inferior estimators with large bias (over-smoothing) but struggles to distinguish slightly overfitting estimators.

When J r /J s stays bounded away from 0, the sufficient stability condition required by the CLT becomes J s j J * , meaning that if the two estimators under comparison are too similar, the stability conditions can only hold for over-smoothing estimators.

Stability of Online Algorithms

In this subsection, we provide examples of the stability of online algorithms, specifically addressing (3.13) in Assumption 3.3. We focus on a simplified setting where the basis functions are given by the individual coordinates of the covariate vector Z, that is, ϕ j (Z i ) = Z ij . Under this simplification, the basis expansion estimator in (3.10) reduces to the setting considered in Section 6.3.

Notation Declaration: In this subsection we will explicitly track the evolution of each coordinate of f over time. Thus, we need to use the following notation.

• f = (f j : j g 1) denotes the true vector of coefficients.

• fi = ( fi,j : j g 1) denotes the estimated vector of coefficients at the end of step i, whose jth coordinate is fi,j . Note that this notation fi has a different meaning from the notation " fj " in Section 6.3 (where fj denotes the jth coordinate of the estimated coefficient f .) Generally speaking, in this subsection the symbol "i" (and similarly "t" and "n") is used to denote a time point or sample size, whereas the symbol "j" is used to denote a coordinate of Z or f .

• f ′ i and f ′ i,j correspond to the perturb-one version of fi and fi,j respectively, with the perturbed sample point clear from context.

• (J n : n g 1) is a sequence of positive numbers indicating the number of basis functions used by sieve SGD at sample size n.

Batch truncated series estimate. To warm up, first consider the truncated estimator f (Jn) given in (6.3). The derivation from the previous subsection implies that E [∇ i ft (X t+1 )] 2 |D i ≲ J t t -2 . Thus if J t ≍ t Ä then (3.13) holds with b = 1 -Ä /2. If the true coefficient f belongs to a Sobolev ellipsoid:

for a pair of positive constants (ϖ, L), then the minimax optimal error rate is achieved by J n ≍ n 1/(1+2ϖ) . Now consider two estimators with basis cardinality sequences J n,1 ≍ n Ä 1 and J n,2 ≍ n Ä 2 respectively, with J n,1 being the better choice. The sufficient stability condition for wRV to consistently select the better sequence J n,1 is quite similar to that in Section 6.3: Ä 2 must be either (i) smaller than Ä * = 1/(1 + 2ϖ) while Ä 1 is equal to or slightly larger than Ä * , or (ii) slightly larger than Ä * while Ä 1 is sufficiently small. In other words, wRV, like CV, is theoretically more reliable when eliminating highly biased estimators (small Ä 2 versus moderate Ä 1 ), or when comparing two very different estimators (moderately large Ä 2 and small Ä 1 ).

Sieve SGD. To illustrate a more challenging and less explored regime of online estimator stability, we consider the sieve SGD estimator defined in (3.10). Due to its relatively recent introduction, theoretical analyses of convergence rates for general sieve SGD algorithms with non-optimal tuning values are not yet available. Thus, the stability analysis of sieve SGD is still an ongoing line of research. Here we provide some initial results in this direction.

In the simplified setting with data generated from (6.2) and basis

where ³ i is the learning rate, Λ i is a diagonal truncation/shrinkage matrix whose jth diagonal entry equals j -2w 1(j f J i ), and J i is the number of basis functions used at sample size i. As discussed in Section 3.2, the tuning parameter of sieve SGD is the entire sequence (J i : i g 1), which typically takes the form J i ≍ i Ä . The multiplier j -2w provides coordinate-wise shrinkage, which was introduced in the original sieve SGD algorithm.

Proposition 6.4. Assume i.i.d. data from (6.2) satisfying Assumption 6.3. Let { fi : i g 1} be the sieve SGD estimate (6.7) with basis size J i = i Ä , learning rate ³ i = ci -a , and shrinkage exponent w > 1/2, such that a + 2wÄ < 1 and c < jg1 j -2w . Suppose |Y i -ï fi-1 , Z i ð| is uniformly bounded by a constant almost surely for all i g 1 then we have, for any t g i g 1

Proposition 6.4 is an example of stability without strong convexity. It should be compared to Theorem 6.1, which shows a similar rate under strong convexity of the objective function. The argument relies on taking the conditional expectation of the SGD objective function given the previous estimate, and exploiting the strong convexity of this conditional expectation. The proof of Proposition 6.4, along with discussion of potential improvements, is presented in Section 6.5. Remark 6.2. Proposition 6.4 involves several technical conditions, mainly for simplicity of presentation. The most restrictive assumption is the uniform boundedness of the residual Y i -ï fi-1 , Z i ð. This can be ensured if we further restrict the estimates fi to a bounded set such as a Sobolev ellipsoid, realizable by a further projection step at the end of each iteration. Alternatively, one may obtain weaker versions of Proposition 6.4, such as moment or sub-Weibull bounds instead of an L ∞ -norm bound. Remark 6.3. In order to illustrate the relevance of Proposition 6.4 in the context of Theorem 3.2, we need to characterize the risk convergence rate of sieve SGD estimates obtained from general combinations of tuning parameters. Such results are not yet directly available in the current literature due to the recent appearance of sieve SGD. However, there is evidence that sieve SGD may have the same convergence rate as the batch mode truncated basis expansion estimate using J n basis functions at sample size n. Proceeding with this hypothesis and assuming the true coefficient f belongs to the Sobolev ellipsoid W (ϖ, L) in (6.6), consider two sieve SGD estimators with J i = i Är , ³ i = ci -a for r = 1, 2. Suppose Ä 1 = Ä * = 1/(1+2ϖ) is a better choice than Ä 2 . Using the batch mode error rate and stability bound in Proposition 6.4, the stability requirements in Theorem 3.2 reduce to 1. a > 1/2 + ϖÄ 2 when Ä 2 < Ä 1 (under-fitting); 2. a > 1 -Ä 2 /2 when Ä 2 > Ä 1 (over-fitting).

A caveat here is that such choices of (a, Ä 1 , Ä 2 ) are incompatible with the requirement of a + 2wÄ < 1 in Proposition 6.4 unless we allow for w < 1/4. We further conjecture that this requirement is an artifact of the current proof and can be improved using a more refined argument. Intuitively, a larger w implies stronger shrinkage, hence greater stability. Evidence supporting this conjecture is that one can establish consistency of sieve SGD with w = 0. See Section 6.6 for relevant literature.

Technical Proofs

Proof of Theorem 6.3. Assume i < j. Let f i t be the estimate at time t with X ′ i as the ith sample point. Define f ij t correspondingly. Then

where the last inequality follows from the Lipschitz property of φ(•; x). For t > j,

(6.9)

Using the mean value theorem, there exist ft-1 ∈ [ ft-1 , f i t-1 ] and

Plugging this into the RHS of (6.9), we obtain

Let ¸t = ∥∇ i ∇ j ft ∥, we get (6.10) where the second term on the RHS uses Lipschitz property of φ (with Lipschitz constant C 2 ).

To control ft-1 -f ′ t-1 , apply Theorem 6.2 directly to ft-1 -f j t-1 and f i t-1 -f ij t-1 , and then use triangle inequality on ft-1 -f ij t-1 = ( ft-1 -f j t-1 )+( f j t-1 -f ij t-1 ) and f i t-1 -f j t-1 = ( f i t-1 -ft-1 )+( ft-1 -f j t-1 ), then we have,

where c = µ (1-a)(´+µ) and ¶ t,j (c) = exp -c((t + 1) 1-a -(j + 1) 1-a ) .

For t g j + 1, define ϱ i,j,t = 8C 2 C 2 0 ´3 t -a ¶ t-1,i (c)i -a + ¶ t-1,j (c)j -a ¶ t-1,i (c)i -a , then (6.10) implies, for t g j + 1, ¸t f (1 -³ t µ)¸t -1 + ϱ i,j,t . (6.11)

We now seek to upper-bound ¸n using the recursion (6.10). To do so, let c ′ = µ/ [(1 -a)´] g c. Thus

¶ n,k (c ′ ) ¶ k-1,i (c) ¶ k-1,j (c)k -a i -a j -a (I)

, (6.12)

where the last inequality appeals to the fact that for positive integers j < t, t k=j+1 (1 -³ k µ) f exp -µ (1-a)´ (t + 1) 1-a -(j + 1) 1-a . Term (I): Using the fact that ¶ k-1,i (c) g ¶ k-1,i (c ′ ), ¶ n,k (c) ¶ k-1,i (c) f ¶ n,i (c)e c(2 1-a -1) , and ¶ k-1,j (c) f 1, we have

´3 e c(2 1-a -1) i -a j -a ¶ n,i (c)

(1 -a)´3 e c(2 1-a -1) i -a j -a ¶ n,i (c)(n 1-a -j 1-a ) . (1 -a)´3 e c(2 1-a -1) i -a j -a n -(3-a) (n 1-a -j 1-a )

(1 -a)´3 e c(2 1-a -1) n -2 . (6.14)

When i g n -c 2 (n + 1) a log n. Suppose n is large enough that n -c 2 (n + 1) a log n g n/2. Use the following upper bounds in the RHS of (6.13): k -a f j -a f i -a f 2 a n -a , ¶ n,i (c) f 1, we have

(1 -a)´3 e c(2 1-a -1) 2 3a n -3a (n -j)

(1 -a)´3 e c(2 1-a -1) 2 3a n -3a (n -i)

(1 -a)´3 e c(2 1-a -1) 2 3a n -3a c 2 (n + 1) a log n

(1 -a)´3 e c(2 1-a -1) 2 3a+1 c 2 n -2a log n . (6.15)

Combining (6.14) and (6.15), we conclude that (I) f c(a, ´, µ, C 0 , C 2 )n -2a log n , (6.16)

for n large enough.

Term (III):

(III) = ¶ n,j (c ′ )¸j (6.18) where the first inequality follows from (6.8), the second from Theorem 6.2, the third from ¶ j-1,i (c) g ¶ j-1,i (c ′ ) and the same argument used in the treatment of term (I), and the fourth one from the fact that j > i.

Finally, a similar argument partitioning i according to i f n/2 and i g n/2 as in the proof of Theorem 6.1 leads to the conclusion where

). Through a simple exercise using Theorem 2.5 and Theorem 2.4, one can show that Y i is (», 1/2)-SW, and ¶ j is (»/ √ n, 1)-SW for a constant » depending on the model only.

Then Proposition 2.4 and Theorem 2.5 further imply the following (with a potentially different value of constant »): (6.23) In order to get concrete bounds for ∥f Jr,Js ∥, we first consider the scenario J r /J s → 0 as n → ∞.

In the case J r /J s → 0 we have ∥f Jr,Js ∥ ≍J -a/2 r .

Plug this into (6.20), (6.21), (6.22), (6.23) we obtain (6.4). Now we work on the risk stability. Under our model, for any truncation value J we have

and

By definition

which is (n -3/2 », 3/2)-SW. Now apply Theorem 2.6 to the sequence ( J j=Jr+1 ∇ i ( ¶ 2 j ) : J r < J f J s ) with respect to the filtration {D n ∪ {X i1 , ..., X i,J } : J r < J f J s } to conclude that

s » , 2 -SW and hence (6.5) follows.

In the other case, J r g cJ s for a constant c > 0. Now ∥f Jr,Js ∥ ≍ J -(1+a)/2 s √ J s -J r and the remaining proof is the same as in the case of J r /J s → 0.

Proof of Proposition 6.4. Let ¸t = ∇ i ft . For t > i, we have, by construction ¸t = ¸t-1 -³ t Λ t Z t Z T t ¸t-1 . (6.24)

Let ˜t = Λ -1/2 t ¸t and Zt = Λ 1/2 t Z t (with 0/0 = 0), then (6.24) becomes ˜t = ˜t-1 -³ t Zt ZT t ˜t-1 . (6.25)

We first control ∥˜i∥ in a similar fashion as in Theorem 6.2: (6.26) because, by assumption, ∥ Zi ∥ f C for w > 1/2 and Y i -ï fi-1 , Z i ð is uniformly bounded. Taking squared norm on both sides of (6.25), we have

Dividing both sides by ∥˜t -1 ∥ 2 and take conditional expectation given D t-1 , we have

t f e -³tJ -2w t , (6.27) where the second inequality follows from that the minimum effective eigenvalue (those multiplied to the non-zero entries of ˜t-1 ) is at least J -2w t . Now using the telescoping version of (6.27) we obtain

≲i -2a exp -c ′ (t + 1) 1-a-2wÄ -(i + 1) 1-a-2wÄ .

The claim follows from the same argument as in Theorem 6.2 by considering i f t/2 and i g t/2. Remark 6.4 (Potential improvements). Inspecting the proof of Proposition 6.4, the main source of potential (and likely) looseness responsible for the counterintuitive dependence on w seems to be bounding ˜T t-1 Λ t ˜t-1 g J -2w t ∥˜t -1 ∥ 2 2 in (6.27). The diagonal entries of Λ t have very different magnitudes. On the other hand, the shrinkage mechanism in the sieve SGD update suggests that | ft,j | decays at a rate of j -2w as j increases. Thus ˜t-1,j = j w ( ft-1,jf ′ t-1,j ) should carry a scaling of j -w . So intuitively, in the product ˜T t-1 Λ t ˜t-1 , large entries in Λ t are aligned with large entries of ˜t-1 and it is overly conservative to use a lower bound based on the smallest positive eigenvalue of Λ t .

Bibliographic Notes

Stability of regularized M-estimators that are suitable for the conditions in Section 2 and Section 4 has been established in Bousquet and Elisseeff (2002) and Austern and Zhou (2024), respectively. These results are concerned with the global optimum of the M-estimation problem, without taking into account the optimization error. First-order stability of SGD under strong convexity and smooth gradients has been established in Hardt et al. (2016), which is covered as a special case of our Theorem 6.1. Our Proposition 6.1 corresponds to Lemma 3.7.3 of Hardt et al. (2016), which provides an elegant proof. The second-order stability of SGD presented in Section 6.2 was first developed in Kissel and Lei (2022).

The stability of (the standardized) difference loss function considered in Section 6.3 is adapted from Kissel and Lei (2022), which contains numerical examples of the first and second-order stability of the difference loss function. The stability of the difference risk function in part (b) of Proposition 6.3 is new and will lead to a deterministic centering CLT for the CV estimate of expected difference risks.

As an example illustrating the necessity of part 3 of Assumption 6.3, Zhang et al. (2023) considered a similar nonparametric regression setting, where the coefficients f j alternate between zero and j -2 according to +log 3 [log 2 (j)], being even or odd. In that paper, the purpose was to study the ability of cross-validation to pick the correct model from two candidates J r and J s with J r = J s -1. If we use this sequence (f j : j g 1) in the setting of Proposition 6.3 and assume J s j n 1/5 , then the proof can be adapted to show that ∇ 1 ℓ (r,s) (X 0 , D n )/Ã (r,s) n = o P (1) if and only if f Js ̸ = 0. Such a dichotomy originates from the variance of the loss difference ℓ (r,s) (X 0 , D n ). When f Js = 0, the difference between the two true models is zero and any variance in the loss difference will be from the random fluctuation in the estimated coefficient fJs . On the other hand, when f Js ̸ = 0, there is a bias term contributing to the loss difference. We do not expect Gaussian approximation to hold for the former case due to the severe degeneracy of the loss difference.

The original idea of sieve SGD first appeared in Dieuleveut and Bach (2016). The form of sieve SGD considered in this article was proposed in Zhang and Simon (2022), which establishes its near optimality in terms of error rates and storage cost (see also Quan and Lin (2024)). A version of sieve SGD implemented without shrinkage (w = 0) has been shown to achieve the optimal convergence rate in Chen and Klusowski (2024), providing evidence for potential improved stability analysis with better dependence on the shrinkage parameter. The stability property in Proposition 6.4 is adapted from Zhang and Lei (2023).

Another important general method for designing stable algorithms is resampling. In Bühlmann and Yu (2002), it has been shown that the bootstrap-aggregating (bagging) method can turn unstable algorithms such as decision trees into stable ones. Similar bagging ideas are also used in stability-based variable selection (Meinshausen and Bühlmann, 2010;Liu et al., 2010). A more systematic study of bagging and resampling methods for algorithmic stability in the context of predictive inference has been initiated in Soloff et al. (2024a) and Soloff et al. (2024b). metric, and in high-dimensional regimes where the quality gap between two models shrinks as the sample size increases; see also Zhan and Yang (2022) for extensions to multiple choices of K.

A corresponding extension of the CV model confidence sets in Section 5.1 would compute a separate confidence set for each split. Aggregating multiple confidence sets in a principled and computationally efficient way is nontrivial; see Gasparin and Ramdas (2024) for a recent attempt.

Choice of K

A practical issue in using cross-validation is the choice of K: which value of K yields the best model selection? At first glance, the question seems somewhat circular: we plan to use K-fold CV to select the best model (or, more generally, a fitting procedure), yet doing so appears to require first selecting K-itself a modeling choice. In other words, to solve the model selection problem, we would first need to solve an even more complicated model selection problem. This observation suggests that there is no single, definitive answer to the choice of K for cross-validation. The aim of this subsection is therefore to collect a few simple facts, known concerns, and practical rules of thumb about selecting K.

For model selection consistency, it is often recommended to use reversed K-fold CV (see Remark 3.1) with a diverging K, so that most of the sample is used for evaluation (Yang, 2007;Zhang and Yang, 2015;Arlot and Celisse, 2010). However, this scheme is effective only in "strong-signal" settings, where the candidate models are well separated (e.g., in the sense of Assumption 3.1). The caveats of this trainingvalidation split are clear: (i) it does not handle the important case in which the difference between models is small; (ii) the quality of models fitted on a much smaller training sample may not reflect their quality at the original sample size. This challenge is described succinctly in Zhang and Yang (2015):

Simply speaking, in terms of comparing the predictive performances of two modeling procedures, a large enough evaluation set is preferred to account for the randomness in the prediction assessment, but at the same time we must make sure that the relative performance of the two model selection procedures at the reduced sample size resembles that at the full sample size. This typically forces the training size to be not too small. Therefore, the choice of splitting ratio needs to balance the above two conflicting directions.

How should we choose K in practice? We suggest the following rule of thumb:

• Use leave-one-out CV (LOO-CV) if computationally feasible;

• Otherwise, use 10-fold CV if feasible.

• Otherwise, use 5-fold CV.

The rationale is simple: we aim to keep the training sample size as close as possible to the full sample size. In many tuning-parameter selection problems, the optimal choice scales as a monomial of the sample size with an exponent between -1 and 1. Thus the optimal choices based on the full sample size n and the training size n tr differ by a factor between 1 -1/K and K/(K -1), which plateaus quickly once K is moderately large. Arlot and Lerasle (2016) provide a more thorough theoretical investigation of the choice of K in least-squares density estimation, reaching a similar conclusion: the variance of a quantity central to model selection depends on K through the factor 1 + 4/(K -1). Hence K = 5 is a good balance between computational cost and statistical accuracy.

Meanwhile, the lack of model selection consistency for these choices of K (K = n, 10, 5) can be addressed via the model confidence set approach in Section 5, or by alternatives such as repeated-CV voting (Zhang and Yang, 2015;Zhan and Yang, 2022).

The Rashomon Effect and Occam Dilemma

The uncertainty of CV risk estimation and the associated model selection confidence set are closely related to the Rashomon effect and the Occam dilemma-two of the three major changes in modern statistical learning highlighted in Breiman's landscape paper Breiman (2001):foot_5 -the Rashomon effect refers to the existence of many models with similar, nearly optimal predictive accuracy;

-the Occam dilemma refers to the tension between model simplicity and predictive accuracy.

See Rudin (2025) for a more detailed discussion and a comprehensive literature review. Here we focus on the connection to the present article: cross-validation and stability.

The Rashomon effect. The presence of multiple nearly optimal models has been documented by many researchers over the past decade. For example, Nevo and Ritov (2017) and Fisher et al. (2019) estimate the Rashomon set-the collection of candidate models that are nearly optimal -by including all models whose empirical risk is within a small tolerance of the empirical risk minimizer, and then use this set to study variable importance. These methods require a user-specified tolerance, whose choice can substantially affect the output. The CV model confidence set developed in Section 5.1, especially the version based on the difference loss, serves a similar purpose but with a more interpretable tuning parameter-the miscoverage level

´. The Rashomon set can also be used to construct more stable predictors with comparable near-optimal accuracy. For instance, one may average the models in the Rashomon set: because each constituent is near-optimal, predictive performance is preserved, and stability often improves through averaging. This idea can be refined via weighted averages, choosing weights to optimize stability or variance, or adapting them to covariates when models perform differently across regions of the feature space.

The Occam dilemma. The potential conflict between simplicity and predictive accuracy is of practical importance, since interpretability and accuracy are both core objectives in statistical learning and model selection. When the Rashomon set contains a simple model, the decision is straightforward: choose the simple model(s) in the set-an instance of Occam's razor, and consistent with the perspective in Section 5.2. Semenova et al. (2022) and Semenova et al. (2023) develop methods to assess whether a simple model exists within the Rashomon set and to find one via a sequence of Rashomon sets induced by progressively reduced hypothesis spaces.

The Occam dilemma is more of a concern when none of the models in the Rashomon set is simple, indicating that the underlying relationship between the response and the covariates cannot be captured by the simple candidate models under consideration. In such cases, it is reasonable to prioritize predictive accuracy over model simplicity and interpretability. Even when no simple model is available, we can still aim for the selected complex model to be stable. One reason to prefer simplicity is its estimation stability, which reduces variance in both estimation and prediction. A complex model that is stable can likewise enjoy a favorable bias-variance trade-off. Random forests and other bagging methods are prominent examples. Thus, stability in a complex model can partially offset the loss of simplicity. A stable and accurate (albeit complex) model is worth preserving: it represents our best approximation to future observations and may even guide the discovery of a simpler, equally accurate model. Hence, stability offers a principled bridge between predictive accuracy and interpretability.