Multi-Variable Conformal Prediction

• Multi-variable conformal prediction is a framework that extends traditional scalar methods to vector-valued and structured outputs for joint prediction regions.
• It employs techniques like coordinate-wise decomposition, copula models, and optimal transport to capture dependencies across multiple output dimensions.
• The approach balances finite-sample coverage guarantees with computational challenges, offering practical calibration strategies for complex predictions.

Multi-variable conformal prediction generalizes conformal prediction from scalar responses to vector-valued and structured outputs. In this setting, the prediction target may be a multivariate response $Y \in \mathbb{R}^d$, a bounding box together with a semantic class, or an entire future trajectory over multiple forecast horizons. The core objective remains the construction of a prediction region $C(x)$ such that $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$, but the classical scalar architecture of “one score, one threshold, one interval” becomes inadequate because the output space is multidimensional, structured, or both (Johnstone et al., 2022, Ries et al., 8 May 2026, Lützow et al., 12 May 2026). A second fundamental difficulty is that conformal prediction is rank-based, whereas $\mathbb{R}^d$ has no canonical total order; recent work therefore develops multi-variable conformal prediction through coordinate-wise decomposition, direct multivariate region construction, optimal-transport-based ranking, density-based scoring, and scenario-theoretic calibration with multiple optimization variables (Klein et al., 5 Feb 2025, Thurin et al., 31 Jan 2025).

1. Formal problem and validity notions

The basic multivariate regression formulation treats the response as a vector $Y_i\in\mathbb{R}^d$, with exchangeable observations $Z_i=(X_i,Y_i)$, and seeks a prediction region $C(x)\subseteq \mathbb{R}^d$ satisfying

$$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$

In full conformal inference, one evaluates a candidate $y\in\mathbb{R}^d$ by augmenting the sample with $(x,y)$, recomputing residual-based scores, and retaining those $C(x)$0 whose conformal $C(x)$1-value exceeds the threshold. Validity remains distribution-free under exchangeability, but the prediction set is now a subset of $C(x)$2 rather than an interval (Johnstone et al., 2022).

Structured outputs sharpen the distinction between marginal and joint validity. In probabilistic object detection, each output is a bounding box

$C(x)$3

together with a class label $C(x)$4. The relevant notions of validity include coordinate-wise marginal coverage, box-level coverage for simultaneous coverage of all four coordinates, class-conditional box coverage, and a joint class-and-box guarantee in a sequential classification-plus-regression pipeline. The full box is correct only if all coordinates are simultaneously covered, and in the joint setting the true class should also lie in a prediction set (Ries et al., 8 May 2026).

Multi-step forecasting provides a third canonical instance. There the target is a future path

$C(x)$5

while the predictor outputs

$C(x)$6

This is a multi-output object rather than a single-step scalar forecast. The literature stresses that naïvely pooling all horizons ignores heterogeneity across horizons, while treating each horizon independently ignores dependence among neighboring forecast steps (Yu et al., 27 Mar 2025).

2. Componentwise and stratified constructions

A common strategy is to decompose the structured output into scalar components, conformalize each component, and then lift the scalar guarantees to a joint guarantee by explicit coverage accounting. For scalar regression, split conformal uses either the absolute residual score

$C(x)$7

or the scaled score

$C(x)$8

with conformal quantile

$C(x)$9

The resulting intervals are

$\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$0

or

$\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$1

In object detection, this is applied separately to the four box coordinates, and Bonferroni yields box-level validity: $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$2 which for $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$3 becomes $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$4. The method is explicitly conservative because it avoids modeling coordinate dependence; class-wise calibration strengthens the target to $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$5 but requires enough calibration samples per class. The same paper also studies a two-step class-plus-box pipeline in which RAPS produces a class prediction set and the localization step uses the largest relevant class-specific localization quantile over that set; the quoted joint class-and-box guarantee factorizes only under an independence assumption between the classification and regression conformal procedures (Ries et al., 8 May 2026).

Copula-based conformal prediction addresses the same joint-coverage problem from a different angle. In multi-target regression, it models dependence among target-wise normalized nonconformity variables

$\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$6

through a copula $\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$7, so that

$\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$8

determines target-wise significance levels consistent with a desired global confidence. The resulting prediction region remains a hyper-rectangle, but its widths are calibrated jointly rather than under independence. Empirically, the independent copula is typically conservative, the Gumbel copula is intermediate, and the empirical copula is usually closest to nominal global validity while producing the smallest hyper-rectangle volumes (Messoudi et al., 2021).

For multi-step forecasting, Dual-Splitting Conformal Prediction (DSCP) also follows a decomposition strategy but adapts the residual pool both across forecast trajectories and across horizons. It uses signed residuals

$\mathbb{P}(Y_{n+1}\in C(X_{n+1}))\ge 1-\alpha$9

clusters calibration forecast vectors with k-means while selecting the number of clusters by silhouette score, and merges adjacent horizon residual pools when a Kolmogorov–Smirnov test satisfies $\mathbb{R}^d$0. Intervals are then constructed from lower and upper residual quantiles of the selected merged residual subsets. The paper states the usual conformal objective

$\mathbb{R}^d$1

but it does not provide a formal theorem establishing finite-sample marginal or joint coverage for DSCP after clustering and adaptive merging (Yu et al., 27 Mar 2025).

3. Exact and geometric joint regions

A more direct line of work constructs multivariate prediction regions in $\mathbb{R}^d$2 rather than products of scalar intervals. In full conformal multi-output regression, the basic score is a residual norm,

$\mathbb{R}^d$3

and the prediction region is

$\mathbb{R}^d$4

When the fitted predictor depends affinely on the candidate response,

$\mathbb{R}^d$5

each residual comparison $\mathbb{R}^d$6 becomes a quadratic inequality

$\mathbb{R}^d$7

The exact conformal region is therefore a rank-thresholded subset of $\mathbb{R}^d$8 induced by quadratic surfaces. Because this geometry is expensive to compute, the same work proposes $\mathbb{R}^d$9, which aggregates many scalar conformal problems defined on projections $Y_i\in\mathbb{R}^d$0, and a multivariate extension of $Y_i\in\mathbb{R}^d$1, which approximates the region boundary by one-dimensional root-finding along rays (Johnstone et al., 2022).

For dependent multivariate time series, ellipsoidal regions provide a different reduction. MultiDimSPCI defines a covariance-normalized quadratic nonconformity score

$Y_i\in\mathbb{R}^d$2

with $Y_i\in\mathbb{R}^d$3, and calibrates future score quantiles sequentially by fitting a quantile regression model to the time series of scalarized residual scores. The resulting prediction region can be written as an ellipsoid or ellipsoidal shell centered at $Y_i\in\mathbb{R}^d$4. The paper gives finite-sample high-probability bounds on the conditional coverage gap under i.i.d. and strong-mixing assumptions, and reports that local covariance estimation can reduce average set size by up to $Y_i\in\mathbb{R}^d$5 (Xu et al., 2024).

A closely related but feature-adaptive construction arises from conformalized Gaussian scoring. If the fitted conditional model is

$Y_i\in\mathbb{R}^d$6

then the highest-density-mass score reduces to the Mahalanobis norm

$Y_i\in\mathbb{R}^d$7

Split conformalization of this score yields feature-adaptive ellipsoids

$Y_i\in\mathbb{R}^d$8

The same Gaussian machinery extends to missing output coordinates, refinement after partial revelation of outputs, and conformal prediction on linear transformations $Y_i\in\mathbb{R}^d$9 of the output (Braun et al., 28 Jul 2025).

4. Optimal transport, center-outward ranking, and credal interpretations

Optimal transport methods address a basic obstacle of multivariate conformal prediction: vector-valued scores have no canonical total order. OTCP resolves this by transporting the multivariate score distribution to the uniform distribution on the unit ball and then ranking scores by the radius of their transported images. With center-outward distribution $Z_i=(X_i,Y_i)$0, the rank of $Z_i=(X_i,Y_i)$1 is

$Z_i=(X_i,Y_i)$2

and the conformal score becomes

$Z_i=(X_i,Y_i)$3

Finite-sample marginal validity is preserved even for approximate transport maps by calibrating the acceptance radius on the empirical distribution of transformed norms rather than relying on exact recovery of the population map (Klein et al., 5 Feb 2025).

An allied construction based on Monge–Kantorovich vector ranks defines reference vectors

$Z_i=(X_i,Y_i)$4

with $Z_i=(X_i,Y_i)$5 uniform on the sphere, and obtains empirical multivariate ranks through an optimal assignment between calibration scores and references. The induced quantile region is

$Z_i=(X_i,Y_i)$6

For split OT-CP, the prediction set is the preimage of this score-space quantile region. The main finite-sample guarantee is stated, in the notation of the paper, as

$Z_i=(X_i,Y_i)$7

while OT-CP+ replaces the global rank map by a $Z_i=(X_i,Y_i)$8-nearest-neighbor conditional rank map and obtains asymptotic conditional coverage (Thurin et al., 31 Jan 2025).

A different conceptual reframing treats conformal prediction as implicitly inducing a credal set of predictive distributions. Under a consonant conformal transducer $Z_i=(X_i,Y_i)$9, split CP defines an upper probability

$C(x)\subseteq \mathbb{R}^d$0

with a corresponding core credal set $C(x)\subseteq \mathbb{R}^d$1. The conformal prediction region is then exactly the $C(x)\subseteq \mathbb{R}^d$2-imprecise highest density region of that credal set. Building on this view, Maximum Mean Imprecision quantifies disagreement inside the induced credal set; for classification it reduces either to the second-largest conformal $C(x)\subseteq \mathbb{R}^d$3-value or to a closed form based on sorted $C(x)\subseteq \mathbb{R}^d$4-values, whereas for common regression scores it becomes constant across test instances, a limitation that the paper argues is likely to be even more pronounced for multivariate continuous outputs (Chau et al., 2 Feb 2026).

5. Dependence, localization, and distribution shift

A general theoretical synthesis interprets conformal methods as procedures that reveal partial information $C(x)\subseteq \mathbb{R}^d$5 about the full data vector and then posit an exact or approximate conditional law $C(x)\subseteq \mathbb{R}^d$6. In this framework, the unified $C(x)\subseteq \mathbb{R}^d$7-value is

$C(x)\subseteq \mathbb{R}^d$8

and exact validity holds if $C(x)\subseteq \mathbb{R}^d$9. More generally,

$$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$0

with a coarser bound in terms of $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$1. The framework recovers standard conformal prediction, weighted conformal prediction, nonexchangeable conformal prediction, and randomly-localized conformal prediction, and it is largely agnostic to the label space $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$2; in principle, it therefore supports vector-valued and structured outputs once a scalar score is defined on the joint output space (Barber et al., 3 Apr 2025).

Under covariate shift, KMM-CP replaces explicit density-ratio estimation by Kernel Mean Matching on the covariates. In the setting

$$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$3

it estimates calibration weights by minimizing RKHS moment discrepancy subject to bounded-weight constraints, then computes a weighted conformal quantile. Its coverage-error bound decomposes into an MMD term and an effective-sample-size term: $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$4 The selective extension SKMM introduces target-side selection variables and restricts conformal correction to a subset of the target covariate space with more reliable support overlap. The paper is explicitly described as only indirectly relevant to multi-variable conformal prediction: the weighting machinery is output-agnostic, but the paper does not develop a multivariate prediction-set construction (Laghuvarapu et al., 27 Mar 2026).

Time-series methods illustrate different responses to nonexchangeability. MultiDimSPCI replaces exchangeability by sequential score calibration and obtains finite-sample high-probability conditional coverage-gap bounds under strong mixing, whereas DSCP uses conformal-style calibration on clustered and merged residual pools but leaves the multi-horizon validity theorem implicit. This contrast suggests a substantive divide between methods whose dependence handling is purely empirical and methods that make the dependence structure part of the theoretical coverage analysis (Xu et al., 2024, Yu et al., 27 Mar 2025).

6. Multi-variable calibration without data splitting

The paper that explicitly names the topic introduces Multi-Variable Conformal Prediction (MCP) as a generalization of split conformal prediction from “one scalar score + one scalar threshold” to vector-valued scores and multiple simultaneous calibration variables. Instead of a scalar thresholded set

$$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$5

MCP defines

$$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$6

where $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$7 and $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$8. The parameters $$\mathbb P\!\left(Y_{n+1}\in C(X_{n+1})\right)\ge 1-\alpha.$$9 can control shape, orientation, position, scale, or multimodal structure, and they are chosen by solving a scenario-theoretic optimization problem over the full calibration set rather than by fitting shape parameters on one split and conformalizing a scalar threshold on another (Lützow et al., 12 May 2026).

RemMCP is the convex variant. It solves a cascade of constrained optimization problems with support-constraint removal, using an outlier budget

$y\in\mathbb{R}^d$0

The framework yields expected-coverage and high-confidence guarantees through scenario theory, and under stronger assumptions the conditional coverage probability satisfies

$y\in\mathbb{R}^d$1

Classical split conformal prediction appears as the special case

$y\in\mathbb{R}^d$2

so SCP is literally embedded in RemMCP. RelMCP replaces hard feasibility by a penalized objective with slack variables,

$y\in\mathbb{R}^d$3

and supports non-convex score functions, including multimodal unions, at the cost of greater conservatism and a weaker theorem for the adaptive penalty-search implementation (Lützow et al., 12 May 2026).

The empirical motivation for MCP is that using all calibration data jointly for shape optimization and calibration can reduce variance across calibration runs while keeping coverage at target levels. In the reported experiments, RemMCP and RelMCP met target coverage with prediction set sizes smaller than or comparable to split-data baselines, and they substantially reduced run-to-run variance. At the same time, the paper stresses that the burden of certification grows with the parameter dimension $y\in\mathbb{R}^d$4: for ellipsoids, increasing output dimension increases $y\in\mathbb{R}^d$5, reduces the allowable removal budget, and can make the method conservative or infeasible unless the calibration set is sufficiently large (Lützow et al., 12 May 2026).

7. Limitations and open directions

Several limitations recur across the literature. Coordinate-wise conformalization with Bonferroni correction is simple and assumption-free, but it is conservative because it ignores dependence among output coordinates. Copula-based alternatives can be tighter, yet they require estimating the dependence structure of target-wise nonconformity scores and remain axis-aligned at the level of the final hyper-rectangle. Class-wise or class-conditional calibration strengthens validity targets but can destabilize quantile estimates when calibration samples are sparse (Ries et al., 8 May 2026, Messoudi et al., 2021).

Direct multivariate methods face geometric and computational trade-offs. Exact full conformal multi-output regression becomes combinatorially difficult because the acceptance region is a subset of $y\in\mathbb{R}^d$6 defined by many quadratic boundaries. OT-based approaches avoid fixed convex templates, but the literature emphasizes the curse of dimensionality in transport-map estimation and the resulting loss of efficiency in higher output dimension. Ellipsoidal methods encode covariance structure but may be misspecified when the residual distribution is strongly non-elliptical, while Gaussian scoring remains most natural when a Gaussian or broader elliptical local approximation is reasonable (Johnstone et al., 2022, Klein et al., 5 Feb 2025, Xu et al., 2024, Braun et al., 28 Jul 2025).

A further limitation is that exact conditional coverage remains unattainable without strong assumptions, so most multi-variable methods trade distribution-free finite-sample marginal validity for approximate conditional behavior, asymptotic guarantees, or empirical calibration. DSCP explicitly illustrates this boundary: it produces horizon-wise intervals informed by trajectory type and horizon similarity, but without a formal finite-sample theorem for the resulting multi-output procedure. KMM-CP addresses covariate-side mismatch rather than output-side multivariate geometry. Credal-set approaches show that conformal prediction can support a principled notion of epistemic predictive uncertainty, but the strongest tractable results currently lie in finite-label classification rather than multivariate continuous regression (Yu et al., 27 Mar 2025, Laghuvarapu et al., 27 Mar 2026, Chau et al., 2 Feb 2026).

Taken together, these works define multi-variable conformal prediction not as a single method but as a design space. One branch decomposes structured outputs into scalar parts and reassembles them with explicit multiplicity correction. Another builds joint geometric regions in $y\in\mathbb{R}^d$7 through quadratic, ellipsoidal, transport-based, or density-based constructions. A third changes the calibration architecture itself, replacing scalar thresholding by joint optimization over multiple calibration variables. The common theme is that once the output is genuinely multivariate, conformal prediction is no longer only a question of threshold calibration; it becomes a question of how to represent geometry, dependence, and heterogeneity while preserving a defensible finite-sample coverage statement.