A Kernel Nonconformity Score for Multivariate Conformal Prediction
arXiv:2604.21595v1 [stat.ML] 23 Apr 2026
Louis Meyer University of Warwick [email protected]
Wenkai Xu University of Warwick [email protected]
Abstract Multivariate conformal prediction requires nonconformity scores that compress residual vectors into scalars while preserving certain implicit geometric structure of the residual distribution. We introduce a Multivariate Kernel Score (MKS) that produces prediction regions that explicitly adapt to this geometry. We show that the proposed score resembles the Gaussian process posterior variance, unifying Bayesian uncertainty quantification with the coverage guarantees of frequentisttype. Moreover, the MKS can be decomposed into an anisotropic Maximum Mean Discrepancy (MMD) that interpolates between kernel density estimation and covariance-weighted distance. We prove finite-sample coverage guarantees and establish convergence rates that depend on the effective rank of the kernelbased covariance operator rather than the ambient dimension, enabling dimensionfree adaptation. On regression tasks, the MKS reduces the volume of prediction region significantly, compared to ellipsoidal baselines while maintaining nominal coverage, with larger gains at higher dimensions and tighter coverage levels.
1
Introduction
Conformal prediction constructs prediction regions with finite-sample coverage guarantees under minimal assumptions on the data distribution [31, 38]. Given a point predictor fˆ and a confidence level 1 − α, the framework produces a region C(X) satisfying P(Y ∈ C(X)) ≥ 1 − α for any model class and any distribution [2]. The univariate case is well-understood, conformalized quantile regression [27] and standardized residual methods [19] provide simple and effective prediction intervals. Exact conditional coverage P(Y ∈ C(X) | X) ≥ 1 − α cannot be achieved distribution-free without structural assumptions [11, 37]. The efficiency of the score, measured by its size, depends on the choice of nonconformity score and the order it induces. For multivariate responses Y ∈ Rd , the nonconformity score determines the underlying “shape” of the prediction region, where the efficiency of the region depends on how faithfully this shape reflects the geometry of the residual distributions. Designing scores that yield data-adaptive regions, without requiring/learning a generative model or imposing a fixed parametric shape, remains an open problem. The difficulty is intrinsic to the multivariate setting. In one dimension, residuals reside in an ordered set R, and admit a natural total order where the quantile is well-defined. In Rd where no canonical ordering exists, the score function S : Rd → R imposes a total order by compressing a d-dimensional residual into a scalar, and this compression unavoidably discards information. Different choices of S discard different structural properties of the residual distribution, and the efficiency of the prediction region depends predominately on the information preserved. Existing approaches split into two families: ellipsoidal approaches and density-based methods. Ellipsoidal approaches score residuals by their Mahalanobis distance under the sample covariance, producing regions aligned with the principal axes of the residual distribution [16, 40]. While still being efficient when the residuals are approximately Gaussian, the ellipsoid must expand into lowPreprint.
density areas to maintain coverage when the residual geometry departs from ellipticity, whether from model misspecification, multimodality, or nonlinear dependences between output coordinates, inflating the region unnecessarily. Density-based methods instead threshold an estimate of the residual density, producing level sets that can take arbitrary shape [9, 15]. In practice, nonparametric density estimation deteriorates as the output dimension grows [29], and even a perfect density estimate ranks residuals by a single scalar that is invariant to directions: two residuals at equal density but in structurally different orientations relative to the calibration data receive identical scores. Ellipsoidal scores capture directional covariance but are restricted to a fixed geometry, while density-based scores capture inherent geometry without preserving the directional information. We introduce a kernel-based nonconformity score that combines desirable properties from both. The score embeds calibration residuals into a reproducing kernel Hilbert space (RKHS) and measures the conformity of a test residual through a regularized quadratic form involving the empirical covariance operator [12] via its feature space. We prove that the proposed score resembles the posterior variance of a Gaussian process (GP) [39] evaluated at the calibration residual locations. The GP model does not need to be correctly specified. The posterior variance serves as a scoring function whose geometric properties we inherit from Bayesian nonparametrics; the coverage guarantee rests on conformal calibration, which requires only exchangeability of the scores. Our approach resembles Bayesian geometry for the design of the prediction region while achieving frequentist-type validity. The MKS identity yields a decomposition into two interpretable components from the kernel literature. The first is a Maximum Mean Discrepancy (MMD) term [13] that measures the discrepancy between the empirical test residual and the calibration distribution, acting as a one-sample typicality test through a kernel density score. The second is a kernel principal component analysis (KPCA) [30] correction term that accounts for the anisotropy of the residual distribution. The prediction regions that result are neither density level sets nor ellipsoids, but kernel-adapted regions shaped by both the density and the directional covariance of the residuals. When the kernel is linear, the score reduces to the regularized Mahalanobis distance, so the framework, in principle, generalizes the ellipsoidal methods of Johnstone et al. [16] and Xu et al. [40]. Our contributions are threefold. (1) We propose MKS: a kernel-based nonconformity score for multivariate conformal prediction and prove that it equals a GP posterior variance (Proposition 4.1). We show that the score decomposes into an MMD term and a KPCA correction (Proposition 4.2); and that the linear kernel recovers the Mahalanobis distance (Proposition 4.4), establishing that our method strictly generalizes ellipsoidal conformal prediction. (2) We provide finite-sample coverage guarantees under exchangeability, with convergence rates depending on the effective rank of the RKHS covariance operator rather than the ambient dimension d (Theorem 3.4); we extend these guarantees to stationary α-mixing sequences (Theorem B.20 in Appendix B.3), following the framework of Xu et al. [40]. (3) We demonstrate on synthetic and real data that the proposed kernel-based score produces smaller prediction regions than Mahalanobis and density-only baselines while maintaining nominal coverage, with the advantage increasing substantially at higher coverage levels and at higher output dimensions.
2
Background
Split conformal prediction. We consider a regression setting with inputs X ∈ X and multivariate responses Y ∈ Rd . Given a dataset {(Xi , Yi )}ni=1 drawn from a distribution PX,Y , the goal is to construct a prediction region C(Xnew ) satisfying P Ynew ∈ C(Xnew ) ≥ 1 − α (1) for a user-specified level α ∈ (0, 1). Following the split conformal framework [19, 22], the data is randomly partitioned into a training set Dtrain of size ntrain and a calibration set Dcal = {(Xi , Yi )}Ti=1 of size T . A model fˆ : X → Rd is fitted on Dtrain , and the residuals εi = Yi − fˆ(Xi ) ∈ Rd are computed on the calibration set Dcal . A nonconformity score function S : Rd → R maps each residual to a scalar measuring its degree of nonconformity, with larger values indicating less typical residuals. The conformal threshold q̂ is defined as the ⌈(1 − α)(T + 1)⌉/T quantile of the calibration scores {S(ε1 ), . . . , S(εT )}, yielding the prediction region C(X) = y ∈ Rd : S y − fˆ(X) ≤ q̂ (2) 2
If the calibration scores and the test score are exchangeable, then C(X) satisfies (1) [38]. This guarantee is distribution-free, i.e. the result holds for any data distribution, any model fˆ, and any score function S. The choice of S does not affect coverage but determines the shape and efficiency of the prediction region, e.g. the length of prediction interval in one dimension. Exchangeability. The coverage guarantee (1) requires the calibration residuals ε1 , . . . , εT and the test residual εT +1 to be exchangeable. Under the split conformal setup with i.i.d. data, this holds by construction. When the data exhibits temporal or distributional dependence, exchangeability generally fails. Barber et al. [3] provide a general framework for conformal prediction beyond exchangeability, bounding the coverage gap in terms of the total variation distance between the true joint distribution and the nearest exchangeable one. Xu et al. [41] derive finite-sample coverage guarantees under α-mixing dependence for ellipsoidal scores, showing that coverage deviations vanish as the mixing coefficients decay. Other forms of non-exchangeability, including covariate shift [35], missing data [42], and contamination [8], have also been studied. Multivariate nonconformity scores. The choice of S determines the geometry of the prediction region (2), since the region is a sublevel set of S. The simplest approach applies univariate conformal prediction independently to each coordinate at level α/d, producing hyper-rectangular regions that ignore all dependence between outputs. Copula-based methods [20, 33] calibrate the marginals jointly but still yield rectangular regions. b + λI)−1 (ε − ε̄) produces ellipsoidal regions aligned The Mahalanobis distance S(ε) = (ε − ε̄)⊤ (Σ with the principal axes of the residual distribution [16]. Xu et al. [40] use this score with a global covariance estimate for multivariate time series, providing coverage guarantees under α-mixing dependence. Braun et al. [5] learn an input-dependent covariance Σ(X), producing locally adapted ellipsoids with improved conditional coverage. Messoudi et al. [21] estimate local covariance structures using nearest neighbors. Density-based approaches use the estimated density p̂(ε) as a nonconformity score, producing regions that are level sets of the density estimate [15, 29]. In practice, nonparametric density estimation becomes unreliable as the output dimension grows, and scalable alternatives require fitting conditional generative models such as normalizing flows [9, 24]. Optimal transport approaches [34, 17] and adaptive norm-based formulations [4] offer more flexible region shapes at the cost of solving expensive optimization problems or training dedicated networks. Dheur et al. [9] provide a comprehensive comparison of multivariate conformal methods within a unified framework. A common limitation of density-based scores is that they rank residuals by a single scalar summary, local density, that is invariant to direction. Two residuals at equal density but different orientations relative to the calibration covariance receive identical scores. Ellipsoidal methods are well suited to residual distributions with approximately elliptical contours, but when the residual geometry is non-elliptical, the ellipsoid must expand into low-density regions to maintain coverage.
3
Multivariate Kernel Score
We introduce a kernel-based nonconformity score that extends the Mahalanobis distance to nonlinear settings. Let {εi }Ti=1 ⊂ Rd be calibration residuals obtained from a split conformal procedure, with PT b = 1 PT (εi − ε̄)(εi − ε̄)⊤ . empirical mean ε̄ = T1 i=1 εi and sample covariance Σ i=1 T −1 Definition 3.1 (Regularized Mahalanobis Distance). The regularized Mahalanobis distance nonconformity score is b + λId −1 (ε − ε̄) MD(ε) = (ε − ε̄)⊤ Σ (3) where λ > 0 is a regularization parameter. Let k : Rd × Rd → R be a positive definite kernel with feature map ϕ : Rd → H, where H is a reproducing kernel Hilbert space satisfying k(x, x′ ) = ⟨ϕ(x), ϕ(x′ )⟩H . Definition 3.2 (Multivariate Kernel Score (MKS)). The kernel nonconformity score is the Mahalanobis distance computed in the feature space H: bϕ + λIH −1 ϕ̃(ε)⟩H êk (ε) = ⟨ϕ̃(ε), C (4) 3
Algorithm 1 Kernel Conformal Prediction Require: Calibration residuals {εi }Ti=1 , kernel k, regularization parameter γ, level α 1: Compute centered Gram matrix K̃ with 1X 1 X 1X K̃ij = k(εi , εj ) − k(εi , εl ) − k(εl , εj ) + 2 k(εl , εm ). T T T l
l
l,m
4: êk (εi ) ← k̃(εi , εi ) − k̃∗ (εi )⊤ (K̃ + γIT )−1 k̃∗ (εi ) 5: end for 6: q̂ ← ⌈(1 − α)(T + 1)⌉/T quantile of {êk (εi )}T i=1
(as in Eq. (5))
2: Precompute (K̃ + γIT )−1 3: for i = 1, . . . , T do
Ensure: Prediction region C(X) = {y : êk (y − fˆ(X)) ≤ q̂}
PT where ϕ̃(ε) = ϕ(ε) − µ̂ϕ is the centered feature map, µ̂ϕ = T1 i=1 ϕ(εi ) is the kernel mean P T bϕ = 1 embedding, and C i=1 ϕ̃(εi ) ⊗ ϕ̃(εi ) is the sample covariance operator on H [12]. T −1 The Tikhonov regularization ensures the score is well-defined. In an RKHS with infinite-dimensional feature space, the covariance operator has eigenvalues decaying to 0 and Cϕ−1 is unbounded, but the regularized operator Cϕ + λIH has operator norm at most 1/λ. In the Gaussian Process interpretation developed in Section 4, λ corresponds to the variance of the observation noise, controlling the balance between interpolating the calibration data and smoothing over it. Using the kernel trick and the Woodbury matrix identity (proof in Appendix A.1), the score can be computed entirely in terms of kernel evaluations. The procedure is summarized in Algorithm 1. Proposition 3.3 (Computable Form). Let K̃ ∈ RT ×T be the centered Gram matrix with entries K̃ij = ⟨ϕ̃(εi ), ϕ̃(εj )⟩H , and let k̃∗ (ε) ∈ RT be the centered kernel vector with [k̃∗ (ε)]i = ⟨ϕ̃(ε), ϕ̃(εi )⟩H . Setting γ = λ(T − 1), the Multivariate Kernel Score is êk (ε) = k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε)
(5)
where k̃(ε, ε) = ∥ϕ̃(ε)∥2H . Given a target miscoverage level α ∈ (0, 1), we define the conformal threshold q̂ as the ⌈(1 − α)(T + 1)⌉/T quantile of {êk (ε1 ), . . . , êk (εT )}. The prediction region for a new input XT +1 with point prediction fˆ(XT +1 ) is then C(XT +1 ) = y ∈ Rd : êk y − fˆ(XT +1 ) ≤ q̂ (6) Theorem 3.4 (Coverage Guarantee). If the residuals ε1 , . . . , εT +1 are exchangeable, then the prediction region (6) satisfies P YT +1 ∈ C(XT +1 ) ≥ 1 − α (7) Under Assumptions B.1–B.4, which include the standard exogeneity condition εT +1 ⊥ XT +1 , the conditional coverage deviation satisfies P YT +1 ∈ C(XT +1 ) | XT +1 − (1 − α) ≤ C · T −r (8) with high probability, where r > 0 depends on the estimation rate ρ and regularization, and C depends on the effective rank r(Σ) = tr(Σ)/∥Σ∥op of the RKHS covariance operator. Since conditional coverage implies marginal coverage, (7) follows as a corollary. Under α-mixing, we establish marginal coverage guarantees in Appendix B.3. Remark. For smooth kernels such as the RBF kernel, the eigenvalues of Σ decay rapidly, so that r(Σ) can be much smaller than both d and T . The full convergence analysis and the proof of Theorem 3.4 are detailed in Appendix B. 4
4
Probabilistic and Geometric Structure
4.1
Gaussian Process Connection
Proposition 4.1 (Gaussian Process Identity). Let f ∼ GP(0, k̃) be a Gaussian process with zero mean and centered kernel k̃. Given observations at locations {εi }Ti=1 with noise variance σn2 = γ, the posterior variance at a test point ε is [39, Eq. 2.26]: VarGP f (ε) | {εi }Ti=1 = k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) (9) Comparing with Proposition 3.3, we have the exact identity êk (ε) = VarGP f (ε) | {εi }Ti=1
(10)
In standard GP regression, the posterior variance quantifies uncertainty about a function at unobserved locations, and depends on the observed function values y through the posterior mean. Our score depends only on the locations {εi }Ti=1 of the calibration residuals, not on any response observed at those locations. The MKS is a purely geometric quantity that measures how well a test residual is explained by the spatial configuration of the calibration data. A second distinction is that we use the centered kernel k̃ rather than the raw kernel k, introducing density sensitivity. For the RBF kernel, k(ε, ε) = 1 for any ε, while after centering, k̃(ε, ε) = 1 + k̄ − 2p̂k (ε) varies with the local density of the calibration distribution. Centering enables the score to distinguish residuals in dense regions from those in sparse ones. The regularization parameter γ enters the GP interpretation as the observation noise variance, and governs how much of the calibration data’s structure the score retains. Small γ corresponds to near-interpolation, where the GP posterior variance drops to nearly zero at each calibration point, and the score surface closely tracks the placement of individual residuals. Large γ smooths over individual observations, making the posterior variance reflect only the broad distributional structure, and the score approaches the MMD, as shown in Proposition 4.2. This interpretation connects Bayesian and frequentist reasoning within a single framework. The GP posterior variance provides geometrically rich machinery for designing the nonconformity score. It captures density, directional covariance, and manifold structure, producing prediction regions that adapt to the shape of the residual distribution. The conformal calibration provides validity, the coverage guarantee (7) holds regardless of whether the GP model is correctly specified, because it relies on exchangeability of the scores, not on Gaussianity of the process. 4.2
MMD-KPCA Decomposition
PT Proposition 4.2 (MMD–KPCA Decomposition). Let QT = T1 i=1 δεi denote the empirical distribution of the calibration residuals, where δεi is the Dirac measure at εi (the probability distribution that assigns all its mass to the point εi ). The kernel nonconformity score satisfies êk (ε) = MMD2 (δε , QT ) −
T X αj (ε)2 j=1
λj + γ
(11)
where MMD2 (δε , QT ) = ∥ϕ̃(ε)∥2H = k(ε, ε) + k̄ − 2p̂k (ε) is the squared Maximum Mean Discrepancy between the Dirac measure at ε and the (empirical) calibration distribution, with PT P p̂k (ε) = T1 i=1 k(ε, εi ) and k̄ = T12 i,j k(εi , εj ). The quantities {λj }Tj=1 are the eigenvalues of K̃, and αj (ε) = vj⊤ k̃∗ (ε) is the projection of ϕ̃(ε) onto the j-th kernel principal component of the calibration data, with {vj }Tj=1 the corresponding eigenvectors. The proof is in Appendix A.3. The first term builds on the Maximum Mean Discrepancy method of Gretton et al. [13], an integral probability metric that measures the distance between two distributions P and Q through their mean embeddings in the RKHS, defined as MMD2 (P, Q) = ∥µP − µQ ∥2H . When the kernel is characteristic, MMD(P, Q) = 0 if and only if P = Q [32]. The MMD tests for distributional equality and is widely used as a two-sample test statistic. In our setting, we apply the same principle with P = δε , a single test residual, and Q = QT , the empirical calibration distribution, so that the term measures the discrepancy between one point and the calibration set, effectively acting 5
as a one-sample typicality test. For the RBF kernel, k(ε, ε) = 1 for all ε and k̄Pdepends only on the calibration set, so the only term that varies with the test point is p̂k (ε) = T1 i k(ε, εi ), a ParzenRosenblatt kernel density estimate up to normalization [23, 28]. The MMD term is monotonically decreasing in estimated density, implying that the residuals in high-density regions receive low scores and those in sparse regions receive high scores. Using this term alone as a nonconformity score would produce prediction regions that are level sets of the kernel density estimate, as in the HPD-split method of Izbicki et al. [15]. The KPCA correction introduces anisotropy. It reduces the score for residuals whose feature-space representation aligns with directions of high calibration variance, where spread is expected and deviations are uninformative. The score is large only where the calibration data provides neither density support nor directional explanation, and small where either is present. The prediction region is therefore neither a density level set nor an ellipsoid, but shaped by both the local density and the directional covariance of the calibration residuals. This decomposition mirrors a familiar progression in Rd . The squared error ∥ε − ε̄∥2 treats all directions equally, and the Mahalanobis distance adds inverse-covariance weighting to penalize deviations along low-variance directions. The MMD term is the RKHS analogue of the squared error, an isotropic distance from the mean embedding that ignores directional structure. The KPCA correction adds inverse-covariance weighting in H, completing the lift from Mahalanobis in Rd to an anisotropic distance in the feature space. The parameter γ controls how much of the directional structure the score uses. Large γ suppresses the correction and recovers isotropic density-only scoring, while small γ lets the full anisotropic structure shape the prediction region. 4.3
Special Cases and Properties
Proposition 4.3 (Non-negativity). The MKS satisfies êk (ε) ≥ 0 for all ε ∈ Rd , with equality if and only if ϕ̃(ε) lies in the span of the calibration data in H and is perfectly interpolated by (K̃ + γIT )−1 . Proof. Immediate from the GP identity (Proposition 4.1): the score equals the posterior variance of a Gaussian process, which is non-negative by construction. Proposition 4.4 (Linear Kernel Recovery). For the linear kernel k(x, x′ ) = x⊤ x′ , the kernel nonconformity score equals λ · MD(ε), where MD(ε) is the regularized Mahalanobis distance of Definition 3.1. Since the constant factor λ does not affect ranking, the two scores produce identical conformal prediction regions. The proof is in Appendix A.4. The Multivariate Kernel Score strictly generalizes the Mahalanobis approach [16, 40], as the linear kernel recovers their score exactly, while nonlinear kernels extend the method beyond ellipsoidal regions. In the linear case, the regularization parameter λ corresponds to ridge regularization of the sample covariance, and the calibration scores reduce to the diagonal of the ridge hat matrix X(X ⊤ X + λI)−1 X ⊤ . The Multivariate Kernel Score admits a decomposition into a component within the span of the calibration data in H and a component orthogonal to it. Proposition 4.5 (Spectral Decomposition). Let {λj , vj }Tj=1 be the eigenvalues and eigenvectors p of K̃, and let ej = wj /∥wj ∥H = wj / λj be the normalized kernel principal components in H, PT PT where wj = i=1 [vj ]i ϕ̃(εi ). Define the orthogonal residual ϕ̃⊥ (ε) = ϕ̃(ε) − j=1 ⟨ϕ̃(ε), ej ⟩H ej . Then T X γ · αj (ε)2 + ∥ϕ̃⊥ (ε)∥2H (12) êk (ε) = λ (λ + γ) j j j=1 P The proof is in Appendix A.5. The within-span component j γ αj (ε)2 /[λj (λj + γ)] measures atypicality along the principal directions of the calibration data. The calibration directions are inversely weighted by their variance, consistent with the Mahalanobis interpretation. The orthogonal component ∥ϕ̃⊥ (ε)∥2H captures the part of the test residual that lies entirely outside the span of the calibration data in H. It dominates for residuals in regions of the output space not explored by the calibration set, where the data provides no structural information and the GP posterior variance reduces to the prior variance. 6
The parameter γ controls the balance between these two components. When γ → 0, the withinspan weights vanish and the score trusts the calibration structure completely and penalizes only the orthogonal component. When γ → ∞, the weights approach 1/λj , the KPCA correction in (11) vanishes, and the score reduces to MMD2 (δε , QT ), recovering isotropic density-only scoring. For calibration residuals, ϕ̃(εi ) lies in the span of the calibration data by construction, so the orthogonal component vanishes and the score admits a closed-form expression in terms of the kernel ridge regression hat matrix. Proposition 4.6 (Leverage Score Identity). Let H = K̃(K̃ + γIT )−1 be the kernel ridge regression hat matrix, with diagonal entries hii = Hii . For each calibration residual εi , the Multivariate Kernel Score satisfies êk (εi ) = γ · hii (13) The proof is in Appendix A.6. The diagonal entries hii are the γ-ridge leverage scores of the centered Gram matrix [1], which measure the influence of each calibration point on the kernel ridge regression fit. The conformal threshold is a quantile q̂ directly determined by the ordering of the leverages {hii }Ti=1 , where high-leverage calibration residuals occupying unusual positions in feature space receive high nonconformity scores, and low-leverage residuals in well-populated regions receive low scores. For a test point ε ∈ / {εi }, the orthogonal component ∥ϕ̃⊥ (ε)∥2H is generally nonzero, and the MKS can be understood as an out-of-sample extension of leverage, penalizing test residuals that lie outside the span of the calibration data in H. High-dimensional geometry. The concentration of norm for sub-Gaussian random vectors [36, Theorem 3.1.1] implies that even √ spherical Gaussian data in d dimensions concentrates on a thin shell of width O(1) around radius d. The highest-density region is therefore an annulus rather than a ball, and any convex prediction region, including ellipsoidal, must enclose the low-density interior at the cost of excess volume that grows with the dimension d. The inefficiency persists even under correct model specification with Gaussian residuals. The MKS, through the density sensitivity of the MMD term (Proposition 4.2), assigns high scores to interior points unsupported by calibration data, avoiding this cost.
5
Experiments
In all experiments, we use the RBF kernel k(ε, ε′ ) = exp(−∥ε − ε′ ∥2 /2ℓ2 ) for demonstration purposes. The primary contribution of this work is conceptual and theoretical; the choice of kernel and its optimization are important orthogonal questions that we defer to future work. The Multivariate Kernel Score introduces two hyperparameters, the lengthscale ℓ and the regularisation γ. Both are in principle tunable via cross-validation on the calibration scores or through the Gaussian process marginal likelihood implied by Proposition 4.1. We use simple heuristics throughout the experiments: γ is set to the 90th percentile eigenvalue of the centered Gram matrix K̃ (see Appendix C.1 for volume estimation details), and ℓ is specified per experiment below. We compare against the Mahalanobis score (Definition 3.1), the Bonferroni method (univariate conformal prediction per coordinate at level α/d), and a Density-only score using MMD2 (δε , QT ) alone (the γ → ∞ limit of our score). Sensitivity to the lengthscale is analyzed in Appendix 3, where the kernel score outperforms the density-only baseline at every value tested. Both the kernel score and the density-only score share the same lengthscale, so differences are attributable to the KPCA correction term alone. For each dataset, we fit a linear model and a Multi-Layer Perceptron (MLP) (Appendix C.3.2) to obtain the residuals and apply all scores to the resulting residuals. 5.1
Synthetic Data
We generate n = 50,000 samples from a bivariate regression task (details in Appendix C.3.1) with non-Gaussian, heteroscedastic noise. We set the kernel lengthscale to l = 0.5. Sensitivity to this choice is analyzed in Appendix 3. We fit both a linear regression and a neural network (2-layer MLP), and apply split conformal prediction with a 50/25/25 train/calibration/test split. Results are averaged over 100 random seeds. Figure 1 displays the prediction regions for both models across four coverage levels. Under the linear model (Figure 1a), the residuals exhibit clear non-elliptical structure due to model misspecification. 7
Linear α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.90016±0.00381 0.90001±0.00394 0.90005±0.00363 0.89886±0.00382
4.13745±0.03328 2.82308±0.02818 3.34529±0.02959 2.48489±0.02346
0.71391±0.01359 0.77994±0.01252 0.66097±0.01517 0.82092±0.01042
0.90007±0.00386 0.90034±0.00362 0.89977±0.00444 0.90052±0.00386
0.51110±0.00833 0.50315±0.00905 0.74175±0.01487 0.45890±0.00774
0.88635±0.00669 0.88902±0.00568 0.88919±0.00623 0.88533±0.00727
0.05
Bonferroni Mahalanobis Density MKS
0.94997±0.00275 0.95007±0.00281 0.95022±0.00285 0.94995±0.00293
4.97647±0.04020 3.54716±0.04695 4.09317±0.03996 3.05467±0.03167
0.85863±0.01030 0.85084±0.01021 0.83089±0.01119 0.90564±0.00796
0.95030±0.00277 0.95036±0.00255 0.95046±0.00281 0.95011±0.00267
0.69037±0.01230 0.71884±0.01362 1.01563±0.02712 0.60802±0.01079
0.94211±0.00431 0.94241±0.00384 0.94337±0.00388 0.94026±0.00456
0.02
Bonferroni Mahalanobis Density MKS
0.97996±0.00163 0.97994±0.00172 0.98014±0.00173 0.97974±0.00169
5.97821±0.06445 4.88194±0.08387 4.96063±0.06658 3.80359±0.04239
0.94382±0.00557 0.92729±0.00743 0.93530±0.00670 0.95963±0.00452
0.98013±0.00178 0.98009±0.00170 0.98018±0.00185 0.98002±0.00179
0.95647±0.02194 1.10661±0.03097 1.38108±0.04531 0.82686±0.01830
0.97543±0.00253 0.97541±0.00260 0.97573±0.00278 0.97490±0.00260
0.01
Bonferroni Mahalanobis Density MKS
0.99008±0.00127 0.99017±0.00105 0.99014±0.00115 0.99002±0.00118
6.71683±0.09442 6.07201±0.13758 5.55893±0.08743 4.41304±0.07076
0.97186±0.00408 0.96284±0.00470 0.96936±0.00439 0.97954±0.00323
0.98998±0.00131 0.98999±0.00127 0.99010±0.00117 0.98999±0.00124
1.19619±0.03910 1.47632±0.05331 1.65518±0.06259 1.02031±0.02455
0.98672±0.00208 0.98673±0.00178 0.98701±0.00180 0.98679±0.00180
2
Table 1: Results on synthetic data across 100 seeds. All values reported as mean ± std. Bold indicates best volume per model and α. Volume is estimated by Monte Carlo (Appendix C.1). WSC denotes worst-slab coverage (Appendix C.2), an approximation of conditional coverage; higher is better. 2.0
1.00
1.5
0.75
1.0
0.50
0.5
0.25
0.0
0.00
0.5
0.25 0.50
1.0
0.75
1.5 2
1
0
1
2
2
1
1
0
1
2
1
1.00
1.0
0.5
0.0
0.5
1.0
1.0
1
(a) Linear model residuals
0.5
0.0
0.5
1.0
1
(b) MLP residuals
Figure 1: Prediction regions on synthetic data at coverage levels 1 − α ∈ {0.90, 0.95, 0.98, 0.99}. In each panel, left: Mahalanobis (blue); right: MKS (red), overlaid on calibration residuals The ellipsoids must expand into low-density regions to achieve the target coverage, wasting volume on regions where residuals are unlikely to occur. The MKS contour follows the geometry of the calibration data, producing tighter regions that exclude these areas. Under the neural network (Figure 1b), the residuals are closer to elliptical but the MKS retains a volume advantage. Table 1 reports coverage, volume, and worst-slab coverage (WSC) across all methods, models, and coverage levels . The Multivariate Kernel Score achieves the smallest volume at every α level under both models. While maintaining nominal marginal coverage and approximate conditional coverage as measured by WSC1 (details Appendix C.2), the kernel score produces prediction regions whose volumes are 9 to 31% smaller than Mahalanobis and 21 to 38% smaller than the density-only score. The advantage grows at tighter coverage levels. As the target coverage increases, the ellipsoid must extend into the tails of the residual distribution where the calibration data provides little support, and the volume cost of this expansion grows faster for the ellipsoid than for the MKS contour, which expands along the support of the data. The effect is most pronounced under the misspecified linear model, where the residual geometry is far from elliptical, but the MKS still improves upon Mahalanobis under the neural network, confirming that it adapts to the residual geometry without penalty when the distribution is closer to elliptical. The worst-slab coverage supports this reading. Under the linear model, the MKS achieves the best WSC at every α level, so the volume reduction reflects a genuinely tighter fit rather than a sacrifice in conditional coverage. Under the MLP, WSC values are comparable across all presented methods, which is consistent with the more symmetric geometry and leaving less room for improvement. 1WSC searches over linear slabs, structurally favoring methods with whose regions align with linear subspaces, e.g. Mahalanobis ellipsoid and isotropic density level sets. The anisotropic KPCA term can produce boundaries less aligned with linear slabs even when the region is tighter. Comparable WSC at smaller volume suggests substantial coverage improvement.
8
Volume Ratio (MKS / Mahalanobis)
MKS Coverage
Dataset
d
1−α=.90
.95
.98
.99
1−α=.90
.95
.98
.99
House
2 3 4
0.88655±0.017 0.95057±0.018 0.60678±0.057
0.91529±0.017 0.95180±0.023 0.59417±0.048
0.84290±0.035 0.73051±0.038 0.41423±0.054
0.68921±0.040 0.59717±0.044 0.26740±0.057
0.89979±0.00602 0.89991±0.00567 0.90038±0.00631
0.94981±0.00421 0.94960±0.00434 0.94959±0.00407
0.97976±0.00215 0.97934±0.00236 0.97894±0.00249
0.99006±0.00165 0.98962±0.00194 0.98868±0.00193
Bio
2 3 4
0.51606±0.096 0.49689±0.099 0.53949±0.095
0.53702±0.099 0.52137±0.100 0.52587±0.090
0.58862±0.103 0.45920±0.082 0.47425±0.080
0.68670±0.106 0.49004±0.091 0.41697±0.044
0.89985±0.00490 0.89948±0.00386 0.89964±0.00428
0.94992±0.00310 0.94977±0.00289 0.94953±0.00315
0.97992±0.00170 0.97971±0.00199 0.97959±0.00182
0.98964±0.00142 0.98964±0.00145 0.98959±0.00131
Blog
2 3 4
0.48138±0.038 0.59225±0.046 0.30738±0.094
0.63779±0.042 0.65759±0.043 0.36040±0.106
0.77891±0.050 0.51207±0.032 0.26011±0.076
0.59402±0.044 0.33584±0.034 0.14225±0.041
0.89959±0.00393 0.89933±0.00348 0.89901±0.00338
0.94961±0.00281 0.94975±0.00252 0.94956±0.00226
0.97970±0.00164 0.97972±0.00188 0.97962±0.00189
0.98958±0.00143 0.98923±0.00144 0.98941±0.00138
Table 2: Real data results under Ridge regression residuals, averaged over 50 seeds. Left: volume ratio of the MKS score to Mahalanobis (values below 1 indicate smaller MKS regions). Right: empirical coverage of the MKS score (target is 1−α).
The density-only score consistently produces larger regions than the Mahalanobis ellipsoid under both models, demonstrating that the improvement of the MKS is not due to density estimation alone but to the covariance adjustment provided by the KPCA correction. A comparison between our MKSand the density score across a range of lengthscale parameters is provided in Table 3 of Appendix C.3.3, showing that our method outperforms regardless of the value of the parameter.
5.2
Real Data
We evaluate the method on three regression datasets previously used in Feldman et al. [10] and Romano et al. [27], where we augment the original univariate responses with additional variables from the same source, producing response dimensions d = 2, 3, and 4. For real data experiments across varying output dimensions, we use an automatic lengthscale selection2 . The lengthscale scales linearly with d to compensate for the concentration of pairwise distances in higher-dimensional residual spaces, ensuring sufficient smoothing for robust coverage across dimensions. All details are provided in Appendix C.4.1. We report results from a Linear model forecast residual. MLP regression results are provided in Appendix C.4.3. Table 2 reports the volume ratio of regions produced using the Multivariate Kernel Score relative to the Mahalanobis baseline, computed over 50 random splits. The MKS produces smaller prediction regions across all datasets and dimensions, with the advantage increasing both with the output dimension d and the coverage level 1 − α. As the number of output coordinates grows, correlations between outputs create residual distributions that an ellipsoid cannot approximate without enclosing large empty regions, and the MKS exploits this structure. The effect is compounded at tighter coverage levels, matching the pattern observed on synthetic data in Section 5.1. Table 2 also confirms that coverage is maintained across all configurations. The empirical coverage of the MKS tracks the nominal alpha within sampling variability at every dataset, dimension and coverage level. Coverage values are comparable across all methods, consistent with the distribution free guarantee of the conformal framework. Moreover, all methods achieve comparable approximate conditional coverage in higher dimensions and tighter coverage levels, as measured empirically through the WSC (Appendix C.4.3). Full results including Bonferroni and density-based non-conformity scores are reported in Appendix C.4.3. The Bonferroni method, ignoring correlations and yielding hyper-rectangles, produces the largest volumes across all settings. The density-only score performs comparably to Mahalanobis, with relative performance depending on the residual geometry of each dataset. The MKS , which adds the KPCA correction to the same density estimate, consistently matches or improves upon both, with the covariance adjustment becoming increasingly important as the output dimension grows and the residual distribution departs further from isotropy. Overall, the Multivariate Kernel Score maintains nominal coverage while reducing region volumes by 5 to 86% relative to the Mahalanobis baseline, with the largest gains at higher dimensions and tighter coverage levels.
2We set ℓ = d ∗ median{∥ε
i − εj ∥}/2, where d is the output dimension and median takes values over all pairs of i ̸= j.
9
6
Conclusion
We introduced the Multivariate Kernel Score (MKS) for multivariate conformal prediction that generalizes the Mahalanobis distance to reproducing kernel Hilbert spaces, producing prediction regions adapted to the residual distribution without imposing a parametric shape. The proposed score can be perceived via two distinct perspectives: (1) the Gaussian process posterior variance; (2) the decomposition into a density term and a directional covariance correction term. The two-sided interpretation of MKS allows us to connect Bayesian uncertainty quantification with distribution-free coverage guarantees of frequentist-type. By compressing multivariate residuals into scalars through the RKHS geometry, the MKS defines a nonparametric center-outward ordering on Rd , with the conformal threshold selecting the quantile level. Finite-sample convergence rates depend on the effective rank of the RKHS covariance operator rather than the ambient dimension. On real data, the MKS significantly reduces prediction region volumes relative to the Mahalanobis baseline while maintaining nominal coverage, with the largest gains at higher dimensions and tighter coverage levels. Future research directions include adaptation to non-exchangeable settings such as incorporating covariate shift, and conformal prediction for function-valued responses.
Acknowledgments and Disclosure of Funding This work was supported by G-Research through the NextGen Scholarship.
References [1] Ahmed El Alaoui and Michael W. Mahoney. Fast randomized kernel methods with statistical guarantees, 2015. [2] Anastasios N. Angelopoulos and Stephen Bates. A gentle introduction to conformal prediction and distribution-free uncertainty quantification, 2022. [3] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. Conformal prediction beyond exchangeability. The Annals of Statistics, 51(2):816–845, 2023. [4] Sacha Braun, Liviu Aolaritei, Michael I Jordan, and Francis Bach. Minimum volume conformal sets for multivariate regression. arXiv preprint arXiv:2503.19068, 2025. [5] Sacha Braun, Eugène Berta, Michael I Jordan, and Francis Bach. Multivariate conformal prediction via conformalized gaussian scoring. arXiv preprint arXiv:2507.20941, 2025. [6] Krisztian Buza. BlogFeedback. https://doi.org/10.24432/C58S3F.
UCI Machine Learning Repository, 2014.
DOI:
[7] Maxime Cauchois, Suyash Gupta, and John C Duchi. Knowing what you know: valid and validated confidence sets in multiclass and multilabel prediction. Journal of machine learning research, 22(81):1–42, 2021. [8] Jase Clarkson, Wenkai Xu, Mihai Cucuringu, Yvik Swan, and Gesine Reinert. Split conformal prediction under data contamination. arXiv preprint arXiv:2407.07700, 2024. [9] Victor Dheur, Matteo Fontana, Yorick Estievenart, Naomi Desobry, and Souhaib Ben Taieb. A unified comparative study with generalized conformity scores for multi-output conformal regression. arXiv preprint arXiv:2501.10533, 2025. [10] Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated multiple-output quantile regression with representation learning, 2022. [11] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. The limits of distribution-free conditional predictive inference. Information and Inference: A Journal of the IMA, 10(2):455–482, 2021. [12] Kenji Fukumizu, Francis R Bach, and Michael I Jordan. Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces. Journal of Machine Learning Research, 5:73– 99, 2004. 10
[13] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. The journal of machine learning research, 13(1):723–773, 2012. [14] house. House sales in king county, USA. https://www.kaggle.com/harlfoxem/ housesalesprediction/metadata. Accessed: July, 2021. [15] Rafael Izbicki, Gilson Shimizu, and Rafael B Stern. Cd-split and hpd-split: Efficient conformal regions in high dimensions. Journal of Machine Learning Research, 23(87):1–32, 2022. [16] Chancellor Johnstone and Bruce Cox. Conformal uncertainty sets for robust optimization. In Conformal and Probabilistic Prediction and Applications, pages 72–90. PMLR, 2021. [17] Michal Klein, Louis Bethune, Eugene Ndiaye, and Marco Cuturi. Multivariate conformal prediction using optimal transport. arXiv preprint arXiv:2502.03609, 2025. [18] Vladimir Koltchinskii and Karim Lounici. Concentration inequalities and moment bounds for sample covariance operators. Bernoulli, pages 110–133, 2017. [19] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman. Distribution-free predictive inference for regression. Journal of the American Statistical Association, 113(523):1094–1111, 2018. [20] Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Copula-based conformal prediction for multi-target regression. Pattern Recognition, 120:108101, 2021. [21] Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Ellipsoidal conformal inference for multi-target regression. In Conformal and Probabilistic Prediction with Applications, pages 294–306. PMLR, 2022. [22] Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In European conference on machine learning, pages 345–356. Springer, 2002. [23] Emanuel Parzen. On estimation of a probability density function and mode. The annals of mathematical statistics, 33(3):1065–1076, 1962. [24] Vincent Plassier, Alexander Fishkov, Mohsen Guizani, Maxim Panov, and Eric Moulines. Probabilistic conformal prediction with approximate conditional validity. arXiv preprint arXiv:2407.01794, 2024. [25] Prashant Rana. Physicochemical Properties of Protein Tertiary Structure. UCI Machine Learning Repository, 2013. DOI: https://doi.org/10.24432/C5QW3H. [26] Emmanuel Rio. Asymptotic Theory of Weakly Dependent Random Processes. Springer, 2017. [27] Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019. [28] Murray Rosenblatt. Remarks on Some Nonparametric Estimates of a Density Function. The Annals of Mathematical Statistics, 27(3):832–837, 1956. [29] Max Sampson and Kung-Sik Chan. Flexible conformal highest predictive conditional density sets, 2025. [30] Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller. Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5):1299–1319, 1998. [31] Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of machine learning research, 9(3), 2008. [32] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert RG Lanckriet. Hilbert space embeddings and metrics on probability measures. The Journal of Machine Learning Research, 11:1517–1561, 2010. 11
[33] Sophia Huiwen Sun and Rose Yu. Copula conformal prediction for multi-step time series prediction. In The Twelfth International Conference on Learning Representations, 2024. [34] Gauthier Thurin, Kimia Nadjahi, and Claire Boyer. Optimal transport-based conformal prediction. arXiv preprint arXiv:2501.18991, 2025. [35] Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas. Conformal prediction under covariate shift. Advances in neural information processing systems, 32, 2019. [36] Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018. [37] Vladimir Vovk. Conditional validity of inductive conformal predictors. In Asian conference on machine learning, pages 475–490. PMLR, 2012. [38] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in a random world. Springer, 2005. [39] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006. [40] Chen Xu, Hanyang Jiang, and Yao Xie. Conformal prediction for multi-dimensional time series by ellipsoidal sets, 2024. [41] Chen Xu and Yao Xie. Conformal prediction for time series. IEEE transactions on pattern analysis and machine intelligence, 45(10):11575–11587, 2023. [42] Margaux Zaffran, Aymeric Dieuleveut, Julie Josse, and Yaniv Romano. Conformal prediction with missing values. In International Conference on Machine Learning, pages 40578–40604, 2023.
12
A
Derivations and Proofs
A.1
Derivation of Proposition 3.3 (Computable Form)
We show that the RKHS Mahalanobis distance in Definition 3.2 reduces to the finite-dimensional expression in Proposition 3.3 via the Woodbury matrix identity. Proof. Let Φ = [ϕ(ε1 ), . . . , ϕ(εT )] denote the collection of feature vectors in H, and let H = IT − T1 11⊤ ∈ RT ×T be the centering matrix. The centered feature vectors are Φ̃ = ΦH, and the bϕ = 1 Φ̃Φ̃⊤ . sample covariance operator is C T −1 We apply the Woodbury matrix identity. For A = λIH , U = √T1−1 Φ̃, and C = IT , the identity (A + U U ⊤ )−1 = A−1 − A−1 U (I + U ⊤ A−1 U )−1 U ⊤ A−1 gives:
bϕ + λIH C
−1
−1 1 1 1 ⊤ = IH − 2 Φ̃ IT + Φ̃ Φ̃ Φ̃⊤ λ λ (T − 1) λ(T − 1)
(14)
The centered Gram matrix is K̃ = Φ̃⊤ Φ̃ = HKH ∈ RT ×T , where Kij = k(εi , εj ). Substituting into (14): −1 −1 1 1 1 bϕ + λIH Φ̃⊤ (15) Φ̃ IT + K̃ C = IH − 2 λ λ (T − 1) λ(T − 1) Let ϕ̃(ε) = ϕ(ε) − µ̂ϕ be the centered feature map of a test residual. Applying Definition 3.2: −1 bϕ + λIH êk (ε) = ϕ̃(ε), C ϕ̃(ε) H
−1 1 1 1 ⊤ = k̃(ε, ε) − 2 k̃∗ (ε) IT + K̃ k̃∗ (ε) λ λ (T − 1) λ(T − 1)
(16)
where k̃∗ (ε) = Φ̃⊤ ϕ̃(ε) ∈ RT is the centered kernel vector with entries [k̃∗ (ε)]i = ⟨ϕ̃(ε), ϕ̃(εi )⟩H . Setting γ = λ(T − 1), we have IT + γ1 K̃ = γ1 (K̃ + γIT ), so (IT + γ1 K̃)−1 = γ(K̃ + γIT )−1 . Substituting into (16): 1 γ k̃(ε, ε) − 2 k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) λ λ (T − 1) i 1h k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) = λ
êk (ε) =
where the second line uses γ/[λ(T − 1)] = 1. Since the factor 1/λ is a positive constant that does not affect the ranking of residuals, we define the Multivariate Kernel Score as in Proposition 3.3: êk (ε) = k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) For completeness, the centered quantities in terms of the raw kernel are: T T 2X 1 X k̃(ε, ε) = k(ε, ε) − k(ε, εi ) + 2 k(εi , εj ), T i=1 T i,j=1 T T T 1X 1 X 1X k(ε, εj ) − k(εj , εi ) + 2 k(εj , εℓ ) [k̃∗ (ε)]i = k(ε, εi ) − T j=1 T j=1 T j,ℓ=1
13
A.2
Proof of Proposition 4.1 (Gaussian Process Identity)
Proof. Let f ∼ GP(0, k̃) and suppose we observe f at locations ε1 , . . . , εT with i.i.d. Gaussian noise of variance γ. By standard GP conditioning [39, Eq. 2.26], the posterior variance at a test point ε is VarGP f (ε) | {εi }Ti=1 = k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) where K̃ij = k̃(εi , εj ) and [k̃∗ (ε)]i = k̃(ε, εi ). This is identical to the expression in Proposition 3.3. Note that the posterior variance depends only on the observation locations and the kernel, not on the observed function values. The identity êk (ε) = VarGP [f (ε) | {εi }Ti=1 ] follows by direct comparison. A.3
Proof of Proposition 4.2 (MMD–KPCA Decomposition)
Proof. We establish two intermediate steps and combine them for the result. The MMD term. The squared maximum mean discrepancy between two distributions P and Q with respect to a kernel k is [13] MMD2 (P, Q) = Ex,x′ ∼P k(x, x′ ) + Ey,y′ ∼Q k(y, y ′ ) − 2 Ex∼P, y∼Q k(x, y) PT For P = δε (the Dirac measure at ε) and Q = QT = T1 i=1 δεi (the empirical measure of the calibration residuals), each expectation evaluates as follows. Since x = x′ = ε almost surely under δε : ′ ′ ′ ′ ∼δ [k(x, x )] = k(ε, ε). Since y, y are drawn independently from QT : Ey,y ′ ∼Q [k(y, y )] = Ex,xP ε T P 1 1 i,j k(εi , εj ) = k̄. The cross term: Ex∼δε ,y∼QT [k(x, y)] = T i k(ε, εi ) = p̂k (ε). Therefore T2 MMD2 (δε , QT ) = k(ε, ε) + k̄ − 2p̂k (ε) Equivalently, MMD2 (P, Q) = ∥µP − µQ ∥2H where µP = Ex∼P [ϕ(x)] is the kernel mean embedding [13]. Since µδε = ϕ(ε) and µQT = µ̂ϕ , we have MMD2 (δε , QT ) = ∥ϕ(ε) − µ̂ϕ ∥2H = ∥ϕ̃(ε)∥2H = k̃(ε, ε) The KPCA correction term. Let K̃ = V ΛV ⊤ be the eigendecomposition, with V = [v1 , . . . , vT ] orthogonal and Λ = diag(λ1 , . . . , λT ). Then (K̃ + γIT )−1 = V (Λ + γIT )−1 V ⊤ =
T X
1 vj vj⊤ λ + γ j=1 j
and therefore k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) =
T X (vj⊤ k̃∗ (ε))2 j=1
λj + γ
=
T X αj (ε)2 j=1
λj + γ
where αj (ε) = vj⊤ k̃∗ (ε). PT Verification that αj (ε) is a KPCA projection. Define the RKHS vectors wj = i=1 [vj ]i ϕ̃(εi ) ∈ bϕ with eigenvalues λj /(T − 1): H. These are eigenfunctions of the sample covariance operator C * + T T X X 1 bϕ wj = C ϕ̃(εℓ ) ϕ̃(εℓ ), [vj ]i ϕ̃(εi ) T −1 i=1 ℓ=1
=
1 T −1
T X ℓ=1
H
T X ϕ̃(εℓ ) [vj ]i K̃ℓi = i=1
T
=
1 T −1
λj λj X [vj ]ℓ ϕ̃(εℓ ) = wj T −1 T −1 ℓ=1
14
T X ℓ=1
ϕ̃(εℓ ) [K̃vj ]ℓ
The projection of ϕ̃(ε) onto wj is ⟨ϕ̃(ε), wj ⟩H =
T X
[vj ]i ⟨ϕ̃(ε), ϕ̃(εi )⟩H = vj⊤ k̃∗ (ε) = αj (ε)
i=1
confirming that αj (ε) is the projection of ϕ̃(ε) onto the j-th kernel principal component. Combining.
From Proposition 3.3:
êk (ε) = k̃(ε, ε) − k̃∗ (ε)⊤ (K̃ + γIT )−1 k̃∗ (ε) = MMD2 (δε , QT ) −
T X αj (ε)2 j=1
A.4
λj + γ
Proof of Proposition 4.4 (Linear Kernel Recovery)
Proof. For the linear kernel k(x, x′ ) = x⊤ x′ , the feature map is ϕ(x) = x, so µ̂ϕ = ε̄, ϕ̃(εi ) = εi −ε̄, and K̃ij = (εi − ε̄)⊤ (εj − ε̄). Let Ẽ = [ε1 − ε̄, . . . , εT − ε̄] ∈ Rd×T . Then K̃ = Ẽ ⊤ Ẽ, k̃∗ (ε) = Ẽ ⊤ (ε− ε̄), and k̃(ε, ε) = ∥ε− ε̄∥2 . The score in Proposition 3.3 becomes: êk (ε) = ∥ε − ε̄∥2 − (ε − ε̄)⊤ Ẽ (Ẽ ⊤ Ẽ + γIT )−1 Ẽ ⊤ (ε − ε̄) h i = (ε − ε̄)⊤ Id − Ẽ (Ẽ ⊤ Ẽ + γIT )−1 Ẽ ⊤ (ε − ε̄)
(17)
By the push-through identity, Ẽ (Ẽ ⊤ Ẽ + γIT )−1 Ẽ ⊤ = (Ẽ Ẽ ⊤ + γId )−1 Ẽ Ẽ ⊤ . Substituting into (17): h i êk (ε) = (ε − ε̄)⊤ Id − (Ẽ Ẽ ⊤ + γId )−1 Ẽ Ẽ ⊤ (ε − ε̄) h i = (ε − ε̄)⊤ (Ẽ Ẽ ⊤ + γId )−1 (Ẽ Ẽ ⊤ + γId ) − Ẽ Ẽ ⊤ (ε − ε̄) = γ · (ε − ε̄)⊤ (Ẽ Ẽ ⊤ + γId )−1 (ε − ε̄) b and γ = λ(T − 1): Since Ẽ Ẽ ⊤ = (T − 1)Σ b + λId ) −1 (ε − ε̄) êk (ε) = λ(T − 1) · (ε − ε̄)⊤ (T − 1)(Σ b + λId )−1 (ε − ε̄) = λ · MD(ε) = λ · (ε − ε̄)⊤ (Σ The constant factor λ > 0 does not affect the ranking of residuals, so the MKS and the regularized Mahalanobis distance produce identical conformal prediction regions.
A.5
Proof of Proposition 4.5 (Spectral Decomposition)
Proof. From the proof of Proposition 4.2, the projection of ϕ̃(ε) onto wj satisfies ⟨ϕ̃(ε), wj ⟩H = αj (ε), so the projection onto the normalized direction is αj (ε) ⟨ϕ̃(ε), ej ⟩H = p λj Since the ej are orthonormal and ϕ̃⊥ (ε) is orthogonal to all of them, Pythagoras’ theorem gives ∥ϕ̃(ε)∥2H =
T X
⟨ϕ̃(ε), ej ⟩2H + ∥ϕ̃⊥ (ε)∥2H =
j=1
T X αj (ε)2 j=1
15
λj
+ ∥ϕ̃⊥ (ε)∥2H
Substituting into the MMD–KPCA decomposition (Proposition 4.2), with ∥ϕ̃(ε)∥2H MMD2 (δε , QT ): êk (ε) = ∥ϕ̃(ε)∥2H −
T X αj (ε)2 j=1
=
T X αj (ε)2 j=1
=
T X
A.6
T X αj (ε)2 j=1
j=1
=
λj + γ
+ ∥ϕ̃⊥ (ε)∥2H −
λj αj (ε)2
=
λj + γ
1 1 − + ∥ϕ̃⊥ (ε)∥2H λj λj + γ
T X γ · αj (ε)2
λ (λj + γ) j=1 j
+ ∥ϕ̃⊥ (ε)∥2H
Proof of Proposition 4.6 (Leverage Score Identity)
Proof. For a calibration residual εi , ϕ̃(εi ) lies in the span of the calibration data in H by construction, so ∥ϕ̃⊥ (εi )∥2H = 0 and the spectral decomposition (Proposition 4.5) reduces to êk (εi ) =
T X γ αj (εi )2 j=1
λj (λj + γ)
.
Since k̃∗ (εi ) = K̃:,i and K̃vj = λj vj , the KPCA projection at a calibration point is αj (εi ) = vj⊤ K̃:,i = λj [vj ]i , Substituting: T X γ λj [vj ]2 i
êk (εi ) =
j=1
λj + γ
=γ
T X
λj [vj ]2i = γ · hii λ + γ j j=1
where hii = [V Λ(Λ+γIT )−1 V ⊤ ]ii is the i-th diagonal entry of the matrix H = K̃(K̃ +γIT )−1 .
B
Theoretical Analysis
This section provides the full theoretical analysis for the coverage guarantees stated in Theorem 3.4. We first state the assumptions and auxiliary results, then prove asymptotic conditional coverage under i.i.d. residuals (Section B.2), where the standard exogeneity condition makes the conditional guarantee valid, and finally establish marginal coverage under α-mixing dependence (Section B.3), where temporal correlation between XT +1 and εT +1 precludes a conditional guarantee. Throughout this section, we work at the RKHS operator level. The true nonconformity score at time t is et = ⟨ϕ(εt ) − µ, Σ−1 λ (ϕ(εt ) − µ)⟩H and the estimated nonconformity score is b −1 (ϕ(ε̂t ) − µ̂)⟩H êt = ⟨ϕ(ε̂t ) − µ̂, Σ λ bλ = Σ b + λIH are the regularized true and empirical cross-covariance where Σλ = Σ + λIH and Σ operators. Here µ and Σ are the population mean embedding and covariance operator (fixed), while b are estimated from the calibration set. The estimated score êt corresponds to êk (ε̂t ) in µ b and Σ Proposition 3.3 after the Woodbury reduction (Appendix A.1). We also define the empirical CDFs of the true and estimated scores: T 1X e FT +1 (x) = 1{et ≤x} T t=1
T 1X b FT +1 (x) = 1{êt ≤x} . T t=1
16
B.1
Assumptions
+1 Assumption B.1 (i.i.d. Residuals and Exogeneity). The residuals {εt }Tt=1 are independent and identically distributed. Furthermore, the CDF Fe of the true nonconformity score et is Lipschitz continuous with constant L > 0. Finally, the test residual εT +1 is independent of the test input XT +1 .
Remark. The Lipschitz condition on Fe ensures that small perturbations of the scores lead to small changes in coverage probability. This assumption is relaxed in Section B.3, where we extend the results to stationary α-mixing sequences. The exogeneity condition εT +1 ⊥ XT +1 is standard in regression and ensures that the true nonconformity score eT +1 is independent of XT +1 , which allows the marginal bound in the proof to hold conditionally on any fixed XT +1 = x. Under α-mixing dependence (Section B.3), XT +1 typically contains lags of Yt and is correlated with εT +1 , so this condition fails and only marginal coverage is established there. Assumption B.2 (Estimation Quality in RKHS). There exists a sequence {δT }T ≥1 satisfying δT = O(T −ρ ) for some ρ > 0 such that: T
1X ∥ϕ(ε̂t ) − ϕ(εt )∥2H ≤ δT2 T t=1
and
∥ϕ(ε̂T +1 ) − ϕ(εT +1 )∥H ≤ δT .
Remark. This requires the prediction error in the RKHS feature space to be small on average over the calibration set and pointwise at the test point. A sufficient condition is that the kernel satisfies ∥ϕ(x) − ϕ(y)∥H ≤ Lϕ ∥x − y∥ for some Lϕ > 0, in which case Assumption B.2 holds whenever T
1X ∥ε̂t − εt ∥2 ≤ (δT /Lϕ )2 . T t=1 For the RBF kernel, the feature map is not globally Lipschitz, but the condition can be verified under additional regularity on the residual distribution (e.g., sub-Gaussian εt ). The polynomial rate ρ at which the error vanishes depends on the estimation method and the problem structure. Faster estimators (larger ρ) yield tighter coverage guarantees. Assumption B.3 (Boundedness). The kernel k is bounded: there exists B > 0 such that k(x, x) ≤ B 2 for all x in the domain. Remark. This is equivalent to ∥ϕ(x)∥H ≤ B for all x. It ensures that all moments of ϕ(ε) exist (E[∥ϕ(ε)∥q ] ≤ B q ), which is required for the concentration inequalities used below. Assumption B.4 (Trace-Class Covariance). The population covariance operator Σ = E[(ϕ(ε) − µ) ⊗ (ϕ(ε) − µ)] is trace-class: tr(Σ) < ∞. B.2
Coverage Guarantee and Convergence Rate (i.i.d.)
We prove the coverage guarantee through a sequence of lemmas, each bounding one source of error. These are combined in the proof of Theorem B.12. B.2.1
Auxiliary Results
We collect standard results used throughout the proofs. Lemma B.5 (Cauchy–Schwarz for Self-Adjoint Operators). For any self-adjoint operator A on H and any x ∈ H, |⟨x, Ax⟩H | ≤ ∥x∥2H ∥A∥op . Lemma B.6 (Resolvent Identity). For invertible operators A and B on H, A−1 − (A + B)−1 = A−1 B(A + B)−1 . Proof. We have A−1 = A−1 (A + B)(A + B)−1 = (I + A−1 B)(A + B)−1 = (A + B)−1 + A−1 B(A + B)−1 . Rearranging gives the result. 17
Lemma B.7 (McDiarmid’s Inequality). Let g : X m → R satisfy the bounded differences condition: for all i ∈ {1, . . . , m}, there exists ci < ∞ such that sup
g(x1 , . . . , xm ) − g(x1 , . . . , xi−1 , x̃i , xi+1 , . . . , xm ) ≤ ci .
x1 ,...,xm ,x̃i
Then for independent random variables X1 , . . . , Xm and any ξ > 0, 2ξ 2 P P g(X1 , . . . , Xm ) − E[g(X1 , . . . , Xm )] > ξ ≤ exp − m 2 . i=1 ci
B.2.2
Empirical CDF Concentration
Lemma B.8 (Convergence of Empirical CDF under i.i.d.). Under Assumption B.1, for any T ≥ 1 and δ ∈ (0, 1), with probability at least 1 − δ, r log(2/δ) e sup FT +1 (x) − Fe (x) ≤ . (18) 2T x∈R Proof. Under Assumption B.1, the true scores e1 , . . . , eT are i.i.d. with CDF Fe . By the Dvoretzky– Kiefer–Wolfowitz (DKW) inequality, for any ξ > 0, P sup FeT +1 (x) − Fe (x) > ξ ≤ 2 exp(−2T ξ 2 ). x
Setting the right-hand side equal to δ and solving: δ = 2 exp(−2T ξ 2 ) gives ξ =
B.2.3
p
log(2/δ)/(2T ).
Score Approximation Error
Lemma B.9 (Difference Between Estimated and True Nonconformity Scores). Under Assumptions B.1, B.2, B.3, and B.4, with probability at least 1 − 2δ, T
1X 4BδT 1 |et − êt | ≤ + 4B∥∆µ ∥ + ∥∆µ ∥2 T t=1 λ λ " r r 4B 2 r(Σ) _ r(Σ) _ log(1/δ) _ log(1/δ) + 2 CK ∥Σ∥op λ T T T T # r 2 2B 2 log(1/δ) 8B 2 T 2 + + 8BδT + 4δT , + √ +B T T −1 T −1 T (19) where ∆µ = µ − µ̂ is the total mean estimation error, satisfying with probability at least 1 − δ: r 2B 2 log(1/δ) ∥∆µ ∥ ≤ δT + √ + B (20) T T and CK > 0 is the constant from Koltchinskii et al. [18, Corollary 2]. Proof. We decompose the score difference via a telescoping argument. For each t = 1, . . . , T , define the intermediate quantities: (t)
I1 = ⟨v̂t , Σ−1 λ v̂t ⟩H (t)
I2 = ⟨ŵt , Σ−1 λ ŵt ⟩H where we use the notation: vt = ϕ(εt ) − µ v̂t = ϕ(ε̂t ) − µ ŵt = ϕ(ε̂t ) − µ̂
(true features, population-centered), (estimated features, population-centered), (estimated features, sample-centered). 18
b −1 Note that et = ⟨vt , Σ−1 λ vt ⟩ and êt = ⟨ŵt , Σλ ŵt ⟩, and that ŵt = v̂t + ∆µ where ∆µ = µ − µ̂. By the triangle inequality: (t)
(t)
(t)
(t)
|et − I1 | + |I1 − I2 | + | {z } | {z }
|et − êt | ≤
(t)
T1 : residual error
(t)
T2 : mean error
|I2 − êt | | {z }
.
(t)
T3 : covariance error
We bound each term separately. (t)
Bounding T1 (residual estimation error) (t)
−1 T1 = ⟨vt , Σ−1 λ vt ⟩ − ⟨v̂t , Σλ v̂t ⟩ −1 = ⟨vt − v̂t , Σ−1 λ vt ⟩ + ⟨v̂t , Σλ (vt − v̂t )⟩
= ⟨vt − v̂t , Σ−1 λ (vt + v̂t )⟩
(self-adjointness of Σ−1 λ )
≤ ∥vt − v̂t ∥ · ∥vt + v̂t ∥ · ∥Σ−1 λ ∥op
(Lemma B.5).
Since ∥Σ−1 λ ∥op ≤ 1/λ and, by the triangle inequality and Assumption B.3, ∥vt + v̂t ∥ = ∥ϕ(εt ) + ϕ(ε̂t ) − 2µ∥ ≤ ∥ϕ(εt )∥ + ∥ϕ(ε̂t )∥ + 2∥µ∥ ≤ 4B, we obtain (t)
T1 ≤
4B ∥ϕ(εt ) − ϕ(ε̂t )∥. λ
Summing over t and applying the Cauchy–Schwarz inequality: T T 1 X (t) 4B 1 X T ≤ · ∥ϕ(εt ) − ϕ(ε̂t )∥ T t=1 1 λ T t=1 v u T X 4B u t1 ≤ ∥ϕ(εt ) − ϕ(ε̂t )∥2 λ T t=1
≤
4BδT , λ
where the last step uses Assumption B.2. (t)
Bounding T2 (mean estimation error).
Since ŵt = v̂t + ∆µ :
(t)
−1 T2 = ⟨v̂t , Σ−1 λ v̂t ⟩ − ⟨v̂t + ∆µ , Σλ (v̂t + ∆µ )⟩ −1 −1 −1 = ⟨v̂t , Σ−1 λ v̂t ⟩ − ⟨v̂t , Σλ v̂t ⟩ − 2⟨v̂t , Σλ ∆µ ⟩ − ⟨∆µ , Σλ ∆µ ⟩ −1 = 2⟨v̂t , Σ−1 λ ∆µ ⟩ + ⟨∆µ , Σλ ∆µ ⟩ −1 ≤ 2 ⟨v̂t , Σ−1 λ ∆µ ⟩ + ⟨∆µ , Σλ ∆µ ⟩ .
By Cauchy–Schwarz and ∥v̂t ∥ = ∥ϕ(ε̂t ) − µ∥ ≤ 2B: −1 2 ⟨v̂t , Σ−1 λ ∆µ ⟩ ≤ 2∥v̂t ∥ · ∥Σλ ∥op · ∥∆µ ∥ ≤ −1 2 ⟨∆µ , Σ−1 λ ∆µ ⟩ ≤ ∥∆µ ∥ · ∥Σλ ∥op ≤
4B ∥∆µ ∥, λ
1 ∥∆µ ∥2 λ
Since this bound is independent of t: T
1 X (t) 1 T ≤ 4B∥∆µ ∥ + ∥∆µ ∥2 . T t=1 2 λ 19
We now bound ∥∆µ ∥ = ∥µ − µ̂∥. Introducing the sample mean of true features µ̃ = PT (1/T ) t=1 ϕ(εt ): ∥∆µ ∥ = ∥µ − µ̂∥ ≤ ∥µ − µ̃∥ + ∥µ̃ − µ̂∥. The second term is bounded by Assumption B.2: T 1X ϕ(εt ) − ϕ(ε̂t ) T t=1 v u T u1 X ≤t ∥ϕ(εt ) − ϕ(ε̂t )∥2 ≤ δT . T t=1
∥µ̃ − µ̂∥ =
For ∥µ − µ̃∥, we apply McDiarmid’s inequality (Lemma B.7). Define T
1X ϕ(εt ) − µ . T t=1
g(ε1 , . . . , εT ) = ∥µ̃ − µ∥ = Replacing any single εi changes g by at most
1 2B ∥ϕ(εi ) − ϕ(ε̃i )∥ ≤ T T so the bounded differences condition holds with ci = 2B/T for all i. By Lemma B.7: 2ξ 2 T ξ2 P ∥µ̃ − µ∥ − E∥µ̃ − µ∥ > ξ ≤ exp − = exp − 2 . T · (2B/T )2 2B |g(. . . , εi , . . .) − g(. . . , ε̃i , . . .)| ≤
We bound the expectation using the i.i.d. assumption: 1 X E ⟨ϕ(εt1 ) − µ, ϕ(εt2 ) − µ⟩ E∥µ̃ − µ∥2 = 2 T t ,t 1
2
1 = E∥ϕ(ε) − µ∥2 (cross terms vanish by independence) T 4B 2 ≤ T where the last step uses ∥ϕ(ε) − µ∥2 ≤ 2(∥ϕ(ε)∥2 + ∥µ∥2 ) ≤ 4B 2 . By Jensen’s inequality: p 2B E∥µ̃ − µ∥ ≤ E∥µ̃ − µ∥2 ≤ √ . T Setting δ = exp(−T ξ 2 /(2B 2 )) and solving for ξ: with probability at least 1 − δ, r 2B 2 log(1/δ) ∥µ̃ − µ∥ ≤ √ + B . T T Combining with ∥µ̃ − µ̂∥ ≤ δT : with probability at least 1 − δ, r 2B 2 log(1/δ) ∥∆µ ∥ ≤ δT + √ + B . T T (t)
Bounding T3 Lemma B.5:
(covariance estimation error).
Recall ŵt = ϕ(ε̂t ) − µ̂, so ∥ŵt ∥ ≤ 2B. By
(t) b −1 T3 = ⟨ŵt , Σ−1 λ ŵt ⟩ − ⟨ŵt , Σλ ŵt ⟩ b −1 )ŵt ⟩ = ⟨ŵt , (Σ−1 − Σ λ
λ
b −1 ≤ ∥ŵt ∥2 · ∥Σ−1 λ − Σλ ∥op b −1 ∥op . ≤ 4B 2 · ∥Σ−1 − Σ λ
20
(21)
λ
b − Σ: By the resolvent identity (Lemma B.6) with A = Σ + λI and B = Σ −1 b b −1 b + λI)−1 ∥op ∥Σ−1 (Σ − Σ)(Σ λ − Σλ ∥op = ∥(Σ + λI) b − Σ∥op · ∥(Σ b + λI)−1 ∥op ≤ ∥(Σ + λI)−1 ∥op · ∥Σ
≤
1 b ∥Σ − Σ∥op . λ2
b − Σ∥op ≤ ∥Σ e − Σ∥op + ∥Σ b − Σ∥ e op and bound each part. We decompose ∥Σ e − Σ∥op . Define vt = ϕ(εt ) − µ, ∆ = µ̃ − µ = (1/T ) PT vt , and CT = Statistical error ∥Σ t=1 PT (1/T ) t=1 vt ⊗ vt . e = 1 PT (ϕ(εt ) − µ̃) ⊗ (ϕ(εt ) − µ̃) and ϕ(εt ) − µ̃ = vt − ∆, a direct computation Since Σ t=1 T −1 gives: T T 1X · (vt − ∆) ⊗ (vt − ∆) T − 1 T t=1 X T 1 T T vt ⊗ v t − ∆ ⊗ ∆ = CT − ∆ ⊗ ∆ . = T − 1 T t=1 T −1
e= Σ
Therefore: T (CT − ∆ ⊗ ∆) − Σ T −1 T 1 CT − (∆ ⊗ ∆). = (CT − Σ) + T −1 T −1
e −Σ= Σ
1 ∥∆∥2 : Taking operator norms and using T T−1 ∥∆∥2 = ∥∆∥2 + T −1
e − Σ∥op ≤ ∥CT − Σ∥op + ∥∆∥2 + ∥Σ
1 1 ∥CT ∥op + ∥∆∥2 . T −1 T −1
Since ∥vt ∥ ≤ 2B (Assumption B.3): T
∥CT ∥op ≤
1X ∥vt ∥2 ≤ 4B 2 T t=1
2
T
∥∆∥ ≤
1X ∥vt ∥ ≤ 2B T t=1
2
1 1 4B 2 giving T −1 ∥CT ∥op ≤ T4B −1 and T −1 ∥∆∥ ≤ T −1 .
For ∥CT − Σ∥op : since ∥vt ∥ ≤ 2B, the random vectors vt are bounded, hence subgaussian. Combined with tr(Σ) < ∞ (Assumption B.4), the vt are pregaussian, satisfying the conditions of [18, Corollary 2]. For any δ ∈ (0, 1), with probability at least 1 − δ: r r r(Σ) _ r(Σ) _ log(1/δ) _ log(1/δ) ∥CT − Σ∥op ≤ CK ∥Σ∥op (22) T T T T where r(Σ) = E∥vt ∥2 /∥Σ∥op ≤ 4B 2 /∥Σ∥op . Using the bound on ∥∆∥2 from the mean estimation analysis (which holds with probability at least 1 − δ, a union bound over the covariance and mean events gives: with probability at least 1 − 2δ, r r r(Σ) _ r(Σ) _ log(1/δ) _ log(1/δ) e ∥Σ − Σ∥op ≤ CK ∥Σ∥op T T T T (23) r 2 2B 2 log(1/δ) 8B 2 + √ +B + . T T −1 T 21
b − Σ∥ e op . Define w̃t = ϕ(εt ) − µ̃ (true features, sample-centered), ∆t = Prediction error ∥Σ ¯ = (1/T ) PT ∆t = µ̂ − µ̃ (the mean prediction ϕ(ε̂t ) − ϕ(εt ) (a pointwise prediction error), and ∆ t=1 error). ¯ = w̃t + θt , where θt = ∆t − ∆ ¯ is the demeaned Then ŵt = ϕ(ε̂t ) − µ̂ = (ϕ(εt ) + ∆t ) − (µ̃ + ∆) prediction error. Expanding: b − Σ∥ e op = ∥Σ =
T 1 X ŵt ⊗ ŵt − w̃t ⊗ w̃t T − 1 t=1 op T 1 X w̃t ⊗ θt + θt ⊗ w̃t + θt ⊗ θt T − 1 t=1 op
T 1 X 2∥w̃t ∥ ∥θt ∥ + ∥θt ∥2 . ≤ T − 1 t=1
¯ ≤ ∥∆t ∥ + ∥∆∥. ¯ By Assumption B.2: We have ∥w̃t ∥ = ∥ϕ(εt ) − µ̃∥ ≤ 2B and ∥θt ∥ = ∥∆t − ∆∥ v u T T u1 X 1X ¯ ∆t ≤ t ∥∆t ∥2 ≤ δT . ∥∆∥ = T t=1 T t=1 Using (1/T )
P
t ∥∆t ∥ ≤ δT and (1/T )
b − Σ∥ e op ≤ ∥Σ ≤
2 2 t ∥∆t ∥ ≤ δT :
P
T 1X T 4B(∥∆t ∥ + δT ) + (∥∆t ∥ + δT )2 · T − 1 T t=1
T 8BδT + 4δT2 . T −1
Combining. Summing and dividing by T : T i 1 X (t) 4B 2 h e b − Σ∥ e op . T3 ≤ 2 ∥Σ − Σ∥op + ∥Σ T t=1 λ
e − Σ∥op (Equation (23)) and ∥Σ b − Σ∥ e op hold simultaneously with probability at The bounds on ∥Σ least 1 − 2δ. Adding T1 + T2 + T3 gives (19). B.2.4
ECDF Coupling
Lemma B.10 (Distance Between Empirical CDFs of True and Estimated Scores). Under Assumptions B.1, B.2, B.3, and B.4, with probability at least 1 − 2δ, sup FbT +1 (x) − FeT +1 (x) ≤ (L + 1) CS + 2 sup FeT +1 (x) − Fe (x) x∈R
x∈R
1/2 PT where CS = (1/T ) t=1 |et − êt | and L is the Lipschitz constant of Fe . Proof. The argument follows [40, Lemma 4.11]. For any threshold x: T T 1 X 1X b e FT +1 (x) − FT +1 (x) = 1{êt ≤x} − 1{et ≤x} ≤ 1{|et −x|≤|et −êt |} . T t=1 T t=1
For any η > 0, we split each indicator:
1{|et −x|≤|et −êt |} ≤ 1{|et −x|≤η} + 1{|et −êt |>η} . 22
(24)
For the second term, by Markov’s inequality: T T 1 1X C2 1X 1{|et −êt |>η} ≤ · |et − êt | = S . T t=1 η T t=1 η
For the first term, using the Lipschitz property of Fe : T
1X 1{|et −x|≤η} ≤ FeT +1 (x + η) − FeT +1 (x − η) T t=1 ≤ Fe (x + η) − Fe (x − η) + 2 sup |FeT +1 (x) − Fe (x)| x
≤ 2Lη + 2 sup |FeT +1 (x) − Fe (x)|. x
Combining: |FbT +1 (x) − FeT +1 (x)| ≤ 2Lη + CS2 /η + 2 supx |FeT +1 (x) − Fe (x)|. The stated bound with constant (L + 1) follows from optimizing over η; see [40, Lemma 4.11] for the sharp constant. B.2.5
Test Point Score Error
Lemma B.11 (Test Point Nonconformity Score Error). Under Assumptions B.1, B.2, B.3, and B.4, with probability at least 1 − 2δ, |eT +1 − êT +1 | ≤ ζT (25) where ζT has the same form as the right-hand side of (19). Proof. By the same telescoping decomposition: (T +1)
|eT +1 − êT +1 | ≤ T1 (T +1)
(T +1)
+ T2
(T +1)
+ T3
.
(T +1)
The bounds on T2 and T3 are identical to those in Lemma B.9, since they depend only on (T +1) b µ̂, Σ, and Σ (shared between calibration and test point). For T1 , the same derivation gives: (T +1)
T1
≤
4BδT 4B ∥ϕ(εT +1 ) − ϕ(ε̂T +1 )∥ ≤ λ λ
using the test-point bound in Assumption B.2. P Summing all three terms gives |eT +1 − êT +1 | ≤ ζT , where ζT has the same expression as (1/T ) t |et − êt | = CS2 from Lemma B.9. B.2.6
Coverage Guarantee
Theorem B.12 (Asymptotic Conditional Coverage Guarantee under i.i.d.). Under Assumptions B.1, B.2, B.3, and B.4, for any T ≥ 1, α ∈ (0, 1), and δ ∈ (0, 1), with probability at least 1 − 4δ: r log(2/δ) α b P(YT +1 ∈ CT +1 | XT +1 ) − (1 − α) ≤ 12 + 4(L + 1)(CS + ζT ) + 2δ. 2T Proof. Under Assumption B.1, εT +1 ⊥ XT +1 , so eT +1 is a function of εT +1 alone and satisfies Fe (eT +1 ) | XT +1 ∼ Uniform(0, 1). Furthermore, conditioned on XT +1 = x, the estimated score êT +1 = êk (f (x) − fˆ(x) + εT +1 ) depends on x only through the deterministic shift f (x) − fˆ(x), with all remaining randomness from εT +1 ⊥ x. All subsequent bounds therefore hold conditionally on any fixed XT +1 = x, and the | XT +1 notation is justified throughout this proof. bα : For any β ∈ [0, α], by construction of C T +1 bTα+1 ⇐⇒ β ≤ FbT +1 (êT +1 ) ≤ 1 − α + β YT +1 ∈ C 23
Since Fe (eT +1 ) ∼ Uniform(0, 1) under Assumption B.1, by the elementary bound |1{a ≤ x ≤ b} − 1{a ≤ y ≤ b}| ≤ |1{a ≤ x} − 1{a ≤ y}| + |1{x ≤ b} − 1{y ≤ b}|: b α | XT +1 ) − (1 − α) P(YT +1 ∈ C T +1 ≤ P |Fe (eT +1 ) − β| ≤ |FbT +1 (êT +1 ) − Fe (eT +1 )| + P |Fe (eT +1 ) − (1 − α + β)| ≤ |FbT +1 (êT +1 ) − Fe (eT +1 )| p Let AT denote the event in Lemma B.8: supx |FeT +1 (x) − Fe (x)| ≤ log(2/δ)/(2T ), with P(AT ) ≥ 1 − δ. For any γ ∈ [0, 1]: P |Fe (eT +1 ) − γ| ≤ |FbT +1 (êT +1 ) − Fe (eT +1 )| ≤ P |Fe (eT +1 ) − γ| ≤ |FbT +1 (êT +1 ) − Fe (êT +1 )| + |Fe (êT +1 ) − Fe (eT +1 )| AT + δ Conditioning on AT , we bound |FbT +1 (êT +1 ) − Fe (êT +1 )| by the triangle inequality: |FbT +1 (êT +1 ) − Fe (êT +1 )| ≤ sup |FbT +1 (x) − FeT +1 (x)| + sup |FeT +1 (x) − Fe (x)| x
x
By Lemma B.10 and the definition of AT : r
log(2/δ) 2T
r
log(2/δ) 2T
sup |FbT +1 (x) − FeT +1 (x)| ≤ (L + 1)CS + 2 x
Adding supx |FeT +1 (x) − Fe (x)| ≤
p
log(2/δ)/(2T ) on AT :
|FbT +1 (êT +1 ) − Fe (êT +1 )| ≤ (L + 1)CS + 3 By the Lipschitz continuity of Fe and Lemma B.11:
|Fe (êT +1 ) − Fe (eT +1 )| ≤ L|êT +1 − eT +1 | ≤ LζT Combining: r
log(2/δ) |FbT +1 (êT +1 ) − Fe (eT +1 )| ≤ (L + 1)CS + 3 + LζT 2T r log(2/δ) + (L + 1)(CS + ζT ) ≤3 2T Since Fe (eT +1 ) ∼ Uniform(0, 1), P(|Fe (eT +1 ) − γ| ≤ η | AT ) ≤ 2η. Therefore: P |Fe (eT +1 ) − γ| ≤ |FbT +1 (êT +1 ) − Fe (eT +1 )| AT r log(2/δ) ≤6 + 2(L + 1)(CS + ζT ) 2T Adding P(AcT ) ≤ δ and applying with γ = β and γ = 1 − α + β: bTα+1 | XT +1 ) − (1 − α) P(YT +1 ∈ C r i h log(2/δ) ≤2 6 + 2(L + 1)(CS + ζT ) + δ 2T r log(2/δ) = 12 + 4(L + 1)(CS + ζT ) + 2δ 2T The probability is at least 1 − 4δ by a union bound over: DKW (1 − δ), mean concentration (1 − δ), covariance concentration (1 − δ), and test-point estimation (1 − δ). 24
B.2.7
Convergence Rate
Theorem B.13 (Conditional Coverage Convergence Rate under i.i.d.). Under Assumptions B.1–B.4, there exists a constant C > 0 depending only on B, CK , ∥Σ∥op , L such that with probability at least 1 − 4T −1 : ( C · T β−1/4 (log T )1/4 if ρ ≥ 1/2, λ = T −β , β ∈ (0, 1/4) bTα+1 | XT +1 )−(1−α) ≤ P(YT +1 ∈ C C · T β−ρ/2 if ρ < 1/2, λ = T −β , β ∈ (0, ρ/2) Proof. Set δ = T −1 (so log(1/δ) = log T ) and λ = T −β . We determine the asymptotic rate of each term in Lemma B.9. Residual error term.
Mean error term.
4BδT = O(T β−ρ ) λ
(26)
From (21):
r p 2 log T 2B = O δT + T −1/2 log T ∥∆µ ∥ ≤ δT + √ + B T T √ √ If ρ ≥ 1/2, the T −1/2 log T term dominates δT = O(T −ρ ), so ∥∆µ ∥ = O(T −1/2 log T ). Since ∥∆µ ∥2 = O(T −1 log T ) is dominated by ∥∆µ ∥ for large T : p 1 4B∥∆µ ∥ + ∥∆µ ∥2 = O T β−1/2 log T λ If ρ < 1/2, then δT dominates and: 1 4B∥∆µ ∥ + ∥∆µ ∥2 = O(T β−ρ ) λ In both cases, this term vanishes for β < min(ρ, 1/2). We write: p 1 4B∥∆µ ∥ + ∥∆µ ∥2 = O T β−min(ρ,1/2) log T λ
(27)
Covariance error term. Define the effective rank bound r0 = 4B 2 /∥Σ∥op ≥ r(Σ). For T ≥ er0 , we have log T ≥ r0 ≥ r(Σ), so: r r r(Σ) log T r(Σ) log T ≤ , ≤ T T T T p Since log T /T dominates log T /T for all T , the Koltchinskii bound (22) simplifies to: r log T ∥CT − Σ∥op = O (28) T The remaining terms in (23):
2B √ +B T
8B 2 = O(T −1 ), T −1
r
2 log T T
2
log T =O T
T 8BδT + 4δT2 = O(δT ) = O(T −ρ ). T −1
(29) (30)
p Since O( log T /T ) dominates both O(log T /T ) and O(T −1 ), and O(T −ρ ) dominates O(T −1 ) for ρ < 1, we denote the full covariance error bracket by A: r log T e − Σ∥op + ∥Σ b − Σ∥ e op = O A = ∥Σ + O(T −ρ ). T 25
Multiplying by 4B 2 /λ2 = O(T 2β ): p 4B 2 2β −1/2 −ρ · A = O(T ) · O T log T + O(T ) . λ2
(31)
For ρ ≥ 1/2, the first term dominates: p 4B 2 · A = O(T 2β−1/2 log T ) 2 λ which vanishes for β < 1/4. For ρ < 1/2, the second term dominates: 4B 2 · A = O(T 2β−ρ ) λ2 which vanishes for β < ρ/2. Dominant term.
Comparing the exponents from (26), (27), and (31): β − ρ ≤ β − min(ρ, 1/2) < 2β − min(ρ, 1/2)
so the covariance term (31) dominates. Thus: √ O(T 2β−1/2 log T ) CS2 = ζT = O(T 2β−ρ )
if ρ ≥ 1/2, if ρ < 1/2
(32)
where we use that CS2 is the calibration average bound from Lemma B.9 and ζT is the test-point bound from Lemma B.11. Taking square roots: CS = O(T β−1/4 (log T )1/4 ) for ρ ≥ 1/2, and CS = O(T β−ρ/2 ) for ρ < 1/2. From Theorem B.12 with δ = T −1 , with probability at least 1 − 4T −1 : r log(2T ) α b P(YT +1 ∈ CT +1 | XT +1 ) − (1 − α) ≤ 12 + 4(L + 1)(CS + ζT ) + |2T{z−1} . {z } 2T | {z } | O(T −1 ) score error
Final bound.
O(T −1/2
√
log T )
Since ζT = CS2 and CS → 0 (for β < min(ρ/2, 1/4)), we have CS + ζT ≤ 2CS for large T . The CS term dominates the DKW and T −1 terms because β − min(ρ/2, 1/4) > −1/2 > −1 for β, ρ > 0, giving: ( C · T β−1/4 (log T )1/4 if ρ ≥ 1/2 α bT +1 | XT +1 ) − (1 − α) ≤ P(YT +1 ∈ C C · T β−ρ/2 if ρ < 1/2 where C depends on B, CK , ∥Σ∥op , L. B.3
Coverage under α-Mixing Dependence
We extend the coverage analysis to sequences with temporal dependence. The key simplification is that the mean and covariance operator are assumed known (Assumption B.15), eliminating the mean and covariance estimation error and isolating the effect of mixing on conformal calibration. B.3.1
Additional Assumptions
Assumption B.14 (Stationary Strong Mixing). P∞ The process {εt }t≥1 is stationary and α-mixing with mixing coefficients {αk }k≥0 , satisfying k=0 αk ≤ M < ∞. Furthermore, the CDF Fe of the true nonconformity score is L-Lipschitz. Assumption B.15 (Known Mean and Covariance). The population mean µ = E[ϕ(ε)] and covariance operator Σ = E[(ϕ(ε) − µ) ⊗ (ϕ(ε) − µ)] are known. 26
b = Σ in the construction of the nonconRemark. Under Assumption B.15, we set µ̂ = µ and Σ formity scores for analysis. This eliminates the mean estimation error (term T2 ) and the covariance estimation error (term T3 ) from the decomposition in Lemma B.9, leaving only the residual estimation error (term T1 ). This assumption isolates the effect of temporal dependence on conformal calibration. Analyzing covariance operator estimation under α-mixing in infinite-dimensional RKHSs would require concentration inequalities considerably more involved than their i.i.d. counterparts and is beyond the scope of this work. In practice, Σ can be estimated from an independent sample or from a sufficiently long stationary pre-training period. B.3.2
Inheritance of Mixing
Lemma B.16 (Inheritance of Mixing by Nonconformity Scores). Let {εt }t≥1 be a stationary αmixing sequence with mixing coefficients {αk }k≥0 . Define et = ⟨ϕ(εt ) − µ, Σ−1 λ (ϕ(εt ) − µ)⟩H . Then {et }t≥1 is stationary and α-mixing with αk ({et }) ≤ αk ({εt }) for all k ≥ 0. Proof. The map ε 7→ e(ε) = ⟨ϕ(ε) − µ, Σ−1 λ (ϕ(ε) − µ)⟩H is measurable as a composition of continuous functions: ϕ is continuous (for RBF, Matérn, polynomial kernels) and v 7→ ⟨v, Σ−1 λ v⟩ is continuous on H since ∥Σ−1 ∥ ≤ 1/λ < ∞. Therefore: op λ σ(et : t ∈ I) ⊆ σ(εt : t ∈ I) for any index set I ⊆ Z. Since α-mixing coefficients are defined as suprema over sub-σ-algebras, the containment gives αk ({et }) ≤ αk ({εt }) for all k ≥ 0. B.3.3
ECDF Concentration under Mixing
Lemma B.17 (Empirical CDF Concentration under Mixing). Under Assumption B.14, for the true nonconformity scores {et }Tt=1 , with probability at least 1 − (M (log T )2 /(2T ))1/3 : sup FeT +1 (x) − Fe (x) ≤ x∈R
(M/2)1/3 (log T )2/3 . T 1/3
(33)
P∞ Proof. By Lemma B.16, the scores {et }Tt=1 are stationary and α-mixing with k=0 αk ≤ M . The result follows from Rio’s inequality [26, Proposition 7.1] applied to the indicator processes {1{et ≤x} }Tt=1 , following the argument of Xu et al. [40, Lemma B.11]. B.3.4
Test Point Error under Known Covariance
Lemma B.18 (Test Point Nonconformity Score Error under Known Mean and Covariance). Under Assumptions B.14, B.15, B.2, and B.3: 4BδT |eT +1 − êT +1 | ≤ ζ̃T = . (34) λ b = Σ. In the telescoping decomposition, T (T +1) = 0 Proof. Under Assumption B.15, µ̂ = µ and Σ 2 (T +1) (no mean error) and T3 = 0 (no covariance error). Only the residual error remains: 4B 4BδT (T +1) |eT +1 − êT +1 | = T1 ≤ ∥ϕ(εT +1 ) − ϕ(ε̂T +1 )∥ ≤ . λ λ B.3.5
Coverage under Mixing
Theorem B.19 (Marginal Coverage Guarantee under α-Mixing). Under Assumptions B.14, B.15, B.2, and B.3, for any T ≥ 2, α ∈ (0, 1), and λ > 0, with probability at least 1 − (M (log T )2 /(2T ))1/3 : 1/3
bTα+1 ) − (1 − α) ≤ 16 · (M/2) (log T ) P(YT +1 ∈ C T 1/3 p where C̃S = 4BδT /λ and ζ̃T = 4BδT /λ = C̃S2 . 27
2/3
+ 4(L + 1)(C̃S + ζ̃T ),
Under α-mixing dependence, XT +1 typically contains lags of Yt and is therefore correlated with εT +1 . The exogeneity condition εT +1 ⊥ XT +1 required for conditional coverage (Assumption B.1) fails in this setting, and the bound is strictly marginal. Exact conditional coverage is known to be unachievable distribution-free [37, 11]. Approximate conditional coverage is assessed empirically via worst-slab coverage in Section 5.1. Proof. Under Assumption B.15, the same argument as Lemma B.18 applies to each calibration point: for t = 1, . . . , T , 4BδT |et − êt | ≤ = C̃S2 , λ PT so (1/T ) t=1 |et − êt | ≤ C̃S2 . By the same argument as Lemma B.10 (which does not rely on independence of the sequence): sup |FbT +1 (x) − FeT +1 (x)| ≤ (L + 1)C̃S + 2 sup |FeT +1 (x) − Fe (x)|. x
x
Define AT = {supx |FeT +1 (x) − Fe (x)| ≤ (M/2)1/3 (log T )2/3 /T 1/3 }. By Lemma B.17: 1/3 M (log T )2 P(AT ) ≥ 1 − . 2T Following the framework of Theorem B.12, conditioning on AT : |FbT +1 (êT +1 ) − Fe (eT +1 )| ≤ sup |FbT +1 (x) − FeT +1 (x)| + sup |FeT +1 (x) − Fe (x)| + Lζ̃T x
x
1/3
≤ (L + 1)C̃S + 3 ·
(M/2)
(log T )2/3
+ Lζ̃T T 1/3 (M/2)1/3 (log T )2/3 . ≤ (L + 1)(C̃S + ζ̃T ) + 3 · T 1/3 Since Fe (eT +1 ) ∼ Uniform(0, 1) under stationarity, P(|Fe (eT +1 ) − γ| ≤ η | AT ) ≤ 2η. Note also that 1/3 M (log T )2 (M/2)1/3 (log T )2/3 c P(AT ) ≤ <2· for T ≥ 2. 2T T 1/3 Therefore, for any γ ∈ [0, 1]: P |Fe (eT +1 ) − γ| ≤ |FbT +1 (êT +1 ) − Fe (eT +1 )| ≤ 2(L + 1)(C̃S + ζ̃T ) + 8 ·
(M/2)1/3 (log T )2/3 . T 1/3
Applying with γ = β and γ = 1 − α + β and summing gives the stated bound. B.3.6
Convergence Rate under Mixing
Theorem B.20 (Marginal Coverage Convergence Rate under α-Mixing). Under Assumptions B.14, B.15, B.2, and B.3, there exists C > 0 depending only on B, M , L such that with probability at least 1 − (M (log T )2 /(2T ))1/3 : ( C · T −1/3 (log T )2/3 if ρ > 2/3, λ = T −β , β ∈ (0, ρ − 2/3) α b P(YT +1 ∈ CT +1 ) − (1 − α) ≤ C · T (β−ρ)/2 if ρ ≤ 2/3, λ = T −β , β ∈ (0, ρ) Proof. With δT = O(T −ρ ) and λ = T −β : r 4BδT C̃S = = O(T (β−ρ)/2 ), λ 28
ζ̃T =
4BδT = O(T β−ρ ) λ
(35)
Since β < ρ, both vanish as T → ∞. For large T where C̃S < 1: C̃S + ζ̃T = C̃S (1 + C̃S ) ≤ 2C̃S so 4(L + 1)(C̃S + ζ̃T ) = O(T (β−ρ)/2 ). From Theorem B.19: b α ) − (1 − α) ≤ O(T −1/3 (log T )2/3 ) + P(YT +1 ∈ C T +1 | {z } mixing ECDF term
O(T (β−ρ)/2 ) | {z }
.
score approximation term
The score term decays faster than the mixing term when: β−ρ 1 <− 2 3
⇐⇒
2 β <ρ− . 3
Case ρ > 2/3. We can choose β ∈ (0, ρ − 2/3). Then: (ρ − 2/3) − ρ 1 β−ρ < =− , 2 2 3 so the mixing ECDF term dominates, giving rate O(T −1/3 (log T )2/3 ). Case ρ ≤ 2/3. We choose β ∈ (0, ρ). Since ρ ≤ 2/3, for any β > 0: β−ρ −ρ 1 ≥ ≥− , 2 2 3 so the score term dominates, giving rate O(T (β−ρ)/2 ).
C
Evaluation Details and Additional Experimental Results
C.1
Volume Estimation
For the Mahalanobis score, the prediction region is an ellipsoid and its volume admits a closed-form expression π d/2 b 1/2 · q̂ d/2 · det(Σ) V= Γ(d/2 + 1) b is the sample covariance. For all other methods, the where q̂ is the conformal threshold and Σ prediction region has no closed form and we estimate its volume by Monte Carlo. We sample uniformly from a bounding box constructed by padding the range of the calibration residuals by 30% in each coordinate, score each sample, and compute the fraction falling inside the prediction region. We use 500,000 MC samples for all experiments. We use the same random samples (controlled by a fixed seed) for all scorers within each experimental seed. C.2
Worst-Slab Coverage
Worst-slab coverage (WSC) [7] approximates conditional coverage by searching for subgroups of the test data with low empirical coverage. Let ci = 1{Yi ∈ C(Xi )} denote the coverage indicator for test point i. We consider slabs of the form {x : a ≤ v ⊤ x ≤ b} where v is a unit direction. For each direction v, we project the test features onto v and compute the empirical quantiles of the projections {v ⊤ Xi }. We construct slabs that contain approximately 20% of the data. Worst-slab coverage is then defined as P ci 1{Qv (q) ≤ v ⊤ Xi ≤ Qv (q + 0.2)} WSC = min Pi . ⊤ v,q i 1{Qv (q) ≤ v Xi ≤ Qv (q + 0.2)} In practice, we approximate this search by sampling 200 random unit directions and evaluating 20 quantile positions per direction. Each slab is required to contain at least 30 test points. The same set of random directions is used for all methods and random seeds. 29
C.3 C.3.1
Additional Details and Results on Synthetic Data Data Generating Process
We consider a bivariate regression task with a single input X ∼ Uniform(0, 1) and response Y = f (X) + ε ∈ R2 , where the true regression function is f1 (x) = 3x2 − 1.5x + sin(4πx), f2 (x) = 2x3 − x + 0.5 cos(3πx). The noise is constructed as ε = (z1 , z2 +0.5 z12 −0.5σ12 ) where z1 ∼ N (0, 0.32 ) and z2 ∼ N (0, 0.12 ) are independent. The noise is also centered by subtracting its sample mean. We generate n = 50,000 samples and use a 50/25/25 train/calibration/test split. C.3.2
Models
We use two fitted models to compare the conformal prediction results: a simple ordinary least squares (LinearRegression, wrapped in MultiOutputRegressor) and a more complex fitting via a two-layer MLP with hidden sizes (64, 64) trained for up to 1000 iterations. C.3.3
Density vs Kernel Ablation
Table 3 compares the kernel-based and density scores across seven lengthscale values under both models. The MKS achieves lower volume than the density score at every lengthscale and every coverage level, confirming that the KPCA correction is the source of the improvement rather than a particular choice of parameter. For instance, under the linear model at α = 0.01, the MKS volume ranges from 3.85 (l = 0.1) to 5.25 (l = 5.0), while the density score ranges from 4.14 to 7.36 over the same values. The gap is consistent across all lengthscales and widens at tighter coverage levels. C.4 C.4.1
Additional Details and Results on Real Data Datasets
We use three datasets previously used in Feldman et al. [10] and Romano et al. [27]: House (King County house sales, Kaggle) [14], Bio (physicochemical properties of protein tertiary structure, UCI #265) [25], and Blog (BlogFeedback, UCI #304) [6]. Following Feldman et al. [10], we augment the original univariate response with one additional variable to obtain a bivariate target (d = 2). We further extend each dataset to d = 3 and d = 4 by including additional response variables from the same source, listed in Table 4. All datasets are split 50/25/25 into training, calibration, and test sets. Features and targets are standardized using the training set statistics. Results are averaged over 50 random seeds controlling the split permutation. Full results including the comparison among all methods, coverages, volume and WSC are provided in Section C.4.3. Dataset
n
Response
Additional response 1
Add. response 2
Add. response 3
House Bio Blog
21,613 45,730 52,397
price RMSD col. 281
lat F7 col. 61
sqft_living15 F8 col. 62
sqft_lot15 F9 col. 280
Table 4: Dataset characteristics and target variables. C.4.2
Models
We also use two models for our real data experiment, with a simpler and a more complex model choice. Ridge regression (regularisation parameter chosen to be 1.0) and a two-layer MLP with hidden sizes (128, 64) trained for up to 500 iterations. For multivariate targets, Ridge is wrapped in MultiOutputRegressor, while the MLP handles multi-output natively. All models are implemented with SciKit-Learn. The two models represent different levels of misspecification: Ridge regression cannot capture the nonlinear structure in the data and produces residuals with pronounced nonelliptical geometry. This allows us to evaluate the MKS under different fitting regimes. 30
31
MLP
Linear
Model
0.90052±0.00582 0.90139±0.00522
0.90070±0.00568 0.90110±0.00541
0.90090±0.00540 0.90090±0.00517
0.90103±0.00539 0.90037±0.00519
Density MKS
Density MKS
Density MKS
Density MKS
0.80
1.00
2.00
5.00
0.80511±0.01986 0.53423±0.01079
0.80086±0.01970 0.49489±0.01089
0.78709±0.01986 0.48445±0.01015
0.77771±0.01979 0.47791±0.00917
0.74302±0.01815 0.46954±0.00888
0.66943±0.01504 0.46516±0.00848
0.48773±0.00952 0.45945±0.00831
5.07019±0.09638 3.04592±0.04627
4.96452±0.08779 2.79988±0.03948
4.53486±0.07131 2.65114±0.04094
4.20237±0.06372 2.60484±0.04225
3.33558±0.04632 2.50061±0.03817
2.74557±0.04104 2.42346±0.03463
2.40666±0.03372 2.32507±0.02593
Volume
α = 0.10
0.88257±0.01086 0.88171±0.01133
0.88248±0.01072 0.87605±0.01227
0.88205±0.01071 0.87395±0.01198
0.88210±0.01096 0.87367±0.01127
0.88233±0.01053 0.87310±0.01210
0.88200±0.01067 0.87229±0.01051
0.87529±0.01077 0.87014±0.01040
0.59986±0.02441 0.77005±0.02118
0.58943±0.02418 0.80586±0.01531
0.58110±0.02373 0.81214±0.01574
0.58948±0.02272 0.81433±0.01681
0.66976±0.02183 0.81824±0.01708
0.81057±0.01539 0.82324±0.01776
0.80686±0.01878 0.81424±0.01750
WSC
0.95092±0.00442 0.95115±0.00387
0.95092±0.00445 0.95112±0.00409
0.95107±0.00436 0.95055±0.00413
0.95098±0.00437 0.95040±0.00389
0.95097±0.00417 0.95041±0.00314
0.95086±0.00428 0.95011±0.00349
0.95015±0.00399 0.94791±0.00382
0.95144±0.00453 0.95166±0.00406
0.95145±0.00475 0.95116±0.00452
0.95122±0.00491 0.95099±0.00459
0.95127±0.00488 0.95046±0.00412
0.95197±0.00445 0.95081±0.00403
0.95109±0.00407 0.95064±0.00358
0.94990±0.00327 0.94345±0.00365
Coverage
1.08088±0.03066 0.76402±0.02121
1.07643±0.03049 0.69110±0.01666
1.06211±0.02985 0.66165±0.01484
1.05129±0.02940 0.64439±0.01236
1.00827±0.02728 0.62425±0.01229
0.91608±0.02466 0.61330±0.01168
0.65108±0.01412 0.59864±0.01196
6.00493±0.14783 3.77140±0.07772
5.92955±0.13990 3.45757±0.06024
5.50945±0.08870 3.27998±0.06195
5.14043±0.06718 3.20406±0.05628
4.08088±0.05998 3.07326±0.05210
3.43340±0.06112 2.98107±0.04969
2.96526±0.04963 2.85422±0.03656
Volume
α = 0.05
0.94048±0.00720 0.94062±0.00597
0.94043±0.00719 0.93871±0.00647
0.94062±0.00718 0.93619±0.00671
0.94057±0.00709 0.93495±0.00617
0.94052±0.00662 0.93329±0.00657
0.94057±0.00664 0.93171±0.00718
0.93538±0.00717 0.92833±0.00786
0.78548±0.02030 0.88700±0.01343
0.78014±0.02203 0.90210±0.01187
0.77838±0.02235 0.90152±0.01184
0.78571±0.02241 0.90048±0.01211
0.83886±0.01876 0.90452±0.01236
0.90014±0.01083 0.90848±0.01016
0.90229±0.01342 0.90400±0.01245
WSC
0.98019±0.00217 0.98009±0.00233
0.98015±0.00217 0.98025±0.00200
0.98012±0.00207 0.98021±0.00226
0.98011±0.00202 0.98010±0.00225
0.98018±0.00205 0.97971±0.00240
0.98039±0.00193 0.97937±0.00250
0.97989±0.00233 0.97816±0.00292
0.98146±0.00327 0.98095±0.00240
0.98143±0.00325 0.98121±0.00226
0.98142±0.00324 0.98101±0.00203
0.98157±0.00322 0.98117±0.00223
0.98135±0.00323 0.98057±0.00195
0.98051±0.00219 0.98010±0.00199
0.97970±0.00170 0.97470±0.00158
Coverage
1.41761±0.05991 1.13080±0.03827
1.41468±0.05860 1.01828±0.03039
1.40416±0.05389 0.95271±0.02952
1.39582±0.05130 0.90393±0.02528
1.35669±0.04473 0.84846±0.02124
1.24713±0.03970 0.81882±0.02065
0.87171±0.02430 0.78432±0.02221
6.88021±0.18977 4.61485±0.08904
6.80431±0.17419 4.31707±0.08659
6.36322±0.13938 4.11646±0.07803
5.99954±0.11421 4.03307±0.07528
4.93295±0.09635 3.83086±0.07402
4.27948±0.08756 3.67712±0.07187
3.65743±0.06676 3.45636±0.06150
Volume
α = 0.02
0.97257±0.00384 0.97324±0.00358
0.97248±0.00389 0.97257±0.00338
0.97238±0.00365 0.97214±0.00390
0.97252±0.00374 0.97162±0.00441
0.97248±0.00366 0.97033±0.00451
0.97276±0.00369 0.96933±0.00553
0.97105±0.00434 0.96614±0.00675
0.91548±0.01560 0.95452±0.00762
0.91400±0.01592 0.96081±0.00595
0.91424±0.01586 0.96148±0.00616
0.91810±0.01576 0.96157±0.00729
0.94090±0.01312 0.96090±0.00806
0.95948±0.00613 0.96133±0.00792
0.96081±0.00635 0.95533±0.00645
WSC
0.99007±0.00137 0.99009±0.00158
0.99004±0.00139 0.98995±0.00140
0.99000±0.00138 0.98993±0.00148
0.99000±0.00137 0.98978±0.00144
0.99000±0.00140 0.98989±0.00156
0.98993±0.00140 0.98983±0.00142
0.98969±0.00157 0.98829±0.00146
0.99081±0.00210 0.99062±0.00196
0.99090±0.00202 0.99028±0.00186
0.99093±0.00183 0.98991±0.00204
0.99099±0.00186 0.99007±0.00186
0.99035±0.00192 0.98990±0.00175
0.99012±0.00172 0.98985±0.00157
0.98909±0.00145 0.98487±0.00194
Coverage
1.63999±0.08516 1.42811±0.06867
1.63769±0.08294 1.28713±0.04928
1.62945±0.07586 1.20139±0.04692
1.62273±0.07164 1.13500±0.04229
1.59159±0.05810 1.05777±0.03602
1.48702±0.05129 1.00377±0.02982
1.04602±0.03084 0.92214±0.03192
7.35608±0.21110 5.24838±0.13560
7.28806±0.20310 4.95440±0.11322
6.84038±0.16236 4.75424±0.10536
6.48079±0.13796 4.66199±0.09676
5.50337±0.10867 4.42099±0.09306
4.92633±0.11617 4.21487±0.08557
4.13983±0.07949 3.84754±0.06354
Volume
α = 0.01
Table 3: Density vs MKS ablation across lengthscale values on 2D synthetic data (21 seeds). Bold indicates lowest volume.
0.90062±0.00573 0.90115±0.00560
5.00
Density MKS
0.90103±0.00585 0.90080±0.00557
Density MKS
2.00
0.50
0.90123±0.00566 0.90029±0.00568
Density MKS
1.00
0.90032±0.00564 0.90093±0.00514
0.90194±0.00548 0.90108±0.00518
Density MKS
0.80
Density MKS
0.90167±0.00545 0.90185±0.00561
Density MKS
0.50
0.30
0.90099±0.00557 0.90083±0.00573
Density MKS
0.30
0.90019±0.00543 0.89840±0.00473
0.90129±0.00527 0.89964±0.00529
Density MKS
0.10
Density MKS
0.89872±0.00522 0.89065±0.00426
Density MKS
0.10
Coverage
Method
ℓ
WSC
0.98462±0.00204 0.98471±0.00299
0.98452±0.00206 0.98414±0.00239
0.98452±0.00194 0.98424±0.00277
0.98448±0.00194 0.98400±0.00313
0.98448±0.00202 0.98357±0.00363
0.98414±0.00210 0.98276±0.00419
0.98329±0.00347 0.97986±0.00565
0.95757±0.00961 0.97686±0.00520
0.95724±0.00948 0.97871±0.00450
0.95786±0.00875 0.97905±0.00490
0.96014±0.00841 0.97943±0.00523
0.97052±0.00691 0.97871±0.00471
0.97895±0.00520 0.97929±0.00457
0.97790±0.00495 0.97267±0.00552
C.4.3
Full Results
In Tables 5–13, we report comprehensive comparison as described in on coverage, volume, and worst-slab-coverage(WSC) for all datasets, dimensions, coverage levels, and models as described in C.4.1. The MKS achieves the smallest volumes across all configurations, confirming the volume ratio results in Table 2. The Bonferroni method produces the largest regions in all cases. The density-only score and Mahalanobis ellipsoid are closer in performance, with the relative ordering depending on the dataset and dimension. Coverage is maintained at the nominal level across all methods, consistent with the distribution-free guarantee. WSC values are comparable across methods, with the MKS showing slightly lower WSC in some configurations at low α, consistent with the non-linear, anisotropic nature of the score as mentioned in Section 5.1.
Bio d = 2
Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89986±0.00385 0.89950±0.00401 0.89990±0.00389 0.89985±0.00490
6.75296±1.71711 5.56503±1.32203 4.71097±0.12988 2.75441±0.19967
0.81252±0.01168 0.80202±0.01310 0.79870±0.01041 0.73552±0.01274
0.89962±0.00312 0.89976±0.00300 0.89997±0.00324 0.89979±0.00359
6.05152±1.59344 4.94924±1.23485 3.75990±0.11716 2.11502±0.13028
0.84597±0.00835 0.83852±0.00989 0.83146±0.00772 0.77817±0.01177
0.05
Bonferroni Mahalanobis Density MKS
0.95008±0.00326 0.95019±0.00324 0.95026±0.00321 0.94992±0.00310
9.57435±2.43364 7.90659±1.85625 6.68404±0.21344 4.07637±0.29597
0.88893±0.00880 0.88266±0.00966 0.88062±0.00822 0.84279±0.00870
0.95015±0.00243 0.95031±0.00245 0.95023±0.00245 0.94956±0.00295
9.29171±2.43219 7.58293±1.89455 5.66443±0.17985 3.37645±0.22497
0.90848±0.00725 0.90357±0.00808 0.89593±0.00618 0.86161±0.00890
0.02
Bonferroni Mahalanobis Density MKS
0.97993±0.00183 0.97980±0.00206 0.98002±0.00177 0.97992±0.00170
14.27852±3.48611 11.79345±2.73223 9.86869±0.40005 6.68368±0.53531
0.94127±0.00637 0.93820±0.00694 0.93739±0.00534 0.92408±0.00584
0.97980±0.00205 0.97969±0.00215 0.97967±0.00208 0.97959±0.00209
14.27309±3.61125 11.63511±2.78044 8.65577±0.36513 6.15941±0.53668
0.95159±0.00629 0.94918±0.00667 0.94352±0.00486 0.92888±0.00639
0.01
Bonferroni Mahalanobis Density MKS
0.98993±0.00153 0.98992±0.00146 0.98987±0.00137 0.98964±0.00142
19.48478±4.20208 16.26672±3.57259 13.69825±0.83963 10.84783±1.19679
0.96331±0.00543 0.96212±0.00510 0.96087±0.00435 0.95571±0.00557
0.98972±0.00130 0.98971±0.00133 0.98985±0.00128 0.98950±0.00143
19.00227±4.48698 15.54972±3.52471 11.89990±0.65938 9.86605±0.98616
0.96973±0.00470 0.96830±0.00495 0.96411±0.00373 0.95859±0.00530
Table 5: Bio dataset (d = 2), Ridge and MLP. Mean ± std over 50 seeds.
Bio d = 3
Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89969±0.00359 0.89976±0.00408 0.89927±0.00428 0.89948±0.00386
30.02751±7.69485 21.64955±5.20752 20.09597±0.46127 10.28282±0.65734
0.67521±0.00981 0.68549±0.01119 0.70991±0.01125 0.67072±0.01096
0.89958±0.00436 0.90004±0.00456 0.89978±0.00375 0.89990±0.00477
22.12253±6.06727 15.76829±4.13603 14.13009±0.63231 6.94000±0.50385
0.76461±0.00971 0.76211±0.01072 0.77441±0.00882 0.73332±0.01151
0.05
Bonferroni Mahalanobis Density MKS
0.94966±0.00285 0.94974±0.00294 0.94968±0.00306 0.94977±0.00289
50.02735±12.56544 34.54247±8.05692 32.14227±0.88929 17.26134±1.16080
0.82107±0.00974 0.82375±0.00997 0.83499±0.00907 0.81059±0.00961
0.95002±0.00281 0.94974±0.00282 0.94996±0.00319 0.94930±0.00336
39.60544±10.42029 27.58456±6.91951 24.52894±1.10924 12.92087±0.99177
0.85878±0.00759 0.85726±0.00830 0.86897±0.00833 0.84026±0.00989
0.02
Bonferroni Mahalanobis Density MKS
0.97968±0.00174 0.97961±0.00188 0.97978±0.00199 0.97971±0.00199
95.17278±19.23166 64.71753±14.46727 60.32933±2.52745 28.61148±1.85680
0.91273±0.00769 0.91248±0.00867 0.91613±0.00825 0.91550±0.00821
0.97956±0.00195 0.97976±0.00198 0.97970±0.00189 0.97928±0.00196
72.06027±17.85519 48.71245±11.87846 43.01599±1.92057 25.45658±2.25113
0.92881±0.00710 0.92910±0.00728 0.93637±0.00599 0.92196±0.00684
0.01
Bonferroni Mahalanobis Density MKS
0.99004±0.00129 0.98991±0.00128 0.98974±0.00142 0.98964±0.00145
145.03113±29.83276 98.86763±23.09573 78.77614±3.28850 46.60653±4.47993
0.95426±0.00582 0.95499±0.00576 0.95464±0.00611 0.95376±0.00657
0.98957±0.00140 0.98989±0.00145 0.98984±0.00131 0.98935±0.00142
104.90286±24.88488 70.63898±15.87836 63.30613±3.88502 41.29964±3.65368
0.95867±0.00562 0.96033±0.00533 0.96347±0.00413 0.95599±0.00579
Table 6: Bio dataset (d = 3), Ridge and MLP. Mean ± std over 50 seeds. 32
Bio d = 4 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.90019±0.00385 0.89963±0.00414 0.89985±0.00383 0.89964±0.00428
70.80023±18.38201 40.55173±9.50937 53.79147±2.03691 21.05211±1.76196
0.69807±0.01069 0.70700±0.01138 0.71142±0.01091 0.68528±0.01304
0.90004±0.00429 0.89924±0.00403 0.89975±0.00400 0.89942±0.00388
50.35894±12.53578 29.76705±7.14846 37.15305±1.66572 14.55809±1.01036
0.76343±0.01010 0.75579±0.01188 0.77225±0.00972 0.73253±0.01157
0.05
Bonferroni Mahalanobis Density MKS
0.95005±0.00293 0.95003±0.00282 0.95000±0.00294 0.94953±0.00315
137.83500±33.72638 80.39881±18.64258 98.82912±4.29079 40.72004±3.35456
0.84677±0.00977 0.84152±0.00974 0.84042±0.00885 0.83254±0.01118
0.94963±0.00267 0.94988±0.00248 0.95011±0.00271 0.94968±0.00259
113.04785±27.63441 61.30482±14.55846 75.22207±3.48944 31.23205±2.41832
0.87175±0.00829 0.86802±0.00810 0.87465±0.00722 0.85159±0.00852
0.02
Bonferroni Mahalanobis Density MKS
0.98017±0.00166 0.98027±0.00187 0.97990±0.00197 0.97959±0.00182
404.97876±85.56561 185.19879±42.95393 209.05947±9.75432 84.66618±7.39257
0.93408±0.00619 0.93498±0.00660 0.92055±0.00739 0.92597±0.00634
0.98003±0.00219 0.98009±0.00198 0.97976±0.00169 0.97950±0.00179
264.83531±58.95916 134.24717±27.35345 151.93896±8.83569 73.15290±6.05987
0.94730±0.00679 0.94626±0.00627 0.93883±0.00584 0.92894±0.00590
0.01
Bonferroni Mahalanobis Density MKS
0.99001±0.00125 0.99004±0.00120 0.98999±0.00134 0.98959±0.00131
699.11904±116.42029 417.42339±72.95149 312.60368±24.06309 171.51322±18.46514
0.96915±0.00438 0.96695±0.00432 0.95750±0.00550 0.95725±0.00512
0.99008±0.00121 0.99001±0.00134 0.98999±0.00131 0.98944±0.00126
580.56035±99.22350 313.47197±53.64378 258.95104±16.59086 136.00975±14.86120
0.97358±0.00402 0.97235±0.00440 0.96384±0.00505 0.95857±0.00443
Table 7: Bio dataset (d = 4), Ridge and MLP. Mean ± std over 50 seeds. House d = 2 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89993±0.00595 0.89952±0.00530 0.89996±0.00538 0.89979±0.00602
6.91293±0.28132 6.07404±0.22147 7.13803±0.15414 5.38247±0.14907
0.78255±0.01437 0.75950±0.01477 0.78087±0.01371 0.74503±0.01512
0.90048±0.00563 0.90037±0.00540 0.89953±0.00570 0.90022±0.00594
4.85886±0.20871 4.23579±0.17488 5.00007±0.18663 3.94353±0.13939
0.78705±0.01187 0.77534±0.01278 0.80346±0.01297 0.76441±0.01475
0.05
Bonferroni Mahalanobis Density MKS
0.95030±0.00383 0.95010±0.00400 0.95010±0.00375 0.94981±0.00421
9.56767±0.36476 8.52586±0.31740 9.51130±0.25101 7.80109±0.25569
0.85745±0.01124 0.85127±0.01139 0.86766±0.00998 0.85082±0.01181
0.95016±0.00444 0.95010±0.00408 0.95013±0.00442 0.94990±0.00385
7.65780±0.37194 6.57688±0.30881 7.59197±0.29353 6.11101±0.25413
0.86862±0.01146 0.86139±0.01146 0.88396±0.01036 0.85647±0.01127
0.02
Bonferroni Mahalanobis Density MKS
0.98026±0.00223 0.98024±0.00238 0.98030±0.00227 0.97976±0.00215
15.72499±0.99432 13.41890±0.88564 13.80108±0.50327 11.28610±0.45260
0.91726±0.00871 0.91451±0.00915 0.92553±0.00724 0.91900±0.00775
0.97997±0.00315 0.97996±0.00319 0.98034±0.00320 0.97936±0.00311
13.11495±0.88412 10.91229±0.71338 12.01604±0.72647 9.74079±0.60296
0.93058±0.01023 0.92711±0.01072 0.94152±0.00978 0.92674±0.01104
0.01
Bonferroni Mahalanobis Density MKS
0.99023±0.00186 0.99019±0.00182 0.99012±0.00186 0.99006±0.00165
23.17162±2.19761 22.64279±2.07488 18.82786±1.58660 15.54018±0.86849
0.95354±0.00817 0.95310±0.00821 0.95399±0.00804 0.95369±0.00742
0.98970±0.00217 0.98985±0.00204 0.99004±0.00200 0.98921±0.00236
19.04634±1.47191 16.09192±1.18692 15.96740±1.10450 13.22862±0.96917
0.95680±0.00854 0.95637±0.00836 0.96340±0.00659 0.95604±0.00929
Table 8: House dataset (d = 2), Ridge and MLP. Mean ± std over 50 seeds. House d = 3 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.90013±0.00490 0.90021±0.00563 0.90089±0.00512 0.89991±0.00567
20.89480±0.90285 16.67193±0.70435 22.53800±0.62636 15.84239±0.60545
0.73449±0.01453 0.72109±0.01534 0.76266±0.01312 0.72481±0.01641
0.90077±0.00589 0.90035±0.00599 0.90030±0.00598 0.90028±0.00629
14.83318±0.67728 11.08563±0.53872 14.52913±0.69473 10.50678±0.50450
0.75021±0.01558 0.73265±0.01596 0.77171±0.01571 0.73014±0.01627
0.05
Bonferroni Mahalanobis Density MKS
0.95054±0.00408 0.94991±0.00441 0.95049±0.00422 0.94960±0.00434
38.14848±2.19057 27.94215±1.59993 35.22317±1.29551 26.57591±1.29750
0.83456±0.01205 0.82241±0.01379 0.85295±0.01165 0.83117±0.01286
0.94946±0.00422 0.94914±0.00501 0.94982±0.00514 0.94875±0.00526
28.64170±1.67484 20.68068±1.31609 25.76468±1.46887 19.24210±1.18160
0.84675±0.01358 0.83469±0.01631 0.86488±0.01494 0.83735±0.01699
0.02
Bonferroni Mahalanobis Density MKS
0.98014±0.00280 0.97994±0.00279 0.97990±0.00280 0.97934±0.00236
81.96302±7.32291 65.37222±5.80240 60.01819±4.23344 47.58553±2.89582
0.91972±0.01047 0.91636±0.01068 0.91883±0.00986 0.91340±0.00921
0.97973±0.00296 0.97963±0.00310 0.97977±0.00317 0.97851±0.00332
62.65107±5.10575 43.88950±3.24665 48.82813±3.55911 36.98856±2.49890
0.92557±0.01088 0.92057±0.01154 0.93380±0.00983 0.91987±0.01173
0.01
Bonferroni Mahalanobis Density MKS
0.99031±0.00185 0.99026±0.00194 0.99037±0.00174 0.98962±0.00194
130.27228±14.46327 122.34194±13.44742 96.07691±10.33696 72.65389±5.81531
0.95624±0.00797 0.95506±0.00890 0.95661±0.00774 0.95291±0.00884
0.98990±0.00212 0.98979±0.00236 0.99006±0.00223 0.98888±0.00251
107.87534±10.68322 75.67968±7.75971 73.40386±6.38579 56.93327±5.32941
0.95752±0.00792 0.95580±0.00912 0.96107±0.00799 0.95430±0.00941
Table 9: House dataset (d = 3), Ridge and MLP. Mean ± std over 50 seeds. 33
House d = 4 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.90069±0.00528 0.90010±0.00570 0.90057±0.00601 0.90038±0.00631
69.97675±9.20081 40.01080±5.05735 57.32679±3.24989 24.07322±2.10074
0.72485±0.01456 0.71051±0.01628 0.74052±0.01626 0.70942±0.01710
0.90026±0.00626 0.90015±0.00576 0.90017±0.00625 0.89955±0.00623
51.56111±7.14190 27.66793±3.69525 36.06235±2.69036 17.23644±1.62225
0.74011±0.01442 0.72300±0.01407 0.74855±0.01378 0.71253±0.01541
0.05
Bonferroni Mahalanobis Density MKS
0.95051±0.00435 0.95010±0.00421 0.95018±0.00414 0.94959±0.00407
193.30028±30.79762 106.26195±16.84056 119.77957±8.66023 62.53671±6.70189
0.82575±0.01439 0.82019±0.01434 0.83249±0.01247 0.81375±0.01331
0.94948±0.00434 0.94918±0.00444 0.95021±0.00434 0.94905±0.00436
139.61219±20.52921 71.80895±9.98681 85.05221±8.07272 48.71377±4.85942
0.83525±0.01353 0.82533±0.01561 0.84324±0.01191 0.82255±0.01389
0.02
Bonferroni Mahalanobis Density MKS
0.98049±0.00285 0.98028±0.00294 0.98030±0.00286 0.97894±0.00249
786.22762±153.66223 508.59046±102.86986 416.29345±87.02325 206.16876±21.16318
0.91915±0.01089 0.91617±0.01186 0.91554±0.01155 0.90918±0.01057
0.97970±0.00248 0.97986±0.00279 0.97986±0.00255 0.97811±0.00278
537.13765±87.67112 280.77555±42.99095 256.68070±39.72057 151.63365±20.04697
0.91950±0.00994 0.91802±0.01097 0.91704±0.00973 0.91029±0.01039
0.01
Bonferroni Mahalanobis Density MKS
0.99039±0.00167 0.99033±0.00162 0.99022±0.00166 0.98868±0.00193
1907.48674±305.65238 1537.53082±277.95149 1159.12814±240.79155 399.57582±53.94181
0.95710±0.00798 0.95693±0.00776 0.95515±0.00781 0.94975±0.00851
0.99014±0.00180 0.99016±0.00199 0.99011±0.00172 0.98807±0.00201
1342.70906±266.10634 784.01829±121.01162 731.24547±171.99074 302.25430±46.71982
0.95632±0.00806 0.95643±0.00875 0.95486±0.00798 0.94786±0.00886
Table 10: House dataset (d = 4), Ridge and MLP. Mean ± std over 50 seeds. Blog d = 2 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89900±0.00384 0.89890±0.00382 0.89921±0.00383 0.89959±0.00393
8.29337±0.71425 6.70800±0.56138 5.78854±0.11370 3.21050±0.11630
0.79786±0.00747 0.79437±0.00773 0.79100±0.00702 0.73308±0.01090
0.90000±0.00356 0.89995±0.00347 0.89956±0.00329 0.89923±0.00344
3.47217±0.91605 2.94403±0.92669 1.85948±0.63569 1.72493±0.46442
0.77852±0.01429 0.77172±0.01355 0.73680±0.01089 0.72938±0.01147
0.05
Bonferroni Mahalanobis Density MKS
0.94936±0.00265 0.94943±0.00277 0.94958±0.00261 0.94961±0.00281
11.21232±0.90834 9.11673±0.72328 8.09491±0.21440 5.79019±0.25816
0.86863±0.00671 0.86659±0.00694 0.86076±0.00642 0.83434±0.00840
0.94981±0.00215 0.94970±0.00235 0.94974±0.00251 0.94918±0.00239
6.12415±1.40650 5.19817±1.35799 3.84685±0.83740 3.57286±0.67226
0.86334±0.00906 0.86041±0.00931 0.84318±0.00901 0.84186±0.00925
0.02
Bonferroni Mahalanobis Density MKS
0.98002±0.00164 0.98003±0.00176 0.98002±0.00179 0.97970±0.00164
21.29593±1.74039 18.95573±1.75739 20.76735±2.26818 14.69500±0.81123
0.91759±0.00638 0.91696±0.00684 0.91713±0.00707 0.91565±0.00689
0.97995±0.00166 0.98003±0.00160 0.97995±0.00151 0.97900±0.00158
14.60027±3.90896 12.36117±3.34895 11.00846±1.63998 8.91418±1.39787
0.93472±0.00679 0.93388±0.00654 0.92591±0.00553 0.92462±0.00599
0.01
Bonferroni Mahalanobis Density MKS
0.99019±0.00139 0.99023±0.00135 0.99022±0.00134 0.98958±0.00143
50.76856±4.32594 43.78691±3.35545 42.43467±3.28301 25.94511±2.10044
0.95818±0.00583 0.95839±0.00575 0.95847±0.00576 0.95537±0.00585
0.98987±0.00137 0.98985±0.00139 0.98972±0.00114 0.98861±0.00119
29.20655±10.05672 24.65639±10.89113 19.10084±3.59817 15.73014±2.25929
0.96254±0.00535 0.96218±0.00535 0.95862±0.00432 0.95597±0.00433
Table 11: Blog dataset (d = 2), Ridge and MLP. Mean ± std over 50 seeds. Blog d = 3 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89984±0.00353 0.89977±0.00372 0.89968±0.00359 0.89933±0.00348
19.63098±1.55911 13.63895±1.07565 16.07511±0.48056 8.03864±0.37649
0.76831±0.00891 0.75466±0.00797 0.78485±0.00787 0.72723±0.00856
0.90055±0.00367 0.90025±0.00370 0.89956±0.00367 0.89938±0.00346
8.28306±3.05518 6.13269±2.67363 4.23297±2.01865 3.64035±1.35910
0.76843±0.00845 0.76148±0.00806 0.74583±0.00907 0.73693±0.00871
0.05
Bonferroni Mahalanobis Density MKS
0.94987±0.00276 0.95008±0.00261 0.94963±0.00247 0.94975±0.00252
38.53108±2.85214 27.62096±1.99581 27.63325±1.06093 18.09991±0.86087
0.84566±0.00824 0.84444±0.00703 0.85206±0.00647 0.83616±0.00710
0.95009±0.00251 0.95009±0.00259 0.95009±0.00252 0.94935±0.00231
20.41023±6.22721 14.55600±5.56613 11.13314±3.44430 8.90112±2.38038
0.86677±0.00689 0.86376±0.00791 0.85273±0.00678 0.84755±0.00688
0.02
Bonferroni Mahalanobis Density MKS
0.98027±0.00171 0.98023±0.00175 0.98012±0.00185 0.97972±0.00188
160.66764±15.42586 101.09682±8.13264 101.30861±10.72005 51.60622±3.15894
0.93024±0.00611 0.92932±0.00661 0.92500±0.00649 0.92344±0.00685
0.98001±0.00176 0.98004±0.00185 0.98015±0.00164 0.97908±0.00183
68.89382±17.21461 45.51245±15.14107 48.51137±8.59343 25.64780±5.57789
0.93707±0.00564 0.93601±0.00625 0.93035±0.00554 0.92722±0.00603
0.01
Bonferroni Mahalanobis Density MKS
0.99002±0.00125 0.98999±0.00130 0.99006±0.00133 0.98923±0.00144
394.07433±57.18158 291.25393±39.44550 243.97121±30.99577 96.68725±6.64635
0.96347±0.00424 0.96338±0.00440 0.96193±0.00466 0.95676±0.00517
0.98984±0.00124 0.98983±0.00113 0.98992±0.00125 0.98833±0.00115
182.68447±52.56609 114.26591±42.88878 138.87714±24.61825 48.83056±8.89393
0.96425±0.00421 0.96383±0.00398 0.96072±0.00415 0.95644±0.00397
Table 12: Blog dataset (d = 3), Ridge and MLP. Mean ± std over 50 seeds. 34
Blog d = 4 Ridge α
Method
MLP
Coverage
Volume
WSC
Coverage
Volume
WSC
0.1
Bonferroni Mahalanobis Density MKS
0.89987±0.00346 0.89963±0.00338 0.89956±0.00339 0.89901±0.00338
36.02848±9.77178 16.59707±3.91954 31.92792±3.14113 4.83457±1.03373
0.76808±0.00700 0.75219±0.00718 0.77937±0.00756 0.73287±0.00914
0.89920±0.00338 0.89892±0.00360 0.89947±0.00356 0.89880±0.00364
21.80796±48.99889 9.39265±12.75505 4.55082±2.35156 2.48198±1.36772
0.76238±0.01124 0.75204±0.01125 0.74327±0.01059 0.73215±0.01123
0.05
Bonferroni Mahalanobis Density MKS
0.94987±0.00258 0.95014±0.00258 0.94955±0.00252 0.94956±0.00226
93.33837±22.31949 45.13031±10.62577 67.15577±5.41781 15.31169±2.50086
0.84815±0.00700 0.84679±0.00697 0.85655±0.00618 0.84429±0.00669
0.94953±0.00310 0.94954±0.00298 0.94985±0.00272 0.94907±0.00253
67.17839±121.40008 29.75976±32.27370 26.04943±45.06430 10.04471±7.00155
0.86725±0.00873 0.86347±0.00866 0.85655±0.00752 0.84979±0.00698
0.02
Bonferroni Mahalanobis Density MKS
0.98017±0.00190 0.98021±0.00191 0.98015±0.00188 0.97962±0.00189
717.27415±204.60465 283.73237±78.54986 400.26628±59.40722 68.65380±8.74056
0.93324±0.00619 0.93305±0.00611 0.92666±0.00653 0.92734±0.00669
0.97970±0.00186 0.97978±0.00181 0.97997±0.00205 0.97868±0.00210
357.88343±423.84731 156.24391±150.82879 194.34181±101.36700 54.43470±66.02764
0.93894±0.00580 0.93751±0.00544 0.93142±0.00614 0.92820±0.00661
0.01
Bonferroni Mahalanobis Density MKS
0.99002±0.00150 0.99004±0.00147 0.99012±0.00150 0.98941±0.00138
2768.70057±841.15401 1358.22839±398.51791 1516.02620±324.13823 179.08912±18.70267
0.96479±0.00497 0.96476±0.00495 0.96269±0.00535 0.95993±0.00460
0.98976±0.00128 0.98989±0.00135 0.98993±0.00136 0.98810±0.00132
1405.88576±1340.47362 603.26515±606.47872 969.01006±397.03194 118.25062±84.42974
0.96556±0.00442 0.96521±0.00467 0.96211±0.00457 0.95776±0.00483
Table 13: Blog dataset (d = 4), Ridge and MLP. Mean ± std over 50 seeds.
35