The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
41,021 characters
Cross-Fitting Under Nonregularity: Normality and Inference via Locality
\maketitle
\begin{abstract}
Cross-fitting is routine in much of applied research.
While conventional confidence intervals that ignore cross-fold dependence are asymptotically valid in several settings, they undercover in many applications that share a common form of nonregularity: from the classic cross-validation problem of testing whether a fitted model outperforms another, to testing for heterogeneous treatment effects with machine learning, to estimating the value of a potentially non-unique optimal treatment regime.
Exploiting a new locality condition, I show that a large class of cross-fitting estimators still satisfies a central limit theorem despite the nonregularity, but with an asymptotic variance that must be adjusted for the cross-fold correlation.
Then, I propose a method for estimating this correlation and construct new confidence intervals that attain asymptotically nominal coverage.
Finally, I show that the proposed confidence intervals attain approximately nominal coverage in a simulation study with random forests and neural networks.
\end{abstract}
\clearpage
\section{Introduction}
\label{sec.intro}
Cross-fitting is widely popular in applied research.
By using every observation both for training and for testing, it improves efficiency over sample-splitting, but doing so introduces statistical dependence across folds that can make inference challenging.
In many settings, including semiparametric estimation and debiased machine learning, this dependence is asymptotically negligible, and confidence intervals that treat folds as independent have approximately nominal coverage \citep{schickAsymptoticallyEfficientEstimation1986,hubbardStatisticalInferenceData2016b,chernozhukovDoubleDebiasedMachine2018b,neweyCrossFittingFastRemainder2018}.
In several important applications, however, a form of nonregularity keeps the dependence from vanishing, and conventional confidence intervals undercover: from the classic cross-validation problem of comparing the predictive performance of fitted models \citep{stoneAsymptoticsCrossvalidation1977,efronImprovementsCrossValidation632+1997,bayleRelativeInstabilityModel2026}, to testing for treatment effect heterogeneity \citep{chernozhukovFisherSchultzLecture2025a,wagerCommentFisherSchultzLecture2025a}, to learning the value of a potentially non-unique optimal treatment regime \citep{robinsOptimalStructuralNested2004a,luedtkeStatisticalInferenceMean2016}.
I revisit the asymptotics of cross-fitting and make two contributions.
First, exploiting a new locality argument, I show that a large class of split-sample estimators still satisfies a central limit theorem despite the nonregularity, but with an asymptotic variance that depends on the cross-fold correlation.
Second, I propose a method for estimating this correlation and constructing confidence intervals based on cross-fitting that attain asymptotically nominal coverage in many nonregular settings.
To introduce the problem I study, and previous results, consider a simple example: estimating the conditional predictive performance with 2-fold cross-validation \citep{efronEstimatingErrorRate1983,picardCrossValidationRegressionModels1984,efronImprovementsCrossValidation632+1997,blumBeatingHoldoutBounds1999,dudoitAsymptoticsCrossvalidatedRisk2005a,yangComparingLearningMethods2006,goldsteinTestingRelativePerformance2016,bayleCrossvalidationConfidenceIntervals2020,leiCrossValidationConfidence2020,austernAsymptoticsCrossvalidation2025,luoLimitsAssumptionfreeTests2025}, a routine exercise in applications such as clinical prediction \citep{jiangCalculatingConfidenceIntervals2008,pirracchioMortalityPredictionIntensive2015}.
Given an i.i.d. sample of an outcome $Y$ and covariates $X$, $(Y_i, X_i) \sim P$ for $i = 1, \dots, n$, an analyst randomly splits $\bc{1,\dots,n}$ into two folds $\mathsf{s}_1$ and $\mathsf{s}_2$ of equal size, and uses a learning algorithm, such as a random forest, to train a model ${\hat{\mu}}_k$ to predict $Y$ from $X$ on each fold $\mathsf{s}_k$, $k = 1, 2$.
The analyst then evaluates each predictor on the other fold, and averages the two estimates:
$$\wh{\mathrm{Err}} = \frac{1}{2} \bp{ \wh{\mathrm{Err}}_1 + \wh{\mathrm{Err}}_2 }, \qquad \wh{\mathrm{Err}}_k = \frac{1}{n/2} \sum_{i \notin \mathsf{s}_k} L\of{Y_i, {\hat{\mu}}_k(X_i)},$$
where $L$ is a loss function, such as $L(y, \hat{y}) = (y - \hat{y})^2$.
A natural and popular target is the out-of-sample loss of the fitted predictors \citep{efronEstimatingErrorRate1983,dudoitAsymptoticsCrossvalidatedRisk2005a},
$$\mathrm{Err} = \frac{1}{2} \bp{ \mathrm{Err}_1 + \mathrm{Err}_2 }, \qquad \mathrm{Err}_k = \E{ L\of{\tilde{Y}, {\hat{\mu}}_k(\tilde{X})} \;\middle|\; \bc{Y_i, X_i}_{i \in \mathsf{s}_k} },$$
where $(\tilde{Y}, \tilde{X}) \sim P$ independent of the sample.
Conventional inference on $\mathrm{Err}$ uses a concentration argument that renders cross-fold dependence negligible.
If $L\of{Y_i, {\hat{\mu}}_k(X_i)} = \bar{L}_n\of{Y_i, X_i} + o_p(1)$, with $\bar{L}_n\of{y, x} = \E{L\of{y, {\hat{\mu}}_k(x)}}$, then, under standard regularity conditions,
\begin{equation} \label{eq.cv_concentration}
\sqrt{n/2} \bp{\wh{\mathrm{Err}}_k - \mathrm{Err}_k} = \frac{1}{\sqrt{n/2}} \sum_{i \notin \mathsf{s}_k} \bc{\bar{L}_n\of{Y_i, X_i} - \E{\bar{L}_n\of{\tilde{Y}, \tilde{X}}}} + o_p(1),
\end{equation}
cross-fold dependence is first-order negligible, $\wh{\mathrm{Err}}$ is approximately an average of i.i.d. terms, and $\sqrt{n} (\wh{\mathrm{Err}} - \mathrm{Err})$ is asymptotically normal.
This is the core mechanism in much of the cross-validation literature with random centering \citep{bayleCrossvalidationConfidenceIntervals2020,austernAsymptoticsCrossvalidation2025,leiModernTheoryCrossValidation2025}, which often verifies concentration around $\bar{L}_n$ through algorithmic stability, and is very close to the arguments in \citet{hubbardStatisticalInferenceData2016b} and \citet{favaTrainingTestingMultiple2025}.
In practice, however, $\mathrm{Err}$ is usually not of interest on its own, but in comparison with the error of another model or a benchmark.
Let ${\hat{\nu}}_k$ denote another model trained on fold $\mathsf{s}_k$ using a different learner, such as a neural network or the sample average ${\hat{\nu}}_k(x) = (n/2)^{-1} \sum_{i \in \mathsf{s}_k} Y_i$.
In this case, the cross-validation estimator is
$$\wh{\Delta} = \frac{1}{2} \bp{ \wh{\Delta}_1 + \wh{\Delta}_2 }, \qquad \wh{\Delta}_k = \frac{1}{n/2} \sum_{i \notin \mathsf{s}_k} \wh{\delta}_k\of{Y_i, X_i},$$
$$\wh{\delta}_k(y, x) = L\of{y, {\hat{\mu}}_k(x)} - L\of{y, {\hat{\nu}}_k(x)},$$
and the estimand is
$$\Delta = \frac{1}{2} \bp{\Delta_1 + \Delta_2}, \qquad \Delta_k = \E{ \wh{\delta}_k\of{\tilde{Y}, \tilde{X}} \;\middle|\; \bc{Y_i, X_i}_{i \in \mathsf{s}_k} }.$$
In this setting, the concentration argument in \cref{eq.cv_concentration} does not immediately apply, especially under the leading null hypothesis $\Delta = 0$.
If both ${\hat{\mu}}_k$ and ${\hat{\nu}}_k$ converge to the same limit, $\wh{\delta}_k(y, x) = o_p(1)$, and hence
$\sqrt{n} ({\wh{\Delta} - \Delta}) = o_p(1)$, which is uninformative for inference.
A more useful result would require a deterministic function $\bar{\delta}_n(y,x) = \mathbb{E}[\wh{\delta}_k(y,x)]$ such that
\begin{equation} \label{eq.cv_degenerate}
\frac{\wh{\delta}_k(Y_i, X_i)}{\sigma_n} = \frac{\bar{\delta}_n(Y_i, X_i)}{\sigma_n} + o_p(1),
\end{equation}
with $\sigma_n^2 = \operatorname{Var}[\wh{\delta}_k(\tilde{Y}, \tilde{X})]$, $\sigma_n^2 \to 0$.
This requirement is counterintuitive: if for example $E[Y|X]$ is constant, the difference in performances $\wh{\delta}_k(Y_i, X_i)$ is pure sampling noise, and it would not be unnatural that $\bar{\delta}_n(y,x) = 0$.
The failure of \cref{eq.cv_degenerate} is a known form of nonregularity \citep{robinsOptimalStructuralNested2004a,laberAdaptiveConfidenceIntervals2011a}, previously addressed with procedures other than cross-fitting.
\citet{yadlowskyEvaluatingTreatmentPrioritization2025b} use a single train/holdout split.
Sequential schemes evaluate each observation with a rule trained only on earlier data \citep{luedtkeStatisticalInferenceMean2016,wagerCommentFisherSchultzLecture2025a}, and hence train and evaluate models on less data than cross-fitting does.
Repeated sample-splitting combined with rules such as the median or twice the median of $p$-values, or majority votes over intervals, controls size \citep{vandewielTestingPredictionError2009,diciccioExactTestsMultiple2020,gasparinMergingUncertaintySets2024,gasparinCombiningExchangeableValues2025,chernozhukovFisherSchultzLecture2025a}, as do confidence intervals that inflate the standard error or add noise \citep{rinaldoBootstrappingSampleSplitting2019,daiSignificanceTestsFeature2024,verdinelliDecorrelatedVariableImportance2024,bayleRelativeInstabilityModel2026} and the multivariate approach of \citet{favaTrainingTestingMultiple2025}, but all of these procedures are conservative and/or sacrifice power.
I focus instead on cross-fitting itself, and construct confidence intervals that attain asymptotically nominal coverage by exploiting a novel locality condition.
My first contribution is a central limit theorem (CLT) that accommodates this nonregularity.
My CLT applies to a broad class of split-sample estimators, covering $K$-fold and leave-one-out cross-fitting.
Unlike in previous cross-validation CLTs, the asymptotic variance in my result depends on the correlation between fold-level estimators, which turns out to be positive in simulations, so that omitting it would lead to undercoverage.
My CLT builds on the intuition that the average of dependent random variables can still be asymptotically normal as long as the dependence structure is ``local''.
In the model comparison example,
$$\frac{\wh{\Delta}}{\sigma_n} = \frac{1}{n} \bp{\sum_{i \in \mathsf{s}_1} \frac{\wh{\delta}_2\of{Y_i, X_i}}{\sigma_n} + \sum_{j \in \mathsf{s}_2} \frac{\wh{\delta}_1\of{Y_j, X_j}}{\sigma_n}},$$
and the threat to normality is the statistical dependence between $\wh{\delta}_2(Y_i, X_i)$ and $\wh{\delta}_1(Y_j, X_j)$ for $i \in \mathsf{s}_1$ and $ j \in \mathsf{s}_2$.
Dependence, however, need not invalidate normality.
Indeed, there is a long tradition of CLTs for dependent data.
In the time series literature, for example, it is typical to control dependence by assuming that two observations are ``nearly independent'' if they are ``far apart'' in time.
I formalize the notion of locality in the cross-fitting setting: changing an observation $(Y_i, X_i)$ with $i \in \mathsf{s}_1$ can affect all of $\{\wh{\delta}_1(Y_j, X_j)\}_{j \in \mathsf{s}_2}$, but this effect must be ``small'' if $(Y_i, X_i)$ is ``far'' from $(Y_j, X_j)$.
This is natural for nonparametric learners: predictions ${\hat{\mu}}_1(X_j)$ are more sensitive to a training point $(Y_i, X_i)$ when $X_j$ is close to $X_i$.
But the locality condition is not universal: for parametric learners, changing one training point can affect all test points equally.
My CLT paves the way for inference, but raises a new challenge: estimating the cross-fold correlation.
My second contribution is to address this problem.
I propose a resampling-based estimator that applies to any black-box learner and many estimators, achieving nominal coverage asymptotically under mild regularity conditions.
\Cref{sec.setup} lays out the setting and notation, and \Cref{sec.clt_2cf} derives the central limit theorem.
\Cref{sec.inference} is devoted to inference and feasible confidence intervals, and \Cref{sec.simulations} contains a Monte Carlo study supporting the theoretical results.
\Cref{sec.conclusion} concludes.
The main text focuses on two-fold cross-fitting to simplify the exposition.
In \Cref{app.clt_general}, I extend the result to $K$-fold cross-fitting, including leave-one-out.
All proofs are deferred to the appendix.
\section{Setting and Notation}
\label{sec.setup}
I define the setting and notation, introduce the cross-fitting estimator and its estimand, and present three applications that fit into this framework, including the cross-validation model comparison problem discussed in the introduction.
Consider an i.i.d. sample $\Dc = \{W_i\}_{i=1}^n$ of even size $n$, with $W_i \sim P$ taking values in $\Wc$.
For example, $W = (Y, X)$ for an outcome $Y$ and covariates $X$.
Further, consider a random split of the indices $\bs{n} = \bc{1,\dots,n}$ into two disjoint folds $\mathsf{s}_1,\mathsf{s}_2$ of equal size, independently of the data.
Without loss of generality, I treat the split as fixed.
For any subset $\mathsf{s} \subseteq \bs{n}$, let $\Dc_\mathsf{s} = \bc{W_i : i \in \mathsf{s}}$ denote the corresponding subsample and ${\tilde{\mathsf{s}}} = \bs{n} \setminus \mathsf{s}$ the complement of $\mathsf{s}$.
Let $\Ac$ denote a training algorithm, a measurable map $\Ac : \bigcup_{m \ge 1} \Wc^m \to \Hc$, invariant under permutations of its arguments, that takes a training sample $\Dc$ of any size and outputs a model $\Ac(\Dc) \in \Hc$.
Moreover, let $\bc{f_\eta : \eta \in \Hc}$ be a known family of functions $f_\eta : \Wc \to \Rb$, and let ${\hat{\eta}}_k = \Ac(\Dc_{\mathsf{s}_k})$ be the model trained on the $k$-th fold.
I study the cross-fitting estimator
$$\hat{\theta}_{\hat{\eta}} = \frac{1}{2} \sum_{k=1}^{2} \hat{\theta}_{{\hat{\eta}}_{k}}, \qquad \hat{\theta}_{{\hat{\eta}}_{k}} = \frac{1}{n/2} \sum_{i \notin \mathsf{s}_{k}} f_{{{\hat{\eta}}}_{k}}(W_i),$$
where ${\hat{\eta}} = ({\hat{\eta}}_1, {\hat{\eta}}_2)$.
The estimand is
$$\theta_{\hat{\eta}} = \frac{1}{2} \sum_{k=1}^{2} \theta_{{\hat{\eta}}_k}, \qquad \theta_{{\hat{\eta}}_k} = P f_{{\hat{\eta}}_{k}} = \int f_{{\hat{\eta}}_{k}}(w) \mathrm{d} P(w).$$
For example, if $f_\eta(w) = (y - \eta(x))^2$, $\theta_{\hat{\eta}}$ is the out-of-sample mean squared error of a predictor that selects one of the two fitted models uniformly at random.
The three examples below fit into this setting and exhibit the nonregularity that motivates this paper.
\begin{example}[Model comparison with cross-validation \citep{blumBeatingHoldoutBounds1999,dudoitAsymptoticsCrossvalidatedRisk2005a,bayleCrossvalidationConfidenceIntervals2020,austernAsymptoticsCrossvalidation2025,favaTrainingTestingMultiple2025,bayleRelativeInstabilityModel2026}] \label{ex.setup.cv}
\,
Let $W = (Y,X)$, with outcome $Y$ and covariates $X$, and let $\eta = (\mu, \nu)$ be a pair of predictors of $Y$ from $X$.
The fitted model ${\hat{\eta}}_k = ({\hat{\mu}}_k, {\hat{\nu}}_k)$ then consists of two predictors trained on fold $\mathsf{s}_k$ by different learners, for example a random forest and a neural network.
The score is the loss difference
$$f_{\eta}(y,x) = L(y, \mu(x)) - L(y, \nu(x)),$$
where $L$ is a loss function, such as the squared error $L(y,\hat{y}) = (y - \hat{y})^2$ or the classification error $L(y,\hat{y}) = \ind{y \neq \hat{y}}$ for discrete outcomes.
Each
$$\theta_{{\hat{\eta}}_k} = \int \bs{L(y, {\hat{\mu}}_k(x)) - L(y, {\hat{\nu}}_k(x))} \mathrm{d} P(y,x)$$
is the difference in out-of-sample loss between the two fitted models ${\hat{\mu}}_k$ and ${\hat{\nu}}_k$, and the estimand
$$\theta_{\hat{\eta}} = \frac{1}{2} \bp{\theta_{{\hat{\eta}}_1} + \theta_{{\hat{\eta}}_2}}$$
is its average over the two folds, or, equivalently, the out-of-sample loss difference when $({\hat{\mu}}_k, {\hat{\nu}}_k)$ is drawn uniformly at random from the two fitted pairs.
\end{example}
\begin{example}[Value of a learned treatment rule \citep{luedtkeStatisticalInferenceMean2016,shiBreakingCurseNonregularity2020,whitehouseInferenceOptimalPolicy2026}]
\label{ex.setup.otr}
\,
Let $W = (Y, A, X)$, with outcome $Y$, binary treatment assignment $A$, and covariates $X$, in a randomized experiment with known treatment probability $\P{A = 1 \mid X} = \pi \in (0,1)$.
Moreover, let $Y(0), Y(1)$ denote potential outcomes with $Y = Y(A)$, and let $\eta : \Xc \to \bc{0,1}$ be a treatment rule, so that ${\hat{\eta}}_k$ is a rule learned on fold $\mathsf{s}_k$, for example by empirical welfare maximization \citep{kitagawaWhoShouldBe2018} or by thresholding an estimate of the conditional average treatment effect $\tau(x) = \E{Y(1) - Y(0) \mid X = x}$.
The value of a rule is $V(\eta) = \E{Y(\eta(X))}$, and $V(\eta) = P f_\eta$ for
$$f_\eta(y, a, x) = \frac{a y}{\pi} \eta(x) + \frac{(1-a) y}{1 - \pi} \bc{1 - \eta(x)}.$$
Each $\theta_{{\hat{\eta}}_k} = V({\hat{\eta}}_k)$ is the value of the rule learned on fold $\mathsf{s}_k$, and the estimand $\theta_{\hat{\eta}}$ is the average value of the two learned rules.
The value of the learned treatment regimes is also informative for the value of the optimal treatment regime.
The learned rules are dominated by the optimal rule $\eta^*(x) = \ind{\tau(x) > 0}$, so $\theta_{\hat{\eta}} \le V(\eta^*)$, and a lower confidence bound for $\theta_{\hat{\eta}}$ is also a lower confidence bound for the optimal value.
Moreover, inference on $\theta_{\hat{\eta}}$ is asymptotically equivalent to inference on $V(\eta^*)$ whenever the regret $V(\eta^*) - \theta_{\hat{\eta}} = o_p(n^{-1/2})$, which \citet{luedtkeStatisticalInferenceMean2016} establish under margin and rate conditions, without requiring ${\hat{\eta}}_k$ to converge.
\end{example}
\begin{example}[Best linear predictor (BLP) test of predictive power \citep{chernozhukovFisherSchultzLecture2025a,wagerCommentFisherSchultzLecture2025a}]
\label{ex.setup.blp}
\,
Let $W = (Y, X)$, and let $\eta : \Xc \to \Rb$ be a predictor of $Y$ from $X$.
A natural way to assess the predictive power and calibration of ${\hat{\eta}}_k$ is to regress $Y$ on ${\hat{\eta}}_k(X)$ and a constant in the test fold ${\tilde{\mathsf{s}}}_k$.
A slope of zero means no predictive power, and a slope of one means good calibration.
Holding ${\hat{\eta}}_k$ fixed, the population slope is
$$\theta_{{\hat{\eta}}_k} = Pf_{{\hat{\eta}}_k}, \qquad f_\eta(y, x) = \frac{(y - \E{Y})\bc{\eta(x) - P\eta}}{P(\eta - P\eta)^2}.$$
The ordinary least squares estimator differs from $\hat{\theta}_{{\hat{\eta}}_k}$ due to the oracle scale $P({\hat{\eta}}_k - P{\hat{\eta}}_k)^2$ and the oracle centerings $\E{Y}$ and $P\eta$, but this distinction is typically first-order negligible since the sample variance of ${\hat{\eta}}_k(X)$ is ratio-consistent for $P({\hat{\eta}}_k - P{\hat{\eta}}_k)^2$ under mild moment conditions.
Under the null that $\E{Y \mid X}$ is constant, if ${\hat{\eta}}_k$ is consistent for this constant, $P({\hat{\eta}}_k - P{\hat{\eta}}_k)^2 \pto 0$, and the score diverges.
The standardized score $(f_{{\hat{\eta}}_k} - P f_{{\hat{\eta}}_k}) / \operatorname{Var}[{f_{{\hat{\eta}}_k}(W)}]^{1/2}$ may then have no deterministic limit, which is the nonregularity my CLT addresses.
The same construction tests for treatment effect heterogeneity in a randomized controlled trial by replacing $Y$ with its Horvitz--Thompson transform, which yields the best linear predictor of \citet[Theorem~3.2]{chernozhukovFisherSchultzLecture2025a}.
\end{example}
\section{A CLT for Two-Fold Cross-Fitting}
\label{sec.clt_2cf}
I establish a central limit theorem for two-fold cross-fitting that accommodates the examples of the previous section.
First, I discuss the required assumptions of the CLT, including a novel locality condition.
Then, I establish a qualitative CLT, followed by a quantitative result that establishes convergence to normality at the i.i.d. rate $n^{-1/2}$ in the Wasserstein metric under stronger conditions.
To simplify exposition, I focus on data-generating processes (DGPs) where
$$P f_{\eta}=0 \qquad \text{and} \qquad P f_{\eta}^2 = 1$$
for every $\eta \in \Hc$.
Alternatively, the general result follows from replacing each score $f_{{\hat{\eta}}_k}(w)$ by $(f_{{\hat{\eta}}_k}(w) - P f_{{\hat{\eta}}_k}) / \sigma_{{\hat{\eta}}_k}$, with $\sigma^2_{\eta} = P(f_\eta - P f_\eta)^2$.
While this normalization is infeasible, it is a natural target for a CLT, since the centering reflects the estimand of interest $\theta_{\hat{\eta}}$, and $\sigma^2_{{\hat{\eta}}_k}$ can often be estimated consistently from the test sample ${\tilde{\mathsf{s}}}_k$, which I address when deriving feasible standard errors.
Importantly, this normalization allows the variance of the unstandardized score $\sigma^2_{{\hat{\eta}}_k}$ to converge to zero in probability, requiring only that it be almost surely positive for every $n$.
This section's main result is a new CLT establishing
$$T_n = \frac{\sqrt{n} \hat{\theta}_{\hat{\eta}}}{\sqrt{1 + \rho_n}} \leadsto \Nc \of{0, 1},$$
with
$$\rho_n = \Corr{\hat{\theta}_{{\hat{\eta}}_1}, \hat{\theta}_{{\hat{\eta}}_2}} > -1.$$
I derive feasible standard errors in \Cref{sec.inference}.
I collect in \Cref{ass.clt_2cf} the conditions I use for the CLT.
I use the convention, without loss of generality, that $1,2 \in \mathsf{s}_1$.
The assumption uses three further definitions.
First, define the mean-square stability coefficient
$$\gamma_n = \E{\bc{f_{{\hat{\eta}}_1}(W) - f_{{\hat{\eta}}'_1}(W)}^2 },$$
where $W \sim P$ independent of the data, $\Dc'_{\mathsf{s}_1} = \bc{W'_1} \cup \bc{W_i : i \in \mathsf{s}_1 \setminus \bc{1}}$ is the training sample with $W_1$ replaced by an independent draw $W'_1 \sim P$, and ${\hat{\eta}}'_1 = \Ac(\Dc'_{\mathsf{s}_1})$ is the corresponding refitted model.
$\gamma_n$ is the mean-square stability coefficient studied by \citet{kaleCrossValidationMeanSquareStability2011} and \citet{bayleCrossvalidationConfidenceIntervals2020}, applied here to the centered and conditionally normalized score.
Second, define the local covariance coefficients
$$r_i = \frac{n}{2} \Cov{\sqrt{\frac{n}{2}} \hat{\theta}_{{\hat{\eta}}_1}, \sqrt{\frac{n}{2}} \hat{\theta}_{{\hat{\eta}}_2} \;\middle|\; \Dc_{-i} }, \qquad i=1,\dots,n,$$
where $\Dc_{-i} = \bc{W_j : j \neq i}$ contains all observations except $W_i$.
Each coefficient measures the covariance between the fold estimators when only observation $i$ varies, isolating the interaction between using $W_i$ to train ${\hat{\eta}}_1$ and to evaluate $\hat{\theta}_{{\hat{\eta}}_2}$, when $i \in \mathsf{s}_1$.
Note that this normalization makes $r_i$ a local analogue of the cross-fold correlation, since
$$\E{r_i} = \rho_n.$$
Moreover, for a random variable $X$, write
$$
\norm{X}_p
=
\begin{cases}
\E{\abs{X}^p}^{1/p}, & 1 \le p < \infty,\\
\inf\bc{M \ge 0 : \P{\abs{X} \le M} = 1}, & p = \infty.
\end{cases}
$$
\begin{assumption}
\label{ass.clt_2cf}
\,
For some $\bar{c} < \infty$, $\ubar{c} > 0$, $p \in (2, \infty]$, with the convention $1/p = 0$ when $p = \infty$, and sufficiently large $n$, the following conditions hold:
\begin{assumptionenum}
\item \label{ass.clt_2cf.moments} (Moments).
$$\norm{f_{{\hat{\eta}}_1}(W)}_p \le \bar{c};$$
\item \label{ass.clt_2cf.degeneracy} (Non-degeneracy).
$$\rho_n \ge -1 + \ubar{c};$$
\item \label{ass.clt_2cf.stability} (Mean-square stability).
$$n^{\frac{1}{2(1-1/p)}} \gamma_n \to 0;$$
\item \label{ass.clt_2cf.locality} (Locality).
$$\Cov{r_1, r_2} \to 0.$$
\end{assumptionenum}
\end{assumption}
\begin{lemma}[Sufficient Conditions]
\label{lem.clt_2cf.sufficient}
When $\Var{r_1} > 0$,
$$\abs{\Cov{r_1, r_2}} \le n \gamma_n \abs{\Corr{r_1, r_2}}.$$
In particular, \Cref{ass.clt_2cf.stability,ass.clt_2cf.locality} hold under any of the following:
\begin{enumerate}[label=\upshape(\alph*)]
\item
$n \gamma_n = O(1), \Corr{r_1, r_2} = o(1)$;
\item
$p = \infty, \sqrt{n} \, \gamma_n = o(1), \sqrt{n} \, \Corr{r_1, r_2} = O(1)$;
\item
$n \gamma_n = o(1)$.
\end{enumerate}
\end{lemma}
\Cref{ass.clt_2cf.moments} is a uniform $L^p$ moment bound for some $p > 2$, with $p = \infty$ for bounded scores.
\Cref{ass.clt_2cf.degeneracy} rules out near-perfect negative correlation between the fold-level estimators.
It holds, for example, whenever $\rho_n \ge 0$, as in all simulations of \Cref{sec.simulations}.
\Cref{ass.clt_2cf.stability} is a mean-square stability condition, restricting the sensitivity of the fitted score to replacing one training observation.
Its exponent $1/[2(1-1/p)]$ decreases from $1$ as $p \downarrow 2$ to $\frac{1}{2}$ at $p = \infty$, so that lighter tails permit a less stable learner.
Related conditions have been used to establish asymptotic normality of cross-validation estimators of test error under the stronger requirement $n \gamma_n = o(1)$ \citep{bayleCrossvalidationConfidenceIntervals2020,austernAsymptoticsCrossvalidation2025,leiModernTheoryCrossValidation2025}, which implies $\rho_n \to 0$ (\Cref{remark.stability}).
This condition can often be verified using algorithm-specific leave-one-out or replace-one stability bounds, available under suitable regularity and tuning conditions \citep[for a discussion, see, e.g.,][Section~6]{leiModernTheoryCrossValidation2025}.
Finally, \Cref{ass.clt_2cf.locality} is a locality condition.
It requires the dependence created by using the same data point $W_i$ to train model ${\hat{\eta}}_1$ and evaluate $f_{{\hat{\eta}}_{2}}(W_i)$ to be local across observations.
Note that $r_i$ measures the dependence between $\hat{\theta}_{{\hat{\eta}}_1}$ and $\hat{\theta}_{{\hat{\eta}}_2}$ as only observation $W_i$ varies, isolating the interaction between using $W_i$ to train ${\hat{\eta}}_1$ and to evaluate $\hat{\theta}_{{\hat{\eta}}_2}$, when $i \in \mathsf{s}_1$.
Under the stability requirement $n \gamma_n = O(1)$, it suffices that $\Corr{r_1, r_2} \to 0$ at any rate (\Cref{lem.clt_2cf.sufficient}(a)).
This condition is plausible for nonparametric learners, including machine learning algorithms, because changing a training observation mostly affects predictions near its own covariate value.
The coefficient $r_i$ is then driven by the test observations whose covariates lie close to those of $W_i$, so that $r_1$ and $r_2$ depend on nearly disjoint parts of the test fold and are approximately uncorrelated.
In parametric models, however, all observations act through the same few parameters, so the contributions $r_i$ can remain correlated, and asymptotic normality may fail.
Note that $\E{r_i} = \rho_n$, and that \Cref{ass.clt_2cf.locality} allows $\rho_n$ to be non-zero asymptotically.
\Cref{lem.clt_2cf.sufficient} gives three sufficient conditions for \Cref{ass.clt_2cf.stability,ass.clt_2cf.locality}, trading off stability against locality.
Since $\Var{r_1} \le n \gamma_n$, \Cref{ass.clt_2cf.locality} holds whenever $n \gamma_n \Corr{r_1, r_2} \to 0$.
Under condition (a), stability alone keeps $n \gamma_n$ bounded, so the correlation may vanish at any rate.
Condition (b) relaxes stability by strengthening locality: once the correlation vanishes at rate $n^{-1/2}$, stability is needed only at the $\sqrt{n}$ rate for bounded scores, $\sqrt{n} \gamma_n \to 0$.
Condition (c) is the stability requirement of previous cross-validation CLTs, and it implies \Cref{ass.clt_2cf.locality} with no restriction on the correlation.
Hence, my CLT nests the regime $n \gamma_n \to 0$, in which $\rho_n \to 0$.
\begin{theorem}[Cross-Fitting Central Limit Theorem]
\label{thm.clt_2cf}
Let \Cref{ass.clt_2cf} hold.
Then,
$$T_n \leadsto \Nc \of{0, 1}.$$
\end{theorem}
\Cref{thm.clt_2cf} is, to my knowledge, the first to establish a central limit theorem for cross-fitting estimators that accommodates the nonregularity that the standardized influence functions may not converge to a fixed limit.
Compared to previous CLTs in the cross-validation \citep{dudoitAsymptoticsCrossvalidatedRisk2005a,bayleCrossvalidationConfidenceIntervals2020,austernAsymptoticsCrossvalidation2025,leiModernTheoryCrossValidation2025} and semiparametric estimation/double machine learning \citep{schickAsymptoticallyEfficientEstimation1986,chernozhukovDoubleDebiasedMachine2018b} literatures, the asymptotic variance of $\sqrt{n} \hat{\theta}_{\hat{\eta}}$ in my result depends on the cross-fold correlation, which may be nonzero.
Note that the stronger stability condition $n \gamma_n = o(1)$, required by previous cross-validation CLTs, implies $\rho_n = o(1)$.
Compared to those results, I weaken the stability requirement (\Cref{ass.clt_2cf.stability}), and introduce the locality \Cref{ass.clt_2cf.locality}.
The next result strengthens \Cref{thm.clt_2cf}.
Under stronger conditions, the normal approximation holds at the rate $n^{-1/2}$, the Berry--Esseen rate for averages of independent random variables with finite third moments.
I state the result in the Wasserstein distance, defined for random variables $X$ and $Y$ by
$$d_W\of{X, Y} = \sup_{h \in \Lip{1}} \, \bigl| \E{h(X)} - \E{h(Y)} \bigr|,$$
where
$$\Lip{1} = \bc{h \colon \Rb \to \Rb \;:\; \abs{h(x) - h(y)} \le \abs{x - y} \text{ for all } x, y \in \Rb}.$$
\begin{theorem}[Quantitative Cross-Fitting Central Limit Theorem]
\label{thm.clt_2cf.quantitative}
Let \Cref{ass.clt_2cf} hold with $p \ge 3$.
Suppose, in addition, that
$$
n^{\frac{1}{1-1/p}}\gamma_n=O(1),
\qquad
n \Cov{r_1, r_2}=O(1).
$$
Then,
$$
d_W\of{T_n, \Nc\of{0,1}}=O(n^{-1/2}).
$$
\end{theorem}
\Cref{thm.clt_2cf.quantitative} shows that the dependence between folds does not slow the normal approximation as long as the locality condition holds at a fast enough rate.
Its conditions trade off tails against stability.
The moment condition requires a uniformly bounded $L^p$ norm of the score with $p \ge 3$, as in the Berry--Esseen theorem.
The stability condition requires $\gamma_n = O(n^{-1/(1-1/p)})$, which ranges from $\gamma_n = O(n^{-3/2})$ at $p = 3$ to $\gamma_n = O(n^{-1})$ for bounded scores, so that heavier tails require a more stable learner.
For bounded scores, the stability condition has the same order as the sufficient condition $n \gamma_n = O(1)$ in \Cref{lem.clt_2cf.sufficient}(a), and the $O(n^{-1/2})$ rate then costs only a strengthening of locality, from $\Cov{r_1, r_2} = o(1)$ to $O(n^{-1})$.
Since $\Cov{r_1, r_2} = \Var{r_1} \Corr{r_1, r_2}$ and $\Var{r_1} \le n \gamma_n$ (\Cref{lem.clt_2cf.stability_bounds}), the latter holds whenever $\Corr{r_1, r_2} = O(n^{-1})$.
\begin{remark}[Comparison with previous stability conditions]
\label{remark.stability}
Previous stability-based CLTs for cross-validation require first-order loss stability at the $o(n^{-1})$ scale relative to the variance of the leading score
\citep{bayleCrossvalidationConfidenceIntervals2020,austernAsymptoticsCrossvalidation2025,leiModernTheoryCrossValidation2025}.
Applied to the centered and conditionally normalized score used here, this condition is equivalent to
\[
n\gamma_n\to0.
\]
Indeed, if $\widetilde\sigma_n^2=\operatorname{Var}[\mathbb{E}[f_{{\hat{\eta}}_1}(W)\mid W]]$, their condition is
$n\gamma_n/\widetilde\sigma_n^2\to0$, while
\[
0\le 1-\widetilde\sigma_n^2\le \frac n4\gamma_n.
\]
Hence, previous results cover the regime $n\gamma_n=o(1)$, which also forces $\rho_n\to0$.
By contrast, \Cref{thm.clt_2cf} allows weaker stability, for example $n\gamma_n=O(1)$ when the locality condition makes the covariance contributions asymptotically uncorrelated.
\end{remark}
\section{Estimation of \texorpdfstring{$\rho$}{rho} and Inference}
\label{sec.inference}
I propose a resampling-based procedure to construct feasible confidence intervals based on \Cref{thm.clt_2cf}, and show that they attain nominal coverage asymptotically under regularity conditions.
Let $n=4m$ be the total sample size, and return to the unnormalized scores of \Cref{sec.setup}, that is, $f_{{\hat{\eta}}_k}(w)$ may have non-zero mean and non-unit variance.
Throughout this section, $\rho_n$ denotes the cross-fold correlation of the centered and normalized fold estimators,
$$\rho_n = \Corr{\frac{\hat{\theta}_{{\hat{\eta}}_1} - \theta_{{\hat{\eta}}_1}}{\sigma_{{\hat{\eta}}_1}}, \frac{\hat{\theta}_{{\hat{\eta}}_2} - \theta_{{\hat{\eta}}_2}}{\sigma_{{\hat{\eta}}_2}}},$$
which is the correlation of \Cref{sec.clt_2cf} computed with each score $f_{{\hat{\eta}}_k}$ replaced by $(f_{{\hat{\eta}}_k} - \theta_{{\hat{\eta}}_k}) / \sigma_{{\hat{\eta}}_k}$.
I consider directly observable scores, which fit \Cref{ex.setup.cv,ex.setup.otr}, and give an implementation for the BLP \Cref{ex.setup.blp} in \Cref{app.inference.blp}.
Keep the original folds fixed.
For $b=1,\dots,B$ and $k=1,2$, independently draw a uniform subset $\mathsf{s}_{k,b}\subset\mathsf{s}_k$ of size $m$ and fit
$$
{\hat{\eta}}_{k,b}=\Ac(\Dc_{\mathsf{s}_{k,b}}).
$$
Define
\begin{align}
\hat{\sigma}_{{\hat{\eta}}_{k,b}}^2
&=
\frac1{2m-1}\sum_{i\in{\tilde{\mathsf{s}}}_k}
\bs{
f_{{\hat{\eta}}_{k,b}}(W_i)
-
\frac1{2m}\sum_{j\in{\tilde{\mathsf{s}}}_k}f_{{\hat{\eta}}_{k,b}}(W_j)
}^2,
\label{eq.inference.pooled_scale}\\
Z_{k,b}
&=
\frac{
\sum_{i\in\mathsf{s}_{3-k,b}}f_{{\hat{\eta}}_{k,b}}(W_i)
-
\sum_{i\in{\tilde{\mathsf{s}}}_k\setminus\mathsf{s}_{3-k,b}}f_{{\hat{\eta}}_{k,b}}(W_i)
}{
\sqrt m\,\hat{\sigma}_{{\hat{\eta}}_{k,b}}
},
\label{eq.inference.contrast}\\
\wh\rho_n
&=
\max\bc{
-1,\,
\min\bc{
1,\,
\frac1B\sum_{b=1}^B Z_{1,b}Z_{2,b}
}
}.
\label{eq.inference.rhohat}
\end{align}
Set $Z_{k,b}=0$ if its denominator is zero.
Note that $Z_{k,b}$ evaluates the same fitted score on two disjoint samples, so its expected value is zero.
Let $\hat{\sigma}_{{\hat{\eta}}_k}^2$ be the test-fold sample variance
\begin{equation}
\label{eq.inference.original_scale}
\hat{\sigma}_{{\hat{\eta}}_k}^2
= \frac{1}{2m-1}\sum_{i\in{\tilde{\mathsf{s}}}_k}
\bs{f_{{\hat{\eta}}_k}(W_i)-\hat{\theta}_{{\hat{\eta}}_k}}^2,
\qquad k=1,2.
\end{equation}
The proposed standard error and confidence interval are
\begin{align}
\wh{\se}\of{\hat{\theta}_{\hat{\eta}}}
&=
\sqrt{
\frac{
\hat{\sigma}_{{\hat{\eta}}_1}^2+\hat{\sigma}_{{\hat{\eta}}_2}^2
+2\wh\rho_n\hat{\sigma}_{{\hat{\eta}}_1}\hat{\sigma}_{{\hat{\eta}}_2}
}{2n}
},
\label{eq.inference.se}\\
\widehat{\mathrm{CI}}_{\alpha}
&=
\bs{
\hat{\theta}_{\hat{\eta}}
\pm z_{1-\alpha/2}\wh{\se}\of{\hat{\theta}_{\hat{\eta}}}
},
\label{eq.inference.ci}
\end{align}
where $z_{1-\alpha/2}$ is the standard normal quantile.
The original estimator and its target are unchanged.
\begin{theorem}[Asymptotic Validity of Confidence Interval]
\label{thm.inference}
Let \Cref{ass.clt_2cf} hold after centering and normalizing the scores.
Suppose $0 < \sigma_\eta^2 < \infty$ for every $\eta\in\Hc$.
If, for some $\rho > - 1$,
\begin{equation}
\Var{\E{Z_{1,1}Z_{2,1}\mid\Dc}}\to0,
\qquad
\rho_{n} \to \rho,
\qquad
B\to\infty,
\label{eq.inference.resampling}
\end{equation}
then $\wh\rho_n-\rho_n\Lto[2]0$.
If, additionally,
\begin{equation}
\frac{\sigma_{{\hat{\eta}}_1}}{\sigma_{{\hat{\eta}}_2}}\pto1,
\label{eq.inference.scale_balance}
\end{equation}
then, for every fixed $\alpha\in(0,1)$,
$$
\P{\theta_{\hat{\eta}}\in\widehat{\mathrm{CI}}_{\alpha}}\to1-\alpha.
$$
\end{theorem}
The first condition in \cref{eq.inference.resampling} requires the covariance averaged over resampling splits to stabilize across datasets.
The second requires the cross-fold correlation to converge to some limit greater than $-1$.
\section{Simulation Study}
\label{sec.simulations}
I examine the finite-sample performance of the confidence intervals derived in \Cref{sec.inference} in a simulation study.
In particular, I consider the cross-validation (CV) model comparison (\Cref{ex.setup.cv}) and the BLP test of predictive power (\Cref{ex.setup.blp}).
I use sample sizes $n \in \bc{1000, 2000}$ and 2,400 Monte Carlo replications for each data-generating process and sample size.
The correlation estimator uses $B=500$ repetitions.
Observations are independent and identically distributed, with $X_i \sim \Nc(0, I_5)$ and $U_i \sim \Nc(0, 1)$ independent of $X_i$.
I consider two no-signal designs.
In the homoskedastic Gaussian design,
\[
Y_i=U_i.
\]
In the heteroskedastic, conditionally skewed design,
\[
Y_i
=
\omega(X_i)
\frac{
\exp{\lambda(X_i)U_i}
-\exp{\lambda(X_i)^2/2}
}{
\sqrt{
\exp{2\lambda(X_i)^2}
-\exp{\lambda(X_i)^2}
}
},
\]
where, writing $\Lambda(t)=\bc{1+\exp{-t}}^{-1}$,
\[
\omega(x)=0.5+1.5\Lambda(x_1+0.5x_2),
\qquad
\lambda(x)=0.25+0.75\Lambda(x_2-x_3).
\]
Both designs satisfy $\E{Y_i \mid X_i}=0$.
The conditional variance is one in the Gaussian design and $\omega(X_i)^2$ in the skewed design.
The latter also permits conditional skewness to vary with the covariates.
I estimate the predictive model with two learners.
The first is a single-hidden-layer neural network from the \texttt{R} package \texttt{nnet}, with 10 hidden units, a linear output, no skip connections, and weight decay equal to $0.1$.
The second is a random forest from the \texttt{R} package \texttt{ranger}, with 200 trees, three candidate covariates per split, and a minimum node size of five.
For CV, I compare the squared-error loss of the fitted learner with that of the training fold outcome mean.
In the notation of \Cref{ex.setup.cv}, $L(y, \hat{y}) = (y - \hat{y})^2$, ${\hat{\mu}}_k$ is the learner trained on $\mathsf{s}_k$, and ${\hat{\nu}}_k$ is the constant predictor equal to the average of $Y_i$ over $i \in \mathsf{s}_k$.
For both applications, I report the coverage rate $\P{\theta_{\hat{\eta}} \in \widehat{\mathrm{CI}}_{\alpha}}$ at the nominal $95\%$ level.
Since $\E{Y_i \mid X_i} = 0$ in both designs, the BLP estimand is $\theta_{\hat{\eta}}=0$.
For CV, the estimand of \Cref{ex.setup.cv} is
$$\theta_{\hat{\eta}} = \frac{1}{2} \sum_{k=1}^{2} \int \bs{\bc{y - {\hat{\mu}}_k(x)}^2 - \bc{y - {\hat{\nu}}_k}^2} \mathrm{d} P(y, x),$$
the average out-of-sample loss difference, holding the fitted models fixed.
This target need not be zero under a no-signal design.
I compute it by integrating analytically over the outcome noise and numerically over 65,000 independent covariate draws per fitted model.
\begin{table}[tbp]
\centering
\small
\caption{Coverage rates of nominal 95\% BLP and CV confidence intervals.}
\label{tab.simulations.coverage}
\begin{tabular}{@{}lrcccc@{}}
\toprule
& & \multicolumn{2}{c}{BLP} & \multicolumn{2}{c}{CV} \\
\cmidrule(lr){3-4}
\cmidrule(l){5-6}
Design & $n$ & Estimated $\wh\rho_n$ & $\wh\rho_n = 0$ & Estimated $\wh\rho_n$ & $\wh\rho_n = 0$ \\
\midrule
\multicolumn{6}{c}{Neural network} \\
\addlinespace
Gaussian & 1000 & 94.38 & 90.42 & 94.21 & 91.67 \\
& 2000 & 95.04 & 91.13 & 94.75 & 91.54 \\
Skewed & 1000 & 94.33 & 91.54 & 94.38 & 92.75 \\
& 2000 & 95.13 & 92.38 & 94.42 & 92.13 \\
\midrule
\multicolumn{6}{c}{Random forest} \\
\addlinespace
Gaussian & 1000 & 94.13 & 84.50 & 94.54 & 85.50 \\
& 2000 & 95.67 & 86.88 & 95.54 & 87.54 \\
Skewed & 1000 & 94.75 & 87.96 & 95.29 & 89.00 \\
& 2000 & 95.33 & 89.33 & 95.79 & 90.38 \\
\bottomrule
\end{tabular}
\medskip
\begin{minipage}{\linewidth}
\footnotesize
\textit{Notes:}
Entries are coverage percentages of $\theta_{{\hat{\eta}}}$,
based on 2,400 replications per design and sample size.
Gaussian denotes the homoskedastic Gaussian design,
Skewed denotes the heteroskedastic, conditionally skewed design.
Estimated $\wh\rho_n$ uses \cref{eq.inference.rhohat} with $B=500$;
$\wh\rho_n = 0$ ignores the cross-fold correlation.
\end{minipage}
\end{table}
\Cref{tab.simulations.coverage} reports the results.
With the estimated correlation, coverage rates range from $94.13\%$ to $95.67\%$ for BLP and from $94.21\%$ to $95.79\%$ for CV, close to the nominal $95\%$ level across the learners, designs, and sample sizes considered.
Setting $\wh\rho_n = 0$, which ignores the cross-fold correlation, yields undercoverage in every case, with coverage rates from $84.50\%$ to $92.38\%$ for BLP and from $85.50\%$ to $92.75\%$ for CV.
The undercoverage is most severe for the random forest, whose coverage falls below $90\%$ in the Gaussian design for both applications and sample sizes.
\section{Conclusion} \label{sec.conclusion}
I revisited the problem of inference with cross-fitting when the standardized fitted score need not converge to a deterministic limit, a nonregularity shared by cross-validation model comparison, tests of predictive power and treatment effect heterogeneity, and the value of learned treatment rules.
Exploiting a new locality condition, I established a central limit theorem for cross-fitting estimators with an asymptotic variance that depends on the cross-fold correlation, and I proposed confidence intervals based on a resampling estimator of this correlation that attain asymptotically nominal coverage.
In simulations with random forests and neural networks, these intervals delivered coverage close to the nominal level, whereas intervals that ignore the cross-fold correlation undercovered.
\clearpage
\bibliography{refs2.bib}
\clearpage