arXiv:2605.12341v1 [stat.ML] 12 May 2026
Multi-Variable Conformal Prediction: Optimizing Prediction Sets without Data Splitting
Laura Lützow TU Munich & MCML, Germany [email protected] Marco C. Campi University of Brescia, Italy [email protected]
Simone Garatti Politecnico di Milano, Italy [email protected]
Lars Lindemann ETH Zürich, Switzerland [email protected]
Matthias Althoff TU Munich & MCML, Germany [email protected]
Abstract Conformal prediction constructs prediction sets with finite-sample coverage guarantees, but its calibration stage is structurally constrained to a scalar score function and a single threshold variable — forcing shapes of prediction sets to be fixed before calibration, typically through data splitting. We introduce multi-variable conformal prediction (MCP), a framework that extends conformal prediction to vector-valued score functions with multiple simultaneous calibration variables. Building on scenario theory as a principled framework for certifying data-driven decisions, MCP unifies prediction set design and calibration into a single optimization problem, eliminating data splitting without sacrificing coverage guarantees. We propose two computationally efficient variants: RemMCP, grounded in constrained optimization with constraint removal, which admits a clean generalization of split conformal prediction; and RelMCP, based on iterative optimization with constraint relaxation, which supports non-convex score functions at the cost of possibly greater conservatism. Through numerical experiments on ellipsoidal and multi-modal prediction sets, we demonstrate that RemMCP and RelMCP consistently meet the target coverage with prediction set sizes smaller than or comparable to those of baselines with data split, while considerably reducing variance across calibration runs — a direct consequence of using all available data for shape optimization and calibration simultaneously.
1
Introduction
Reliable uncertainty quantification is essential for deploying machine learning models in high-stakes settings, such as scientific inference, autonomous systems, and safety-critical decision-making. Beyond point predictions, these settings demand prediction sets — regions guaranteed to contain the true output with a user-specified probability — that are valid under minimal distributional assumptions. Conformal prediction has emerged as the leading framework for this purpose: given any pre-trained model and a held-out calibration set, it constructs prediction sets with finite-sample marginal coverage guarantees under the mild assumption of exchangeable data [Vovk et al., 1999, Lei et al., 2013, Angelopoulos and Bates, 2022]. Despite its appeal, classical conformal prediction is structurally limited. The calibration stage involves a scalar nonconformity score function and a single scalar nonconformity threshold. This architectural choice forces prediction sets into rigid, often overly conservative shapes that cannot adapt to the geometry of complex outputs or error structures. The straightforward extension to multi-dimensional outputs calibrates a separate threshold per output dimension and corrects the marginal coverage Preprint.
Standard Split Conformal Prediction
Data
Optimization of the prediction set
Subset 1
Multi-Variable Conformal Prediction Joint optimization and calibration over parameters q ∗ ∈ Rnq
Calibration of threshold q̃ ∗ ∈ R
Subset 2
Prediction set Ỹ(x, q̃ ∗ ) = {y | s̃(x, y) ≤ q̃ ∗ }
Target coverage 1 − ε
High variance
Prediction set Y(x, q ∗ ) = {y | s(x, y, q ∗ ) ≤ 0}
Low variance
Coverage probability
Figure 1: Comparison of standard split conformal prediction (SCP) and multi-variable conformal prediction (MCP). SCP decouples the optimization of the prediction-set shape and threshold calibration by splitting the available data into two subsets, thereby reducing the effective sample size available to each stage. MCP eliminates this split by jointly optimizing and calibrating multiple shape parameters over the full dataset. Both approaches achieve the target coverage 1 − ϵ in expectation, but MCP yields a narrower coverage distribution by using all available data for calibration.
using the Bonferroni method [Stankevičiūtė et al., 2021], resulting in hyperrectangular prediction sets that are often highly conservative [Feldman et al., 2023]. Conservatism can be partially reduced by introducing dimension-specific normalization factors fitted from additional data, which adapt the side lengths of the hyperrectangle to the scale of the residuals [Cleaveland et al., 2024]. Other methods reduce conservatism by modeling output dependencies through copulas [Messoudi et al., 2021, Sun and Yu, 2024] or graph-based architectures [Cini et al., 2025], or by adopting more flexible shapes of prediction sets such as ellipsoids fitted from additional data [Johnstone and Cox, 2021, Messoudi et al., 2022, Xu et al., 2024], zonotopes [Gray et al., 2025], norm balls [Braun et al., 2025], multi-modal regions constructed as unions of parametrized score functions [Tumu et al., 2024], and nonconvex sets obtained via normalizing flows [Fang et al., 2025, Luo and Zhou, 2025] or variational autoencoders [Feldman et al., 2023]. However, all of these approaches share a common limitation: the shape of the prediction set is fixed before calibration, decoupling the two stages and requiring separate data splits, which reduces data efficiency and inflates variance across calibration runs. In parallel, scenario theory has emerged in control and optimization as a framework for certifying data-driven decisions with finite-sample guarantees [Calafiore and Campi, 2006, Garatti and Campi, 2022]. While its guarantees are structurally reminiscent of conformal coverage bounds, they arise from a fundamentally different viewpoint based on randomized constraint satisfaction. Scenariobased methods have been used to construct prediction sets in specific settings, including interval predictors [Campi et al., 2009, 2015, Garatti et al., 2019, Sadeghi et al., 2019], rectangular sets for multivariate data [de Angelis et al., 2021], and nonconvex reachable sets in data-driven reachability analysis [Dietrich et al., 2024]. Recent work has begun to formalize connections between scenario theory and conformal prediction [Coppola et al., 2024, Lin and Bansal, 2024, O’Sullivan et al., 2025], which was exploited by zono-conformal prediction [Lützow et al., 2025], where multiple variables are calibrated simultaneously to construct zonotopic prediction sets. Yet this remains a specific instantiation: a general framework that provides coverage guarantees and admits arbitrary vector-valued score functions and multiple calibration variables has not been developed. We close this gap with multi-variable conformal prediction (MCP), a framework that extends the conformal calibration stage to accommodate vector-valued score functions and multiple simultaneous optimization variables. Unlike existing approaches, MCP couples prediction set design and calibration into a single optimization problem, as illustrated in Fig. 1 — allowing the shape of sets to adapt directly to the available data and enabling joint calibration of multiple parameters without sacrificing distribution-free coverage guarantees or requiring data splitting. Our main contributions are: • We introduce MCP, a general framework that extends conformal calibration to vector-valued score functions and multiple simultaneous optimization variables, eliminating the need for data splitting for optimizing the shape of prediction sets (Sec. 3). 2
• We develop RemMCP (Sec. 3.1), a variant based on scenario optimization with constraint removal that recovers split conformal prediction as a special case and admits finite-sample expected-coverage and high-confidence coverage guarantees. • We develop RelMCP (Sec. 3.2), a variant based on unconstrained optimization with constraint relaxation that supports non-convex score functions and achieves target high-confidence coverage guarantees via an iterative calibration algorithm grounded in non-convex scenario theory. • We provide guidance on the design of score and cost functions to trade off expressiveness, computational complexity, and conservatism (Sec. 3.3). • Through numerical experiments, we show that MCP yields prediction sets that are smaller than or comparable to split conformal prediction baselines while consistently meeting the target coverage and with considerably lower variance across calibration runs (Sec. 4). Together, these results establish MCP as a flexible and theoretically grounded extension of conformal inference, bridging ideas from split conformal prediction and scenario-based optimization.
2
Preliminaries
This section introduces the two methodological pillars of the proposed framework. We first review scenario theory and its finite-sample guarantees for data-driven optimization problems, which form the theoretical backbone of our approach. We then revisit split conformal prediction, presenting it in a formulation that extends naturally to the multi-variable setting developed in the subsequent section. 2.1
Scenario Theory for Data-Driven Optimization
Scenario theory provides a principled framework for optimization under uncertainty. Originally introduced in Calafiore and Campi [2005, 2006], Campi and Garatti [2008] as an algorithmic tool to robust design, it has since evolved into a comprehensive methodology for data-driven decisionmaking; see Campi et al. [2021] for a recent survey. We focus on two schemes most relevant to MCP: optimization with constraint removal and optimization with constraint relaxation. Let q ∈ Rnq denote cal the decision variable and let Mcal = {δ (m) }nm=1 ⊂ ∆ncal be i.i.d. samples drawn from an unknown distribution over the sample space ∆. Each sample δ induces a constraint q ∈ Qδ , while preferences ∗ over feasible solutions are encoded through a cost function J(q). Since the optimizer qM depends cal on random samples, it is itself a random variable whosereliability is characterized by the violation ∗ ∗ probability V (qM ) = P δ ∈ ∆ : qM ̸∈ Qδ | Mcal , which is the probability that a previously cal cal unseen sample violates the learned constraint. Optimization with constraint removal. To balance solution quality and conservatism, we allow the optimizer to discard a subset R ⊂ Mcal of sampled constraints, yielding the scenario program ∗ qM = argmin J(q) cal q∈Rnq
s.t. q ∈ Qδ , ∀δ (m) ∈ Mcal \ R.
(1)
This class of schemes is studied in [Campi and Garatti, 2011]. Among the many admissible removal schemes, we adopt the one from [Romao et al., 2023a,b], which is especially well-suited to MCP ∗ and admits particularly tight guarantees. Let supp(qM ) denote the support set of the optimizer, cal ∗ which contains the samples δ ∈ Mcal whose removal would change qM . Given a removal budget cal ρ ∈ N0 , the set R = Rρ is constructed iteratively via [Romao et al., 2023a,b] ∅ if k = 0, Rk = (2) ∗ Rk−1 ∪ supp(qk−1 ) ∪ Zk−1 otherwise, ∗ where qk∗ is the optimizer at stage k for Mcal \Rk , and Zk−1 consists of nq −|supp(qk−1 )| additional ∗ scenarios from Mcal \ (Rk−1 ∪ supp(qk−1 )), selected via a lexicographic tie-breaking rule to ensure |Rk | = knq at each stage. The guarantees for this scheme rely on three assumptions.
Assumption 1. For any ncal ∈ N and Mcal ⊂ ∆ncal , the optimal solution of (1) exists and is unique. Assumption 2. For each k = 0, . . . , ρ, the scenario program in (1) over Mcal \ Rk is nondegenerate with probability one: solving (1) using only supp(qk∗ ) yields the same solution as using Mcal \ Rk . Assumption 3. The objective J is convex in q and Qδ is a convex set for all δ ∈ ∆. 3
Uniqueness in Assumption 1 can be relaxed via a convex tie-breaking rule. Existence and convexity (Assumptions 1 and 3) can be enforced by construction through appropriate choices of the cost function and constraint sets. Under convexity, the nondegeneracy condition in Assumption 2 is mild: it simply requires that the constraints do not have accumulation points [Romao et al., 2023b]. In practice, these assumptions are straightforward to satisfy. Under Assumptions 1 to 3, the optimizer qρ∗ satisfies [Romao et al., 2023b, Thm. 3] P V (qρ∗ ) ≤ ε ≥ 1 −
nq (ρ+1)−1
X j=0
ncal j ε (1 − ε)ncal −j , j
(3)
bounding the confidence with which the violation probability does not exceed ε. The expected violation probability satisfies the complementary bound [O’Sullivan et al., 2025, Thm. 3] nq (ρ + 1) E V (qρ∗ ) ≤ , (4) ncal + 1 which shows that the average violation probability over random draws of Mcal decreases as the calibration set grows or the removal budget shrinks. Notably, (4) continues to hold when the samples in Mcal are not i.i.d. but merely exchangeable, since the proof relies only on the invariance of the distribution to permutations of the observations. The following two additional assumptions are needed for the bound in (3) to hold with equality [Romao et al., 2023b, Thm. 5]: Assumption 4. For each k = 0, . . . , ρ, the scenario program in (1) over Mcal \ Rk is fully ∗ supported with probability one: for any Mcal ⊂ ∆ncal with ncal > nq , the support set supp(qM ) cal has cardinality nq with probability one. Assumption 5. For any Mcal ⊂ ∆ncal and M̃cal ⊂ Mcal , for any k = 0, . . . , ρ and δ ∈ supp(q̃k ), Sk−1 we have q ∗ (M) ∈ / Qδ for any M ⊂ M̃cal \ ( j=0 supp(q̃j ) ∪ {δ}) with |M| = nq , where q̃k is the optimizer at stage k for M̃cal . The violation condition (Assumption 5) requires that each removed sample is genuinely violated by the solution obtained without it. Optimization with constraint relaxation. As an alternative to constraint removal, constraint relaxation [Garatti and Campi, 2025, Sec. 4.1] replaces hard constraints with penalized violations. This approach is also applicable to non-convex problems and can, in some cases, yield less conservative solutions, albeit with coverage guarantees that are evaluated a posteriori rather than specified a priori. We adopt this approach in Sec. 3.2 to derive the RelMCP formulation. 2.2
Split Conformal Prediction
Split conformal prediction (SCP), also known as inductive conformal prediction, is a distribution-free method for constructing prediction sets with finite-sample marginal coverage guarantees [Angelopoulos and Bates, 2022]. Let M be a dataset consisting of samples (x, y) from an unknown distribution over the data space ∆, which we randomly split into a training dataset Mtrain and a calibration cal dataset Mcal = {(x(m) , y (m) )}nm=1 ⊂ ∆ncal . SCP is based on the following procedure: 1. Training: (a) We fit a point predictor f : Rnx → Rny using Mtrain . (b) We define a score function s̃ : Rnx × Rny → R, which measures how poorly a candidate output y conforms to the prediction f (x). The structure of the score function determines the geometry of the prediction set, while the scale of the prediction sets is calibrated in the subsequent step to achieve the desired coverage. Optionally, the score function may be fitted using additional data [Johnstone and Cox, 2021, Feldman et al., 2023, Tumu et al., 2024]. 2. Calibration: (a) We compute the outlier budget for a target miscoverage level ε ∈ (0, 1) as ρ = ⌊ε(ncal + 1)⌋ − 1. ∗
(5)
(b) We set the conformal threshold q̃ ∈ R to the (ncal − ρ)-smallest nonconformity score in cal {s̃(x(m) , y (m) )}nm=1 , which is obtained by applying s̃ to the calibration set Mcal . This is equivalent to setting q̃ ∗ to the empirical (1 − ε)-quantile of the scores [Tumu et al., 2024]. 4
3. Prediction: The prediction set for a new input x is constructed via Ỹ(q̃ ∗ , x) = {y ∈ Rny | s̃(x, y) ≤ q̃ ∗ } . (6) Letting η̃(q̃ ∗ ) = P (x, y) ∈ ∆ : y ∈ Ỹ(q̃ ∗ , x) | Mcal denote the coverage probability of the prediction set for a new test point (x, y). Under exchangeability of the data points from Mcal and the test point, this construction guarantees finite-sample coverage [Angelopoulos et al., 2026, Sec. 4.1] E [η̃(q̃ ∗ )] ≥ 1 − ε,
(7)
where the expectation is over the calibration data. If the data points satisfy the slightly stronger condition of being i.i.d. and if the score distribution is continuous, the random variable η̃(q̃ ∗ ) follows a Beta distribution [Angelopoulos et al., 2026, Sec. 4.1]: η̃(q̃ ∗ ) ∼ Beta(ncal − ρ, ρ + 1),
(8)
which allows computing the confidence level 1 − β in the high-confidence coverage guarantee P(η̃(q̃ ∗ ) ≥ 1 − ε) ≥ 1 − β. In regression settings, a standard choice of nonconformity score is the p-norm of the residual, s̃(x, y) = ∥y − f (x)∥p , which induces prediction regions given by p-norm balls centered at f (x) [Braun et al., 2025]. Alternatively, one may employ dimension-wise scores s̃j (x, y) = y(j) − f(j) (x) ,
(9)
and calibrate a separate threshold for each dimension j = 1, . . . , ny . This yields axis-aligned rectangular prediction sets. Joint coverage can be ensured via a Bonferroni correction by calibrating each marginal model at level 1 − ε/ny [Stankevičiūtė et al., 2021]. Multi-modal prediction regions can be constructed using s̃(x, y) = minj s̃θj (x, y) , where s̃θi are score functions parameterized by θj and fitted on additional data [Tumu et al., 2024]. For example, unions of ellipsoids arise by defining s̃θj (x, y) = (y − f (x) − ĉj )⊤ Σ̂−1 j (y − f (x) − ĉj ), with fitted parameters θi = (Σ̂j , ĉj ), ny ×ny ny Σ̂j ≻ 0 ∈ R , ĉj ∈ R .
3
Multi-Variable Generalization
Classical SCP relies on a scalar score function s̃ and calibrates a single scalar threshold q̃ ∗ . We extend this framework by introducing a vector-valued score function s(x, y, q) ∈ Rns and a parameter vector q ∈ Rnq that determine whether an output y conforms to the prediction f (x) of a point predictor f . cal Given Mcal = {(x(m) , y (m) )}nm=1 ⊂ ∆ncal , the goal is to find q ∗ such that the prediction set Y(x, q ∗ ) = {y ∈ Rny | s(x, y, q ∗ ) ≤ 0}
(10)
satisfies a prescribed coverage guarantee on η(q ∗ ) = P ((x, y) ∈ ∆ : y ∈ Y(q ∗ , x) | Mcal ) while optimizing a user-defined objective J(q) that encodes preferences, such as minimal prediction set volume. For a target miscoverage level ε ∈ (0, 1), we consider the expected coverage guarantee E[η(q ∗ )] ≥ 1 − ε
(11)
P (η(q ∗ ) ≥ 1 − ε) ≥ 1 − β.
(12)
and the high-confidence guarantee
We develop two tractable formulations. The first is based on a constrained scenario program with constraint removal, which recovers SCP as a special case but requires the score function to be convex in q. The second relies on unconstrained optimization with constraint relaxation; this approach applies more broadly, including non-convex settings. Design considerations for the score and cost function are discussed in Section 3.3. 3.1
MCP via Constraint Removal
MCP based on constrained optimization with constraint removal (RemMCP) constructs prediction sets that satisfy the expected coverage guarantee in (11). It proceeds as follows: 1. Training: 5
(a) We fit a point predictor f : Rnx → Rny using Mtrain . (b) We define a convex score function s(x, y, q) ∈ Rns and a convex cost function J(q) ∈ R, which can optionally be fitted using additional data. 2. Calibration: (a) We compute the outlier budget as j ε(n + 1) k cal ρ= − 1. (13) nq If the right-hand side of (13) is smaller than 0, no feasible ρ ∈ N0 exists. In this case, we need to increase the size of Mcal or, alternatively, reduce the parameter dimension nq , which requires redefining the score and cost functions. (b) We compute the parameter q ∗ = qρ∗ by solving the cascading optimization problems qk∗ = argmin J(q) q∈Rnq
s.t. s(x(m) , y (m) , q) ≤ 0, ∀(x(m) , y (m) ) ∈ Mcal \ Rk ,
(14)
for k = 0, . . . , ρ, with Rk constructed via (2). 3. Prediction: The prediction set for a new input x is constructed via (10). This prediction set construction guarantees that every calibration data point y (m) in Mcal \ Rρ is covered by its respective prediction set Y(x(m) , q ∗ ). We obtain the following coverage guarantees: Theorem 1 (Probabilistic Guarantees of RemMCP). Let Assumptions 1 and 2 hold for the optimization problems in (14) and let Mcal consist of ncal exchangeable samples from an unknown distribution. Then the coverage of RemMCP satisfies the expected coverage guarantee in (11). If the samples in Mcal are i.i.d., the high-confidence guarantee in (12) holds with nq (ρ+1)−1
β=
X j=0
ncal j ε (1 − ε)ncal −j , j
(15)
which correspondsto the cumulative distribution function of a beta distribution I1−ε ncal − nq (ρ + 1) + 1, nq (ρ + 1) . Under the additional Assumptions 4 and 5, the coverage follows exactly η(q ∗ ) ∼ Beta ncal − nq (ρ + 1) + 1, nq (ρ + 1) . (16) We next show that SCP is a special cases of RemMCP and that RemMCP is less conservative than dimension-wise SCP. All proofs are provided in Appendix A. Corollary 1. The prediction sets (6) and coverage guarantees (7) and (8) of SCP with scalar score s̃(x, y) and continuous score distribution are recovered from RemMCP by setting s(x, y, q) = s̃(x, y) − q, J(q) = q. Corollary 2. The prediction sets obtained by RemMCP with −y + f (x) − q s(x, y, q) = , y − f (x) − q J(q) = 1⊤ q,
(17a) (17b)
(18a) (18b)
assuming a continuous score distribution, are contained in the prediction sets produced by the dimension-wise SCP method in (9) with Bonferroni correction [Stankevičiūtė et al., 2021] for any target miscoverage level ϵ ∈ (0, 1). 3.2
MCP via Constraint Relaxation
An alternative approach for MCP, which we name RelMCP, bases the calibration phase on unconstrained optimization via constraint relaxation while targeting the high-confidence guarantee in (12) for a prescribed confidence level 1 − β. RelMCP proceeds as follows. 1. Training: As in Sec. 3.1 but the score and cost function do not need to be convex in q. 6
2. Calibration: For a grid of penalty parameters ϕ1 , . . . , ϕneval > 0, which control the trade-off between objective optimality and empirical constraint satisfaction, we solve qi∗ = argmin J(q) + ϕi q∈Rnq
ncal X
ξm (q),
(19)
m=1
with ξm (q) = max 0, maxj s(j) (x(m) , y (m) , q) and i = 1, . . . , neval , where qi∗ need not be a global optimum. The achievable miscoverage εi for each qi∗ is evaluated a posteriori as the solution to j nX cal −1 1 β di = 1, (20) ncal neval ncal (1 − ε )ncal −j i di j=di
where di is the solution complexity [Garatti and Campi, 2025]; namely, the number of samples satisfying ξm (qi∗ ) > 0, plus the cardinality of the smallest additional subset of samples that, together with them, uniquely determines qi∗ . We then select as q ∗ the solution qi∗ corresponding to the smallest penalty ϕi for which the target coverage is achieved, i.e., εi ≤ ε. 3. Prediction: The prediction set for a new input x is constructed via (10) using the calibrated parameters q ∗ = qi∗ obtained from the calibration step. Theorem 2 (Probabilistic Guarantees of RelMCP). Let Mcal consist of ncal i.i.d. samples drawn from an unknown distribution and let q ∗ be a valid solution of RelMCP. Then RelMCP satisfies the high-confidence bound in (12). The achievable miscoverage εi in (20) depends on neval : a larger grid reduces the per-evaluation confidence budget β/neval , tightening the achievable coverage (1 − εi ) at each grid point. Since predefining a grid that achieves the target coverage without unnecessary conservatism is difficult in practice, we propose the adaptive procedure in Alg. 1, which iteratively searches for a suitable penalty parameter; a detailed description is provided in Appendix B. Because the penalty is adapted to the calibration data, the coverage guarantee of Thm. 2 does not apply strictly to this algorithm. Nevertheless, our experiments show that it consistently exceeds the target coverage in practice. Algorithm 1 RelMCP Calibration Inputs: Calibration data Mcal , score function s(x, y, q), cost function J(q), target miscoverage ε ∈ (0, 1), target misconfidence β ∈ (0, 1), initial penalty weight ϕ1 > 0, increase factor f+ > 1, decrease factor f− ∈ (0, 1), maximum number of iterations imax ∈ N 1: i ← 1; ival ← 0; iinv ← 0 # Initialize iteration counters 2: while i > 0 do 3: qi∗ , di ← solve (19) with penalty ϕi # Solve optimization problem 4: εi ← solve (20) with di and neval = i # Compute miscoverage level 5: if εi ≤ ε then # Current solution meets target coverage 6: ival ← i # Record as the best valid solution 7: else if ξm (qi∗ ) = 0 for all m = 1, . . . , ncal then # Coverage insufficient but all constraints satisfied 8: return None # Return since no valid solution was found 9: else # Target coverage not met 10: iinv ← i # Record as the best invalid solution 11: 12: 13: 14: 15:
16:
if ival ̸= 0 then # Valid solution was found previously ε′ival ← (20) using dival and neval = i + 1 # Check coverage for next iteration if ε′ival > ε or i ≥ imax then # Will not meet target coverage or maximum number of iterations reached return qi∗val , ϕival , dival , εival , ival # Return best valid solution max(1, i−imax ) if ival = 0 # Increase best invalid penalty weight ϕiinv · f+ ϕi+1 ← ϕival · f− if iinv = 0 # Decrease best valid penalty weight 1 otherwise # Set to midpoint between best penalty weights 2 (ϕival + ϕiinv ) i ← i + 1 # Next iteration
7
3.3
Practical Design Guidelines
The choice of score and cost functions governs the computational complexity of calibration, the tractability of test-time prediction, and the conservatism of the resulting prediction sets. We discuss the main design strategies here; concrete instantiations for rectangular, ellipsoidal, zonotopic, sublevel, and multi-modal prediction sets are provided in Appendix C. Score function design. The score function can be designed by specifying a parametrized family of prediction sets Y(x, q) for which set membership can be efficiently evaluated, and define the score function implicitly via s(x, y, q) = I{y − f (x) ∈ / Y(x, q)}, where I(·) is the indicator function. This yields a closed-form description of the set boundary. Alternatively, we can define the prediction set as the sublevel set Y(x, q) = {y | g(x, y, q1 ) ≤ q2 } of an arbitrary function g, with score function s(x, y, q) = g(x, y, q1 ) − q2 , simultaneously calibrating the sublevel threshold q2 and additional shape parameters q1 . This yields highly flexible prediction set shapes but does not produce a closedform set boundary. Unions of prediction sets can be obtained by combining score functions via the min-operator analogously to [Tumu et al., 2024], with the advantage that all component parameters are optimized and calibrated jointly without requiring an additional data split. Cost function design. The cost function typically penalizes prediction set volume, though other objectives are equally admissible. When a closed-form proxy is unavailable—as is the case for volume over general sublevel set functions—Monte Carlo estimates offer a practical alternative. Concretely, we can optimize theP empirical quality of the prediction set averaged over a representative set of inputs X , i.e., J(q) = x∈X ℓ Y(x, q) , where ℓ(·) denotes an appropriate performance measure such as the volume, diameter, or a surrogate thereof. Connection to existing frameworks. SCP [Vovk, 2013], interval predictor models [Campi et al., 2009], zono-conformal prediction [Lützow et al., 2025], reachset-conformant identification [Liu et al., 2023, Lützow and Althoff, 2026], and data-driven reachability analysis [Dietrich et al., 2024] are all special cases of MCP. The score or constraint functions and the cost functions developed independently across these communities therefore serve as a valuable source of inspiration for novel MCP instantiations.
4
Numerical Experiments
We empirically evaluate MCP against standard SCP baselines using the vehicle prediction benchmark from [Tumu et al., 2024, Sec. 4.2]. The point predictor f (x) maps the vehicle state over the past 5 time steps, x, to the vehicle position T time steps into the future, y (time-step size 0.1 s). For each run we randomly generate ncal = 2,000 calibration points for a vehicle at an intersection and target expected coverage of at least 1 − ε with ε = 0.05. Additional experimental details and results are provided in Appendices D and E. We consider the following methods for generating ellipsoidal and multi-modal prediction sets: SCP-1: Standard SCP with score s̃(x, y) = ∥y − f (x)∥2 . This score yields spherical prediction sets and serves as the geometry-agnostic baseline. SCP-2: SCP in which the ncal calibration points are split: a held-out subset of n′cal = 14 ncal points is reserved for threshold calibration, while the remainder is used to estimate A) the residual covariance Σ̂ of an ellipsoidal prediction set, via the score s̃(x, y) = (y − f (x))⊤ Σ̂−1 (y − f (x)) [Johnstone and Cox, 2021, Messoudi et al., 2022], or B) ellipsoid centers ĉi via k-means and per-cluster sample covariances Σ̂i of a union-of-ellipsoids prediction region, via the score s̃(x, y) = mini∈{1,2,3} (y − f (x) − ĉi )⊤ Σ̂−1 i (y − f (x) − ĉi ). SCP-3: As SCP-2, but reserving n′cal = 21 ncal points for threshold calibration. RemMCP: The MCP approach via constraint removal (Sec. 3.1), which jointly optimizes the covariance Σ of an ellipsoidal prediction set via the score s(x, y) = (y − f (x))⊤ Σ−1 (y − f (x)) over the full calibration set while minimizing ellipsoid volume (see Appendix C). RelMCP: The MCP approach via constraint relaxation (Sec. 3.2), targeting a confidence level 1 − β at least as large as that of SCP-1 (given by (15)), using A) the same convex single-ellipsoid score and cost as RemMCP for a direct comparison, or B) the same cost with the non-convex union-of-ellipsoids score s(x, y) = mini∈{1,2,3} (y −f (x)− ci )⊤ Σ−1 i (y − f (x) − ci ), where the centers ci and covariances Σi are jointly optimized during calibration. 8
SCP-1
SCP-2A
SCP-3A
RemMCP
RelMCP-A
5
SCP-2B
SCP-3B
RelMCP-B
5
e2
2
0
0
0
−2
−5 −4
−2
0
2
−5 −1
4
0
e1
−4
1
−2
e1
(a) T = 2 (residuals scaled by 103 ).
0
2
4
e1
(b) T = 50.
(c) T = 50.
Figure 2: Vehicle prediction sets and calibration residuals (∗) for run 1. SCP-2A
SCP-3A
RemMCP
RelMCP-A
SCP-2B
0.97
0.97
0.97
0.95
0.95
0.95
0.93
0.93
0.93
Volume
Coverage
SCP-1
50
6
40
4
30
2
(a) T = 2 (Volume was scaled by 106 ).
SCP-3B
RelMCP-B
60
30
(b) T = 30.
(c) T = 50.
Figure 3: Empirical coverage (top) and prediction set volume (bottom) over 10,000 test points and 100 calibration runs. The black diamond marks the mean; the SCP-2B and SCP-3B mean volumes at T = 2 (80,139 and 87,022) lie outside the plot range.
Fig. 2 visualizes prediction sets and calibration residuals for one representative run at T = 2, T = 30, and T = 50. At T = 2 (Fig. 2a), the model predictions remain close to the true positions, yielding compact, nearly unimodal residuals well approximated by a single ellipsoid. This regime allows us to assess whether joint optimization produces tighter prediction sets than decoupled alternatives even when the geometry is simple. At T = 30 and T = 50 (Figures 2b and 2c), the residuals cluster around three distinct modes, providing a more demanding test of whether joint calibration remains competitive when modeling complex multi-modal structure. Fig. 3 reports empirical coverage and prediction set volume across 100 independent calibration runs. All methods achieve the target coverage on average. RelMCP consistently exceeds it by a small margin, a consequence of the conservatism introduced by its iterative penalty search and a conservative heuristic, which is used to compute the solution complexity di . Across all three prediction horizons, the multi-modal methods SCP-2B, SCP-3B, and RelMCPB produce substantially smaller prediction sets than the spherical SCP-1 and the single-ellipsoid approaches SCP2A, SCP-3A, RemMCP, and RelMCP-A. Among the multimodal methods, RelMCP-B is the most robust: whereas the mean volumes of SCP-2B and SCP-3B at T = 2 are 80,139 and 87,022 respectively (inflated by a small number of degenerate runs), RelMCP-B achieves a stable mean volume below 40. Among single-ellipsoid methods, RemMCP yields the tightest sets. All SCP-2 and SCP-3 9
SCP-1 SCP-2 SCP-3 RemMCP
0.9
1−ε
0.95
1
Coverage Probability η(q)
Figure 4: Coverage distributions for ε = 0.05 obtained from (16). The displayed curve for RemMCP is a lower bound in general, and holds exactly under Assumptions 4 and 5.
variants exhibit a higher variance, especially in coverage, a direct consequence of the data split reducing the number of points available for threshold calibration. For T = 30, the standard deviation in coverage ranges from 0.0066 to 0.0101 for SCP-2 and SCP-3 variants, compared with 0.0058 for SCP-1, 0.0056 for RemMCP, 0.0046 for RelMCP-A, and 0.0039 for RelMCP-B. This behavior is accurately predicted by the theoretical coverage distributions derived via (16) and visualized in Fig. 4: with only n′cal = 41 ncal points for threshold estimation, the coverage distribution of SCP-2 has substantially heavier tails than those of SCP-1 or RemMCP, meaning SCP-2 is more likely to fall well short of (or well above) the target coverage on any given run. Allocating half the data to calibration (SCP-3) narrows the distribution but correspondingly reduces the data available for shape estimation. Achieving low variance without sacrificing data efficiency requires using all calibration data for shape optimization and threshold calibration simultaneously — precisely what MCP is designed to do. Regarding computation, SCP-1, SCP-2, and SCP-3 complete in under 10 ms. RemMCP requires approximately 0.2 s and RelMCP-A around 0.8 s. RelMCP-B incurs the highest cost at roughly 21 s, driven by repeated optimization over a high-dimensional parameter vector.
5
Conclusions
We have introduced multi-variable conformal prediction (MCP), a framework that extends conformal calibration to vector-valued score functions and multiple simultaneous optimization variables, eliminating the need for data splitting between shape optimization and threshold calibration. By grounding MCP in scenario theory, we derived two tractable variants — RemMCP, based on constrained optimization with constraint removal, and RelMCP, based on constraint relaxation — each equipped with finite-sample coverage guarantees under mild assumptions. RemMCP recovers split conformal prediction as a special case and is provably less conservative than dimension-wise Bonferroni-corrected conformal prediction. Numerical experiments demonstrate that both variants consistently meet the target coverage while producing prediction sets smaller than or comparable to split conformal baselines — with lower variance across calibration runs, a direct consequence of using all available calibration data simultaneously for shape optimization and threshold calibration. Limitations of MCP include the increased computational cost relative to split conformal prediction and the need for a sufficiently large calibration set, as MCP does not leverage training data for shape estimation. Future work includes developing richer score and cost function designs for structured outputs and classification tasks, tightening finite-sample bounds for RelMCP, and extending MCP to settings with distribution shift—directions that we expect will further broaden the practical impact of the framework.
Acknowledgments and Disclosure of Funding This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – SFB 1608 – 501798263.
References Anastasios N. Angelopoulos and Stephen Bates. A gentle introduction to conformal prediction and distribution-free uncertainty quantification. Preprint, 2022. https://arxiv.org/abs/2107.07511. Anastasios N. Angelopoulos, Rina Foygel Barber, and Stephen Bates. Theoretical Foundations of Conformal Prediction. Cambridge University Press, 2026. Preprint, https://arxiv.org/abs/2411.11824. Nathaniel Bowman and Michael T. Heath. Computing minimum-volume enclosing ellipsoids. Mathematical Programming Computation, 15(4):621–650, 2023. Sacha Braun, Liviu Aolaritei, Michael I. Jordan, and Francis Bach. Minimum volume conformal sets for multivariate regression. Preprint, 2025. https://arxiv.org/abs/2503.19068. Giuseppe C. Calafiore and Marco C. Campi. Uncertain convex programs: randomized solutions and confidence levels. Mathematical Programming, 102(1):25–46, 2005. Giuseppe C. Calafiore and Marco C. Campi. The scenario approach to robust control design. IEEE Transactions on Automatic Control, 51(5):742–753, 2006. 10
Marco C. Campi and Simone Garatti. The exact feasibility of randomized solutions of uncertain convex programs. SIAM Journal on Optimization, 19(3):1211–1230, 2008. Marco C. Campi and Simone Garatti. A sampling-and-discarding approach to chance-constrained optimization: Feasibility and optimality. Journal of Optimization Theory and Applications, 148: 257–280, 2011. Marco C. Campi, Giuseppe C. Calafiore, and Simone Garatti. Interval predictor models: Identification and reliability. Automatica, 45(2):382–392, 2009. Marco C Campi, Simone Garatti, and Federico A Ramponi. Non-convex scenario optimization with application to system identification. In IEEE Conference on Decision and Control, pages 4023–4028, 2015. Marco C. Campi, Algo Carè, and Simone Garatti. The scenario approach: a tool at the service of data-driven decision making. Annual Reviews in Control, 52:1–17, 2021. Andrea Cini, Alexander Jenkins, Danilo Mandic, Cesare Alippi, and Filippo Maria Bianchi. Relational conformal prediction for correlated time series. In International Conference on Machine Learning, 2025. Matthew Cleaveland, Insup Lee, George J. Pappas, and Lars Lindemann. Conformal prediction regions for time series using linear complementarity programming. In Conference on Artificial Intelligence, 2024. Rudi Coppola, Andrea Peruffo, Lars Lindemann, and Manuel Mazo. Scenario approach and conformal prediction for verification of unknown systems via data-driven abstractions. In European Control Conference, pages 558–563, 2024. Marco de Angelis, Roberto Rocchetta, Ander Gray, and Scott Ferson. Constructing consonant predictive beliefs from data with scenario theory. In International Symposium on Imprecise Probabilities and Their Applications, 2021. Elizabeth Dietrich, Alex Devonport, and Murat Arcak. Nonconvex scenario optimization for datadriven reachability. In Conference on Learning for Dynamics and Control, pages 514–527, 2024. Zhenhan Fang, Aixin Tan, and Jian Huang. CONTRA: Conformal prediction region via normalizing flow transformation. In International Conference on Learning Representations, 2025. Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated multiple-output quantile regression with representation learning. Journal of Machine Learning Research, 24(24):1–48, 2023. Simone Garatti and Marco C. Campi. Risk and complexity in scenario optimization. Mathematical Programming, 191(1):243–279, 2022. Simone Garatti and Marco C. Campi. Non-convex scenario optimization. Mathematical Programming, 209(1):557–608, 2025. Simone Garatti, Marco C. Campi, and Algo Carè. On a class of interval predictor models with universal reliability. Automatica, 110(C), 2019. article no. 108542. Ander Gray, Vignesh Gopakumar, Sylvain Rousseau, and Sebastien Destercke. Guaranteed prediction sets for functional surrogate models. In Conference on Uncertainty in Artificial Intelligence, 2025. Chancellor Johnstone and Bruce Cox. Conformal uncertainty sets for robust optimization. In Proceedings of the 10th Symposium on Conformal and Probabilistic Prediction and Applications, volume 152, pages 72–90, 2021. Jing Lei, James Robins, and Larry Wasserman and. Distribution-free prediction sets. Journal of the American Statistical Association, 108(501):278–287, 2013. Albert Lin and Somil Bansal. Verification of neural reachable tubes via scenario optimization and conformal prediction. In Alessandro Abate, Mark Cannon, Kostas Margellos, and Antonis Papachristodoulou, editors, Conference on Learning for Dynamics and Control, volume 242, pages 719–731, 2024. 11
Stefan B. Liu, Bastian Schürmann, and Matthias Althoff. Guarantees for real robotic systems: Unifying formal controller synthesis and reachset-conformant identification. IEEE Transactions on Robotics, 39(5):3776–3790, 2023. Rui Luo and Zhixin Zhou. Volume-sorted prediction set: Efficient conformal prediction for multitarget regression. Preprint, 2025. https://arxiv.org/abs/2503.02205. Laura Lützow and Matthias Althoff. Scalable reachset-conformant identification of linear systems. IEEE Control Systems Letters, 8:520–525, 2024. Laura Lützow and Matthias Althoff. Reachset-conformant system identification. IEEE Transactions on Automatic Control, 71(5):3014–3029, 2026. Laura Lützow, Michael Eichelbeck, Mykel Kochenderfer, and Matthias Althoff. Zono-conformal prediction: Zonotope-based uncertainty quantification for regression and classification tasks. Journal of Machine Learning Research, 26(294):1–37, 2025. Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Copula-based conformal prediction for multi-target regression. Pattern Recognition, 120, 2021. article no. 108101. Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Ellipsoidal conformal inference for multi-target regression. In Symposium on Conformal and Probabilistic Prediction with Applications, volume 179, pages 294–306, 2022. Niall O’Sullivan, Licio Romao, and Kostas Margellos. Bridging conformal prediction and scenario optimization. In IEEE Conference on Decision and Control, pages 6114–6121, 2025. Licio Romao, Kostas Margellos, and Antonis Papachristodoulou. Probabilistic feasibility guarantees for convex scenario programs with an arbitrary number of discarded constraints. Automatica, 149, 2023a. article no. 110601. Licio Romao, Antonis Papachristodoulou, and Kostas Margellos. On the exact feasibility of convex scenario programs with discarded constraints. IEEE Transactions on Automatic Control, 68(4): 1986–2001, 2023b. Jonathan Sadeghi, Marco de Angelis, and Edoardo Patelli. Efficient training of interval neural networks for imprecise training data. Neural Networks, 118:338–351, 2019. Eleftherios Spyromitros-Xioufis, Grigorios Tsoumakas, William Groves, and Ioannis Vlahavas. Multi-target regression via input space expansion: Treating targets as inputs. Machine Learning, 104:55–98, 2016. Kamilė Stankevičiūtė, Ahmed Alaa, and Mihaela van der Schaar. Conformal time-series forecasting. In International Conference on Neural Information Processing Systems, 2021. article no. 475. Sophia Huiwen Sun and Rose Yu. Copula conformal prediction for multi-step time series prediction. In International Conference on Learning Representations, 2024. Renukanandan Tumu, Matthew Cleaveland, Rahul Mangharam, George J. Pappas, and Lars Lindemann. Multi-modal conformal prediction regions by optimizing convex shape templates. In Conference on Learning for Dynamics and Control, 2024. Vladimir Vovk. Conditional validity of inductive conformal predictors. Machine Learning, 92(2–3): 349–376, 2013. Volodya Vovk, Alexander Gammerman, and Craig Saunders. Machine-learning applications of algorithmic randomness. In International Conference on Machine Learning, pages 444–453, 1999. Chen Xu, Hanyang Jiang, and Yao Xie. Conformal prediction for multi-dimensional time series by ellipsoidal sets. In International Conference on Machine Learning, 2024. article no. 2268. Botong Zhang, Shuo Li, and Osbert Bastani. Conformal structured prediction. In International Conference on Learning Representations, pages 8563–8583, 2025.
12
A
Proofs
This section provides the proofs of the theorems and the corollary of the main paper. A.1
Proof of Thm. 1
The programs in (14) are scenario programs of the form (1) with δ = (x, y), convex sets Qδ = {q : s(x, y, q) ≤ 0}, and convex objective J(q). Thus, Assumption 3 holds. If additionally Assumptions 1 and 2 hold and the calibration data points are exchangeable (as required by Thm. 1), (4) applies. Since the coverage satisfies η(q) = 1 − V (q) for any q, where V (q) denotes the violation probability of the corresponding scenario program, this yields E[η(q ∗ )] = 1 − E[V (q ∗ )] ≥ 1 −
nq (ρ + 1) ≥ 1 − ε, ncal + 1
where the last inequality follows from substituting (13), establishing E[η(q ∗ )] ≥ 1 − ε. If the samples are i.i.d., the high-confidence bound in (12) with β given by (15) follows directly from (3). Under the additional Assumptions 4 and 5, (3) holds with equality [Romao et al., 2023b, Thm. 5], and hence so does the high-confidence bound in (12) for RemMCP, which yields (16). A.2
Proof of Corollary 1
We first show that SCP and RemMCP yield identical prediction sets, and then derive the coverage guarantees of SCP from those of RemMCP. Prediction Sets. For an arbitrary SCP score function s̃(x, y), the SCP prediction set defined in (6) coincides with the MCP prediction set in (10) with score function (17a) if q ∗ = q̃ ∗ . Hence, it suffices to show that the parameter q ∗ obtained from the RemMCP procedure with (17) equals the SCP threshold q̃ ∗ . As in SCP, the calibration procedure of RemMCP begins by computing the outlier budget for a given target miscoverage level ε. For nq = 1, the outlier budget in (13) reduces to ρ = ⌊ε(ncal + 1)⌋ − 1, which coincides with the SCP budget in (5). The optimal parameter q ∗ in RemMCP is obtained by solving the cascading optimization problem in (14). For the score and cost functions in (17), which are convex in q and therefore admissible for RemMCP, this reduces to qk∗ = argmin q q∈Rnq
s.t. s̃(x(m) , y (m) ) ≤ q, ∀(x, y) ∈ Mcal \ Rk .
(21)
At each stage k = 0, . . . , ρ, the optimizer qk∗ of (21) equals the maximum nonconformity score s̃(x(m) , y (m) ) over the samples in Mcal \ Rk . By (2) with nq = 1, the set of discarded samples Rk at stage k is constructed by adding the current support scenario to the previously discarded samples in Rk−1 . The support sample at stage k is the data point whose removal would alter the optimizer qk∗ , and under a continuous score distribution, the maximum nonconformity score is attained uniquely with probability one; hence, exactly one data point, the one corresponding to the maximum nonconformity score, serves as the support sample at each stage. Consequently, after ρ stages, the set Rρ contains the ρ samples with the largest nonconformity scores, and q ∗ = qρ∗ equals the (ncal − ρ)-smallest nonconformity score. This coincides exactly with the SCP threshold q̃ ∗ . Therefore, the prediction sets of SCP and RemMCP are identical. Coverage guarantees. Next, we derive the coverage guarantees of SCP in (7) and (8) from those of RemMCP. Since the score is strictly monotone in q, the optimizer of (14) over any finite sample set exists and is unique (Assumption 1). Moreover, when the score distribution is continuous, the program is nondegenerate with probability one (Assumption 2). Assuming exchangeability of the data points, the conditions of Thm. 1 are satisfied, yielding the expected coverage guarantee in (11) for RemMCP, which matches (7). When nq = 1, the distributional result for the coverage of RemMCP in (16) reduces to η(q ∗ ) ∼ Beta(ncal − ρ, ρ + 1), which coincides with (8). Both results require i.i.d. calibration data. The additional assumptions for (16) (Assumptions 4 and 5) hold almost surely in this setting, since the maximum is attained uniquely and each removed score strictly exceeds the remaining threshold with probability one under a continuous score distribution. The continuity condition is identical to that required for (8). 13
A.3
Proof of Corollary 2
With the score function in (18a), the MCP prediction set in (10) is an axis-aligned hyperrectangle centered at f (x) with half-width vector q ∗ . Similarly, dimension-wise SCP with the score function in (9) produces hyperrectangular prediction sets centered at f (x), where the half-width q̃j∗ in each dimension j is calibrated independently using standard SCP. The MCP prediction set is contained in ∗ (or equal to) the dimension-wise SCP set if and only if q(j) ≤ q̃j∗ for all j = 1, . . . , ny . We show that this relation holds while both methods target the same expected coverage level 1 − ε. To guarantee joint coverage of 1 − ε, dimension-wise SCP calibrates each marginal half-width q̃j∗ at level ε/ny via Bonferroni correction. By (5), this yields an outlier budget of ρ = ⌊(ε/ny )(ncal + 1)⌋ − 1 per dimension, which coincides with the RemMCP budget ρ in (13) for target miscoverage ε and nq = ny . Dimension-wise SCP discards, for each dimension j, the ρ samples with the largest deviations |y(j) − f(j) (x)|. Let R̃ρ,j denote the set of samples discarded in dimension j, and define S ny R̃ρ = j=1 R̃ρ,j . Since dimensions are treated independently, the same sample may be discarded in multiple dimensions, implying |R̃ρ | ≤ ρny . In contrast, the RemMCP discard set Rρ contains exactly ρny samples. We now show that Rρ ⊇ R̃ρ . With the score and cost functions in (18), the RemMCP optimization problem becomes qk∗ = argmin 1⊤ q q∈Rnq
s.t. |y − f (x)| ≤ q ∀(x, y) ∈ Mcal \ Rk .
∗ Thus, for each dimension j, the component qk(j) equals the maximum deviation |y(j) − f(j) (x)| over the samples in Mcal \ Rk . Under a continuous score distribution, this maximum is attained uniquely in each dimension with probability one. The corresponding samples form the support set supp(qk∗ ) and are added to the discard set Rk+1 at the next stage (cf. (2)). After ρ stages, Rρ therefore contains, for each dimension j, the ρ samples with the largest values of |y(j) − f(j) (x)|. Hence, Rρ ⊇ R̃ρ .
To illustrate this, consider a sample m̃ that attains the largest deviation in dimension j and the second-largest deviation in dimension j ′ ̸= j. If ρ ≥ 2, dimension-wise SCP discards this sample in both dimensions, resulting in |R̃ρ | < ρny . In contrast, RemMCP removes this sample already at stage k = 1 (as the support constraint for dimension j), and can subsequently discard an additional sample corresponding to the next-largest deviation in dimension j ′ . This would lead to a strictly larger set of discarded samples. Since Rρ ⊇ R̃ρ , the overall solution q̃ ∗ from dimension-wise SCP is feasible for (14). RemMCP can ∗ therefore exploit the remaining outlier budget to further reduce the half-widths, yielding q(j) ≤ q̃j∗ for all j = 1, . . . , ny . Consequently, Y(x, q ∗ ) ⊆ Ỹ(x, q̃1 ) × · · · × Ỹ(x, q̃ny ) for every x, i.e., the prediction sets produced by RemMCP are contained in (or equal to) those from dimension-wise SCP, while both achieve the target coverage level 1 − ε. A.4
Proof of Thm. 2
By introducing additional optimization variables ξ¯m ≥ 0, which compute the slack ξm (q), problem (19) is equivalent to qi∗ =
argmin n
J(q) + ϕi
cal q∈Rnq ,ξ̄∈R≥0
ncal X
ξ¯m
(22a)
m=1
s.t. s(x(m) , y (m) , q) ≤ ξ¯m 1, ∀m = 1, . . . , ncal ,
(22b)
with ξ¯ = [ξ¯1 . . . ξ¯ncal ]⊤ . This is a non-convex scenario program with constraint relaxation, as introduced in Garatti and Campi [2025, Sec. 4.1]. Applying Garatti and Campi [2025, Thm. 10] to the i-th program yields P (η(qi∗ ) ≥ 1 − εi ) ≥ 1 −
β neval
,
where εi is the unique solution of (20). We evaluate neval distinct penalty parameters, and since any evaluation may be incorrect with probability β/neval , all evaluations and, thus, also the selection are 14
simultaneously guaranteed with confidence 1 − β; see also Garatti and Campi [2025, Footnote 29]. Thus, the solution must satisfy P (η(qi∗ ) ≥ 1 − εi ) ≥ 1 − β. For εi ≤ ε, we obtain P (η(qi∗ ) ≥ 1 − ε) ≥ P (η(qi∗ ) ≥ 1 − εi ) ≥ 1 − β. This establishes (12).
B
RelMCP Calibration Algorithm
Alg. 1 describes an iterative procedure for finding a feasible penalty parameter ϕi for the unconstrained optimization problem in (19) of RelMCP. The penalty parameter governs a fundamental trade-off: a large ϕi penalizes constraint violations heavily, producing solutions that satisfy more calibration constraints (lower ξm ), which in turn decreases the solution complexity di and tightens the achievable coverage guarantee; a small ϕi permits more violations, yielding a higher complexity di and looser coverage guarantees, but also less conservative solutions. The goal of the algorithm is to find the value of ϕi that yields the best (least conservative) prediction sets for the target miscoverage level ε and the confidence level 1 − β. We provide a detailed step-by-step explanation of the proposed algorithm below. Initialization. The algorithm maintains two iteration counters: ival , indexing the best valid solution found so far — the one with the smallest penalty ϕi such that the miscoverage εi ≤ ε — and iinv , indexing the best invalid solution — the one with the largest penalty ϕi for which εi > ε still holds. Both counters are initialized to zero in Line 1, indicating that no solution has yet been evaluated. Main loop. At each iteration i, the algorithm solves the penalized program (19) with the current penalty ϕi in Line 3, obtaining a calibrated parameter qi∗ and solution complexity di . The miscoverage εi is then computed in Line 4 from (20) using di and neval = i, where the dependence on i reflects a union bound over all iterations evaluated so far (cf. Thm. 2). Three outcomes are possible: 1. Valid solution (εi ≤ ε): The current solution meets the target coverage guarantee and is recorded as the new best valid solution in Line 6. By the penalty update rule in Line 15, all subsequent penalties satisfy ϕj < ϕi for j > i, so future iterations can only yield a less conservative (i.e., better) valid solution or become invalid. 2. All constraints satisfied, coverage insufficient (εi > ε and ξm = 0 for all m): All calibration constraints are already satisfied with zero slack, so no further increase in ϕ can reduce the solution complexity di or improve the coverage guarantee. The algorithm returns None in Line 8, indicating that no valid solution was found for the given calibration set size and target level ε. 3. Invalid solution (εi > ε with some ξm > 0): The coverage guarantee is not yet met, but active slack variables indicate that increasing ϕ may reduce violations and lower di . The iteration is recorded as the new best invalid solution in Line 10. By the update rule in Line 15, all subsequent penalties satisfy ϕj > ϕi for j > i, so future iterations can only yield a more conservative (i.e., better) invalid solution or become valid. Before continuing, the algorithm checks whether the current best valid solution would remain valid in the next iteration by computing the coverage ε′ival via (20) with neval = i + 1 in Line 12. If ε′ival > ε, returning immediately is necessary, since the confidence budget is too tight to validate this solution in any future iteration. If the maximum number of iterations imax is reached, the algorithm likewise returns. Otherwise, it continues searching for a smaller ϕ, retaining ival as a fallback if the next iteration is invalid. After each iteration, ϕ is updated in Line 15 according to a bisection-like scheme, guided by the history of valid and invalid solutions: • No valid solution yet (ival = 0): The penalty is increased by multiplying the last invalid max(1, i−imax ) penalty ϕiinv by f+ > 1, with the exponent growing beyond one once i exceeds imax to accelerate the search. • No invalid solution yet (iinv = 0): The penalty is decreased by multiplying the last valid penalty ϕival by f− ∈ (0, 1), probing whether a smaller ϕ can still certify the target coverage. 15
• Both valid and invalid solutions exist: The penalty is set to the midpoint 12 (ϕival + ϕiinv ), performing binary search between the best known valid and invalid penalties. Termination. The algorithm returns the best valid solution (qi∗val , ϕival , dival , εival , ival ) as soon as one of two conditions holds: (i) the best valid solution would no longer be certifiable in the next iteration due to the shrinking per-iteration confidence budget, or (ii) the maximum number of iterations imax is reached and a valid solution is found. If no valid solution is found but all constraints are already satisfied, the algorithm returns None. Practical considerations. The initial penalty ϕ1 can strongly influence convergence: a value that is too small may cause many iterations to be spent in the invalid regime, while a value that is too large may over-penalize violations and yield unnecessarily conservative prediction sets. As shown in the experiments in Appendix E, a poor initialization can prevent the algorithm from finding any valid solution within imax iterations, particularly at small miscoverage levels. In practice, ϕ1 can be tuned heuristically by warm-starting from a related problem instance. The increase and decrease factors f+ and f− control the speed of the search: larger f+ and smaller f− lead to faster but coarser exploration, while values closer to 1 yield finer but slower convergence. Coverage guarantee. Thm. 2 provides a finite-sample, high-confidence coverage guarantee for RelMCP under the assumption that a valid solution q ∗ — one satisfying εi ≤ ε — is selected from a predetermined grid of penalty parameters ϕ1 , . . . , ϕneval . Alg. 1 instead constructs this grid adaptively: each penalty ϕi depends on the outcomes of all preceding solver calls, which themselves depend on the calibration data. This data-dependence breaks the independence assumption underlying Thm. 2, so the theorem’s guarantee does not apply strictly to the algorithm’s output. Empirically, Alg. 1 consistently achieves the target coverage (see Sec. 4 and Appendix E), suggesting that the adaptive search does not meaningfully inflate the effective miscoverage. A rigorous treatment of the data-adaptive case is left for future work.
C
Example Instantiations for the Score and Cost Functions
We present representative score and cost function instantiations for regression tasks to illustrate the flexibility of MCP. Although the examples focus on regression, the framework extends straightforwardly to other settings such as classification or structured prediction [Zhang et al., 2025]. Throughout, we assume that a point predictor f (x) has been pretrained on a separate training dataset. The input⊤ ⊤ ⊤ ⊤ ⊤ output pair (x, y) may represent time series, x = [x⊤ 1 · · · xnk ] , y = [y1 · · · ynk ] , with prediction sets Y = Y1 × · · · × Ynk over a time horizon nk . C.1
Convex Prediction Set Shapes
For convex prediction set families, the score function follows directly from the set membership condition. Intervals. A natural and widely applicable choice is the multi-dimensional interval Y(x, q) = {y | h(x, q1 ) ≤ y − f (x) ≤ h(x, q2 )} with lower and upper bounds h(x, q1 ), h(x, q2 ) ∈ Rny . The corresponding score and cost functions are −y + f (x) + h(x, q1 ) s(x, y, q) = , (23a) y − f (x) − h(x, q2 ) 1 X J(q) = ∥h(x, q2 ) − h(x, q1 )∥, (23b) nX x∈X
where X is a set of representative inputs (for example, derived from the training set Mtrain ) with cardinality |X | = nX . A sample (x, y) is conformal if and only if y lies within the interval; the cost penalizes total interval width across the training inputs. Ellipsoids. MCP naturally accommodates ellipsoidal prediction sets of the form Y(x, q) = {y | (y − f (x) − h(x, q))⊤ Σ(x, q)−1 (y − f (x) − h(x, q)) ≤ 1}, 16
where Σ(x, q) ≻ 0 is the covariance matrix of the ellipsoid, determining its shape, and f (x) + h(x, q) is the ellipsoid center. The score and cost functions are s(x, y, q) = (y − f (x) − h(x, q))⊤ Σ(x, q)−1 (y − f (x) − h(x, q)) − 1, 1 X J(q) = − log det Σ(x, q)−1 , nX
(24a) (24b)
x∈X
where the cost minimizes ellipsoid volume [Bowman and Heath, 2023]. When the ellipsoid shape is constant across inputs, Σ can be parameterized by q ∈ Rnq with nq = 12 ny (ny + 1) free parameters. If h does not depend on q, the score function is convex in q, making RemMCP directly applicable. Zonotopes. Zonotopes are centrally symmetric polytopes well-suited to high-dimensional uncertainty quantification due to their compact parameterization and flexible geometry. A zonotopic prediction set takes the form Y(x, q) = {f (x) + c(x, q) + G(x, q)λ | λ ∈ [−1, 1]}, where f (x) + c(x, q) ∈ Rny is the zonotope center and G(x, q) ∈ Rny ×ngen is the generator matrix. The corresponding score and cost functions are s(x, y, q) = 1 − I{∃λ ∈ [−1, 1] : y = f (x) + c(x, q) + G(x, q)λ}, nrot 1 XX 1⊤ |Ri G(x, q)|1, J(q) = nX i=1
(25a) (25b)
x∈X
where Ri ∈ Rny ×ny are random rotation matrices and | · | denotes element-wise absolute value. The cost is a tractable proxy for zonotope volume [Lützow and Althoff, 2024, Lützow et al., 2025]. When both c(x, q) and G(x, q) are linear in q, the calibration program (14) reduces to a linear program via standard auxiliary variable reformulations. C.2
Sublevel Sets
When the convex shapes above are insufficiently expressive, prediction sets can be defined as sublevel sets of arbitrary functions, enabling nonconvex geometries. Following Dietrich et al. [2024], we construct the score function as a mixture of N radial basis functions with centers µi ∈ Rny , widths σi > 0, and threshold γ: s(x, y, q) =
N (y−f (x)−µi )2 X −1 σ2 i e 2 − γ,
(26a)
i=1
J(q) =
N X
σi2 ,
(26b)
i=1 ⊤ ⊤ with parameter vector q = [µ⊤ 1 · · · µN σ1 · · · σN γ] . The cost penalizes the total spread of the basis functions, promoting compact prediction sets.
C.3
Unions of Prediction Sets
Multi-modal prediction sets can be obtained by combining component score functions via the minoperator, following the approach of Tumu et al. [2024], but without requiring an additional data split. Given partial score functions sqi (x, y) and volume proxies Jqi , i = 1, . . . , N —each defined by any ⊤ of the instantiations above—the composite score and cost with q = [q1⊤ . . . qN ] are s(x, y) = min (sqi (x, y)) , i
J(q) =
N X
Jqi .
(27a) (27b)
i=1
The resulting prediction set is the union of the N component sets, with all shape parameters optimized and calibrated jointly over the full calibration dataset. 17
D
Experimental Details
This section provides additional details on the setups of the experiments from Sec. 4. Additional results are provided in Appendix E. All computations are carried out in Python on an i9-12900HK processor (2.5GHz) with 64GB of memory. SCP-2. The calibration set is split into a fitting subset of size ncal − n′cal and a calibration subset of size n′cal = 14 ncal . For SCP-2A, a single ellipsoid is fitted by computing the sample covariance Σ̂ of the fitting residuals. For SCP-2B, k-means clustering with K = 3 is applied to the fitting residuals, and the sample covariance of the fitting residuals of each cluster is computed. Each covariance matrix is normalized so that the (1 − ε)-quantile of its per-cluster Mahalanobis distance equals 1 as done in Tumu et al. [2024]. We note that computing the minimum-volume enclosing ellipsoid per cluster via CMA-ES, as proposed in Tumu et al. [2024], resulted in higher variance and larger prediction sets in our experiments than the sample covariance approach, and was therefore not adopted. In both cases, the fitted shape is then held fixed and a single non-negative scaling threshold q̃ ∗ is calibrated on the reserved split. SCP-3. As SCP-2, but using n′cal = 12 ncal data points for calibration. RemMCP. The cascading scenario program is solved using SLSQP at each of the ρ + 1 stages, with a maximum of 2,000 iterations per stage. The inverse covariance matrix is parameterized as Σ−1 = LL⊤ , with L being a lower triangular matrix and initialized to the identity. The solution of each stage is used as the warm start for the next. RelMCP. Each penalized subproblem is solved with L-BFGS-B, using 3 random restarts and a maximum of 2,000 inner optimizer iterations (RelMCP-A) or 10 restarts and a maximum of 1,000 inner optimizer iterations (RelMCP-B). We set imax = 15 for the outer penalty loop with initial penalty ϕ1 = 0.1, increase factor f+ = 2.0, and decrease factor f− = 0.3. RelMCP-A is initialized with the identity covariance matrix; RelMCP-B is warm-started from a k-means fit over the calibration set with per-cluster sample covariance estimates. In both cases, half of the restarts are initialized close to the current parameter guess and the other half near the best solution from the previous penalty iteration. The solution complexity di is computed via a heuristic: we set di to the number of violated constraints plus min(nq , nnv ), where nnv is the number of non-violated data points. While this is not guaranteed to equal the true solution complexity, we observe that the target coverage is consistently exceeded in practice; for high-dimensional parameter vectors, the true solution complexity is expected to be substantially smaller than nq , making the heuristic conservative. Volume and coverage estimation. Prediction set area is estimated via Monte Carlo sampling using 105 uniform draws inside an adaptive bounding box, applied uniformly across all methods. Empirical coverage is evaluated on ntest = 10,000 held-out test points, and results are aggregated over 100 independent calibration runs, each with a freshly sampled calibration set of size ncal = 2,000.
E
Additional Experiments
We present additional experimental results complementing those of Sec. 4, where we evaluate all methods across varying target miscoverage levels and output dimensions over 100 independent runs with randomly drawn calibration and test data. E.1
Varying Target Coverages
We evaluate all methods across a range of target miscoverage levels ε ∈ {0.01, 0.02, 0.05, 0.15} on two tasks: • Vehicle prediction: We use the trajectory predictor from Tumu et al. [2024] to predict the vehicle position at T = 30 time steps (ncal = 2,000, ntest = 10,000), as in Sec. 4. • CASP dataset: We predict the 2-dimensional output of the CASP regression benchmark [Feldman et al., 2023], training a two-layer network with 64 neurons per layer on 5,000 data points. From the remaining 40,730 points we draw ncal = 2,000 for calibration and 10,000 for testing; all features and targets are normalized to [0, 1]. Because the CASP 18
residuals are unimodal, we restrict comparison to the unimodal methods (SCP-1, SCP-2A, SCP-3A, RemMCP, RelMCP-A) for this task. Fig. 5 visualizes prediction sets and calibration residuals for vehicle prediction across ε ∈ {0.01, 0.02, 0.15} (results for ε = 0.05 appear in Sec. 4), and Fig. 6 shows the corresponding results for the CASP dataset. Box plots of empirical coverage and volume are reported in Figures 7 and 8. The relative ordering of methods is mostly consistent across all coverage levels and both tasks, confirming the findings of Sec. 4. For the vehicle benchmark, the multi-modal methods (SCP-2B, SCP-3B, RelMCP-B) consistently produce substantially smaller prediction sets than their unimodal counterparts, reflecting the multi-modal residual structure. Among unimodal methods, RemMCP achieves the smallest prediction sets across both tasks, while the relative performance of RelMCP-A and the SCP variants depends on the task: RelMCP-A outperforms SCP-3A and SCP-2A on the vehicle benchmark for ε ≥ 0.02 but is more conservative on the CASP dataset. As ε increases, the performance of RemMCP and RelMCP improves, since larger outlier budgets allow more aggressive shape optimization. All SCP-2 and SCP-3 variants exhibit higher coverage variance across runs, a direct consequence of the data split reducing the effective calibration set size. RelMCP consistently exceeds the target coverage by a small margin due to the conservatism of its iterative penalty search. At ε = 0.01, RelMCP-B fails to return a certified valid solution for vehicle prediction and is therefore omitted from Fig. 7, since the per-iteration confidence budget β/neval becomes too tight at very low miscoverage levels and a larger calibration set would be required to certify coverage via (20). Computation times and guaranteed confidence levels 1 − β are reported in Tables 1 and 2. RemMCP and RelMCP computation times mostly grow with ε as the outlier budget ρ increases, requiring more stages of the cascading program or more penalty iterations, respectively. Table 1: Computation time [s] for vehicle position prediction at T = 30 across target miscoverage levels (mean ± standard deviation over 100 runs). Method ε = 0.01 ε = 0.02 ε = 0.05 ε = 0.15 SCP-1 SCP-3A SCP-4A RemMCP RelMCP-A SCP-3B SCP-4B RelMCP-B
0.0002 ± 0.0000 0.0008 ± 0.0001 0.0007 ± 0.0002 0.0513 ± 0.0037 0.5728 ± 0.3719 0.0059 ± 0.0038 0.0049 ± 0.0007 21.4239 ± 4.0192
0.0002 ± 0.0000 0.0008 ± 0.0001 0.0007 ± 0.0002 0.1010 ± 0.0059 0.7274 ± 0.2182 0.0056 ± 0.0007 0.0048 ± 0.0007 22.6351 ± 11.9164
0.0002 ± 0.0000 0.0008 ± 0.0003 0.0007 ± 0.0002 0.2281 ± 0.0123 0.8490 ± 0.1962 0.0062 ± 0.0042 0.0049 ± 0.0006 20.9342 ± 10.7426
0.0002 ± 0.0000 0.0008 ± 0.0001 0.0007 ± 0.0002 0.6073 ± 0.0153 0.8110 ± 0.1128 0.0059 ± 0.0013 0.0048 ± 0.0005 30.3584 ± 5.7495
Table 2: Guaranteed confidence 1−β of achieving target coverage 1−ε for vehicle position prediction at T = 30, obtained from (15) and (20). RelMCP values are reported for an example run. A dash (–) indicates that no certified valid solution was found. Method ε = 0.01 ε = 0.02 ε = 0.05 ε = 0.15 SCP-1 SCP-2A SCP-3A RemMCP RelMCP-A SCP-2B SCP-3B RelMCP-B
0.5302 0.5604 0.5427 0.7041 0.7159 0.5604 0.5427 –
0.5217 0.5433 0.5306 0.5852 0.7081 0.5433 0.5306 0.5234
19
0.5143 0.5286 0.5286 0.5552 0.6490 0.5286 0.5286 0.6065
0.5096 0.5191 0.5135 0.5096 0.5433 0.5191 0.5135 0.6236
SCP-1
SCP-2A
SCP-3A
RelMCP-A
SCP-2B
0
−2
−2
−1
0
1
RelMCP-B
−1
1
0
0
−2
SCP-3B
2
2
2
e2
RemMCP
−2 −2
2
−1
0
e1
1
−2
2
e1
(a) ε = 0.01.
0
2
e1
(b) ε = 0.02.
(c) ε = 0.15.
Figure 5: Vehicle position prediction sets and calibration residuals (∗) at T = 30 across varying ε.
SCP-1
SCP-2A
SCP-3A
RemMCP
RelMCP-A ·10−2
·10−2 0.2 5
e2
0.1
0
5
0
0
−0.1
−5
0
−5
−0.2 −0.5
0
−0.4 −0.2
0.5
0
0.2
0.4
e1
e1
0
0.2
−0.2 −0.1
e1
(b) ε = 0.02.
(a) ε = 0.01.
−0.2
0
0.1
0.2
e1
(c) ε = 0.05.
(d) ε = 0.15.
Figure 6: CASP dataset prediction sets and calibration residuals (∗) across varying ε.
SCP-1
SCP-2A
SCP-3A
RemMCP
RelMCP-A
SCP-2B
SCP-3B
RelMCP-B
Coverage
0.90 0.99 0.99
0.98
0.85
0.97 0.98
Volume
10
8 4
5
4 2
(a) ϵ = 0.01
(b) ϵ = 0.02
(c) ϵ = 0.15
Figure 7: Empirical coverage (top) and prediction set volume (bottom) for vehicle position prediction at T = 30 across varying ε over 10,000 test points and 100 calibration runs. The black diamond marks the mean.
20
SCP-1
SCP-2A
SCP-3A
RemMCP
0.89
0.97
Coverage
RelMCP-A
0.99 0.99
Volume
0.98
0.85
0.95
0.98 0.97
0.93
0.1
0.03
0.81
0.01
0.4 0.02 0.2
0.05 0.01 0
(a) ϵ = 0.01
(b) ϵ = 0.02
(c) ϵ = 0.05
(d) ϵ = 0.15
Figure 8: Empirical coverage (top) and prediction set volume (bottom) for the CASP dataset across varying ε over 10,000 test points and 100 calibration runs. The black diamond marks the mean. E.2
Varying Output Dimensions
We evaluate the effect of increasing output dimension on the SCM20D supply-chain dataset [Spyromitros-Xioufis et al., 2016], which contains 8,966 samples with up to 16 output dimensions. We train two-layer networks with 64 neurons per layer on 3,000 data points and consider output dimensions ny ∈ {2, 4, 6, 8}; the remaining data serve as the test set. Because the residuals are unimodal, we restrict comparison to the unimodal methods (SCP-1, SCP-2A, SCP-3A, RemMCP, RelMCP-A), all calibrated at ε = 0.05 with ncal = 1,000 calibration points. We do not consider ny ≥ 10: the ellipsoid covariance parameterization yields nq = 12 ny (ny + 1) = 55 free parameters at ny = 10, for which (13) produces a negative outlier budget at ε = 0.05 and ncal = 1,000, so no RemMCP solution can be certified without increasing ncal or reducing nq . Prediction sets over the first two output dimensions and box plots of empirical coverage and volume are shown in Figures 9 and 10; computation times are in Tab. 3. As ny increases, RemMCP becomes progressively more conservative, substantially overshooting the target coverage at ny = 8. This has a clear theoretical explanation: the ellipsoid covariance is parameterized by nq = 12 ny (ny + 1) free parameters (e.g., nq = 36 at ny = 8), but the number of support constraints at the optimum is likely far smaller, so that when the full-support assumption (Assumption 4) fails to hold in practice, the coverage bound in Thm. 1 becomes loose. This suggests that lower-dimensional parameterizations of the score function would be better suited to this setting.
21
SCP-1
SCP-2A
SCP-3A
RemMCP
RelMCP-A
1
1
0.5
0.5
0
0
−0.5
−0.5
0.5
e2
0.5
0
0
−0.5 −0.5
−1 −0.5
0
−0.5
0.5
0
0.5
e1
e1
(a) ny = 2.
(b) ny = 4.
−1 −1 −0.5
0
0.5
−1 −0.5
1
0
e1
e1
(c) ny = 6.
(d) ny = 8.
0.5
1
Figure 9: SCM20D prediction sets and calibration residuals (∗) for varying output dimension ny (visualized over the first two dimensions), calibrated at ε = 0.05.
Coverage
SCP-1
SCP-2A
SCP-3A
RemMCP
RelMCP-A
0.98
0.98
0.98
0.98
0.95
0.95
0.95
0.95
0.92
0.92
0.92
0.92
0.6
1
0.3
0.5
Volume
2
(a) ny = 2
2 1
(b) ny = 4
1
(c) ny = 6
(d) ny = 8
Figure 10: Empirical coverage (top) and prediction set volume (bottom) for the SCM20D dataset across output dimensions ny ∈ {2, 4, 6, 8} over 100 calibration runs, calibrated at ε = 0.05. The black diamond marks the mean.
Table 3: Computation time [s] for the SCM20D dataset across output dimensions (mean ± standard deviation over 100 runs). Method ny = 2 ny = 4 ny = 6 ny = 8 SCP-1 SCP-2A SCP-3A RemMCP RelMCP-A
0.0001 ± 0.0000 0.0004 ± 0.0001 0.0003 ± 0.0000 0.0751 ± 0.0023 0.5514 ± 0.1792
0.0001 ± 0.0000 0.0008 ± 0.0001 0.0009 ± 0.0001 0.1061 ± 0.0080 1.3431 ± 1.0503
22
0.0001 ± 0.0000 0.0008 ± 0.0003 0.0008 ± 0.0002 0.1186 ± 0.0156 6.7078 ± 3.5835
0.0001 ± 0.0000 0.0009 ± 0.0002 0.0009 ± 0.0003 0.1096 ± 0.0131 32.4974 ± 14.5166
F
Limitations and Future Work
While MCP offers a principled and flexible framework for joint prediction set design and calibration, several limitations merit discussion. The most immediate practical limitation is computational: by jointly optimizing shape and threshold over the full calibration set, MCP incurs substantially higher calibration cost than SCP, which calibrates a single scalar threshold in closed form. This cost grows with the number of optimization variables nq and the size of the calibration set ncal , and is most pronounced for RelMCP, which requires repeated solving of penalized subproblems during the iterative penalty search. In settings where calibration time is critical, SCP remains the more practical choice. A related limitation concerns data requirements. SCP baselines can leverage the training dataset for shape parameter estimation. MCP, by contrast, relies entirely on the calibration set for both stages simultaneously. While this is precisely what eliminates the variance-inflating data split, it also means that MCP requires a sufficiently large calibration set to certify coverage — a requirement that becomes more demanding as nq grows. Within the two proposed variants, RemMCP is restricted to convex score and cost functions, which limits the prediction set shapes it can directly optimize. While convexity covers a broad and practically relevant class of shapes — including ellipsoids, zonotopes, and hyperrectangles — it excludes more expressive non-convex parameterizations. Moreover, the outlier budget ρ scales inversely with nq , so RemMCP becomes increasingly conservative as the number of optimization variables grows, since more calibration data is consumed to certify each additional degree of freedom. RelMCP supports arbitrary score functions, but introduces its own limitations: the iterative penalty search is not guaranteed to find a feasible certified solution and the inherent conservatism of the procedure depends sensitively on the initialization of the penalty parameter. Finally, while we provide a general framework and establish its theoretical guarantees, we do not offer an extensive treatment of score and cost function design beyond a representative set of instantiations. The space of admissible score and cost functions is large, and the choice of parameterization has a significant effect on the expressiveness, computational tractability, and conservatism of the resulting prediction sets. These limitations point to several directions for future work. On the theoretical side, tightening the finite-sample bounds for RelMCP and developing initialization strategies with stronger convergence guarantees would strengthen its practical reliability. On the applied side, a systematic investigation of score and cost function designs for specific output structures — including time series, classification, and structured prediction tasks — would substantially broaden the applicability of MCP. Beyond the two variants proposed here, other calibration schemes could be developed that may offer better trade-offs between conservatism and computational complexity. Extending MCP to settings with distribution shift or covariate shift is another natural and practically important direction.
23