arXiv CSopen access

Accelerating Conformal Prediction via Approximate Leave-One-Out

Accelerating Conformal Prediction via Approximate Leave-One-Out Jiachen Cong, Jingbo Liu∗1 1

Department of Statistics, University of Illinois Urbana-Champaign, Champaign, IL, 61820, USA

arXiv:2606.31915v1 [stat.ML] 30 Jun 2026

July 1, 2026

Abstract While conformal prediction provides a general framework for uncertainty quantification in predictive inference, its application is often limited by computational cost. Recent methods, including Jackknife+ and Jackknife-minmax, achieve faster computation by trading a slight loss of efficiency relative to full conformal prediction, but still requires computing leave-one-out refits for all observations. In this paper, we further accelerate conformal prediction by incorporating approximate leave-one-out (ALO) estimators, and establish asymptotic coverage and efficiency. While our proof draws on methods developed for analyzing the consistency of ALO cross-validation risk estimators in high-dimensional statistics, it requires adaptations to handle conformal prediction, where leave-i-out residuals are needed for predictions at xn+1 rather than just at the training covariate xi . Simulation results validate our theoretical findings, showing that the ALO-based methods achieve coverage and efficiency comparable to the exact methods, while significantly reducing the runtime.

1

Introduction

Quantifying uncertainty is a central problem in statistics, often addressed by constructing prediction intervals for a new response yn+1 given training data D = {(xi , yi )}ni=1 and xn+1 . Conformal prediction methods provide distribution-free guarantees, with Full Conformal Prediction (FCP) Vovk et al. [2005] being a canonical approach. However, FCP is computationally expensive due to its requirement of searching over candidate response values. In contrast, Jackknife+ Barber et al. [2021] directly output prediction intervals and are more suitable for continuous-response settings, but still incur a computational cost scaling linearly with n. In this work, we focus on Jackknife+ and Jackknife-minmax as practical alternatives to FCP. To address this limitation, recent work has explored ALO techniques. For example, Clarté and Zdeborová [2024] proposed using AMP-based estimator to accelerate FCP; Bhatt et al. [2021], Abad Martinez et al. [2023], Prinster et al. [2023] proposed influence function-based approach for constructing ALO estimators. Another line of work Karoui et al. [2013], Karoui [2018] proves ALO error bounds and use them to show central limit theorems of the estimator, but the latter requires stronger assumptions. However, existing approaches lack rigorous error control for ALO in this setting, which may lead to suboptimal efficiency in practice. Moreover, Clarté and Zdeborová [2024] assumes i.i.d. covariates with diagonal covariance, whereas we consider a more general Gaussian model N (0, Σ) with Σ (not necessarily diagonal). In this paper, we develop accelerated versions of Jackknife+ and Jackknife-minmax by incorporating ALO estimators derived via Newton updates and the Woodbury identity. Our approach significantly reduces computational cost while preserving the statistical guarantees of the original methods. Contribution. ∗

Emails: [email protected], [email protected]. This research was supported in part by NSF Grant DMS-2515510.

1

• We establish uniform error bounds for ALO-based predictions with a new covariate xn+1 , showing that the approximation error ⊤ x⊤ n+1 θ̃/i − xn+1 θ̂/i

(1.1)

remains uniformly controlled, whereas prior work [Wang et al., 2018, Rad and Maleki, 2020, Auddy et al., 2023] focused only on approximation for xi , i ∈ {1, 2, . . . , n}. Our proof leverages some decomposition techniques to decompose (1.1) to terms can be controlled, thereby reducing the new-covariate case to the setting studied in [Rad and Maleki, 2020]. Compared with [Clarté and Zdeborová, 2024], we weaken the assumption on the covariance structure of the distribution of xi ; in particular, we no longer require the covariance matrix to be diagonal. • We propose the accelerated Jackknife+ and Jackknife-minmax methods for constructing prediction intervals based on ALO estimator θ̃/i , and prove they asymptotically retain the same coverage probability and efficiency (interval length) as their original counterparts. • In a special high-dimensional linear model, we establish that Full Conformal, Split Conformal, Jackknife+, and Jackknife-minmax prediction intervals are asymptotically equivalent. Related work. A substantial body of work has examined the construction and error bounds of ALO estimators. [Beirami et al., 2017, Wang et al., 2018, Rad et al., 2020, Rad and Maleki, 2020, Xu et al., 2021, Auddy et al., 2023, Zou et al., 2024] develop explicit forms of ALO estimators using the Woodbury identity and analyze their theoretical and empirical performance in estimating out-of-sample risk. In particular, [Rad and Maleki, 2020, Wang et al., 2018, Zou et al., 2024] investigate ALO estimators under smoothed regularizers, covering ℓ1 penalties, non-smooth but piecewise twice-differentiable regularizers, and more general non-smooth convex Lipschitz regularizers. Another line of work [Karoui et al., 2013, Karoui, 2018] prove ALO error bounds and use them to show central limit theorems of the estimator, but the latter requires stronger assumptions. Notations. Consider a model parametrized by θ ∈ Rp trained on data D. Let r(θ) = (1 − η)r0 (θ) + ηθ⊤ θ 1×p stand for the i-th row of where η ∈ (0, 1) and r0 is a smooth non-negative convex regularizer. Let x⊤ i ∈R ⊤ n×p ⊤ (n−1)×1 X := (x1 , x2 , ..., xn ) ∈ R ; and we define y := (y1 , y2 , ..., yn ) . y/i ∈ R and X/i ∈ R(n−1)×p stand for y and X, excluding the ith entry yi and the ith row x⊤ i , respectively. We assume that   n X  θ̂ := arg minp ℓ(yj , x⊤ θ) + λr(θ) (1.2) j  θ∈R  j=1

is the empirical risk minimizer (ERM), where the loss function ℓ(yj , x⊤ j θ) and regularizer r(θ) are twice⊤ differentiable in θ, and ℓ(yj , xj θ) is convex in θ. Similarly, we define the leave-i-out ERM estimator by θ̂/i := arg minp θ∈R

Define ℓ̇i (θ) :=

 X 

j̸=i

  ⊤ ℓ(yj , xj θ) + λr(θ) . 

(1.3)

∂ℓ(yi , z) ∂ 2 ℓ(yi , z) , ℓ̈i (θ) := . ∂z ∂z 2 z=x⊤ θ z=x⊤ θ i

i

Similarly, we define r̈(θ) = ∇2θ r(θ). Define f (x⊤ θ) as the prediction model f : R → R, with parameter θ. Polynomials of log(n) are denoted + {d } as the ⌈(1 − α)(n + 1)⌉-th smallest value of d , . . . , d , and q̂ − {d } as by polylog(n). We define q̂n,α i 1 n n,α i the ⌊α(n + 1)⌋-th smallest value of d1 , . . . , dn . We use o(1) to denote a deterministic sequence that vanishes as n → ∞, and oP (1) to denote a sequence of random variables that converges to zero in probability.

2

2

Preliminaries

2.1

Full conformal prediction

Let Dn = {(xi , yi )}ni=1 be the training data, and let xn+1 denote a test covariate for which the response yn+1 is unobserved. Full conformal prediction (FCP) outputs a prediction set ĈαFCP (xn+1 ) satisfying  P yn+1 ∈ ĈαFCP (xn+1 ) ≥ 1 − α, (2.4) for any user-specified miscoverage level α ∈ (0, 1), under the exchangeability of (x1 , y1 ), . . . , (xn , yn ), (xn+1 , yn+1 ). To construct the prediction set, FCP considers each candidate value y ∈ R as a hypothetical response for the test point and refits the model on the augmented dataset Dn(y) = Dn ∪ {(xn+1 , y)}. (y)

Let θ̂(y) denote the estimator obtained from Dn . Define the corresponding nonconformity scores by σi (y) := f (x⊤ i θ̂(y)) − yi , and

i = 1, . . . , n,

σn+1 (y) := f (x⊤ n+1 θ̂(y)) − y .

Let σ(1) (y) ≤ · · · ≤ σ(n+1) (y) denote the order statistics of {σi (y)}n+1 i=1 . Then the full conformal prediction set is given by  Ĉ(xn+1 ) = y ∈ R : σn+1 (y) ≤ σ(⌈(1−α)(n+1)⌉) (y) . (2.5) Equivalently, FCP includes a candidate value y if the residual of the hypothetical test point conforms with the residuals computed from the augmented sample. This method is distribution-free and enjoys finite-sample coverage, but it is computationally expensive because θ̂(y) must be recomputed for each candidate value y, or for each value on a fine grid in practice; see [Vovk et al., 2005, Angelopoulos and Bates, 2022].

2.2

Jackknife+

Specifically, we consider two variants of FCP in our paper, called Jackknife+ and Jackknife-minmax Barber et al. [2021]. The construction of prediction intervals is introduced in (2.6) and (2.7), respectively. Let RiLOO := |yi − f (x⊤ i θ̂/i )|, and define h i   jackknife+ − LOO + ⊤ LOO Ĉn,α := q̂n,α f (x⊤ θ̂ ) − R , q̂ f (x θ̂ ) + R , (2.6) n+1 /i i n,α n+1 /i i jack-mm Ĉn,α :=

h

 LOO  LOO i + ⊤ + min f (x⊤ θ̂ ) − q̂ R , max f (x θ̂ ) + q̂ , /i /i n+1 n,α i n+1 n,α Ri

i=1,...,n

i=1,...,n

(2.7)

In Barber et al. [2021], the authors established following results:

2.3

jackknife+ P{yn+1 ∈ Ĉn,α } ≥ 1 − 2α,

(2.8)

jack-mm P{yn+1 ∈ Ĉn,α } ≥ 1 − α.

(2.9)

Leave-one-out approximation

Following Rad and Maleki [2020], using the Newton step and the Woodbury identity, the approximate leave-iout estimator θ̃/i ∈ Rp×1 is defined as follows θ̃/i := θ̂ + where J :=

P

n ⊤ j=1 xj xj ℓ̈j (θ̂) + λ diag

h

J−1 xi ℓ̇i (θ̂) −1 1 − x⊤ i J xi ℓ̈i (θ̂)

i r̈(θ̂) . 3

,

(2.10)

2.4

Assumptions in [Rad and Maleki, 2020] and [Zou et al., 2024]

Assumption 1. (a) r0 is a non-negative, convex and twice differentiable function. (b) Loss function ℓ(y, z) is non-negative, convex and continuously differentiable with respect to z. Remark 1. In our main analysis, we assume that r0 is twice differentiable so that the Hessian is well-defined. For certain non-smooth regularizers, such as the ℓ1 penalty, smoothing approximations can be employed as in [Zou et al., 2024, Rad and Maleki, 2020] to construct ALO estimators with smoothed regularizers. While extending our results to this setting is possible, it lies beyond the main scope of the present paper. Assumption 2. (a) X = (x1 , · · · , xn )⊤ where xi ∈ Rp are i.i.d. N (0, Σ) samples. Let ρmax denote the largest eigenvalue of Σ. Moreover, there exist constants 0 < cX ≤ CX such that p−1 cX ≤ σmin (Σ) ≤ σmax (Σ) ≤ p−1 CX . (b) We assume observations are independent and identically distributed draws from some unknown joint distribution ∗ ∗ p p(yi |x⊤ i θ )q(xi ), where θ ∈ R represents the true parameter. Remark 2. We consider the regime where n, p → ∞, n/p = δ0 ∈ (0, ∞), and elements of θ∗ are O(1); we have √ ∥θ∗ ∥ = O( p), and hence E[(x⊤ θ∗ )2 ] = O(1) [Zou et al., 2024]. Assumption 3. (a) We define ℓ̇(·) := (ℓ̇1 (·), . . . , ℓ̇n (·))⊤ ; ℓ̇/i (·) := (ℓ̇1 (·), . . . , ℓ̇i−1 (·), ℓ̇i+1 (·), . . . , ℓ̇n (·))⊤ ; ℓ̈/i (·) := (ℓ̈1 (·), . . . , ℓ̈i−1 (·), ℓ̈i+1 (·), . . . , ℓ̈n (·))⊤ . We assume that c1 (n) = O(polylog(n)) and c2 (n) = O(polylog(n)), and qn → 0 are all functions of n, such that with probability at least 1 − qn for all i = 1, . . . , n c1 (n) > ∥ℓ̇(θ̂)∥∞ ; c2 (n) > sup

∥ℓ̈/i {(1 − t)θ̂/i + tθ̂} − ℓ̈/i (θ̂)∥2 ∥θ̂/i − θ̂∥2

t∈[0,1]

c2 (n) > sup t∈[0,1]

∥r̈{(1 − t)θ̂/i + tθ̂} − r̈(θ̂)∥2 ∥θ̂/i − θ̂∥2

;

.

(b) There exists a constant ν > 0 and a sequence q̃n → 0 such that for all i = 1, ..., n !     ⊤ inf σmin λ diag r̈(tθ̂ + (1 − t)θ̂/i ) + X/i diag ℓ̈/i (tθ̂ + (1 − t)θ̂/i ) X/i

≥ ν

t∈[0,1]

with probability at least 1 − q̃n . Remark 3. For part (a) of the assumption, [Rad and Maleki, 2020, Section 4] provided its justification for the ridge, elastic net, logistic regression, robust regression, and Poisson regression by choosing c1 (n) and c2 (n) to be polynomials of n. For part (b) of the assumption, since we chose r(θ) = (1 − η)r0 (θ) + ηθ⊤ θ, we have ν = 2λη.

2.5

Approximation error bound in [Rad and Maleki, 2020]

Under Assumption 1 to Assumption 3, [Rad and Maleki, 2020] proved that !  ℓ̇i (θ̂) Hii C0 ⊤ ⊤ max xi θ̂/i − xi θ̂ − ≤√ 1≤i≤n 1 − H p ii ℓ̈i (θ̂) 8n holds with probability at least 1 − 4ne−p − 8n − (n−1) 3 − qn − q̃n , where C0 is defined as same as (3.18). p3 This result however is not sufficient for our application to conformal prediction, and we will adapt the proof technique to bound the approximation error of x⊤ n+1 θ̂/i .

4

2.6

Lemma 25 in [Rad and Maleki, 2020]

Lemma 1 (Lemma 25 in [Rad and Maleki, 2020]). Assume that X ⊤ (D + Γ)X and X ⊤ DX are positive definite, and define: Γ := diag(γ),

ω̄max := σmax (XX ⊤ ),

νmin := σmin (X ⊤ (D + Γ)X),

A := X ⊤ DX.

Then, z ⊤ X ⊤ (D + Γ)X

3

−1

z − z ⊤ X ⊤ DX

−1

z

≤



∥γ∥2 +

 ω̄max 2 ∥γ∥24 XA−1 z 4 . νmin

(2.11)

Main results

Assumption 4. Under assumptions 1 to 3, there exists a large enough constant L < ∞ such that P{y ∈ [a, b]|x⊤ θ} ≤ L · |b − a| and |f (b) − f (a)| ≤ L · |b − a| for ∀a, b ∈ R. Compared with Assumptions 1-3, the newly added Assumption 4 is used to control the effect of the deviation of approximate intervals from the exact ones on the resulting coverage guarantee. For simplicity, under Assumption 3, we define   16n 16n −p p1 (n) := + 3 + 8ne + qn + q̃n , (3.12) (n − 1)3 p   8n 8n −p + 3 + 4ne + qn + q̃n , (3.13) p2 (n) := (n − 1)3 p and when n/p = δ0 ∈ (0, ∞) it is clear that n→∞

n→∞

p1 (n) −−−−−→ 0,

3.1

(3.14)

p2 (n) −−−−−→ 0.

Algorithm

Motivated by (2.10), we introduce R̃iLOO := |yi −f (x⊤ i θ̃/i )|, i = 1, 2, . . . , n. Then the jackknife+ prediction interval based on estimators in (2.10) is constructed as in (3.15), h i   jackknife+ − LOO + ⊤ LOO C̃n,α := q̂n,α f (x⊤ θ̃ ) − R̃ , q̂ f (x θ̃ ) + R̃ . (3.15) /i /i n+1 i n,α n+1 i and the jackknife-minmax prediction interval based on estimators in (2.10) is constructed as in (3.16) h  LOO  LOO i jack-mm + ⊤ + C̃n,α := min f (x⊤ θ̃ ) − q̂ R̃ , max f (x θ̃ ) + q̂ . (3.16) n+1 /i n,α i n+1 /i n,α R̃i i=1,...,n

jackknife+

For simplicity, let C̃n,α jack-mm Ĉn,α = [Ljkm , Ujkm ].

i=1,...,n

jackknife+

= [L̃jk+ , Ũjk+ ], Ĉn,α

jack-mm

= [Ljk+ , Ujk+ ], C̃n,α

= [L̃jkm , Ũjkm ], and

Algorithm 1 Accelerating Jackknife+ with ALO Estimators Input: (x1 , y1 ), (x2 , y2 ), . . . , (xn , yn ), xn+1 , α Output: Prediction interval for yn+1 1. Calculate the full-sample estimator θ̂ (see Eq. (1.2)). 2. Calculate the approximate leave-i-out estimator θ̃/i (see Eq. (2.10)) for i ∈ {1, 2, ..., n}. 3. Obtain the approximate i-th LO residual R̃iLOO := |yi − f (x⊤ i θ̃/i )| and the approximate ⊤ leave-i-out estimation f (xn+1 θ̃/i ) with xn+1 for i ∈ {1, 2, ..., n}. 4. Obtain the prediction interval for yn+1 , defined as [L̃jk+ , Ũjk+ ].

5

Algorithm 2 Accelerating Jackknife-minmax with ALO Estimators Input: (x1 , y1 ), (x2 , y2 ), . . . , (xn , yn ), xn+1 , α Output: Prediction interval for yn+1 1. Calculate the full-sample estimator θ̂ (see Eq. (1.2)). 2. Calculate the approximate leave-i-out estimator θ̃/i (see Eq. (2.10)) for i ∈ {1, 2, ..., n}. 3. Obtain the approximate i-th LO residual R̃iLOO := |yi − f (x⊤ i θ̃/i )| and the approximate ⊤ leave-i-out estimation f (xn+1 θ̃/i ) with xn+1 for i ∈ {1, 2, ..., n}. 4. Obtain the prediction interval for yn+1 , defined as [L̃jkm , Ũjkm ].

3.2

Approximation guarantees

Theorem 3.1. Let Assumption 1 to Assumption 3 hold with ρmax = c/p and n/p = δ0 ∈ (0, ∞). Moreover, suppose that n is large enough such that qn + q̃n < 0.5. Then with probability at least 1 − p1 (n), where p1 (n) is defined in (3.12), the following bound is valid: C0 ⊤ max |x⊤ n+1 θ̃/i − xn+1 θ̂/i | ≤ √ , p

1≤i≤n

(3.17)

where C0 :=

216c3/2 ν3

!

! √  1/2 + c3/2 ( δ + 3)2 p p 5 c c log n 0 c21 (n)c2 (n) + c31 (n)c22 (n) 1 + δ0 ( δ0 + 3)2 . log p ν2 (3.18)

To prove Theorem 3.2, we introduce the following lemma. Lemma 2. Let {ai }ni=1 be a sorted sequence such that a1 ≤ a2 ≤ · · · ≤ an , and suppose there exist sequences {bi }ni=1 satisfying |bi − ai | ≤ ε, for some ε > 0 and all i = 1, . . . , n. Then for any q ∈ {1, 2, ..., n} we have |b(q) − aq | ≤ ϵ, where b(q) is the q-th smallest value of {bi }ni=1 . Proof of the lemma 2 is deferred to Section S.5 in the appendix. Theorem 3.2. Under Assumption 1 to Assumption 4, when n, p → ∞, n/p = δ0 ∈ (0, ∞), with probability at least (1 − p1 (n) − p2 (n)), |Ũjk+ − Ujk+ | = o(1), |L̃jk+ − Ljk+ | = o(1); |Ũjkm − Ujkm | = o(1), |L̃jkm − Ljkm | = o(1). Theorem 3.3. Under Assumption 1 to Assumption 4, when n, p → ∞, n/p = δ0 ∈ (0, ∞), we have jackknif e+ P{yn+1 ∈ C̃n,α } ≥ 1 − 2α − o(1).

(3.19)

Theorem 3.4. Under Assumption 1 to Assumption 4, when n, p → ∞, n/p = δ0 ∈ (0, ∞), we have jack−mm P{yn+1 ∈ C̃n,α } ≥ 1 − α − o(1).

6

(3.20)

4

Numerical experiments

4.1

Synthetic data

In this section, we conduct simulations to support the results in Theorem 3.3-3.4. We report the mean of coverage, mean of operation time, mean of interval length and mean of Jaccard index of original methods (labeled as “JK+” and “JK-minmax”, where “JK” stands for "Jackknife") and accelerated methods (labeled as “Fast JK+” and “Fast JK-minmax”). Recall that the Jaccard index between two sets S1 , S2 is defined as J (S1 , S2 ) =

|S1 ∩ S2 | ∈ [0, 1]. |S1 ∪ S2 |

Values closer to 1 indicate more precise approximations. In Table 3 we report the Jaccard index between Jackknife+ and Fast Jackknife+, Jackknife-minmax and Fast Jackknife-minmax. In our simulation, we set α = 0.1, which means our target coverage level is 1 − α = 0.9 (in the jackknife+ case, as shown in [Barber et al., 2021], the coverage level is 1 − 2α = 0.8). We use training sample size n = 100, test sample size ntest = 100, and repeat the experiment at each dimension p = 50, 100, 200, with √ i.i.d. data points (xi , yi ) generated as xi ∼ N (0, Ip / p) and yi |xi ∼ N (x⊤ i θ, 1). The true coefficient vector θ is randomly generated from a standard normal distribution . For the Ridge regression model, we define the loss function as ℓ(y, x⊤ θ) = (y − x⊤ θ)2 /2 and the regularization term as r(θ) = 12 r0 (β) + 21 θ⊤ θ, where Pseudo-Huber regularizer (this setting allows our results to be compared with those reported in q simulation Pp βj2 [Clarté and Zdeborová, 2024]) r0 (β) = j=1 4( 1 + 4 − 1), with the ridge parameter λ = 0.1, 1, separately. To obtain stable results, we repeat the procedures above for 50 iterations and report the averaged outcomes. Table 1: Comparison of JK+ and Fast JK+. Parameters

Model

Coverage Time(s) Length

(n=100, p=50, λ=1) (n=100, p=100, λ=1) (n=100, p=200, λ=1) (n=100, p=50, λ=0.1) (n=100, p=100, λ=0.1) (n=100, p=200, λ=0.1)

JK+ JK+ JK+ JK+ JK+ JK+

0.880 0.892 0.897 0.872 0.888 0.889

0.075 0.109 0.185 0.091 0.172 0.354

3.808 4.189 4.409 4.135 4.729 4.662

(n=100, p=50, λ=1) (n=100, p=100, λ=1) (n=100, p=200, λ=1) (n=100, p=50, λ=0.1) (n=100, p=100, λ=0.1) (n=100, p=200, λ=0.1)

Fast JK+ Fast JK+ Fast JK+ Fast JK+ Fast JK+ Fast JK+

0.880 0.893 0.897 0.872 0.888 0.890

0.004 0.008 0.043 0.005 0.009 0.044

3.809 4.190 4.410 4.135 4.734 4.672

Table 2: Comparison of JK-minmax and Fast JK-minmax. Parameters

Model

Coverage Time(s) Length

(n=100, p=50, λ=1) (n=100, p=100, λ=1) (n=100, p=200, λ=1) (n=100, p=50, λ=0.1) (n=100, p=100, λ=0.1) (n=100, p=200, λ=0.1)

JK-minmax JK-minmax JK-minmax JK-minmax JK-minmax JK-minmax

0.911 0.920 0.918 0.926 0.942 0.937

0.086 0.103 0.189 0.094 0.164 0.388

4.187 4.549 4.718 4.925 5.574 5.387

(n=100, p=50, λ=1) (n=100, p=100, λ=1) (n=100, p=200, λ=1) (n=100, p=50, λ=0.1) (n=100, p=100, λ=0.1) (n=100, p=200, λ=0.1)

Fast JK-minmax Fast JK-minmax Fast JK-minmax Fast JK-minmax Fast JK-minmax Fast JK-minmax

0.911 0.920 0.918 0.926 0.943 0.937

0.004 0.008 0.042 0.004 0.008 0.047

4.188 4.550 4.719 4.925 5.750 5.400

7

Table 3: Prediction-interval Overlap (Jaccard Index) of Fast JK+ & JK+ and Fast JK-minmax & JK-minmax. Parameters (n=100, p=50, λ=1) (n=100, p=100, λ=1) (n=100, p=200, λ=1) (n=100, p=50, λ=0.1) (n=100, p=100, λ=0.1) (n=100, p=200, λ=0.1)

Fast JK+ & JK+

Fast JK-minmax & JK-minmax

0.9997 0.9997 0.9998 0.9996 0.9983 0.9978

0.9997 0.9996 0.9997 0.9996 0.9982 0.9976

Table 1 and Table 2 present the mean of coverage, mean of operation time and mean of interval length of Jackknife+, Fast Jackknife+, Jackknife-minmax and Fast Jackknife-minmax, respectively. Table 3 shows the similarity, measured by Jaccard Index, between the prediction intervals constructed by the accelerated methods and their corresponding original methods. We find that our accelerated methods substantially reduce the average computational time while maintaining coverage, and without significantly altering the length (efficiency) of the prediction intervals. The prediction intervals constructed by the accelerated methods exhibit a high degree of similarity to those from the original methods. We can observe that, in most cases, a smaller value of λ tends to result in a wider interval, and the interval length also increases as p grows. Moreover, our methods exhibit more significant acceleration when the dimensionality of the covariates is higher. Compared with the synthetic-data simulation results reported in Table 1 of Clarté and Zdeborová [2024], both of our acceleration methods achieve higher coverage than Taylor-AMP, while also providing more efficient prediction (in terms of shorter average interval length) compared with Taylor-AMP, SCP [Vovk et al., 2005] and CQP [Romano et al., 2019]. With respect to Table 2 of Clarté and Zdeborová [2024], fast Jackknife+ and fast Jackknife-minmax exhibit higher Jaccard similarity to exact LOO than Taylor-AMP and SCP. Under Gaussian settings (Table 3 of Clarté and Zdeborová [2024]), our methods deliver more efficient prediction than Bayes posterior and FCP combined with Taylor-AMP. In comparison with Table 4 of Clarté and Zdeborová [2024], both of our accelerated procedures achieve higher prediction efficiency and higher coverage than Taylor-AMP and approximate homotopy [Ndiaye and Takeuchi, 2019], while exhibiting comparable computation time. Furthermore, our framework systematically explores multiple dimensions (p = 50, 100, 200), corresponding to different n/p ratios, and two regularization strengths (λ = 0.1, 1), demonstrating consistent performance across regimes. These improvements can be attributed to two main factors. First, unlike Clarté and Zdeborová [2024], which assumes i.i.d. features with a diagonal covariance matrix, our framework accommodates a general covariance structure Σ. This relaxation is more realistic in practice and ensures that the resulting ALO characterization remains accurate even when the features are correlated, thereby preventing the deterioration observed for AMP and Taylor-AMP under non-isotropic designs. Second, while AMP-based methods rely on an ALO heuristic whose approximation error is not rigorously controlled, our approach benefits from a explicit and tighter ALO error bound derived via the Newton step and Woodbury identity. This leads to more accurate leave-one-out predictions, which in turn improves the quality of the prediction intervals and enhances efficiency.

4.2

Application to Real Data

We evaluate the proposed methods on the Concrete Compressive Strength Dataset [Yeh, 1998] and the Energy Efficiency Dataset [Tsanas and Xifara, 2012]. For both datasets, we consider a linear prediction model f (x⊤ θ) = x⊤ θ, where the full sample estimator is obtained by minimizing the regularized empirical risk objective  s    p n  2 1 X  2 X θ 1 j 2    θ̂ = arg minp yi − x⊤ θ + λ 2 1 + − 1 + ∥θ∥ , i 2  θ∈R  2 4 2 j=1

i=1

8

with the regularization parameter fixed at λ = 0.1. For each dataset, the first 80% of the observations are used for training, and the remaining 20% are used for testing. As shown in Table 4 and Table 5, the accelerated methods achieve nearly identical empirical coverage and comparable average interval lengths to the original Jackknife+ and Jackknife-minmax methods, while substantially reducing the computational time. Table 4: Performance on the Concrete Compressive Strength Dataset Method Jackknife+ Fast Jackknife+ Jackknife-minmax Fast Jackknife-minmax

Coverage Rate

Operation Time (s)

Average Interval Length

0.9660 0.9660 0.9709 0.9709

1.7275 0.0651 1.5735 0.0414

37.0917 37.0851 37.8552 37.8671

Table 5: Performance on the Energy Efficiency Dataset Method Jackknife+ Fast Jackknife+ Jackknife-minmax Fast Jackknife-minmax

4.3

Coverage Rate

Operation Time (s)

Average Interval Length

0.9351 0.9351 0.9481 0.9481

1.6567 0.0509 1.0640 0.0371

12.0404 12.0403 12.2182 12.2181

Code

The code used to produce the results in Table 1, Table 2, Table 3, Table 4 and Table 5 can be found in the following github repository: https://github.com/JiachenCong/Accelerating_Conformal_Prediction_ via_Approximate_Leave-One-Out.git. All experiments were run on an Apple M4 laptop with 16 GB of memory.

5

Discussion

In this section, we address two natural questions regarding our algorithm’s asymptotic performance in high-dimensional settings Karoui [2018]: Q1: We applied the leave-one-out approximation formula to accelerate jackknife+ and jackknife-mm, rather than full conformal prediction, so that rigorous results in [Rad and Maleki, 2020, Auddy et al., 2023, Zou et al., 2024] can be used to control the approximation error. Does this lead to worse efficiency (longer intervals)? Q2: [Karoui et al., 2013] used leave-one-out to prove the central limit theorem M-estimators. Is it possible to further accelerate our algorithm by directly using the central limit theorem result? In this section, we address these questions by studying the asymptotic performances of full conformal prediction, split conformal prediction, jackknife+, and jackknife-mm in a high-dimensional setting rigorously studied in [Karoui, 2018], although some of the conclusions are expected to hold in broader settings. Following [Karoui, 2018], we consider the model yi = Xi⊤ β0 + ϵi , and n  τ 1X  β̂ = arg minp ρi yi − Xi⊤ β + ∥β∥2 . β∈R n 2 i=1

9

(5.21)

Similarly, we define β̂/i as the estimator calculated without the i-th sample. Note that, following [Karoui, 2018], √ we assume the scaling ∥Xi ∥2 = Θ( n) and ∥β∥2 = Θ(1), which is consistent with the scaling in the preceding √ sections if we take Xi := nxi and β := √1n θ. Within the general framework in the preceding sections, (5.21) can be understood as the special case where the regularizer is ridge, the log-likelihood ℓ(y, z) = ρ(y − z) for some ρ, and yi = Xi⊤ β + ϵi under p(yi |Xi⊤ β), where ϵ1 , . . . , ϵn are i.i.d. noise. The main contribution of [Karoui, 2018] was to use leave-one-out (both observation and predictor) to prove the asymptotic distribution of β̂. Though more general sufficient conditions can be found in [Karoui, 2018], for simplicity let us assume the following condition (see [Karoui, 2018, Section 2.1]): C1 p/n has a finite nonzero limit. C2 ρ is convex, whose first and second derivatives are uniformly bounded: ∥ρ′ ∥∞ < ∞, ∥ρ′′ ∥∞ < ∞. Furthermore, sign(ρ′ (t)) = sign(t) and ρ(t) ≥ ρ(0) = 0 for all t ∈ R. C3 X1 , . . . , Xn are i.i.d. following N (0, Ip ), or with i.i.d. entries with bounded support, symmetric density, and unit variance. C4 ϵ1 , . . . , ϵn are i.i.d. following a differentiable, symmetric, unimodal distribution with variance σϵ2 , whose density f satisfies limt→∞ tf (t) = 0. C5 ∥β0 ∥2 remains bounded as n → ∞, and ∥β0 ∥∞ = O(n−ϵ ) for some ϵ > 1/4. We denote by proxρ the proximal map: 1 proxρ (x) := arg min{ρ(y) + (x − y)2 }. 2 y∈R

(5.22)

Under the conditions above, the general result in Theorem 2.1 in [Karoui, 2018] implies the following: Theorem 5.5. [Karoui et al., 2013] Let δ := limn→∞ n/p ∈ (0, ∞), and b0 := limn→∞ ∥β0 ∥22 . Let c, r > 0 be the solution to the following equations:   E (prox′cρ (Z) = 1 − δ −1 + τ, (5.23)   −1 2 2 2 −2 2 δ E (Z − proxcρ (Z)) + τ c b0 = δ r (5.24) where we defined the scalar ϵ ∼ N (0, 1) and Z ∼ N (0, r2 ) + ϵ. Then we have limp→∞ ∥β̂ − β0 ∥ = r and limp→∞ var(∥β̂ − β0 ∥2 ) = 0 in probability. Denote by UJ and Umm the upper limits of the conformal prediction intervals for the (n + 1)-th observation produced by the jackknife+ and jackknife-minmax algorithms. For any given M > 0, let Uf (M ) be the upper 1 2 limit of the discretized interval Ĉ(Xn+1 ) ∩ [−M, M ] ∩ {0, ± M ,M , . . . }, where Ĉ(Xn+1 ) is the full conformal interval. Similarly, define Lf (M ), LJ , Lmm as the lower limits. Under our assumptions, the differences of these intervals vanish as n, M → ∞: p ⊤ β̂ + Theorem 5.6. Define U := Xn+1 σϵ2 + r2 zα/2 , where zα/2 denotes the (1 − α/2) quantile of a standard normal random variable. Let an := max{|U − UJ |, |U − Umm |} and bn,M := maxi=1,...,n |Uf − U |. Then in the sense of weak convergence: lim an = 0

(5.25)

lim lim bn,M = 0

(5.26)

n→∞ M →∞ n→∞

⊤ β̂ − The same holds for the lower bounds, with L := Xn+1

10

p σϵ2 + r2 zα/2 .

The proof of Theorem 5.6 can be found in the appendix. This theorem demonstrates that under the idealized assumptions in Karoui [2018], all considered conformal prediction algorithms are asymptotically equivalent, partially answering Q1. Regarding Q2, we note that the asymptotic approximation in Theorem 5.6 relies on the CLT result, which requires independence of X1 , . . . , Xn and is much stronger than leave-observation-out approximation. Conversely, jackknife+ (and our proposed ALO version) remains valid in the general exchangeable case, which partially answers Q2. It is an interesting direction for future research to identify settings in which leave-one-out approximations remain asymptotically valid while asymptotic normality fails. A related work Liu [2024] shows that leave-one-feature-out approximations remain asymptotically accurate for general correlated sub-Gaussian feature vectors, and the techniques developed there may be extendable to the leave-one-observation-out setting.

11

References Javier Abad Martinez, Umang Bhatt, Adrian Weller, and Giovanni Cherubin. Approximating full conformal prediction at scale via influence functions. Proceedings of the AAAI Conference on Artificial Intelligence, 37 (6):6631–6639, Jun. 2023. doi:10.1609/aaai.v37i6.25814. URL https://ojs.aaai.org/index.php/AAAI/ article/view/25814. Anastasios N. Angelopoulos and Stephen Bates. A gentle introduction to conformal prediction and distributionfree uncertainty quantification, 2022. URL https://arxiv.org/abs/2107.07511. Arnab Auddy, Haolin Zou, Kamiar Rahnama Rad, and Arian Maleki. Approximate leave-one-out cross validation for regression with ℓ1 regularizers (extended version), 2023. URL https://arxiv.org/abs/2310.17629. Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. Predictive inference with the jackknife+. The Annals of Statistics, 49(1):486–507, 2021. https://doi.org/10.1214/20-AOS1965. Ahmad Beirami, Meisam Razaviyayn, Shahin Shahrampour, and Vahid Tarokh. On optimal generalizability in parametric learning. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 3458–3468, Red Hook, NY, USA, 2017. Curran Associates Inc. ISBN 9781510860964. Umang Bhatt, Adrian Weller, and Giovanni Cherubin. Fast conformal classification using influence functions. In Lars Carlsson, Zhiyuan Luo, Giovanni Cherubin, and Khuong An Nguyen, editors, Proceedings of the Tenth Symposium on Conformal and Probabilistic Prediction and Applications, volume 152 of Proceedings of Machine Learning Research, pages 303–305. PMLR, 08–10 Sep 2021. URL https://proceedings.mlr. press/v152/bhatt21a.html. Lucas Clarté and Lenka Zdeborová. Building conformal prediction intervals with approximate message passing, 2024. URL https://arxiv.org/abs/2410.16493. Noureddine El Karoui. On the impact of predictor geometry on the performance of high-dimensional ridgeregularized generalized robust regression estimators. Probability Theory and Related Fields, 170(1-2):95–175, 2018. doi:10.1007/s00440-016-0732-0. Noureddine El Karoui, Derek Bean, Peter J. Bickel, Chinghway Lim, and Bin Yu. On robust regression with high-dimensional predictors. Proceedings of the National Academy of Sciences, 110(36):14557–14562, 2013. doi:10.1073/pnas.1307842110. Jingbo Liu. Stability of a generalized debiased lasso with applications to resampling-based variable selection. arXiv preprint arXiv:2405.03063, 2024. Eugene Ndiaye and Ichiro Takeuchi. Computing full conformal prediction set with approximate homotopy, 2019. URL https://arxiv.org/abs/1909.09365. Drew Prinster, Suchi Saria, and Anqi Liu. Efficient approximate predictive inference under feedback covariate shift with influence functions. In International Symposium on Conformal and Probabilistic Prediction with Applications, 2023. URL https://api.semanticscholar.org/CorpusID:262070652. Kamiar Rahnama Rad and Arian Maleki. A scalable estimate of the out-of-sample prediction error via approximate leave-one-out cross-validation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82(4):965–996, 2020. Kamiar Rahnama Rad, Wenda Zhou, and Arian Maleki. Error bounds in estimating the out-of-sample prediction error using leave-one-out cross validation in high-dimensions. In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 4067–4077. PMLR, 26–28 Aug 2020. URL https://proceedings.mlr.press/v108/rad20a.html. 12

Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Athanasios Tsanas and Angeliki Xifara. Energy Efficiency. UCI Machine Learning Repository, 2012. DOI: https://doi.org/10.24432/C51307. Vladimir Vovk, Alex Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. 01 2005. doi:10.1007/b106715. Shuaiwen Wang, Wenda Zhou, Arian Maleki, Haihao Lu, and Vahab Mirrokni. Approximate leave-one-out for high-dimensional non-differentiable learning problems, 2018. URL https://arxiv.org/abs/1810.02716. Ji Xu, Arian Maleki, Kamiar Rahnama Rad, and Daniel Hsu. Consistent risk estimation in moderately high-dimensional linear regression. IEEE Transactions on Information Theory, 67(9):5997–6030, 2021. doi:10.1109/TIT.2021.3095375. I-Cheng Yeh. Concrete Compressive Strength. https://doi.org/10.24432/C5PK67.

UCI Machine Learning Repository, 1998.

DOI:

Haolin Zou, Arnab Auddy, Kamiar Rahnama Rad, and Arian Maleki. Theoretical analysis of leave-one-out cross validation for non-differentiable penalties under high-dimensional settings, 2024. https://arxiv.org/abs/2402.08543.

13

Appendix In this appendix, we provide detailed proofs of the theorems presented in the main paper.

Organization Table 6 presents the organization of the appendix. Table 6: Organization of Supplementary Materials

S.1

Section

Purpose

Appendix S.1 Appendix S.2 Appendix S.3 Appendix S.4 Appendix S.5 Appendix S.6

Presents existing results adapted from [Rad and Maleki, 2020] for our proof. Provides the complete proof of Theorem 3.1. Provides the complete proof of Theorem 3.2. Provides the complete proof of Theorem 3.3 and Theorem 3.4. Provides the complete proof of Lemma 2. Provides the complete proof of Theorem 5.6 and relative lemmas.

Existing Results

In this section, we present existing results used in our proof, which are adapted from [Rad and Maleki, 2020].

S.1.1

[Rad and Maleki, 2020, Eq.136]

⊤ ℓ̇ (θ). The leave-one-out estimate, θ̂ = θ̂ + ∆∗ , satisfies g (∆∗ + θ̂) = 0. Define g/i (θ) = λṙ(θ) + X/i /i /i /i /i /i The multivariate mean-value Theorem yields Z 1  ∗ ∗ 0 = g/i (θ̂ + ∆/i ) = g/i (θ̂) + J /i (θ̂ + t∆/i )dt ∆∗/i , (S.1.27) 0

where the Jacobian is Moreover, β̂ satisfies

⊤ J /i (θ) = λ diag[r̈(θ)] + X/i diag[ℓ̈/i (θ)] X/i .

(S.1.28)

0 = λṙ(θ̂) + X ⊤ ℓ̇(θ̂) = g/i (θ̂) + ℓ̇i (θ̂)xi .

(S.1.29)

We obtain

Z 1 ℓ̇i (θ̂)xi = − 0

so that ∆∗/i = −ℓ̇i (θ̂)

Z 1 0

 J /i (θ̂ + t∆∗/i )dt ∆∗/i , −1

J /i (θ̂ + t∆∗/i )dt

leading to the following inequality ∥∆∗/i ∥2 ≤

|ℓ̇i (θ̂)| ν

xi ,

(S.1.30)

(S.1.31)

! ∥xi ∥2 .

(S.1.32)

By Assumption 3.(b), we obtain the Eq.136 in [Rad and Maleki, 2020] ∥∆∗/i ∥2 ≤

c1 (n) ∥xi ∥2 ν

14

(S.1.33)

S.1.2

[Rad and Maleki, 2020, Lemma 27]

Lemma 3 (Lemma 27 in [Rad and Maleki, 2020]). Let x ∼ N (0, Σ) with ρmax := σmax (Σ), where Σ ∈ Rp×p , then i h p 2 P ∥x∥24 > 2(1 + c)ρmax p log p ≤ c . (S.1.34) p Moreover, if ωmax := σmax (XX ⊤ ), (S.1.35)

νmin := σmin (J), where x is independent of the symmetric matrix J ∈ Rp×p and X ∈ Rm×p , then     2 ρmax p −1 2 p log p < c , P ∥J x∥4 > 2(1 + c) 2 p νmin     p ρmax ωmax 2 P ∥XJ −1 x∥24 > 2(1 + c) m log m < c . 2 m νmin

(S.1.36) (S.1.37)

The proof of this lemma can be found on pages 60 and 61 of [Rad and Maleki, 2020].

S.1.3

[Rad and Maleki, 2020, Lemma 11]

Lemma 4 (Lemma 11 in [Rad and Maleki, 2020]). Let x ∼ N (0, Σ) with ρmax ≜ σmax (Σ), where Σ ∈ Rp×p , then   P ∥x∥22 > 5pρmax ≤ e−p . (S.1.38) Furthermore, if X ∈ Rn×p is composed of independently distributed N (0, n1 ) entries, then r  q nt2 p σmax (X ⊤ X) ≥ 1 + + t ≤ e− 2 . P n

S.1.4

(S.1.39)

[Rad and Maleki, 2020, Lemma 12]

Lemma 5 (Lemma 12 in [Rad and Maleki, 2020]). If X ∈ Rn×p is composed of independently distributed N (0, Σ) rows, with ρmax := σmax (Σ), where Σ ∈ Rp×p , then h i √ √ P σmax (XX ⊤ ) ≥ ( n + 3 p)2 ρmax ≤ e−p . The proof of this lemma can be found on pages 34 and 35 of [Rad and Maleki, 2020].

S.1.5

[Rad and Maleki, 2020, Theorem 3]

Theorem S.1.7 (Theorem 3 in [Rad and Maleki, 2020]). Let n/p = δ0 and [Rad and Maleki, 2020, Assumption 5] hold with ρmax = c/p. Moreover, suppose that [Rad and Maleki, 2020, Assumption 6 and 7] are satisfied, and that n is large enough such that qn + q̃n < 0.5. Then with probability at least 1 − 4ne−p −

8n 8n − − qn − q̃n , p3 (n − 1)3

(S.1.40)

the following bound is valid: max

1≤i≤n

⊤ x⊤ i θ̂/i − xi θ̂ −

ℓ̇i (θ̂) ℓ̈i (θ̂)

!

Hii 1 − Hii



Q0 ≤√ , p

(S.1.41)

where Q0 :=

72c3/2 ν3

! √  p p 5 c1/2 + c3/2 ( δ0 + 3)2 2 c log n 2 3 2 1 + δ0 ( δ0 + 3) c1 (n)c2 (n) + c1 (n)c2 (n) . log p ν2

!

(S.1.42) 15

S.2

Proof of Theorem 3.1

We introduce the following definitions:     ⊤ diag ℓ̈ (θ̂ ) X , where X = (x , ..., x ⊤ • (1) ∆∗/i := θ̂/i −θ̂, and J /i := λ diag r̈s (θ̂/i ) +X/i 1 i−1 , xi+1 , ..., xn ) . /i /i /i /i     2 2 log p −1 −1 √ • (2) Events Ẽi := 3C̄i X/i J /i (θ̂/i )xi + J /i (θ̂/i )xi >C p , 4 4     2 2 log p −1 −1 > C √p . and Ẽn+1,i := 3C̄i X/i J /i (θ̂/i )xn+1 + J /i (θ̂/i )xn+1 4

4

• (3) C̄i := 4∥xi ∥2

 2

c1 (n)c2 (n) ν

  2 (n) (1 + ω )∥x ∥ . 1 + 2c1 (n)c max i 2 2 ν

 hR i −1 1 ⊤ diag ∗ • (4) Define Φ1 := J −1 ( θ̂ ) − J ( θ̂ ) + X̄ γ ( θ̂ )dt X̄ and /i /i /i /i /i 0 −(1−t)∆/i /i /i /i  h i −1 ⊤ ∗ Φ2 := J −1 , where for δ, θ ∈ Rp /i (θ̂/i ) − J /i (θ̂/i ) + X̄/i diag γ−∆/i /i (θ̂/i ) X̄/i 

 ℓ̈/i (θ + δ) − ℓ̈/i (θ)
γδ/i (θ) := , λ r̈(θ + δ) − r̈(θ)



 X/i X̄/i := ∈ R(n−1+p)×p , I

  ℓ̈/i (θ) D/i (θ) := diag ∈ R(n−1+p)×(n−1+p) . λr̈(θ) ⊤ D (θ)X̄ . Define J /i (θ) := X̄/i /i /i

Consider ⊤ x⊤ n+1 θ̃/i − xn+1 θ̂/i

= ℓ̇i (θ̂)x⊤ n+1 M xi 1 1 1 ⊤ = ℓ̇i (θ̂)(xi + xn+1 )⊤ M (xi + xn+1 ) − ℓ̇i (θ̂)x⊤ n+1 M xn+1 − ℓ̇i (θ̂)xi M xi 2 2 2 = A1 − A2 − A3 ,

(S.2.43) (S.2.44) (S.2.45)

where we defined Z 1 M := 0

−1

J /i (θ̂/i − (1 − t)∆∗/i )dt

∗ − J −1 /i (θ̂/i − ∆/i ).

(S.2.46)

Clearly we have ⊤ |x⊤ n+1 θ̃/i − xn+1 θ̂/i | ≤ |A1 | + |A2 | + |A3 |.

By Assumption 3.(a), we have

(S.2.47)

ℓ̇i (θ̂) = O(polylog(n)).

(S.2.48)

  ⊤ J /i (θ + δ) = J /i (θ) + X̄/i diag γδ/i (θ) X̄/i .

(S.2.49)

By the definition of γδ/i (θ), we have

16

For the new A1 term, we have (xi + xn+1 )⊤ M (xi + xn+1 ) # "Z −1 1 ∗ − J −1 = (xi + xn+1 )⊤ J /i (θ̂/i − (1 − t)∆∗/i )dt /i (θ̂/i − ∆/i ) (xi + xn+1 ) 0

= (xi + xn+1 )

⊤

"Z 1 0

1

−1

J /i (θ̂/i − (1 − t)∆∗/i )dt

# − J −1 /i (θ̂/i )

(xi + xn+1 )

i h −1 ∗ ) + (xi + xn+1 )⊤ J −1 ( θ̂ ) − J ( θ̂ − ∆ /i /i /i (xi + xn+1 ) /i /i " Z 

≤ (xi + xn+1 )

⊤

J −1 /i (θ̂/i ) −

⊤ diag J /i (θ̂/i ) + X̄/i

1

0

 −1 # γ−(1−t)∆∗/i /i (θ̂/i )dt X̄/i (xi + xn+1 )

 h i  −1  −1 ⊤ ⊤ (xi + xn+1 ) + (xi + xn+1 ) J /i (θ̂/i ) − J /i (θ̂/i ) + X̄/i diag γ−∆∗/i /i (θ̂/i ) X̄/i = (xi + xn+1 )⊤ Φ1 (xi + xn+1 ) + (xi + xn+1 )⊤ Φ2 (xi + xn+1 ) ,

(S.2.50)

where ≤1 is due to (S.2.49). For the (xi + xn+1 )⊤ Φ2 (xi + xn+1 ) term, using [Rad and Maleki, 2020, Lemma 25] we have (xi + xn+1 )⊤ Φ2 (xi + xn+1 )     2 2 ω̄max,i ∗ ∗ ≤ γ−∆/i /i (θ̂/i ) + ( θ̂ )(x + x ) γ−∆/i /i (θ̂/i ) X̄/i J −1 i n+1 /i /i ν 2 4 2     2 2 ω̄max,i −1 ≤ γ−∆∗/i /i (θ̂/i ) + γ−∆∗/i /i (θ̂/i ) X̄/i J −1 ( θ̂ )x + X̄ J ( θ̂ )x /i i /i /i /i n+1 /i ν 2 2 4 4     2 2 ω̄max,i γ−∆∗/i /i (θ̂/i ) X̄/i J −1 ≤ 2 γ−∆∗/i /i (θ̂/i ) + /i (θ̂/i )xi 4 ν 2 2     2 2 ω̄max,i + 2 γ−∆∗/i /i (θ̂/i ) + γ−∆∗/i /i (θ̂/i ) X̄/i J −1 (θ̂/i )xn+1 . (S.2.51) /i ν 2 2 4 Similarly, for the term (xi + xn+1 )⊤ Φ1 (xi + xn+1 ) , we can prove the following bound (xi + xn+1 )⊤ Φ1 (xi + xn+1 ) !  Z 1  Z 1 2 2 ω̄max,i γ−(1−t)∆∗/i /i (θ̂/i ) dt X̄/i J −1 ( θ̂ )x ≤2 γ−(1−t)∆∗/i /i (θ̂/i ) dt + i /i /i ν 4 0 0 2 2 !   Z 1 Z 1 2 2 ω̄max,i +2 γ−(1−t)∆∗/i /i (θ̂/i ) dt + γ−(1−t)∆∗/i /i (θ̂/i ) dt X̄/i J −1 (θ̂/i )xn+1 . /i ν 4 0 0 2 2 (S.2.52) Under Assumption 4, we have γ−∆∗/i /i (θ̂/i )

2

≤ ℓ̈/i (θ̂/i − ∆∗/i ) − ℓ̈/i (θ̂/i ) 2 + λ r̈s (θ̂/i − ∆∗/i ) − r̈s (θ̂/i ) 2 ≤ 2c2 (n) ∆∗/i 2 .

(S.2.53)

17

Likewise, Z 1

Z 1 0

≤

γ−(1−t)∆∗/i /i (θ̂/i ) dt

0

2

Z 1 ≤ 0

Z 1 + 0

γ−(1−t)∆∗/i /i (θ̂/i ) 2 dt ℓ̈/i (θ̂/i − (1 − t)∆∗/i ) − ℓ̈/i (θ̂/i ) 2 dt
λ r̈s (θ̂/i − (1 − t)∆∗/i ) − r̈s (θ̂/i ) 2 dt

≤ 2c2 (n) ∆∗/i 2 .

(S.2.54)

By the result in Section S.1.1, we have ∆∗/i 2 ≤

ℓ̇i (θ̂) ν

! xi 2 ≤

c1 (n) ∥xi ∥2 . ν

(S.2.55)

Using results in (S.2.53), (S.2.54) and (S.2.55), we obtain 1 1 |A1 | ≤ |ℓ̇i (θ̂)| · (xi + xn+1 )⊤ Φ1 (xi + xn+1 ) + |ℓ̇i (θ̂)| · (xi + xn+1 )⊤ Φ2 (xi + xn+1 ) 2 2     2 2 ω̄max,i ≤|ℓ̇i (θ̂)| · γ−∆∗/i /i (θ̂/i ) + γ−∆∗/i /i (θ̂/i ) X̄/i J −1 ( θ̂ )x i /i /i ν 2 2 4     2 2 ω̄max,i + |ℓ̇i (θ̂)| · γ−∆∗/i /i (θ̂/i ) + γ−∆∗/i /i (θ̂/i ) X̄/i J −1 (θ̂/i )xn+1 /i ν 2 2 4 !  Z 1  Z 1 2 2 ω̄max,i γ−(1−t)∆∗/i /i (θ̂/i ) dt X̄/i J −1 ( θ̂ )x + |ℓ̇i (θ̂)| · γ−(1−t)∆∗/i /i (θ̂/i ) dt + i /i /i ν 4 0 0 2 2 !   Z 1 Z 1 2 2 ω̄max,i X̄/i J −1 (θ̂/i )xn+1 γ−(1−t)∆∗/i /i (θ̂/i ) dt + |ℓ̇i (θ̂)| · γ−(1−t)∆∗/i /i (θ̂/i ) dt + /i ν 4 0 0 2 2  2   2 2c1 (n)c2 (n) c (n)c2 (n) −1 1+ (1 + ω )∥x ∥ X̄ J ( θ̂ )x ≤4∥xi ∥2 1 max,i i 2 n+1 /i /i /i ν ν2 4    2 2 2c1 (n)c2 (n) c (n)c2 (n) 1+ (1 + ωmax,i )∥xi ∥2 X̄/i J −1 (θ̂/i )xi + 4∥xi ∥2 1 /i 2 ν ν 4  q  2 2c1 (n)c2 (n) c (n)c2 (n) 4 4 ≤4∥xi ∥2 1 1+ (1 + ωmax,i )∥xi ∥2 X/i J −1 (θ̂/i )xn+1 4 + J −1 (θ̂/i )xn+1 4 /i /i 2 ν ν  q  2 2c1 (n)c2 (n) c1 (n)c2 (n) 4 4 −1 1+ (1 + ωmax,i )∥xi ∥2 X/i J −1 + 4∥xi ∥2 /i (θ̂/i )xi 4 + J /i (θ̂/i )xi 4 ν ν2 h i h i 2 2 2 2 −1 −1 −1 + J ( θ̂ )x ≤C̄i X/i J −1 ( θ̂ )x + J ( θ̂ )x + C̄ X J ( θ̂ )x i /i n+1 4 /i /i /i i 4 /i i 4 , /i n+1 4 /i /i /i (S.2.56) 



⊤ ), ω ⊤ ⊤ where ωmax,i := σmax (X/i X/i max := σmax (XX ) and ω̄max,i := σmax X̄/i X̄/i . By this definition, we can derive that  ⊤  !     X X /i /i ⊤ ⊤ ω̄max,i = σmax X̄/i X̄/i = σmax X̄/i X̄/i = σmax I I     ⊤ ⊤ = σmax I + X/i X/i ≤ 1 + σmax X/i X/i = 1 + ωmax,i , (S.2.57)

For |A2 | and |A3 |, we can use a similar method to show that i 1 h 2 2 −1 ( θ̂ )x + J ( θ̂ )x |A2 | ≤ C̄i X/i J −1 /i n+1 4 /i n+1 4 /i /i 2 18

(S.2.58)

and i 1 h 2 2 −1 −1 |A3 | ≤ C̄i X/i J /i (θ̂/i )xi 4 + J /i (θ̂/i )xi 4 . 2

(S.2.59)

Combining the results in (S.2.56), (S.2.58) and (S.2.59), we obtain i 3 h 2 2 −1 −1 ⊤ |x⊤ θ̃ − x θ̂ | ≤ J ( θ̂ )x + X J ( θ̂ )x C̄ i /i i 4 /i /i /i i 4 n+1 /i n+1 /i /i 2 h i 3 2 2 −1 + C̄i X/i J −1 ( θ̂ )x + J ( θ̂ )x n+1 n+1 /i /i /i /i 4 4 2 h ≤3C̄i ·

max

j∈{i,n+1}

i 2 2 −1 + J ( θ̂ )x X/i J −1 ( θ̂ )x . j j /i /i /i /i 4 4

(S.2.60)

Next, we define the event  G :=

log(p) ⊤ max |x⊤ n+1 θ̂/i − xn+1 θ̃/i | ≥ C √ 1≤i≤n p

 .

Adopting the reasoning on pages 54–55 of [Rad and Maleki, 2020], where S denotes the event that Assumption 3 in our paper holds (trivially we have P(S c ) ≤ qn + q̃n ), and qn and q̃n follow the definitions in Assumption 3, we have P[G] ≤ P[G|S] + P[S c ] log(p) ⊤ |S] + qn + q̃n |x⊤ j θ̂/i − xj θ̃/i | > C √ 1≤i≤n j∈{i,n+1} p

≤ P[ max ( max ≤ ≤

log(p) log(p) ⊤ √ ] + P[max1≤i≤n |x⊤ θ̂/i − x⊤ θ̃/i | > C √ ] P[max1≤i≤n |x⊤ n+1 n+1 i θ̂/i − xi θ̃/i | > C p p

1 − qn − q̃n h i 2 2 −1 √ ] P[3C̄i · X/i J −1 ( θ̂ )x J ( θ̂ )x > C log(p) + i i /i /i p /i /i 4 4

+

P[3C̄i ·

h

1 − qn − q̃n 2 −1 X/i J /i (θ̂/i )xn+1 4 +

i 2 log(p) J −1 ( θ̂ )x /i n+1 4 > C √p ] /i

1 − qn − q̃n

+ qn + q̃n .

+ qn + q̃n

(S.2.61)

h i 2 2 −1 √ ] and For P[3C̄i · X/i J −1 ( θ̂ )x + J ( θ̂ )x > C log(p) i i /i /i p /i /i 4 4 h i 2 2 log(p) −1 −1 P[3C̄i · X/i J /i (θ̂/i )xn+1 4 + J /i (θ̂/i )xn+1 4 > C √p ], we have P[3C̄i ·

h

n i X log(p) 2 2 −1 X/i J −1 ( θ̂ )x + J ( θ̂ )x > C P[Ẽi ]; ] ≤ √ i i /i /i /i /i 4 4 p

(S.2.62)

n i X log(p) 2 2 −1 X/i J −1 ( θ̂ )x + J ( θ̂ )x > C ] ≤ P[Ẽn+1,i ]. √ n+1 n+1 /i /i /i /i 4 4 p

(S.2.63)

i=1

P[3C̄i ·

h

i=1

For simplicity, we define  q  n−1 log(n−1) 1 + ωmax , Ci := 6(1 + c) p log p     2 2 √ −1 −1 > C̄i Ci p log p , Fn+1,i := 3 C̄i X/i J /i xn+1 + J /i xn+1 4 4   C √ Ki := √ ≥ C̄i Ci p , p   √ √ Wn+1 := ∥xn+1 ∥22 > 5pρmax ∪ ωmax > ( n + 3 p)2 ρmax . ρ

max ν2

19

Hence, we now obtain an upper bound for P[Ẽn+1,i ] P[Ẽn+1,i ] ≤ P[Ẽn+1,i | Ki ] + P[Kic ] P(Fn+1,i ) ≤ P[Fn+1,i | Ki ] + P[Kic ] ≤ + P[Kic ] P(Ki ) h i
√ 2 ρmax ω P X/i J −1 x > 2(1 + c) n − 1 log(n − 1) max n+1 2 /i ν 4 ≤ P(Ki ) h i 2 ρmax
√ > 2(1 + c) x P J −1 p log p /i n+1 4 ν2 + + P[Kic ] P(Ki )   2 1 2 1 + + P[Kic ], ≤ (n − 1)c pc P(Ki )

(S.2.64)

where ≤1 is due to (S.1.37). To bound P[Kic ] we define !  q q 2 √ c21 (n)c2 (n) (pρmax )3/2 n n−1 log(n−1) 1+ pρmax C := 96 5 p +3 p log p ν3 √

   2 q √ n 2c (n)c (n) + 3 5 pρ pρ 1 + 1 2 max max  p  . × 1 +   ν2 

p √ √ obtained by setting c = 3, and computing pC̄i Ci after putting 5pρmax and ( n/p + 3 p)2 ρmax bounds in event Wn+1 , into ∥xn+1 ∥2 and ωmax , respectively. Next,     C c c < C̄i Ci ≤ P C < pC̄i Ci Wn+1 + P[Wn+1 ] P[Ki ] = P p = P[C < C] + P[Wn+1 ] = P[Wn+1 ]. (S.2.65) The term P[Wn+1 ] is exponentially small because xn+1 ∼ N (0, Σ) with ρmax = σmax (Σ), leading to i h   √ √ P[Wn+1 ] ≤ P ∥xn+1 ∥22 > 5pρmax + P σmax (XX ⊤ ) > ( n + 3 p)2 ρmax ≤ 2e−p , due to result in Section S.1.3 and Section S.1.4. In summary, since for p ≥ 1 we have 1−e1−p < 2, for c = 3 we obtain 4 4 P[Ẽn+1,i ] ≤ + 3 + 2e−p . 3 (n − 1) p Since xi and xn+1 follow same distribution, that is xi ∼ N (0, Σ), the same argument implies that P[Ẽi ] ≤

4 4 + 3 + 2e−p . 3 (n − 1) p

The detailed proof for the bound of P[Ẽi ] can also be found in [Rad and Maleki, 2020]. As a result, we have   1 8n 8n −p P[G] ≤ + 3 + 4ne + qn + q̃n 1 − qn − q̃n (n − 1)3 p   16n 16n ≤ + 3 + 8ne−p + qn + q̃n . 3 (n − 1) p h  i 16n 16n −p + q + q̃ It is clear that (n−1) n n → 0 when n → ∞. 3 + p3 + 8ne 20

(S.2.66)

similar to the result in Section S.1.5. That is, with probability at least 1 −  Thus, we obtain a result  16n 16n −p + p3 + 8ne − qn − q̃n , the following bound is valid (n−1)3 C0 ⊤ max |x⊤ n+1 θ̃/i − xn+1 θ̂/i | ≤ √ , p

1≤i≤n

where C0 is defined as follows !
! √  p p 5 c1/2 + c3/2 ( δ0 + 3)2 216c3/2 2 3 2 2 c log n 1 + δ0 ( δ0 + 3) c1 (n)c2 (n) + c1 (n)c2 (n) . C0 := ν3 log p ν2 (S.2.67) √ √ Note that in the presentation ofqTheorem 3.1, we replaced respectively q 96 5 and 2 3 with the upper-

bounds 216 and 5, and we replaced n/p.

S.3

n−1 log(n−1) with the upper bound p log p

n log n p log p . We also used δ0 to denote

Proof of Theorem 3.2

Using result Section S.1.5, with probability at least 1 − p2 (n), the following holds Q0 ⊤ max |x⊤ i θ̃/i − xi θ̂/i | ≤ √ . 1≤i≤n p

(S.3.68)

Then with probability at least 1 − p2 (n), following holds ⊤ max |x⊤ i θ̃/i − xi θ̂/i | = o(1).

1≤i≤n

(S.3.69)

Plugging (S.3.69) into RiLOO − R̃iLOO , using Assumption 4, with probability at least 1 − p2 (n) the following result holds ⊤ RiLOO − R̃iLOO = (|yi − f (x⊤ i θ̂/i )| − |yi − f (xi θ̃/i )|) ⊤ ≤ f (x⊤ i θ̂/i ) − f (xi θ̃/i ) ⊤ ≤L · max |x⊤ i θ̂/i − xi θ̃/i | 1≤i≤n

(S.3.70)

=o(1). Using Theorem 3.1 and Assumption 4, with probability at least 1 − p1 (n) the following result holds ⊤ ⊤ ⊤ f (x⊤ n+1 θ̂/i ) − f (xn+1 θ̃/i ) ≤ L · xn+1 θ̃/i − xn+1 θ̂/i = o(1).

(S.3.71)

Then with probability at least (1 − p1 (n) − p2 (n)), the following results hold for i = 1, 2, ..., n: LOO LOO (f (x⊤ ) − (f (x⊤ ) = o(1); n+1 θ̃/i ) − R̃i n+1 θ̂/i ) − Ri

(S.3.72)

LOO LOO (f (x⊤ ) − (f (x⊤ ) = o(1). n+1 θ̃/i ) + R̃i n+1 θ̂/i ) + Ri

(S.3.73)

With (S.3.72) and (S.3.73), using Lemma 2, we can easily derive that − LOO − LOO q̂n,α {f (x⊤ } = q̂n,α {f (x⊤ } + o(1); n+1 θ̃/i ) − R̃i n+1 θ̂/i ) − Ri

21

(S.3.74)

+ LOO + LOO q̂n,α {f (x⊤ } = q̂n,α {f (x⊤ } + o(1). n+1 θ̃/i ) + R̃i n+1 θ̂/i ) + Ri

(S.3.75)

That is, with probability at least (1 − p1 (n) − p2 (n)), the following results hold Ũjk+ = Ujk+ + o(1), L̃jk+ = Ljk+ + o(1).

(S.3.76)

We can utilize similar method to conclude that, with probability at least (1 − p1 (n) − p2 (n)), the following results hold Ũjkm = Ujkm + o(1), L̃jkm = Ljkm + o(1).

(S.3.77)

Theorem 3.2 has been proved.

S.4

Proof of Theorem 3.3 & 3.4

Recall that in Assumption 4 we defined that   16n 16n −p + + 8ne + qn + q̃n p1 (n) = (n − 1)3 p3 and

 p2 (n) =

8n 8n + 3 + 4ne−p 3 (n − 1) p

 + qn + q̃n .

Under Assumption 1 to Assumption 4, and the condition that n, p → ∞, n/p = δ0 ∈ (0, ∞), it is clear that n→∞

n→∞

p1 (n) −−−−−→ 0,

p2 (n) −−−−−→ 0.

Using (S.3.74), (S.3.75) and Assumption 4, we obtain that, with probability at least (1 − p1 (n) − p2 (n)), the following holds Jackknif e+ Jackknif e+ Jackknif e+ P{yn+1 ∈ C̃n,α } = P{yn+1 ∈ Ĉn,α } + L · o(1) = P{yn+1 ∈ Ĉn,α } + o(1). (S.4.78)

According to [Barber et al., 2021], we have Jackknif e+ − LOO + LOO Ĉn,α = [q̂n,α {f (x⊤ }, q̂n,α {f (x⊤ }] n+1 θ̂/i ) − Ri n+1 θ̂/i ) + Ri

(S.4.79)

Jackknif e+ P{yn+1 ∈ Ĉn,α } ≥ 1 − 2α.

(S.4.80)

and As s result, using (S.4.80), we have Jackknif e+ P{yn+1 ∈ C̃n,α } ≥ (1 − 2α − o(1)) · (1 − p1 (n) − p2 (n)) = 1 − 2α − o(1).

(S.4.81)

Theorem 3.3 has been proved. Proof of Theorem 3.4 is similar to Theorem 3.3.

S.5

Proof of lemmas

S.5.1

Proof of Lemma 2

Lemma 2 follows from the following more precise version: Lemma 6. Let {ai }ni=1 be a sorted sequence such that a1 ≤ a2 ≤ · · · ≤ an , and suppose there exist sequences {bi }ni=1 satisfying |bi − ai | ≤ ε, for some ε > 0 and all i = 1, . . . , n. Then for any q ∈ {1, 2, ..., n} we have |b(q) − aq | ≤ ϵ, where b(q) is the q-th smallest value of {bi }ni=1 . Proof. Under the constraint of bi − ai ≤ ϵ, the maximum of b(q) is achieved when bi = ai + ϵ for each i, in which case b(q) = aq + ϵ. Hence b(q) ≤ aq + ϵ. Similarly, b(q) ≥ aq − ϵ, and the proof is completed. 22

S.6

Proof of Theorem 5.6

The proof of Theorem 5.6 relies on some auxiliary results (for simplicity, we define µ̂/i (X) = X ⊤ β̂/i and µ̂(X) = X ⊤ β̂): Lemma 7. Under assumptions C1-C5, we have the following (in the sense of weak convergence) (S.6.82)

max |µ̂/i (Xn+1 ) − µ̂(Xn+1 )| = o(1).

1≤i≤n

Proof. Recall that [Karoui, 2018, Theorem 2.2] showed that  max ∥β̂ − β̂/i − ηi ∥2 = OL2

1≤i≤n

polylog(n) n

 ,

(S.6.83)

where the notations are explained as follows: OL2 (·) denotes a bound on a random variable in the square root of the second moment, and β̂/i + ηi can be thought of as an approximate formula for β̂ using the leave-one-out estimator β̂/i . The correction term ηi :=

1 (Si + τ I)−1 Xi ρ′ proxci ρi (ri,(i) ) , n

(S.6.84)

1 ⊤ X (Si + τ I)−1 Xi ; n i 1 X ′′ Si := ρ (rj,(i) )Xj Xj⊤ , n

(S.6.85)

where we defined ci :=

(S.6.86)

j̸=i

where rj,(i) := yj − Xj⊤ β̂/j

(S.6.87)

denotes the leave-one-out residual (in particular, ri,(i) = RiLOO ). Noting that Si , Xi , and Xn+1 are independent, and Xi and Xn+1 are i.i.d. with zero mean and identity ⊤ (S + τ I)−1 X has zero mean and O( 1 ) second moment. Under the assumption covariance, we see that n1 Xn+1 i i n ⊤ η = o (1) (vanishing in the sense of second of ∥ρ′ ∥∞ < ∞, this in turn implies, by (S.6.84), that Xn+1 i L2 √ ⊤ (β̂−β̂ )∥ = moment). Then using the fact that ∥Xn+1 ∥2 = OL2 ( n) and (S.6.83) we obtain max1≤i≤n ∥Xn+1 /i 2 oL2 (1), which is equivalent to (S.6.82). Next, recall that µ̂y (Xi ) := Xi⊤ β̂ y , i = 1, . . . , n + 1, where β̂ y denotes the output of the regression algorithm trained on the augmented dataset {(X1 , y1 ), . . . , (Xn , yn ), (Xn+1 , y)}. Also, adopt the notations Ri := yi − Xi⊤ β̂, ′

ψ(t) := ρ (t),

i = 1, . . . , n,

t ∈ R.

(S.6.88) (S.6.89)

We then have: Lemma 8. For any given M > 0, we have the following (in the sense of convergence in probability): max |µ̂y (Xi ) − µ̂(Xi )| = o(1);

(S.6.90)

max µ̂/i (Xi ) + ci ψ(proxci ρi (ri,(i) )) − µ̂(Xi ) = o(1).

(S.6.91)

max

1 2 y∈[−M,M ]∩{0,± M ,± M ,... } 1≤i≤n

max

1 2 y∈[−M,M ]∩{0,± M ,± M ,... } 1≤i≤n

23

Proof. The proof of (S.6.90) is almost identical to Lemma 7 (by viewing µ̂ as the estimator leaving out the 1 2 (n + 1)-th observation, and then taking the union bound over all y in the set [−M, M ] ∩ {0, ± M , ±M , . . . }), which relies on the error bound on the leave-one-out approximation formula (S.6.83). It still remains, however, to check whether the proof depends on how yn+1 = y is selected. Below, we provide details of justification that (S.6.82) continues to hold if yi is a fixed value in R (rather than randomly generated from the stochastic model), which is equivalent to (S.6.90) (by switching the role of n + 1 and i). This, in turn, relies on two properties: a) (S.6.83) continues to hold when yi is an arbitrary deterministic quantity. ⊤ η = o (1) continues to hold when y is an arbitrary deterministic quantity. b) Xn+1 i i L2 ⊤ (S + It is easy to see that b) is true, under our assumption of ∥ρ′ ∥∞ < ∞, since it is still true that n1 Xn+1 i 1 τ I)−1 Xi has zero mean and O( n ) second moment. For a), we revisit the original proof [Karoui, 2018, p122p125], which is mostly based on several deterministic non-asymptotic bounds therein, except [Karoui, 2018, Lemma 3.7], which relies on convergence in high probability. However, [Karoui, 2018, Lemma 3.7] uses nothing ∥(S +τ I)−1 X ∥ ∥(Si +τ I)−1 Xj ∥ √ but a high probability bound that supj̸=i i τ √n j = O(1), the independence of and Xi , τ n X ⊤ (Si +τ I)−1 Xi

and sub-Gaussianity of Xi , which in turn implies by union bound that supj̸=i j = √1n polylog(n) n with high probability. Thus a) also holds, and we have verified (S.6.90). Next, we observe that (S.6.91) because of (S.6.83) and (S.6.84), where again we used the fact that (S.6.83) continues to hold for arbitrary deterministic yi . Lemma 9. For any given M > 0, we have the following (in the sense of convergence in probability): max

1 2 y∈[−M,M ]∩{0,± M ,± M ,... }

|µ̂(Xn+1 ) + cψ(proxcρ (y − µ̂(Xn+1 ))) − µ̂y (Xn+1 )| = o(1),

(S.6.92)

1 X ⊤ (S + τ I)−1 Xn+1 ; n + 1 n+1 n 1 X ′′ S := ρ (y − µ̂(Xn+1 ))Xj Xj⊤ . n+1

(S.6.93)

where we defined c :=

(S.6.94)

j=1

Proof. The proof is identical to the proof of (S.6.91), again noting that the approximation error bound ((S.6.83), but leaving out n + 1 instead of i) continues to hold with fixed deterministic yn+1 = y. Proof of Theorem 5.6. Lemma 7 and Lemma 6 imply that |UJ − Umm | = o(1), and + |µ̂(Xn+1 ) + q̂n,α {RiLOO } − UJ | = o(1).

(S.6.95)

Next, recall that [Karoui, 2018, Appendix 5] shows that ri,(i) = RiLOO behaves like (in the sense of moment convergence) r   ϵi + E ∥β̂ − β0 ∥2 Zi , where Zi ∼ N (0, 1) independent of ϵi ; furthermore, if i ̸= j, ri,(i) and rj,(j) are asymptotically (pairwise) independent. Since limp→∞ ∥β̂ − β0 ∥ = r, this implies that the empirical distribution (ri,(i) )ni=1 converges to N (0, σϵ2 + r2 ). Thus p + |q̂n,α {RiLOO } − σϵ2 + r2 zα/2 | = o(1) (S.6.96) and we have proved (5.25).

24

1 2 It remains to prove (5.26). Note that by definition, y ∈ [−M, M ] ∩ {0, ± M , ±M , . . . } is included in the full conformal interval if

|y − µ̂y (Xn+1 )| ≤ q̂ + {Ri }.

(S.6.97)

Suppose that ymax and ymin denote the maximum and minimum y under this criterion. By Lemma 9, we see that ymax satisfies ymax − µ̂(Xn+1 ) − cψ(proxcρ (ymax − µ̂(Xn+1 ))) = q̂ + {Ri } + en,M ,

(S.6.98)

where en,M is an error term satisfying limM →+∞ limn→∞ en,M = 0. By (S.6.83), we have Ri − RiLOO = −ci ψ(Ri ) + o(1),

(S.6.99)

implying Ri = proxci ρ (RiLOO + o(1)). Furthermore, using the Sherman–Morrison formula we can show that ci = c + o(1), so that Ri = proxcρ (RiLOO + o(1)). Since proxcρ is an increasing function, we have q̂ + {Ri } = q̂ + {proxcρ (RiLOO + o(1))} = proxcρ (q̂

+

{RiLOO + o(1)}).

(S.6.100) (S.6.101)

Then (S.6.98) shows that proxcρ (ymax − µ̂(Xn+1 )) = proxcρ (q̂ + {RiLOO + o(1)}) + en,M .

(S.6.102)

Since we assumed ∥ρ′′ ∥∞ < ∞, and c = O(1) with high probability, prox−1 cρ must be Lipschitz with high probability. Hence ymax − µ̂(Xn+1 ) = q̂ + {RiLOO + o(1)} + ẽn,M for some ẽn,M satisfying limM →∞ limn→∞ ẽn,M = 0. Then (5.26) follows by (S.6.96).

25

(S.6.103)