A Unified Comparative Study with Generalized Conformity Scores for Multi-Output Conformal Regression
Authors: Victor Dheur
Matteo Fontana
Yorick Estievenart
Naomi Desobry
Souhaib Ben Taieb
Venue: ICML 2025 poster (2025)
Keywords: Uncertainty Quantification
Multivariate regression
Conformal Prediction
Generative Models
Primary area: probabilistic_methods->everything_else
TL;DR: We present a unified study comparing multi-output conformal regression methods and introducing two generalizations of univariate conformity scores.

Abstract: Conformal prediction provides a powerful framework for constructing distribution-free prediction regions with finite-sample coverage guarantees. While extensively studied in univariate settings, its extension to multi-output problems presents additional challenges, including complex output dependencies and high computational costs, and remains relatively underexplored. In this work, we present a unified comparative study of nine conformal methods with different multivariate base models for constructing multivariate prediction regions within the same framework. This study highlights their key properties while also exploring the connections between them. Additionally, we introduce two novel classes of conformity scores for multi-output regression that generalize their univariate counterparts. These scores ensure asymptotic conditional coverage while maintaining exact finite-sample marginal coverage. One class is compatible with any generative model, offering broad applicability, while the other is computationally efficient, leveraging the properties of invertible generative models. Finally, we conduct a comprehensive empirical evaluation across 13 tabular datasets, comparing all the multi-output conformal methods explored in this work. To ensure a fair and consistent comparison, all methods are implemented within a unified code base.

A Unified Comparative Study with Generalized Conformity Scores for Multi-Output Conformal Regression
Victor Dheur' Matteo Fontana Yorick Estievenart'
Naomi Desobry'
Souhaib Ben Taieb 13
Abstract
3
Conformal prediction provides a powerful frame- work for constructing distribution-free prediction regions with finite-sample coverage guarantees. While extensively studied in univariate settings, its extension to multi-output problems presents additional challenges, including complex output dependencies and high computational costs, and remains relatively underexplored. In this work. we present a unified comparative study of nine conformal methods with different multivariate base models for constructing multivariate predic- tion regions within the same framework. This study highlights their key properties while also ex- ploring the connections between them. Addition- ally, we introduce two novel classes of conformity scores for multi-output regression that generalize their univariate counterparts. These scores ensure asymptotic conditional coverage while maintain- ing exact finite-sample marginal coverage. One class is compatible with any generative model, of- fering broad applicability, while the other is com- putationally efficient, leveraging the properties invertible generative models. Finally, we conduct comprehensive empirical evaluation across 13 tabular datasets, comparing all the multi-output conformal methods explored in this work. To en- sure fair and consistent comparison, all methods are implemented within a unified code basd")
PCP
C-PCP
DR-CP
x
0
C-HDR
M-CP
-1
L-CP
-2
-3
0
Y1
Figure I: Examples of bivariate prediction regions with an 80% coverage level for a toy example.
1. Introduction
Quantifying uncertainty in model predictions is crucial in many real-world applications, often involving prediction problems with multiple output variables and complex statis- tical dependencies. For example, in medical diagnostics, the progression of discase can be studied by analysing mul- tiple health indicators that exhibit nonlinear dependencies, such as blood pressure and cholesterol levels of a patient (Rajkomar et al., 2018). Although modern probabilistic AI models can model complex relationships between variables, they may produce unreliable or overly confident predictions (Nalisnick et al., 2018).
Conformal prediction (CP) offers robust framework for improving model reliability by generating distribution-fig prediction regions with a finite-sample coverage guaran- tee (Vovk et al., 1999). Although substantial research has focused on univariate prediction problems (Romano et al., 2019 Sesia and Romano, 202T, Rossellini et al., 2024). multivariate settings have received less attention. Among existing work, Zhou et al. (2024) achieves marginal cover- age by combining univariate prediction regions, but fails to capture dependencies between variables. Other methods, such as density- based approaches (Izbicki et al., 2022) or sample-based techniques (Wang et al., 2023b] Plassier et al., 2025), suffer from high computational costs. An alterna-
'Department of Computer Science, University of Mons, Mons, Belgium Department of Computer Science, Royal Holloway, Uni- versity of London, Egham, United Kingdom DDpaarmeett of Statis- and Science, Mohamed bin Zayed University of Artificial Intelligence, Abu Dhabi, United Arab Emirates. Correspondence to: Victor Dheur evictor acmons.ac.be..
Proceedings of the ndnd International Conference on Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
httpp://github cow/Vekteur mu multi-outpute conformal seappptttoccoooommaeeeeeeesioooi

---

Multi-Oulput Conformal Regression
tive method (Sadinle al., 2019) optimises the size of the region, but does not achieve asymptotic conditional cover- age. For toy bivariate example, Figure D illlstrates the diversity of prediction regions obtained using a selection of conformal methods considered in this paper.
2. Background
Consider a multivariate regression problem where the objec- tive is to predict a d-dimensional response vector y E y Fd based on feature vector E X RP. We assume there exists a true joint distribution Fxy over X x V, and we have access to dataset D ((X() } where
Our first contribution a unified comparative study nine confornal methods with different multivariate base models for constructing multivariate prediction regions within the same framework. This study highlights their key properties while also exploring the connections between them. We ex- amine different conformity scores with different multivariate base predictors, discussing prediction regions derived from the marginal distributions of individual output variables, their joint PDF, or sampling procedures (e.g., generative models).
(XO), yú) iid. Fxy- Given feature vector x, denote the conditional distribution of Y given X = x as Fy|X=a> and the associated probability density function (PDF) as frxx=r.
Using the dataset D, for any x E X, CP allows us to trans- form base predictors, denoted h. into calibrated, distribution- free prediction regions R(x)) S y for the true output y with finite-sample coverage guarantees.
Our second contribution introduces two novel classes of conformity scores for multi-output regression that gener- alize their univariate counterparts, These scores ensure asymptotic conditional coverage while maintaining exact finite-sample marginal coverage.
2.1. Split-conformal prediction
Split oplit-conformal prediction (SCP, Papadopoulos et al., 22022) is computationally efficient variant of conformal predic- tion that divides the dataset D into two disjoint subsets: training set Dimin and calibration set Deat- A model is first trained on Durain to obtain a base predictor h. Based on h, a conformity score (function) 8 X' x y - R is defined, where lower scores indicate better fit between the feature vector x and the response y. The calibration scores S (s y) (e.))) U {+x0} are then computed, from which the (1 - a empirical quantile is calculated as:
The first, CDF-based scores, leverage the cumulative distri- bution function (CDF) of any conformity score to achieve asymptotic conditional coverage. This approach general- izes the univariate HPD-split score, based on univariate highest-density region from Izbicki et al., [020 to multivari- ate prediction regions derived from any conformity score. Additionally, we propose a specific instance of CDF-based scores that builds on PCP from Wang al.. 2023b This method avoids the estimation of a predictive density, instead relying solely samples from any generative model.
q Quantile (s:
ka [Deal
(1)
where ka = [[[Paal + 1)(1 - a)]. This quantile serves as the threshold for constructing prediction regions. For an input x. the (random) prediction region is given by:
The second, latent-based scores is inspired by Feldman et al., 2023 and can be interpreted as an extension of distri- butional conformal prediction (Chernozhukov et al., 2021 to multivariate outputs. Compared to Feldman et al., 2023 does not require directional quantile regression, and the conformalization is performed directly in the latent space, eliminating the need to construct grid. This enhances both computational efficiency and scalability.
A(()) {y E y s(x, v) S 9)
(2)
If the random pair (X, exchangeable with Deals SCP guarantees marginal coverage:
(Y. E R(X)) = Y(s(X,Y) S ②) 2 a, (3)
Finally, as our third contribution, we conduct a large-scale empirical study comparing the different multi-output con- formal methods across 13 tabular datasets with multivariate outputs, evaluating several performance metrics. We con- sider variety of multivariate regression models, namely Multivariate Quantile Function Forecaster (Kan et al., 2022). Distributional Random Forests (Cevid et al., 2022). and multivariate Gaussian Mixture Model parameterized by hypernetwork (Ha et al., 2022, Bishop, 1994).
where the probability taken over (X, and D As- suming no ties in scores, the marginal coverage is exactly |Pal+1. yielding P(Y E R(X)) < + p
Ke
Ideally, the prediction region should achieve conditional coverage at the level 1 - a, i.c.:
I(Y e x(X) xx
(4)
holds almost surely. This is stronger requirement than marginal coverage in (3). However, as Barber et al. (2019) demonstrate, achieving conditional coverage is generally impossible without making additional assumptions about the underlying geteraanng process.
2

---

Multi-Output Conformal Regression
2.2. Multi-output conformal methods
superlevel set of the probability mass function with proba- bility content at least q-
Many conformal prediction methods have been proposed in the literature and implemented within the SCP framework for various base predictors and conformity scores, with specific focus on univariate prediction problems. In this sec- tion, we survey several conformal methods for constructing multivariate prediction regions, using different multivari- ate base predictors and corresponding conformity scores. Specifically, we discuss density- -based, and sample-based methods, which are based on their joint PDF, or a sampling procedure (c.g.. a generative model), respectively, In the following, we describe the conformity scores s for different methods. The methods M-CP and CopulaCPTS, which pro- duce hyperrectangular regions, are detailed in Appendix [ Once a conformity score is defined, the corresponding pre- diction region R can be computed using (2). We detail this relationship for each method in Appendix C Furthermore, in Section 5 analyze the properties and relationships between these methods and provide illustrative examples of the resulting prediction regions.
PCP. Let y(u), 8112) (L) denote sample with L.
points from the (estimated) conditional distribution Fylx==- Probabilistic Conformal Prediction (PCP, Wang et al.,
2023b defines a conformity score as the closest distance to
y:
Spep(x, v) min ly [E[L]
y00 1|2
(9)
where y Fyxx=s LE [L].
(10)
The corresponding region is a union of balls centered at each sampled point y i.e. RFer(z) Uletz (y e y y 101 S a).
HD-PCP.
When predictive density is available alongside sample of L points, Wang et al., 2023b proposed ex- tension to PCP, called PCP. This method uses the same conformity score as in 9 but only retains the [(1 a L samples with the highest density, ensuring that the predic- tion region is concentrated on high-density points.
DR-CP.
Given a predictive density frxxxmrr a natural con- formity score is the negative density:
(5) The corresponding prediction region is density superlevel set, ( {y E y [y I x)) 2 -Q)- Sadinle
SDR-CP(I: (z. y) = =fy z).
ST-DQR. Motivated by the limitation that existing mul- tivariate quantile regression methods do not allow the con- struction of regions with arbitrary shapes, Feldman et al., 2023 proposed to construct convex regions in a latent space 2 using directional quantile regression (Paindaveine and Siman, 2011). These regions are then mapped to the out- put space y using conditional variational autoencoder (CVAE), allowing a non-linear mapping between the two spaces. Specifically, they apply a conformalization step by creating a grid of points within the region in 2, map the points to the output space D, and construct d-balls around the mapped samples, similarly to PCP.
et al. (2019) use this conformity score in the context of classification.
C-HDR. Izbicki et al., 2022 proposed the HPD-split method, which defines a conformity score based on the Highest Predictive Density (HPD):
HPD (v x) S
f(v xx) dy
(6)
= P (f0 x) 2 fly 12) x
3. Generalized Conformity Scores for Multi-Output Regression
where y ~ frix==. The corresponding prediction region highest density region (HDR, Hyndman, 1996) with respect ff level ::
In this section, we introduce two new classes of confor- mity scores: CDF-based and latent-based scores. These scores generalize existing conformity scores for univariate regression to accommodate any conformity score for multi- variate outputs. The former generalizes HPD-split (Izbicki et al., 2020) to any conformity score, allowing to apply this method multivariate outputs. We further propose a specific instance that builds on PCP (Wang et al., 2023b). The latter is inspired by Feldman al., 2023 and can be interpreted as an extension of distributional conformal pre- diction (Chernozhukov et al., 202T) for multivariate outputs. Section S will present a comparative study of the conformity scores introduced in Section 2.2 alongside those introduced in this section.
RCHHDR (x) = {u E y: f(y x) 2 tg].
(8)
where tq sup(t: (f xx) 2 t x x= 2) 4)
Compared to DR-CP, where the threshold -ÿ is independent of x. C-HDR allows the threshold tq vary with x. To compute the HPD in ( Izbicki et al., 2022 use numerical integration, whereas in our experiments, we approximate using Monte Carlo sampling, as described in (3)..
In the context of classification, Adaptive Prediction Sets (Romano et al., 2020) follows a similar principle by con- structing "highest mass region",, which corresponds a
3

---

Multi-Output Conformal Regression
3.1. CDF-based conformity scores
Interestingly, C-PCP shares similarities with the recently proposed CP2-PCP method by Plassier et al., 2025| For given z E X, both methods adapt the radius of the balls based on a second sample from the conditional distribution composed of K points, requiring a total of I. + K samples. A detailed discussion can be found in Appendix
Consider a conformity score sw, and define the random variable W Sw 5w(X,X (X,Y) for random pair (X, For an observation (z. y), we introduce a new conformity score based on the conditional CDF of W given X = x, evaluated at sw (z. Specifically, the score is given by
8CDF (x, y) P(sw (X, Y) S sw(a,y) X = x) FwixxaxSsp (x, y)).
(11)
3.2. Latent-based conformity scores
(12)
Inspired by Feldman et al. (2023), we propose a latent-based conformity score with key distinctions. First, our method does not require the use of directional quantile regression. Additionally, the conformalization step is performed in the latent space, eliminating the need construct grid, which improves both computational efficiency and scalability.
This new conformity score measures the rank of sw (z)
relative to the conditional distribution of W given X x.
This method applies to any conformity score sw and gen- eralizes the (oracle) HPD-split introduced in Izbicki et al., 2020 in the context of univariate regression. Specifically, when sw (I. y) SpR-cp (I. y) used (12) we recover the C-HDR method. Additionally, by the probability integral U(001) 1) for E X, meaning that the conformity score's 's distribution is indepen- dent of x. This property ensures that conditional coverage is achieved as |Deall ~ oo (see Appendix E.2 Lemma 2) A similar observation was made by Izbicki et al., 2020 for C-HDR.
Our base predictor is conditional invertible generative model Q : 2 x X + Y. which maps a latent random variable E z (e.g., drawn from a standard multivariate normal distribution) to the output space V, conditional on X E X' (e.g., using normalizing flows). The model is both conditional and invertible, meaning that
transform, $CDP(X,,)) Y) X
Q Q (v: y, Ve E X,yEy.
We propose the following conformity score, called L-CP (Latent-based Conformal Prediction), defined as:
However, in practice, since the distribution of Y X is unknown, we approximate SCDF using Monte Carlo sampling:
$1.CP(1, (z. y) (Q;
(14)
1 SECDE(2, v) K
"EI(snn (2) sw(2)))
where dz : 2 - R is conformity function in the latent space 2, independent of In our experiments, we Z NN Ia) and dz ||=||-
kE[K]
where (R) FrIx=zI k [K].
(13)
The corresponding prediction region is obtained by mapping a region in the latent space, (2( = (2 E 2 : dz(=) & q), to a region in the output space, (L = {((2;)) : R2 (q)).
Dheur et al., 2024 considered a particular case of this empir- ical CDF-based approach with the SDR-Cr score for a bivari- ate prediction problem involving temporal point processes, where the HDR is estimated via Monte Carlo sampling.
L-CP generalizes Distributional Conformal Prediction (Chernozhukov et al., 2021). which is special case when Y is univariate (d = 1), Z ~ (0, 1), dz(2) = |2- 11. and ((: is the quantile function of Y given x.
We introduce a special case of our new score, called C- PCP (ased Probabilistic Conformal Predic- tion), by setting sw (x, v) spCp(x, y) in (13 which gives:
C-PCP.
Concurrent work by Fang al., 2025 introduces CONTRA sharing the same algorithm as our latent-based methods. While related, the papers diverge in their primary focus. Fang et al., 2023 emphasizes the smaller prediction regions achieved by CONTRA, whereas our work concentrates on the computational complexity and conditional coverage guarantees of the latent-based methods while obtaining re- gion sizes that are small but not smaller than censity-based. methods.
Sc.Cr((x))) 8)
Σ kE[K]
T
min ||[(8]_ 101 s min Ily (011 SE[L]
K
[EELL
Compared to the methods in Izbicki et al., 2020 and Dheur et al.. 2024, this score has the advantage of not requiring the estimation of a predictive density, relying instead on samples from the conditional distribution. Consequently, this score can be applied with any generative model that does not have an explicit density, while still retaining the desirable properties of our based score.
4. Related Work
Conformal Prediction (CP), introduced by Vovk et al., 1999 forms the foundation of our work by providing prediction
4

---

Multi-Oulput Conformal Regression
Oracle
M-CP
CopulaCPTS
DR-CP
C-HDR
2
x
0
2
2
2
2
2
0
Ys
Y>
Y1
Yss
-2
15
PCP
HD-PCP
STDQR
C-PCP
L-CP
2
xx
x
-2
-2
2
Ys
2
Y>
2 Y5
2
2
Figure 2: Prediction regions for bivariate unimodal dataset, conditional on a unidimensional input. The black, green, and yellow contours represent regions with nominal coverage levels of 20%, 40%, and 80%, respectively.
regions with sample coverage guarantees. CP meth- ods are well established for regression with univariate out- puts (Papadopoulos al., 2008; Lei and Wasserman, [014] Romano et al., 2019); Sesia and Romano, 2021) and classi- fication (Romano et al., 2020, Angelopoulos et al.. 2202)). In the multi-output regression setting, we need to capture dependencies between output dimensions, represent more complex prediction regions and handle a larger computa- tional demand.
scores proposed in Section 3.1
5.1. Illustrative examples
We provide illustrative examples of bivariate prediction re- gions for different conformal methods on simulated data, covering both unimodal (Figure 2) and bimodal distribu- tions (Figure I0 in Appendix D.2p The data- aetta-genrrating pro- cesses are given in Appendix D.2 Additionally, we present bivariate prediction regions for a real-world application, pre- dicting a taxi passenger's drop-off location based on the passenger's information (Figures /
To address multivariate prediction challenges, optimal trans- port methods such cyclically monotone mappings (Car- lier et al.. 2016) define multivariate quantile regions with desirable properties such as existence and uniqueness of mappings. Hallin and Šiman (2017), Hallin et al. (2021), ), and Barrio et al. (2044) have proposed extensions of these approaches. Neural network-based techniques leverage nor-
and in Appendix D.I).
In both Figures 2 and [0] the black, green, and yellow con- tours represent prediction regions with nominal coverage levels of 20%, 40%, and 80%, respectively. The top-left panel illustrates the density level sets of the oracle distri- bution FYix. The remaining panels display the prediction regions generated by various conformal methods, all utiliz- ing the MQF base predictor, as explained in Appendix F.2
malizing flows (Kan et al., 22222) Huang et al., 2020) or variational autoencoders (Feldman et al., 2023) to learn flexible quantile regions. Additionally, highest density regions (HDRs) (Hyndman, 1996) handle multimodality and have been applicd in various contexts (Camchl et al., 2024, Izbicki et al., 2022] Dheur et al., (024). Recently, Wang al., proposed constructing prediction regions as hyperballs centered on generated samples, with exten- sions by Plassier et al., (202 improving conditional validity. Other methods use copulas (Messoudi et al., 2021a| Sun and Yu, 2024b) to model the dependency between vari- ables. A provides a more detailed discussion on related work.
We observe the following for the unimodal case in Figure M- CP and CopulaCPTS capture heteroscedasticity but pro- duce rectangular prediction regions, which do not align with the circular level sets of the oracle conditional distribution, resulting in a lack of sharpness. DR- CP fails to maintain conditional coverage, and for X = 1, the absence of black and green contours indicates that the predictive density does not reach the threshold - defined in (5) for coverage levels of 0.2 and 0.4. C-HDR generates prediction regions that closely resemble the oracle level sets. PCP generates highly discontinuous regions, especially at lower coverage levels, where the regions appear as balls centered on individual samples. In contrast, PCP and STDQR yield smoother, more continuous regions but require the estimation of a pre- dictive PDF or the identification of map from the latent space to the output space, respectively.
5. Comparison of Multi-Output Conformal Methods
In this section, we present unified comparison of the con- formity scores introduced in Section 2.2 and the generalized
5

---

In this section, we compare conformal methods based on several key properties. In the following, use - denote Multi-Output Conformal Regression
Table 1: Properties of different multivariate conformal methods. (*) M-CP achieves ACC under certain assumptions
(Appendix p. (**) STDQR and L-CP require a conditional invertible generative model Q
Z x d' + y (+)
CopulaCPTS has pre-training cost of O(C).
Method
Type of region
Asymptotic conditional coverage
Computational complexity
Predictive density not required
Sampling procedure not required
M-CP
Hypernectangle Hyperrectangle
X(*)
O(dM)
CopulaCPTS
x
O(dM)
DR-CP
Density superlevel set
X
O(D)
X
X
C-HDR
Density superlevel set
K +
O(((++5)
PCP
Union of d-balls
X
O(LS)
X
HD-PCP
Union of d-balls
O(L(D +))
X
X
X
STDQR
Union of d-balls
X
O(LS)
((#s)
C-PCP
Union of d-balls
K K->0
O((K L)S) O(Q)
L-CP
Quantile region
((*8)
For our methods, unlike PCP, C-PCP adjusts the radius of the prediction regions improve conditional coverage. This is evident in the example, where the radius of the balls for X = -1 is smaller than for X = 1, as indicated by the tighter regions around the samples. L-CP generates prediction regions that closely align with the oracle level sets, demonstrating good conditional coverage.
the assumptions that Deall - oc and the base predictor corresponds the oracle distribution Fyix.
While the assumption of oracle base predictor is strong, it is crucial to demonstrate that the conformal procedure pre- serves the performance of the base model. Specifically, given I E X, for M-CP and CopulaCPTS, we assume (x(x) Qx1x==((a) and (x) Qx,x=== (a.) with 1....,d; d: for DR-CP, CDRD, and HD-PCP, fr|x== JYIx=zi for L-CP, L- Q(Z;+) d. Y x;; and for PCP and C-PCP, Fyix = Frix
For the bimodal distribution in Figure 10 (Appendix D.2 the prediction regions generated by M-CP and L- CP are connected, failing to capture the bimodal nature of the dis- tribution. For the real-world application, Figures X and 9 (Appendix D.I illustrate predictions under low and high uncertainty, respectively. Our methods, L-CP and C-PC, alongside M-CP and C-HDR, demonstrate the best adaptabil- ity to outputs with varying levels of uncertainty.
Our empirical results (Section 6 demonstrate that meth- ods achieving ACC under these assumptions also exhibit superior approximate conditional coverage across diverse datasets and base predictors. L-CP is the only method that achieves ACC without additional assumptions. C-HD and C-PCP achieve ACC with K + oo. Finally, M- CP achieves ACC under specific assumptions. Assuming that Yi,... .......Y. Ya are conditionally independent given X, M-CP achieves ACC if a. a 4/1 a. Furthermore, under the unrealistic assumption that Y1 X
5.2. Properties
equality distribution and a.s. to denote almost sure equality.
u.s. as Yu X, M- achieves AC a 1 a. The true dependence typically lies between these two extremes. We provide detailed proofs of these statements in Appendix E.2
Marginal coverage. All the conformal methods presented achieve the classical finite-sapll marginal coverage. But, as noted by Wang et al., 2023b (Theorem 1), the marginal coverage of methods such as C-HDR, PCP, HD PCP, and C-PCP also depends on the randomness of the generated samples. Appendix E. we demonstrate that the marginal coverage, conditional on the calibration dataset Deal and the samples drawn from it, follows beta distribution, using standard arguments. CopulaCPTS is the only method that does not enter into the standard split-conformal algorithm and who does not satisfy the above property.
As discussed in Section 5.1 DR-CP fails to achieve ACC. Likewise, PCP, PCP and STDQR do not achieve ACC, as they are constrained bounded volume for any I E X. Assuming each ball has volume of V. PCP generates regions with total volume of at most LV. For a given instance x E X' with high uncertainty, it may be impossible to capture sufficient probability mass to achieve conditional coverage.
producing regions with upper
Asymptotic conditional coverage (ACC), We examine the asymptotic conditional coverage property, which cor- responds conditional coverage as defined in under
Region size. Among the methods that achieve ACC, CDHDR is expected to perform best, it converges to the highest
6

---

Multi-Output Conformal Regression
density regions, which correspond the smallest volume re- gions (Hyndman, (996). Prediction regions from C-PCP are expected to have larger volume since they are constrained to union of L d-balls. Similarly, prediction regions from L-CP are less flexible than those from C-HDR, as they are connected when the region (A) in the latent space is connected for all À E I and Q is continuous. This con- straint may be desirable when more interpretable regions are preferred (Sesia and Romano, (021].
quire the model be invertible. We summarize the different properties in Table
5.3. Connection between sample-hased and density-based methods
Interestingly, the sample-based methods (PCP, HD-PCP, C-PCP) can be viewed as special cases of density-based methods (DR- CP, C-HDR). Let us assume common pre- dictive PDF f is used for the base predictor of these conformal methods. While PCP and C-PCP do not re- quire a PDE, we assume that frix== and Fylx== cor- respond to the same distribution. Let ÿ FY|x== for l E [L] and fe(:; (()) be a PDF with spherical level sets, centered at Y such as a standard multivariate Gaus- sian (1; y Ia)- For E X, we define a new PDF fmax(y x) ma fs(y; / where C' is a nor- malizing constant ensuring that fmax(- | x)) integrates to 1. The following proposition establishes the relationship between these methods.
Among the remaining methods, DR- CP minimizes the mean region size B| [[] under the oracle PDF as (Peal + oo, as shown in Theorem by Sadinle et al., 2019 In con- trast, M-CP and CopulaCPTS are expected to yield larger prediction regions, as they do not explicitly account for dependencies between outputs. While PCP, HD-PCP, and C-PCP can capture multimodality, they are susceptible to the randomness of the sampling procedure, as evidenced by the shape of the regions in F Furthermore, since they rely on finite union of L d-balls, they are subject to the curse of dimensionality in high- sigh-dimensional spaces, where data sparsity necessitates larger balls to maintain marginal coverage.
Proposition 1. PCP is equivalent to DR-CP with fmax. Similarly, HD-PCP is equivalent to DR-CP with f fmux where only [(1 - L] samples with the highest density among (Y [e(1]] are kept. Finally, C-PCP is equivalent to C-HDR with f fmax.
A potential weakness of the mean region size is that it can be disproportionnaelly skewed by inputs with high uncertainty. To mitigate this sensitivity, we also report the median region size as more robust alternative.
Computational complexity. Table reports the compu- tational complexity of each conformity score. For M-CP and CopulaCPTS, let M represent the compute time of the univariate conformity score for a single dimension and C' the optimization time for CopulaCPTS. Let D, S, and Q denote the time required for density evaluation, sam- pling, and calculating the inverse of the quantile function Q-1, respectively. In many cases, M and C' are relatively low, while D. S, and Q are comparable. C-HDR, PCP, HD-PCP. STDQR and C-PCP are significantly slower than M-CP. L- LPC and DR- CP since they need to generate large number of samples to compute the conformity score (we used K- = L = 100 in our experiments).
We provide a proof in Appendix E3 sample-based methods are special cases of density-based ap- proaches, the key advantage of PCP and C-PCP that they rely solely a sampling procedure, without requiring predictive density f as base predictor. Figure summarizes the connections between the main conformal methods.
Although these
An interesting practical takeaway is that DR- CP and C-HDR are linked in the same way as PCP and C-PCP Since DR- CP under the oracle has the smallest expected region size while CDHDR empirically has a smaller median region size, similar observations are expected for PCP and C-PCP. This is veri- fied empirically: PCP has smaller mean region size across all base predictors, while PCP has smaller median re- gion size.
Base predictor. Some conformal methods stand out be- cause they do not need to evaluate the predictive density f or generate samples. M-CP and Copul only require univariate model for each dimension, without needing model for the joint distribution of Y..YY. DR- CP does not require sampling from the model, which is beneficial when using normalizing flows that are slower to invert (e.g., Masked Autoregressive Flows (MAF, Papamakarios, Pavlakou, et al., 2017) or Convex Potential Flows (Huang et al., 2020)). PCP and C-PCP do not require evaluating the predictive density f. making them compatible with any generative model, including diffusion models and GANs. L-CP and STDQR do not require predictive density evaluation but re-
f fmxx
DR-CP
PCP
SECDP
SECDF
C-HDR
C-PCP
f
Figure 3: Connections between different methods.
7

---

Multi-Oulput Conformal Regression
M-CP
CopUlaCPTS
DR-CP
C-HDR
PCP
HD-PCP
STDQR
C-PCP
L-CP
10a
0.9
CEC- (x200)
10a
CEC-Z (X100)
0.8
109
WSC
10-)
0.6
households scm20d
mi
scmia
meps 21
maje 19 mips house
ble blog. data calco
taxt
hauseholds semand
semid
maps 21
maps 19 meps 20 house
bio blog data cal calcon taxt
hauseholds
scm2nd
n12
meps 21 maps 19
meps 20
house
bia blog data calcoft Laxt
Figure 4: Conditional coverage metrics across datasets sorted by size. CEC-X and CEC-Z should be minimized while WSC
should approach a.
6. A Large-Scale Study of Multi-Output Conformal Methods
and CopulaCPTS attain intermediate conditional coverage, with M- CP performing slightly better. These results align with our analysis in Section [.2 where showed that C-PCP, L-CP, and C-HDR achieve ACC, while HD-PCE STDQR, PCP, and CR-CP do not, and M- CP achieves it only under specific conditions. Finally, Figure 12 shows that all methods achieve marginal coverage, as expected.
In this section, we present a large-scale study of multi-output conformal methods using 13 tabular datasets from previous studies (Tsoumakas et al., 2011 Feldman et al., 2023 Wang et 2023b Barrio et al., 2024) Camehl et al., 2024). To ensure sufficient data for training, calibration, and testing, we include only datasets with at least 2,000 instances. The selected datasets contain between 7,207 and 50,000 data points, with the number of input features p ranging from to 279 and the number of output variables d ranging from 2 to 16.
Region size. Figure presents critical difference (CD) diagram (Demšar, 2006) comparing the median region size of all methods across datasets. Higher-ranked methods (further right) perform better. Thick horizontal lines indicate models with no statistically significant difference the 0.05 level (see Appendix [.55 for details).
We consider three base predictors: the Multivariate Quantile Function Forecaster (MQF2), a normalizing flow (Kan et al., 2022). Distributional Random Forests (Cevid et al., 2022),
Among the methods that achieve ACC, C- yields the smallest median region size, as expected, since its regions converge to the highest density regions (Izbicki et al., 2022). C-PCP and L-CP produce slightly larger regions, though the difference is not significant for these datasets. Among the remaining methods, DR-CP RR-C yields the smallest median region. In contrast, M- M-CP and CopulaCPTS generate larger regions, which is expected given their less flexible hyper- rectangular shape. PCP tends to obtain the largest region sizes as it includes samples from low-density areas, whereas STQDR and HD-PCP mitigate this by removing samples from low-density areas, resulting in more compact regions. Fi- nally, Figure [3 in Appendix G provides results for the mean
and multivariate Gaussian mixture model (Bishop. 1994).
We present results for MQF² in the main text, while similar results for the other models are provided in Appendix G We compare the methods using several metrics, including conditional coverage (WSC, CEC-X, and CEC-V), marginal coverage (MC), region size, and computational time. A detailed description of the experimental setup is provided in Appendix
Conditional coverage. Figure presents the results for all datasets, ordered by increasing dataset size. On most datasets, C-PCP, L-CP, and C-HDR obtain the best condi- tional coverage. contrast, PCP, STDQR, PCP, and DR- CP are the least conditionally calibrated. Finally, M-CP
M-CP
C-HDR
STDQR
CopulaCPTS
PCP
C-PCP
DR-CP
HD-PCP
L-CP
10
9
8
6
5
2
102
PCP STDQR CopulaCPTS M-CP HD-PCP
C-HDR
DR-CP C-PCP L-CP
Time
10a
Rank of Median Region Size
households scm20d
rf2
r
scmid
meps 21
20
house
bio blog. data calcofi
blog_data
taxi
meps
meps
Figure 5: CD diagrams with the base predictor MQF based on 10 runs per dataset and method.
Figure 6: Total time in seconds for calibration and test.
8

---

Multi-Oulput Conformal Regression
region size, where DR-CP consistently performs best, under the oracle setting, it minimizes the expected region size, explained in Section 5.2
structure, modeling the joint density, leveraging gener- ative models. These approaches produce expressive prediction regions but are generally more computationally demanding.
Computation time. F shows the total computation time for each method. M-CP and CopulaCPTS have the shortest computation times, as they do not require learning a complex model for the output joint distribution. L-CP and DR-CP follow, benefiting from the absence of per-instance sampling. In contrast, sampling-bar -based methods typically require 100 200 times more computation time.
While conformal prediction (CP) always guarantees marginal coverage, conformity scores whose thresholds do not vary instance-wist fail achieve the desirable property of asymptotic conditional coverage (ACC). In contrast, our proposed scores enable ACC but require estimating the con- ditional distribution of the conformity score-an an inherently challenging task in low-data regimes. Similarly, CP methods based on generative models introduce additional sampling variability. Finally, our large-scale empirical study system- atically compares these conformal methods across multiple multi-output regression datasets, using various evaluation metries, including conditional coverage and prediction re- gion volume.
Application to image dataset. Finally, we extend our com- parison to a regression problem where the output is an im- age, which has a higher dimensionality than the previously considered tabular datasets. Specifically, we use the CIFAR- 10 dataset (Krizhevsky et al.. 2014), consisting of 32x32 RGB images, each labeled with one of 10 possible classes. We train the base predictor Glow (Kingma and Dhariwal, 2018), conditioned on the image label, where the output space is y = [0, 1|3x22x32 (d- = 3072) and the input space is X' = (0.....) 9) (p = 1). The results, detailed in Ap- pendix J lead to similar conclusions regarding conditional coverage, region size, and computational time.
Future work will focus on replacing region prediction with a recalibrated multivariate distribution, equipped with an explicit density function and conformal coverage guaran- tees. We also plan to extend our approach to more complex outputs, including semi-structured and unstructured data (e.g., images, text, and graphs), and to further investigate the theoretical connection between multivariate calibration and conformal prediction
7. Conclusion
We studied the problem of constructing conformal predic- tion regions for multi-output regression, which remains rel- atively underexplored compared to the univariate case. We presented a unified comparative study of several conformal methods along with their associated conformity scores, high- lighting their properties and interconnections. In addition, introduced two new classes of conformity scores: CDF- based scores, including variant compatible with generative models, and latent- based scores, which exploit invertible generative models for improved computational efficiency. Both classes generalize existing conformity scores from the univariate setting.
The choice of conformity score directly influences the ge- ometry and ftexibility of the resulting prediction regions. In the univariate setting, the most flexible regions are typically unions of intervals. In contrast, the multivariate case allows for wider variety of geometries, ranging from hyperrectan- gles and ellipsoids to highly flexible, nonconvex regions that can be disconnected and capture distributional bimodality. simple and computationally efficient approach to con- struct separate univariate prediction regions for each output dimension and apply a correction for joint coverage. How- ever, these methods do not capture dependencies between output dimensions and typically result in rigid, (unions of) hyperrectangular regions with limited flexibility. In con- trast, more flexible methods account for correlations and dependencies across outputs by incorporating the covariance
9

---

Multi-Output Conformal Regression
Impact statement
[12] Victor Chemoohukov, Kaspar Wüthrich, and Yinchu Zhu. "Distributional conformal prediction". Proceed- ings of the National Academy of Sciences of the United States of America 118 (48 2021).
This work contributes to the development of statistically reliable and interpretable machine learning algorithms. En- hancing trust and transparency in predictive modeling helps designing more practical and accessible models for real- world applications.
[13] Janez Demšar. "Statistical Comparisons of Classi- fiers over Multiple Data Sets". Journal of machine learning research: JMLR (2006), pp. 1-30.
[14] Victor Dheur et al. "Distribution- Free Conformal Joint Prediction Regions for Neural Marked Tempo- ral Point Processes". Machine Learning 1 113 (2024), pp. 7055-7102. arXiv: 2401 04612 [C]
References
[1] Brandon Amos, Lei Xu, and J Zico Kolter. "In- put Convex Neural Networks". Proceedings of the 34th International Conference on Machine Learning.
[15] Jacopo Diquigiovanni, Matteo Fontana, and Simone Vantini. "A Conformal approach for functional data prediction". Book of short papers SIS 2021, Pisa, 21-25 Giugno 2021. 2021, pp. 907-910.
Vol. 70. Proceedings of Machine Learning Research. PMLR, 2017, pp. 146-155.
[2] Anastasios Angelopoulos and Stephen Bates. "Con- formal prediction: A gentle introduction". Founda- tions and Trends@ in Machine Learning 16 (4 2023),
[16] Jacopo Diquigiovanni, Matteo Fontana, and Simone Vantini. Distribution Prediction Bands for Mul- tivariate Functional Time Series: an Application to the Italian Gas Market. 1xXiv:2107.00527 [stat]. 2024.
pp. 494-591.
[3] Anastasios N Angelopoulos et al. "Uncertainty Sets for Image Classifiers using Conformal Prediction". International Conference on Learning Representa-
[17] Jacopo Diquigiovanni, Matteo Fontana, and Simone Vantini. "The importance of being a band: Finite- sample exact distribution prediction sets for func- tional data". arXiv [stat.ME] (2021). arXiv: 2102
tions. 2020.
[4] Rina Foygel Barber et al. "The limits of distribution- free conditional predictive inference". arXiv [math.ST] (2019). arXiv: 11603.0664 04684 [math.ST]]
06746 MM
[5] Eustasio del Barrio, Alberto González Sanz, and Marc Hallin. "Nonparametric multiple-output center- outward quantile regression". Journal of the Ameri- can Statistical Association (2024), pp. 1-43.
[18] Eshant English et al. "JANET: Joint Adaptive predictioN-region Estimation for time-serics". arXiv [stat.ML] (2024). arXiv: 2407 06390 [stat.ML] MM]
[19] Zhenhan Fang. Aixin Tan, and Jian Huang. "CON- TRA: Conformal prediction region via normalizing flow transformation". The Thirteenth International Conference on Learning Representations. 2025.
[6] Alessio Benavoli, Giorgio Corani, and Francesca Mangili. "Should We Really Use Post-Hoc Tests Based on Hean-Ranks?? Journal of machine learn- ing research: JMLR 17 (5 2016). pp. 1-10.
[20] Shai Feldman, Stephen Bates, and Yaniv Romano. "Calibrated Multiple-Output Quantile Regression with Representation Learning".. Journal of machine learning research: JMLR 24 (24 2023), pp. 1-48.
[7] Christopher M Bishop. Mixture density networks. Tech. rep. Birmingham, 1994.
[8] Annika Camehl, Dennis Fok, and Kathrin Gruber. "On superlevel sets of conditional densities and multi- variate quantile regression". Journal of Econometrics (105807 2024), p. 105807.
[21] Milton Friedman. "A Comparison of Alternative Tests of Significance for the Problem of m Rank- ings". The Annals of Mathematical Statistics 11 (1 1940), pp. 86-92.
[9] Guillaume Carlier, Victor Chermozhukov, and Alfred Galichon. "Vector quantile regression: An optimal approachh. Annals of statistics 44 (3 2016). pp. 1165-1192.
[22] Léo Grinsztajn, Edouard Oyallon, and G Varoquaux. "Why do tree-based models still outperform deep learning on typical tabular data?" Neural Information Processing Systems 35 (2022), pp. 507-520.
[10] Maxime Cauchois, Suyash Gupta, and John C Duchi. "Knowing what you know: valid and validated con- fidence sets in multiclass and multilabel prediction".
[23] David Ha, Andrew M Dai. and Quoe V Le. "Hyper- Networks". International Conference on Learning Representations. 2022.
Journal of machine learning research: JMLR 22 (1
2021), pp. 3681-3722.
[24] Marc Hallin and Miroslav Šiman. "Multiple-output quantile regression". Handbook of quantile regres- sion (2017), pp. 185-207.
[11] Domagoj Cevid et al. "Distributional Random Forests: Heterogeneity Adjustment and Multivariate Distributional Regression".. Journal of machine learn- ing research: JMLR 23 (333 2022), pp. 1-79.
10

---

Multi-Oulput Conformal Regression
[25] Marc Hallin al. "Distribution and quantile fune- tions, ranks and signs in dimension d: A measure transportation approach". The Annals of Statistics 49 (2 2021), pp. 1139-1165.
[38] Soundouss Messoudi, Sébastien Destercke, and Syl- vain Rousseau. "Copula-based conformal prediction for multi-target regression". Pattern Recognition 120
(2021), p. 108101.
[26] Sture Holm. "A Simple Sequentially Rejective Multi- ple Test Procedure". Scandinavian journal of statis- tics, theory and applications 6 (2 1979), pp. 65-70.
[39] Soundouss Messoudi, Sébastien Destercke, and Syl- vain Rousseau. "Ellipsoidal conformal inference for Multi- Marlit Regression". Conformal and Probabilis- tic Prediction with Applications. Conformal and Prob- abilistic Prediction with Applications. PMLR, 2022,
[27] Eliahu Horwitz and Yedid Hoshen. "Conffusion: Confidence Intervals for Diffusion Models". arXiv
[cs.CV] (2022). arXiv: 2211 09795 [cs.cV]]
pp. 294-306.
[28] Chin- Wei Huang et al. "Convex Potential Flows: Uni- versal Probability Distributions with Optimal Trans- port and Convex Optimization" (2020).
[40] Bric Nalisnick al. "Do Deep Generative Models Know What They Don't 't Know?" International Con- ference on Learning Representations. 2018.
[29] Rob J Hyndman. "Computing and Graphing Highest Density Regions". The American statistician 50 (2 1996), pp. 120-126.
[41] Davy Paindaveine and Miroslav Šiman. "On direc- tional multiple-output quantile regression". Journal of multivariate analysis 102 (2 2011), pp. 193-212
[30] Rafael Izbicki, Gilson Shimizu, and Rafael Stern. "Flexible distribution-free conditional predictive bands using density estimators", Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics. Vol. 108. Proceedings of Machine Learning Research. PMLR, 2020, pp. 3068- 3077.
[42] Harris Papadopoulos, Alex Gammerman, and Volodya Vovk. "Normalized nonconformity measures for regression Conformal Prediction". Proceedings of the 26th IASTED International Conference on Arti- ficial Intelligence and Applications (Innsbruck, Aus- tria). AIA '08. USA: ACTA Press, 2008, pp. 64-69.
[43] Harris Papadopoulos et al. "Inductive Confidence Ma- chines for Regression", Machine Learning: ECML 2002. Springer Berlin Heidelberg, 2002, pp- 345- 356.
[31] Rafael Izbicki, Gilson Shimizu, and Rafael B Stern. "CD-split and HPD-split: Efficient Conformal Re- gions in High Dimensions". Journal of machine learning research: JMLR 23 (87 2022). pp. 1-32.
[44] George Papamakarios, Theo pavlakou, et al. "Masked autoregressive flow for density estimation". Advances in neural information processing systems (2017).
[32] Chancellor Johnstone and Eugene Ndiaye. "Exact and approximate conformal inference for multi- output regression". arXiv [stat.ML] (2022). arXiv: 2210.17405 [sL
[45] Ji Won Park, Robert Tibshirani, and Kyunghyun Cho. Semiparametric conformal prediction.
[33] Kelvin Kan et al. "Multivariate Quantile Function Forecaster". Proceedings of The 25th International Conference on Artificial Intelligence and Statistics. Vol. 151. Proceedings of Machine Learning Research.
arXiv:241 1.02114. 2024.
[46] Vincent Plassier et al. "Probabilistic Conformal Pre- diction with Approximate Conditional Validity". The Thirteenth International Conference on Learning Representations. 2025.
PMLR, 2022, pp. 10603-10621.
[34] Diederik P Kingma and Prafulla Dhariwal. "Glow: Generative flow with invertible Ixl convolu- tions". Neural Information Processing Systems abs/1807.03039 (2018).
[47] Alvin Rajkomar et al. "Scalable and accurate deep learning with electronic health records". .nj digital medicine 1 (1 2018), p. 18.
[48] Yaniv Romano, Evan Patterson, and Emmanuel Can- des. "Conformalized quantile regression". Advances in neural information processing systems (2019).
[35] Alex Krizhevsky, Vinod Nair, Geoffrey Hinton, et al. "The CIFAR-10 dataset". online: http://www. cs. toronto. edu/kriz/cifar. himl 55 (5 2014), p. 2.
[49]
Yaniv Romano, Mattco Sesia, and E Candès. "Clas- sification with valid and adaptive coverage". Neural Information Processing Systems 33 (2020), pp. 3581- 3591.
[36] Jing Lei and Larry Wasserman. "Distribution-free free prediction bands for non narammetet regression". Journal of the Royal Statistical Society. Series B, Sta- tistical methodology 76 (1 2014), pp. 71-96.
[50]
Raphael Rossellini, R Barber, and R Willett. "In- tegrating uncertainty awareness into conformalized Quantile Regression". International Conference on Artificial Intelligence and Statistics 238 (2024). pp. 1548.
[37] Soundouss Mcssoudi, Sébastien Destercke, and Syl- vain Rousseau. Cropula-based conformal prediction for multi-target regression". Pattern Recognition 120
(2021), p. 108101.

---

Multi-Output Conformal Regression
[51] Mauricio Sadinle, Jing Lei, and Larry Wasserman. "Least ambiguous set-valued classifiers with bounded error levels". Journal of the American Statistical As-
[64] Xiaofan Zhou et al. "Conformal prediction: A data perspective". ACM computing surveys (2025).
[65]
Yanfei Zhou, Lars Lindemann, and Matteo Sesia. "Conformalized adaptive forecasting of heteroge- neous trajectories". Proceedings of the 41st Inter- national Conference on Machine Learning. ICML'24. Vienna, Austria: JMLR.org. 2024.
sociation 114 (525 2019), pp. 223-234.
[52] Matteo Sesia and Yaniv Romano. "Conformal pre- diction using conditional histograms". Advances in neural information processing systems (2021).
[53] Kamile Stankeviciute, Ahmed M Alaa, and Mihaela van der Schaar. "Conformal Time-series Forecasting". Advances neural information processing systems 34 (2021).
[54]
Vincent Stimper, Bernhard Schölkopf, and Jose Miguel Hermandez-Lobato. "Resampling Base Dis- tributions of Normalizing Flows". International Con- ference on Artificial Intelligence and Statistics. In- ternational Conference Artificial Intelligence and Statistics. PMLR, 2022, pp. 4915-4936.
[55] Sophia Sun and Rose Yu. Copula Conformal Pre- diction for Multi-step Time Series Forecasting. arXiv:2212.03281 [cs, stat]. 2024.
[56] Sophia Huiwen Sun and Rose Yu. "Copula Confor- mal prediction for multi-step time series prediction". The Twelfth International Conference on Learning Representations. 2024.
[57] Jacopo Teneggi et al. "How to Trust Your Diffusion Model: A Convex Optimization Approach to Confor- mal Risk Control" (2023).
[58] Alexander Timans et al. "Max-Rank: Efficient Mul-
tiple Testing for Conformal Prediction". 258 (2025),
pp. 3898-3906.
[59] Grigorios Tsoumakas et al. "MULAN: Java Li- brary for Multi-Label Learning".. Journal of machine learning research: JMLR 12 (71 2011), pp. 2411-
2414.
[60]
Vladimir Vovk, Alexander Gammerman, and Craig Saunders. "Machine-Learning Applications of Al- gorithmic Randomness". Proceedings of the Six- teenth International Conference on Machine Learn- ing. ICML '99. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1999, pp. 444-453.
[61] Jun Wang al. Conformal Temporal Logic Planning using Large Language Models; Knowing When to Do What and When to Ask Help. 1cXiv:2309.1009
[cs]. 2023.
Zhendong Wang al. "Probabilistic Conformal Pre- diction Using Conditional Random Samples". Inter- national Conference on Artificial Intelligence and Statisties. International Conference on Artificial Intel- ligence and Statistics. PMLR, 2023, pp. 8814-8836. [63] Frank Wilcoxon. "Individual Comparisons by Rank- ing Methods". Biometries Bulletin (6 1945), pp. 80- 83.
[62]
12

---

Multi-Output Conformal Regression
A. Related Work
Our research builds on a broad body of literature that spans several closely related themes. This supplementary section provides a concise overview of these topics.
In the realm of multivariate functional data, Diquigiovanni al., 2021a introduced conformal predictors that create adaptive, finite-sample valid prediction bands, with extensions into time series applications, particularly in energy markets (Diquigiovanni et al., 2024). In image processing, recent applications (Horwitz and Hoshen, 2022 Teneggi et al., 2023) apply CP in a pixel-wise manner, resulting in hyperrectangular regions that may not capture pixel dependencies effectively.
For multi-step steli--stp-aaee step-ahead or multi-horizon forecasting, predictions can be made across multiple outputs simultaneously rather than sequentially, aligning with a multi-output forecasting framework. Stankeviciute et al., [L.T] explored multi-horizon time series forecasting using recurrent neural networks (RNNs), incorporating univariate conformal techniques with nominal coverage adjustments via Bonferroni correction. Similarly, English et al., 2024 adapted the Amplitude- Beglitude-Modulateddd L-inf norm method from Diquigiovanni et al., 2021b for multi-output, multi-step forecasting.
In multi-target regression, Messoudi et al., 20216 applied copula functions in deep neural networks to provide multivariate predictions with guaranteed coverage. Their findings suggest that simple parametric copulas can work for certain datasets, but more complex copulas may be required for well-calibrated predictions, which introduces challenges, as complex copulas typically require significant calibration data. Building on this, Sun and Yu, [2024b) proposed a copula-based method for multi-step time series forecasting, optimizing the calibration and efficiency of confidence intervals. However, this method requires two calibration phases and is primarily feasible with large calibration datasets. Moreover, its validity relies on the empirical copula, limiting applicability to other learnable copula classes. One very recent advancement on the subject, following ideas expressed by Messoudi et al., 20021 in their conclusions, is Park et al., 2024 where the dependence structure between marginal distributions is recovered via the use of vine copulas.
Another set of methodologies that tackle multi-output problems are based on multiplicity approaches for multiple testing. Timans al., 2205 improves over Bonferroni correction using permutation tests, and obtain a tighter and globally valid prediction. Methods based multiplicity correction such controlling the Family- Wiie Error Rate (FWER) valuable for providing error control guarantees across multiple outputs. In contrast, the methods we survey and propose aim for potentially tighter prediction regions by directly modeling the multivariate structure.
In the context of conformal prediction, the flexibility in configuring the prediction region is a degree of freedom for the modeler. To overcome the limitations of traditional hyper-rectangular prediction regions, Messoudi et al., 2022 introduced ellipsoidal uncertainty sets that enable instance-specific adaptation of confidence regions. Johnstone and Ndiaye, 2022 advanced multi-output regression by developing efficient techniques for approximating conformal prediction sets without retraining the model, although their approach relies heavily the predictive model being a linear function of Y. Sun and Yu, 2024b constructed ellipsoidal prediction regions for time series, capable of modeling dependencies between outputs, though this method does not handle multimodality. Our work closely connects with the multivariate conformal prediction literature, where multi-horizon prediction is viewed as a prediction across multiple outputs.
Overall, as this study demonstrates, the flexibility of conformal prediction allows for coherent handling of diverse data types. Multi-output problems represent one facet of a broader taxonomy, as explored by Zhou et al., 2025 who discuss further developments in multi-output conformal prediction.
B. Additional multi-output conformal methods
In this section, we describe the prediction regions ft for M-CP and CopulaCPTS, which both produce hyperrectangular regions.
M-CP. Zhou et al., 2024 applied univariate conformal method to each output i E [d] of the multivariate response. Specifically, given a conformity score 8; for the i-th dimension, joint coverage across all dimensions can be achieved using the following conformity score:
SM-CP(r, y) = -mxx 8 (2., y:)). iE(d]
(15)
A similar score has been considered by Diquigiovanni et al. (202b)) in the context of functional regression.
In this work, we use Conformalized Quantile Regression (CQR, Romano et al., 2019) for each output i E [d], where the
13

---

Multi-Output Conformal Regression
conformity score is given by:
s, M4) max max(4(()) V/s -ū,(x))
(16)
with (x) and (a) representing the lower and upper conditional quantiles of Y = x: at levels aj and au. respectively. In our experiments, we consider equal-tailed prediction intervals, where or = 2. Og - 2 and denotes the miscoverage
level. The corresponding prediction region RM-CP()) = X [[(2) - 9. (x) + q] is hyperrectangle.
CoPUlaCPTS. CopulaCPTS (Sun and Yu, 2024a) is originally designed for time-series but the calibration procedure is valid for any multi-dimensional variable. models the joint probability of uncertainty for each output with a copula. The calibration set divided into two subsets 1 and Deal-2. D serves for estimation of a CDF on the conformity score of each output and D allows calibrate the copula. CopulaCPTS can combine any univariate multivariate conformity scores. In this paper, we use the CQR score 9, (16) for each dimension i [d].
Denote F. the empirical CDF of the conformity scores (5, 1;) (a.y) for i E [d]. and p:11 the corresponding empirical quantile function. In practice, to minimize region sizes while achieving marginal validity, CopulaCPTS computes the optimal si....,5 that minimize the following loss using stochastic gradient descent:
[( 8d)
I (Dat) i=1
Σ
((((),)))<<F F. -(1 -(1-a) -aa) a))
(17)
2eal-2l
Then, the prediction region is defined as:
RCQquaaCPT(()) (y Ey Wi (1,,)( < 8
(18)
Sun and Yu, 2024a proved that this prediction set achieves marginal coverage. However, since CopulaCPTS does not follow the SCP algorithm, it does not achieve properties on the marginal coverage from Appendix E.T
C. Relationship between conformity scores and regions
Section 2.2 and Section 3 in the main text presented several multi-output conformal methods through their conformity score
As explained in Section 2.1 their corresponding prediction region it can be computed as follows:
R(x)) {y E y s(x, v) s g).
In this section, we explicitly derive the prediction region associated with these methods.
M-CP. Following Zhou et al., 2024 the prediction region ÅMEr can be derived from sycr as follows:
SMCr v) s + max (z, Mi) S
(19)
rE(d]
Vi E [d], (1. yi) E
(20)
Vi E [d], [d)(m)x(((()) Vs, Vi 0,(x)) 2
(21)
Vi E [d),4(x) - S ys 0,(x) <
(22)
Wi E [4)[[<x) - - y Ar S (x) +
(23)
Vi E [d]. E [((xz) - (.0,(z) +
(24)
E ((((x) - (x) +
(25)
i=1
E RMC
(26)
DR-CP
The equivalence straightforward.
14

---

Multi-Output Conformal Regression
C-HDRR. Given y fy and U /(8 X =), for any y E V write
SC-HDR(F, y)
(27)
HPD j(1 [x)
(28)
P(f(( | 2 fly 1 X-
(29)
IP(U 2 fl (x) X = xX=
(30)
-PU - fly x)
(31)
F (f(y x)),
(32)
where Xur the conditional CDF of given X = x.
Recall that the prediction region for CDR is given by
RC-HDR (x) = ty E y f(v x) 2 tq}
where = sup(t: PP x 2
(33)
The threshold 4g
33 can equivalently written as follows:
Lg sup(t: P(Í() x) 2 x
(34)
sup(t: P(U 2 t x)
P(U
(35)
= -uuptt 10:::11PPC
(36)
= sup((l 1- Foix=r(())
(37)
[1 i))
(38)
where we use the definition of the upper quantile function in the last step.
Using
33 and 38
we can write
SC-HDR(2, v) S + HPD ((1 x) S
(39)
1 FuXx=o (f(y y) x)) S q
(40)
FuIx=r (f x)) 2 1
(41)
fly 2) 2 Fu - -)
(42)
fly 2t
(43)
y RC-HDR(*).
(44)
PCP. Let ((,,r) represent ball with center and radius r. Following Wang et al., 2023b Rpcr(z) corresponds to union of balls:
show that, for any X,
spCP(r, 5 - min by 8-00 1|2 59
(45)
(E[T]
31 E [L], l1 80 و
(46)
E1 E [L] y E (()
(47)
y E U B(D .0)
(48)
[E(A]
y E Rece(x),
(49)
where (1)) Frlx=z,le e [L].
HD-PCP. For HD-PCP, the reasoning is similar to PCP with the difference that only the [[1 a) samples with the highest density are kept.
15

---

Multi-Output Conformal Regression
CDF-based conformity scores. We note that the region (c) has similar form (x) (v E y sw (w(2))) S $))
except that the threshold sw (I. y) is different and depends In fact, write
ReDF(w) = {y E E sCnF v) S 4}
(50)
{y E y: (W w)) - d}
(51)
{v E y sw v) - F
(52)
In the special case where sw SpCP, since PCP always generates predictions as a union of balls, we can conclude that C-PCP will do the same.
Latent-based conformity scores. Since ((;x)) is bijective, for every E V, there exists a unique E Z such that y ((2;)) Therefore, the condition dz d2 x)y; I q is equivalent to (z( S q. where Q (1)1z).). This gives the prediction region:
RLcr(z) ): (v: x)) < (Q(2;)) : and dz(2) sg).
(53)
(54)
D. Additional illustrative examples
D.I. A real-world application
Following Wang et al., 2023a we apply the multi-output conformal methods to the taxi dataset, where the goal is to predict the drop-off location of a New York taxi passenger based on the passenger's information.
Figures (a) and S(a) displays five randomly selected samples from the dataset, showing the pick-up (red pin) and drop-off
(blue pin) locations of taxi passengers. The remaining panels show a specific input-output pair (ar, v) and the corresponding prediction regions generated by the conformal methods discussed in this paper. The coverage level for these regions set to 0.8, with MQF² as the base predictor, introduced in Section F.2
Figure 7 corresponss to the same data and predictions regions Figure except that it is zoomed for better comparison with Figure 9
Each region is labeled with its size, calculated using the estimator from Section F44| displayed in the bottom left corner. Notably, C-PCP generates regions similar in shape to PCP but with an input-adaptive radius, resulting in smaller region sizes (8.2 compared to 8.67) in this case. Additionally, HD-PCP produces more clustered regions, while PCP and C-PCP show more dispersed regions.
16

---

Multi-Output Conformal Regression
8 8
(a) Sample Data
(b) M-CP
(c) CopulaCPTS
(d) DR-CP
(c) C-HDD
(I) PCP
(g) HD-PCP
(h) STDQR
(i) C-PCP
ii) L-CP
Figure 7: Conformal methods applied on the NYC Taxi dataset for an input with low uncertainty.
2.82
2.18
2.24
0.81
(a) Sample Data
(b) M-CP
(c) CopulaCPTS
(d) DR-CP
(c) C-HDR
4.12
2.41
2.42
1.48
1.05
PCP
(g) HD-PCP
(h) STDQR
(i) C-PCP
(i) L-CP
Figure 8: Conformal methods applied the NYC Taxi dataset for input with low uncertainty.
Figure presents the example for an onput-output pair where the input is associated with higher uncertainty, resulting in larger region sizes. As in the first figure, the shapes of the regions (e.g., unions of hyperrectangles, quantile regions, etc.) remain consistent but expand to cover a larger area. Conformal methods with the best region sizes differ between the two figures, with C-HDR producing the smallest region in the first figure and DR-CP in the second. In this case, C-PCP selects a larger radius than PCP, resulting in larger regions than PCP. The observation that PCP and C-PCP produce more dispersed regions, while HD-PCP generates more clustered regions, also holds true for this higher uncertainty case.
17

---

Multi-Oulput Conformal Regression
58.92
41.49
2.32
36.00
(a) Sample Data
(b) M-CP
(c) CopulaCPTS
(d) DR-CP
C-HDR
14.12
13.96
14.32
45.01
40.18
(f) PCP
(g) HD-PCP
(h) STDQR
(i) C-PCP
()) L-CP
Figure Conformal methods applied on the NYC Taxi dataset for an input with high uncertainty.
D.2. Toy examples
We define two data-generating processes to evaluate the models compared to a known distribution: unimodal heteroscedastic distribution and a bimodal heteroscedastic distribution. The input variable X unidimensional (p 1) and the output variable Y E 2 is bidimensional (d 2). The variables X and Y are scaled linearly such that the mean and variances each dimension are 0 and 1. The figures are inspired by Barrio et al., 2024
Unimodal heteroscedastic process. The first process is illustrated in Figure 2 in the main text. The data generating
process is as follows:
X ~ U(0,1),
(55)
Y X = (((13-z))001 0²
(56)
(57)
where k 200, 0.2, Tn is the 2 identity, and, for I....,k.
(1) COS Oj
p1
(58)
(1) (0.5 sin a,)
M1
(59)
aj (j - 1) π
(60)
1
Detailed metrics for this dataset, supporting Section 5.1 are provided in Table 2
18

---

Multi-Output Conformal Regression
MC
Median Size
CEC-X
CEC-Z
WSC
Test time
Method
(x100)
(x100)
M-CP
0.805
8.470.14
0.110
0.08120.018
0.803
0.0336v.022
CopulaCPT'
0.815 (.011
8.78
0.191 0,049
0.163
4g.15
1.040.021
DR-CP
0.808
7.03p.071
0.613u.045
0.560
0. 0.710uu116
0.0209
C-HDR
0.81 O0.0038
p.80p
0.06370.016
0.0825p.o12
0.798
3.52p.085
PCP
0.805
9.169.089
0.6688.052
0.5870.046
0.7130.0080
1.690.021
HD-PCP
0.804..0037
7.440.056
0.287
0.2566.034
0.758
3.38p.043
STDQR
0.806
7.870.070
0.3433.025
0.305
0. 0.7466.011
1. 1.770.022
C-PCP
0.808
9.14p.12
0.04644.013
0.0484
0.822,.0085
3.44p.056
L-CP
0.803
8.24o.11
0.05440.0073
0.0654p.012
0.811a.011
0.0217a.0005
Table 2: Detailed metrics for the unimodal heteroscedastic process from Figure 2
Bimodal heteroscedastic process. Figure
T0 similar to Figure 2 but with a bimodal distribution for the output, is
introduced in Section
The data generating process is follows:
X ~ (0..,2))
(61)
X x x~ 0.5 NV4 x1a 0.5 -N La/2).
(62)
Oracle
CopulaCPTS
C-HDR
2
88
2
2
0
11
Yi
Y1
PCP
HD-PCP
STDQR
C-PCP
L-CP
Y1
Y
X
Figure 10: Examples of prediction regions on bivariate bimodal dataset, conditional on a unidimensional input.
E. Proofs
E.1. Distribution of the marginal coverage conditional on calibration data
In contrast to M-CP, L-CP, and DR-CP the methods C-HDR, PCP, PCP, and C- CCPCP rely on a non-deterministic conformity score. For each calibration and test point, C-HDR, PCP, -CP, and C-PCP require sampling K, L, L, and + K points, respectively.
Let Dea ((X(), y represent the calibration dataset and (X,Y) be the test instance. Let Scal (S()) } represent samples from the calibration dataset where s is generated based on input X() and Stest the samples generated based on input X. Despite the added sampling uncertainty, these methods still provide marginal coverage guarantee:
P Px.Y.Suu. (Y Ye E R(X))
(63)
19

---

Multi-Output Conformal Regression
Compared 3 the probability is additionally Seal and Sest- This result, specifically for PCP and HD PCP, was demonstrated by Wang et al., 2023b
In Lemma we further show that the marginal coverage conditional on the calibration dataset Deal and the samples Scal follows a beta distribution, using standard arguments. Assuming no ties among the scores, this lemma applies to any conformity score 8.
Lemma 1. Assuming no ties among the scores and i..d. inputs, outputs and samples, the distribution of the coverage, conditional on the calibration dataset and its samples, is given by:
P(Y E R(X) Deat, Seat) ~ (ka, Deall + - ka).
(64)
where ka --1-(aa))) (Paal 1)7. Moreover, P(Y R(X))
Ru
1
Dea
Proof. For the methods C-HDR, PCP, HD-PCP, and C- CCPC, the conformity score is non-deterministic due to sampling uncertainty. To clarify, we define a deterministic conformity score 8: X' y x S, where represents the space of samples for a given method.
For j = I,...., | let S, = s(x()) Y y denote the conformity score on the calibration dataset, and let = s(X,Y, Sea) represent the conformity score for the test instance. Since & deterministic and the tuples (X(1)) y()) s()))..... )..,)., (x(())) y ( ((PPll)) (X,Y, are i.i.d random variables, S are also i..d random variables.
Since
Sipwal S are indentically distributed, they share the same CDF Using the probability integral transform, (S) UO0 1). Thus, (s())).... Fs(S)pul) correspond to uniform variates U1...., UPDal- Since there are no ties among the scores, Fs is strictly increasing, and (S))) = Ug) for j I,..., [Deal], where S/y) and U() are the j-th order statistics. Hence:
P(Y R(X) Deal. Seal) P(S < 8(k.) St....,Spp.l FF((S(.))
(65)
(66)
U(ka) ~ Beta(kka (Deal + 1 - -a).
(67)
(68)
The final step results from the distribution of uniform order statistics. Taking the expectation of the Beta distribution gives:
ka 2 - a. P(Y E R(X)) = E(P(Y e f((x) Deat, Seat) = DDeal +
(69)
E.2. Proofs of asymptotic conditional coverage
E.2.1. L-CP
Proposition 2.
Assuming Deall - oo and Q(Z; - Y|X, L-CP achieves conditional coverage.
Proof. We first show that the conditional coverage of L-CP equal to the CDF of the random variable dz(2) q. i.e.
20

---

Multi-Output Conformal Regression
(zz (0). Given x E X, have:
P(Y E ( XX
(70)
=P(Y e (Q(z;)): E Rz(()) xx x)
(71)
(P(Q-'(Y:)) E Rz(() x=
(Invertibility of ())
(72)
=P(Z E Rz(1))
(Q(7; YIX)
(73)
=P(dz(Z) s) (2)
(74)
(75)
Marginalizing over X, we obtain that the marginal coverage is also equal to Faz (i):
P(Y R Rxx((X) 1x) xx]
(76)
-Ex (P(( P(Y
(77)
=Ex =Ex(F(x(2)(())
(78)
= (Z))
(79)
In the limit of (Deal + c0, thanks to the Glivenko-Canteli theorem, P(Y E Rlcr(X)) 1- -a and the quantile ( obtained by SCP is thus F-1
dz) (1 - a).
Finally, we obtain that the conditional coverage equal to a:
P(Y E R xx=
(80)
= =Fz((z) (F (z( (1-a))
(8
=1 - a.
(82)
E.2.2. C-HDR AND C-PCP
Lemma 2. Assuming Deal + oc, any conformal method with conformity score SCDF 12 achieves conditional coverage,
independently from the conformity score sw of the base method. With the additional assumption that K - Do and f = f. SECDF 13) achieves conditional coverage.
Proof. Let W = sw (X, Y) and consider x E X' and E y. By the probability integral transform, SCDF(2; Y) FwIx==(W X = x) ~ ((,
Marginalizing over X, we obtain:
IPYY E RCDr(X)) = P(SCDE (X,Y) S Q)
(83)
=Ex[P(scDr (X,Y) s 4 X)]
(84)
Ex[IP(U s 4)
(85)
Ex(g]]
(86)
4.
(87)
where U ~ -0, 1). In the limit of [Deall + oc, thanks to the Glivenko-Cantelli theorem, P(Y e RXDF(X)) = 1 a and the quantile & obtained by SCP is thus 1 - a.
Finally, we note that:
IP E ReDr(X) x=x) = P(SCDE < X
(88)
P(U <1-a) a)
(89)
1 a.
(90)
21

---

Multi-Output Conformal Regression
Assuming f S. observe that, for any E X' and y E V SECDF y) - sCDF(2,)) y) as K + by the law of large numbers. Thus, under these conditions, any conformal method with conformity SECDF achieves conditional coverage.
Proposition 3. Assuming Deail -> oo and K oo, both CDR and CCP with the oracle base predictor f- achieve conditional coverage.
Proof. The proof is direct by Lemma 2 with sw (2, SDR-CP(Z, y) for C-HDR and sw(z,)) = =spCP(E, y) for C- PCP.
E.2.3. M-CP
Consider M-CP with exact quantile estimates (x) = Qy, (ax x)) and (.(x) = Qy (a z) where Qy, (a x)) is the quantile function of Y; conditional X evaluated a. This section introduces two propositions where M- requests two different nominal coverage levels au on, namely -I- a and a. The propositions show that M-CP can achieve conditional coverage under two contrasting scenarios: independence total dependence between the dimensions of the output.
Proposition 4. Assuming YI,..., Ya are conditionally independent given X, M-CP achieves conditional coverage if Oa 11
Proof. For any E X' and E [d], we first establish that the 41 ath quantile of the distribution of sX Yi) given X equals O:
P(s 00 X = z) (()) - ux)) x x
(91)
IP((((X) S Y A Y S (x() X =)
(92)
- ((((()) > VY > u(X) | z)
(93)
P(4(() x =x) - xx(X)) (X XX=
(94)
a.)
(95)
auu aq
(96)
1/1 a.
(97)
Using 97 we show that the 1 ath quantile of the distribution of s(X, Y) given is 0:
P(sMMcp(X, Y) s 0 X x) F(s,(X, X)) 0, vi [d] X =x)
(98)
F(s:(X,Y)) < 0 ... Asa S 20|x= =x)
(99)
P(s:(X, Yi) 0 X = x). ((s,(XXYa) 0 X = =x)
(100)
=V1-aa
(101)
=1-a,
(102)
where 1100p is obtained by conditional independence of Yi.... Ya given X. Marginalizing over X, we obtain that the
1 - othh quantile of s(X, Y) is O:
P(sMcr(X,Y) < = Ex [x( (X, Y) 0 X))
(103)
x - a)
(104)
1 a.
(105)
In the limit of (Deal + oo, thanks to the Glivenko-Cantelli theorem, P(Y e RM)r(X)) = 1 - a and the quantile ş obtained by SCP is thus 0. Finally, using
102 and : = O, we obtain that MPCP achieves conditional coverage:
PPYY e PMMc((X = x>) P(sM.c((XX)) < X = x) =1-a.
(106)
22

---

Multi-Output Conformal Regression
Proposition 5. Assuming X
a.5. Ya M-CP achieves conditional coverage if az ai 1 a.
Proof. Let E X. Using 96), we first show that the athh conditional quantile of the distribution of sx Yi), for any i e [d]. O:
(X, YJ S 0 X = x) = =0 at
(107)
-1-60
(108)
Using 1080, we show that the oth quantile of the distribution of xXx, Y) given X is 0:
P(s(X, =) x) P(s;(XX < 0, e [d] x=
(109)
(X,) ......... s Y.) 50
(110)
(X, 0 X
(111)
= a1
(112)
where 111 is due to X Y which implies that, conditional
x, (x() (X) (X(X) and thus s1 11) = sa(X, Ya). Using (105). we obtain that q = 0, Finally, using 112 we
(x(X) and
obtain that M-CP achieves conditional coverage:
P(Y R(X) Xx= =x) P(s 5O X =x) =
(113)
E.3. Connection between sample-base basee and density-based methods
This section proves the connections between sample-based and censity-bascer methods introduced in Section 5.3 We start by restating a known lemma of conformal prediction.
Lemma 3. Consider a conformal prediction method with conformity score s. If g R - IR is a strictly increasing function, then the method with conformity score g s will produce the same prediction regions.
Proof. For any X, consider the prediction region created with as in Section 2.1
R(x) {y E y ((r, s Quantile (fst)selPPall U [oo); k:n)
(114)
Since g is strictly increasing,
f((x)) = {u E v: g(((2,,)) v)) < 9 (Quantile ((8)) U {xo); ka))} f E y: g(s(z, 8)) < Quantile
(115)
U (oo); ka)
(116)
Since (116 corresponds to the prediction region with conformity score g s, this shows that the two methods create the same regions.
Proposition 1. PCP is equivalent to DR-CP with f f= fmax- Similarly, HD -CPP is equivalent to DR- CP with == where only a] samples with the highest density among (1r() he(r are kept. Finally, C-PCP is equivalent to CDR with max-
Proof. In the following proof, note a to signify that there exists strictly increasing function such that g(b).
Consider DR- CP with f fmax. We have:
SDR SDR crR,C)(
-mmaly
(117)
1 max fs (v: (1)) TE(L]
[ ( max (v; 80//C [E[T]]
(118)
min - fu(y; (1)) rE(A]]
(119)
t min [ly 1/21 rE[Z]
(Ja(y; ())) has spherical level sets centered f-(0)
(120)
spCP(2, v).
(121)
23

---

Multi-Output Conformal Regression
We obtain the equivalence between the methods by Lemma 3 The proof for HD-PCP follows the arguments.
We now consider C-HDR with f =max We have:
SC.HDRR (z. v)
Σ
1mfmax (V-(k) x)) fmax() (y 2)) where y frix ls E [K].
(122)
kE(K]
Developing the inequality for lo E [K], we obtain:
Faxx (8 xx)= 2 fmax(t) (y) x)
(123)
max fs (M)) ) 2 max fs 8 (1) [E]L]
fmaxl (x) max fs(y; 800 [EIL]
(124)
[]
min (Y ())) < min - (y; (()) [E[L]
(125)
|E(4]
min ()) (0 ~ minlly 101
(Js(y; (() has spherical level sets centered at (
(126)
[c(
rE(r]
(127)
Noting that
122 with
127 corresponds to the conformity score of C-PCP, we obtain the equivalence.
F. Experimental setup
This section describes experimental setup in more details. Computations performed based on 2 workstations, with with 2 A6000 GPUs and 64 CPU threads, and one with A5000 GPUs and 64 CPU threads, running for 48 hours.
F.1. Datasets
We consider total of 13 datasets that have been used in previous studies. Since focus is on multivariate prediction regions, we select only datasets with output that is at least two-dimensional. Specifically, we include 6 datasets from Feldman et al., 2023 4 datasets from Tsoumakas et al., 2011 (MULAN benchmark), dataset from Wang et al., 2023b datasets from Barrio et al., 2024] and dataset from Camehl et al., 2024
Each dataset split into training, validation, calibration, and test sets with 2048 points reserved for calibration. The remaining data split into 55% for training, 15% for validation and 30% for testing. The preprocessing follows the setup
described in Grinsztajn al., 2022 Table 3 provides the detailed characteristies of each dataset.
Table 3: Characteristics each dataset considered in this study.
Nb instances
Nb features p
Nb targets
Source
Dataset
Camehl Mulan
houscholds
7207
scm20d
8966
60
16
9005
64
rf2
9005
64
8
sem
9803
279
16
Feldman
meps 21
15656
137
meps_19
15785
137
meps 20
17541
137
house
21613
14
2
bio
45730
blog data
50000
55
Del Barrio
calcofi
50000
Wang
taxi
50000
F.2. Base predictors
We consider multiple base predictors and focus on MQF² for our main experiments (Section
24

---

Multi-Output Conformal Regression
The Multivariate Quantile Function Forecaster (MQF2, Kan al., 2022) is normalizing flow that is directly compatible with most of the methods presented since it is invertible, has an explicit PDF, and can be sampled from. M- CPC CopulaCPTS and STDQR require small adaptations from the original methods, discussed below. The quantile function Q and distribution function Q-1 of MQF exhibit cyclical monotonicity, meaning they are the gradient of a convex function (Hallin et al., 2021).
MQF
The main idea behind MQF is to interpret Convex Potential Flows (Huang et al., 2020) as multivariate (vector) quantile functions, in the sense that the representation propent (128) and cyclical monotonicity property (129) are satisfied (Carlier al., (016))
Y Y= ((;3))
Vo E X where 7 ~~ 2(1, 1)",
(128)
T (Q(1):))-Q(0)) - (222
(21 - 22) 20
Va E 22 Z.
(129)
When d = 1, this reduces to the classical univariate quantile function. In practice, we follow Kan et al., 2022 and use a quantile vector that follows a normal distribution ~ No(,, 7), allowing better training.
The underlying model of MQF is partially input-convex neural network (PINN, Amos al., 2017 with two hidden layers, each containing 30 units. Increasing the number of parameters did not significantly improve performance, which is partly due to the efficiency of Convex Potential Flows compared to other normalizing flows (Huang et al., 2022). While hyperparameter tuning for each dataset could enhance performance, it is not the primary focus of this paper.
MQF is trained using maximum likelihood estimation with early stopping, with a patience of 15 epochs, where validation loss is measured every two epochs.
Distributional Random Forests. Distributional Random Forest (Cevid et al., 2022) a model built upon the Random Forest algorithm, which adaptively identifies the relevant training data points for any given test point. More specifically, given test point X, Distributional Random Forest outputs weight w(x(') x)) for each training point x with (+) E Dminn This approach enables accurate estimation of any quantity of interest conditional on x E X. In our experiments, we estimate the conditional distribution X as a Gaussian mixture, with each component centered on a training point and weighted by the Distributional Random Forest.
The PDF at y E y given x X is expressed as:
Draaal
fly (1)zz)= Σ w(x(1) x)) NNy g ola)).
where a is tuned by minimizing the NLL on a grid search. For the Distributional Random Forest, the minimum node size is set to 15, the forest consists of 2000 trees, and the splitting criterion the maximum mean discrepancy (MMD).
Since this method does not operate in latent space, we do not consider L-CP in combination with this base predictor. CD diagrams for this predictor are presented in Appendix G.2
Multivariate Gaussian Mixture Model parameterized by a hypernetwork. As another base predictor, we consider a multivariate Gaussian Mixture Model parameterized by a hypernetork. The hypernetwork is a multilayer perceptron (MLP) that outputs the parameters of a mixture of M multivariate Gaussian distributions. Given x E X,, for each mixture component m E [M], the hypernetwork outputs the logit zm 2x)(x) (for the categorical distribution over the mixture components), the mean Mm (z) (component location), and the lower triangular Cholesky factor Lm (x) (representing the scale of the covariance matrix).
The mixture weights (x: are obtained by applying the softmax function to the logits 2m (a.), ensuring they sum to 1. The covariance matrices (x) for each component are constructed by taking the product Lm (x) (x)) guarantecing that they are positive semi-definite.
The PDF evaluated in y E y conditional to x X is given by:
M
fly x)) Σ Tm())) N(y (m(xx) Em (x)).
m=m
25

---

Multi-Output Conformal Regression
The model is trained using maximum likelihood estimation with M 10.
Similarly to Distributional Random Forests, this method does not operate in a latent space, and thus do not consider L-CP. CD diagrams for this predictor are presented in Appendix G.3
E3. Adaptation of conformal methods into common framework.
To ensure fair comparison among conformal methods, apply the calibration step using the same base predictors. Only M- CP, CopulaCPTS, and STDQR require slight modifications from their original formulations.
For M-CP and CopulaCPTS, direct estimation of marginal distributions for each output Yi, E [d] is infeasible with MQF Instead, we estimate the lower and upper quantiles by first sampling ( he(z) from fy|x== given x E X, and then
computing the empirical quantiles ((14 and ?. ((141-2))) 2>1), Sampling time is not accounted in time computations for these methods. While a more computationally efficient base predictor could be used, this approach ensures direct comparison with other conformal methods by maintaining consistency in the base predictor.
For STDQR, we modify the original method by replacing the conditional variational autoencoder (CVAE) with normalizing flow. Following recommendations for future work from Feldman al., 2023 exploit the property that the output is normally distributed in the latent space and replace the base predictor by a normalizing flow. This adaptation leverages the assumption that the output is normally distributed in the latent space, allowing for an exact inverse transformation and eliminating a potential source of noise. To construct a region Rz with coverage 1 a in the latent space, we select the 1 a proportion of samples closest to the origin, ensuring correct coverage without the need for directional quantile regression. The calibration procedure remains unchanged.
F44 Metrics
Marginal coverage.
Marginal coverage is measured using
MC
Σ
1 e R(z)).
[Drest
(D
Region size.
We report the mean region size
Mean Size
Σ [R(z)).
Drest
(D.
To avoid large regions disproportionately affecting the result, we also report the median of the region sizes
Median Size Quantie((RRs)]ttteee......
Computing the size of the region is challenging in high dimensions. Hence, we propose an unbiased estimator of the region size using importance sampling:
[((x) = J, 1(y R(x))dy Bxns(a)
1(L E R(xx)) f(Y x))
K
(((())) f((x)) f(Y xx)
(130)
K
where (4) ~ fy X=r. & E [K]. This estimator compatible with all base predictors in Appendix F.2 since it both possible to sample from their predictive distribution and evaluate the PDF. Appendix F7 estimator.
we discuss the efficiency of this
Conditional coverage. To ensure robust evaluation of conditional coverage, consider three different conditional coverage metrics, detailed in Appendix [.6] The Worst Slab Coverage (WSC, Cauchois et al., (021) groups inputs into "slabs" and evaluates the worst obtained coverage. The coverage eиror conditional to X (CEC-X) partitions the input space X' and evaluates coverage on each subset. The coverage error conditional to = f(L X), where Ý ~ fyix, (CEC-V, Izbicki et al., 2022) Dheur et al., 2024), creates a partition based on the distribution of V, which is more robust to high-dimensional inputs.
26

---

Multi-Output Conformal Regression
Computing time. We report the total time required for calibration and testing the marginal coverage. Specifically, this requires evaluating conformity scores Deal followed by evaluating conformity scores Diea-
F.5. Multi-Model, Multi-Datasse Comparison
In order determine whether there are significant differences in model performance, we first apply the Friedman test (Friedman, 1940). Following the recommendations of Benavoli al. (2016), then conduct pairwise post-hoc analysis using the Wilcoxon signed-rank test (Wilcoxon, 1945). coupled with Holm's m alpha correction (Holm, 1979) to adjust for multiple comparisons.
The results are visualized using critical difference (CD) diagrams (Demšar, 2006). In these diagrams, models are ranked, with a lower rank (positioned further to the right) indicating better performance. A thick horizontal line connects models whose performances are not statistically different at the 0.05 significance level.
For MC and WSC, the CD diagrams report |MC- (1 a)] and |WSC- (1 ai], both of which should be minimized.
F.6. Metries of Conditional Coverage
Worst Slab Coverage. Introduced Cauchois al., 2021 the Worst Slab Coverage (WSC) metric quantifies the minimal coverage over all possible slabs in R.d, where each slab contains least fraction 5 of the total mass, with 1. For given vector u Rd, the WSC associated with denoted as WSC,, is defined by:
WSC, inf asch {Pom (u R(xz)) a a Suta b) s.t. PD (a S vTx, S b) 2)
(131)
where a, b, E IR. This metric assesses conditional coverage by focusing on inputs that lie within a slab defined by using the inner product vTx, to measure similarity.
To estimate the worst-case slab, we follow the method from Cauchois et al., 2021 uniformly sampling 1,000 vectors 1/ from the unit sphere gd-1 and calculating:
WSC = min WSC,,-
(132)
ey ES
To mitigate overfitting on the test dataset, partition the test set into two subsets, Diess ((2)) U as in Romano et al., 2020 Sesia and Romano, 2021 We identify the worst combination of a, b, and on ( by minimizing the WSC metric with 8 = 0.2, and then evaluate conditional coverage on the separate subset D test
CEC-X. CEC-X approximates conditional coverage by partitioning the input space X E X' C RP. We apply the k- +eans++ clustering algorithm the inputs (+) in the validation dataset D creating partition A = A U...UU U A, over X. The Coverage Error Conditional to X is defined as:
Deol
CEC-X = 1
DDesu PPD E (c(()) (
200 E Ai) - (1-a)) (1-)))(11)) -(1-a))²
(133)
CEC-V. CEC-V similar to CEC-X, but the conditioning is on the distribution of log V = log f(Y X), where T~frie Unlike CEC-X, CEC-V is more robust to high- high-dimensional inputs. This approach originates from the CD-split* method (Izbicki et al., 2022) and has been adapted to multivariate outputs in Dheur et al., 2024
In practice given an input x, a new feature U, is created. First, samples v, from V | X x are generated by sampling .i...., .m ~ fyxx== and evaluating U; = f(ys z). The resulting vector Us = (o(1). m(m) consists of the order statistics ((4) from 01,..., Uma
27

---

Multi-Output Conformal Regression
1000
11"
180
142
341
31"
132
14250
112
122
300
14219
K
K
K
Figure 11: Panels to 4: Trajectories of the log volume estimator with increasing K compared to the true log volume (dashed line) for different output dimensions d. Panel 5: Log volume estimator with K 100 compared the true log volume (dashed line).
The k-means++ clustering algorithm is applied the vectors log ux() in the validation dataset D and partition Av A U U... U.A. A, over R= is obtained. The Coverage Error Conditional to the distribution of is then computed according to 113 using the partition Av.
Dheur et al.. 2024 notes that the distance function corresponding to this partitioning approach is the 2- Wasseessteinn distance with respect to the distribution of V.
E7. Estimator for the region size
In this section, we discuss the efficiency of the region size estimator introduced in Appendix F44 This estimator is based on a density estimator fyxx=z and a sample (*),, E [K], drawn i.i.d. from the conditional distribution
X x for any
X. Specifically, the estimator is given by:
1()(()) R(x)) ((x) = K f(P()) xx)
1
K
While the estimator is unbiased. i... [[V( [R(=)) we want study its variance, Let 1( R(x)) represent the indicator that sample Y lies within the prediction region R(x), and let p IP(Y R(z)) denote the coverage probability obtained from the samples based our density estimator. Using the law of total variance, obtain the following expression for the variance of (x)):
(2))
V
K
x)
T 111
I
E
V
1
V
E
K
f(L
x))
f(Y xx)
1
pV
+ + - P)L
K
f(Y
x))
f(C x))
Assuming that the density estimate corresponds to the true density, i.e. frix=z fr x)) and that R achieves conditional coverage, then p = 1 a, and we obtain:
1 V((()) K (11-000 (1-o)V
(1
+
a
(Y x))
(Y xx)
This indicates that the variance of our estimator only depends on the variance and expectation of the random variable a(Y[a) In this case, the variance does not directly depend on the output dimension d.
Figure II shows how the estimator behaves in a scenario with a specific density estimator and prediction region with varying output dimension and an 80% coverage level. Since there is no dependence on X, we abbreviate the notation as follows:
28

---

Multi-Output Conformal Regression
R R(x), (1) f(u x). and V V(x) for any E X. The density estimator is standard normal distribution
{ E y llvila S Fxa - -a) a)) where x2 is the chi-squared
f(v) Ny: Ia) and the prediction region is ball R =
distribution with d degrees of freedom and x3
its quantile function. can be shown that Py Py~J() (& R) = 1 - a. In
this case, the volume V of R can be computed exactly.
Each of the first four panels in Figure shows five trajectories for log V as K increases from 1 to 100. The true volume, log V. of the prediction region is indicated by a dashed line. We observe that the estimator converges within a reasonable range of the true volume for varying output dimensions d. The last panel illustrates the value of log V function of with log V again marked by dashed line. From this, we observe that the estimator remains close to the true volume across different output dimensions d.
G. Additional results
This section presents additional results for MQF² (Appendix G.1 Distributional Random Forests (Appendix G.2 and the Multivariate Gaussian Mixture Model (Appendix G.3). The experimental setup is described Appendix F
G.1. MQF
Figure presents the marginal coverage and median region size across datasets of increasing size for MQF In Panel 1, all methods except CopulaCPTS attain precise marginal coverage. This is expected since these methods follow the SCP algorithm (Section 2.1 and their marginal coverage conditional on the calibration dataset and samples from the calibration dataset follows Beta distribution whose parameters only depend on the size of the calibration dataset (Appendix E.T While CopulaCPTS attains marginal coverage, the larger variance in its marginal coverage arises because it does not follow the SCP algorithm.
In Panel 2, the median region size is normalized between and for each dataset in order to facilitate comparison. We observe that C- HDR often obtains the smallest median region size. The performance of the other methods can vary highly across datasets for the median region size and is better visualized a CD diagram (see Figure]
M-CP
CopulaCPTS
DR-CP
C-HDR
PCP
HD-PCP
STDQR
C-PCP
L-CP
Marginal coverage
0.90
Median Region Size
1.00
0.85
0.75
0.80
0.50
0.75
0.25
0.70
0.00
households scm2od
n2
m
scmla
meps 21
meps 19
meps
house
bio blog data calcoft
taxi
househalds scm20d
f12
ril
scmic
meps 21
meps 19
mops 20
house
bia blog data calcofi taxi
Figure 12: Marginal coverage and median region size with the base predictor MQF across datasets sorted by size.
Figure 13 presents critical difference diagrams for three conditional coverage metrics (CEC-X, CEC-Z and WSC), the mean region size and median region size, the total calibration and test time. Results are consistent with the results from the main text.
29

---

Multi-Output Conformal Regression
M-CP
CopulaCPTS
DR-CP
C-HDR
PCP
HD-PCP
STDQR
C-PCP
L-CP
houscholds
14.2a48
12.3087
13.2cas
10.6g3s
20.5
17.8041
15.5
18.6
sem20d
67.6
1.12e+02, teall
2.33c+02;
42,0rg
1.05e402 tew01
44,494
09.4
26,0xz
72.01
r2
100054700s
0.0055500033
0.00215unne
0.0006900.002
0.00700g 0036
0.0061 7o
0.006244xxet
0.00262ammz
0.00104,a1e4s
0.0054770w07
0.00555uxes
0.00215u0000
0.000
1700000
4.00624u.x01
0.00262e 0012
0.001044,0004
scmld
0.5 0.528
0.323
0.867amm
0.239%ums
0.68440052
0.67 Ince
0.21 6uuu4
0.197
meps,
0.171
0.2277.013
0.1322u144
0.359um21
0.246
0.283uo15
0.2200mm
0.244oa52
meps__1
19
0.21 4.u22
0.595p42
0.175mm1
0.119hmms
0,3966
0.2666u013
0. 0.307ac4s
0,238un25
0.232au4s
meps__2
0.362acsw
0.223
0.114umz
0.535
0.436
0.472aes2
0,34
0.280komt
house
1.17aems
1.22
0.6644uzt
0.651lum
0,882
0.680k.018
0. 0.799
0.858our
1.19
hin
0. 50006
0.296xcco22
0.257
0.218
0,34300076
0. 0.259
D. 0.209oco077
0.3022
0.2677.001
bloy
0.0948amus
0.03740m266
0.0155emea
0.141ou023
0.125cun2
0.1633
0.1066m2
0.06766017
culcofi
2. Bucos
2. 2.38
1.670023
1.99 1.99
2.33us0
1. 1.89u
1. 1.97xm21
2.8 loooe
2.70
taxi
4,260
4.720.11
2.624a22
2.62
4.03g040
18
3. 3.630uss
4 4.02004
4.94u.12
Table 4: Median region size with the base predictor MQF?
Dataset
M-CP
CopulaCPTS
DR-CP
C-HDR
PCP
HD-PCP
STDQR
C.PCP
I.CP
bouseholds
36.9usa
3
15.7am
40.11a
33.8g,
30.1a1
28.81
62.61.
50.61,4
scm20d
7.03e+062.e000
3.82e+07
5.30e+092.ei0t
1.61e=04
1.56e+0441e138
1. 1.59c+04s
1.37e+09
2,20 2.20e=10p1ee00 10
1.86e+02.eww Dev02
8304022
15.1a9
1.29c+022 teves
1.09e+0243es01
4.05c+0521ee.6
9,85e+0243:.0z
4.01c+02
rf2
0.86e+02, Uswoe
1.83e4020
15.199
3.51e+02
1.29e+02z
1.09++02atea)
4.05c+052
9.87e+024s=-0
4.02c+02)
somld
2. 2.37e405s moees
1.81c405 lowos
78.4; Am-i1
2.73e+08; towin
$7.31
43.01
43.61
1.8c=081sH00
1.52c
meps_21 21
1.21asss
1.17asss
0.315
1.44a44
0.617axz
0.558
0.553
2.07a.1
1.96
meps__1 19
1.14
1. 1.1mm
0.2933..18
1.29
0.58lomm
0.55966m1
0.5377aam1
1.960m7
1.52
meps 20
1.20hees
1.17am
0.30%
1.30pos7
0.6066uZ2
0.562
0.546po1
2.01a.1
1.59
house
18870227
18loms
0.8870.63
1.09gass
1.230004
0.9646630
1. 140.000
1.40000
1.710013
1.40
1.42a39
0.269%
0.486um7
0.390gm14
0.2976x000
0.31
1.4laa
2.57
calcofi
2.04uss
2.22am:
1.42unss
1.950mm
2.22a051
1.70
1.78amss
2.830015
2.34p04
blog_data
0. 0.390amt
0.390ac16
0.0852
0.375
0.17300997
0. 163
0. 18800007
0.613ures
0,520mmr
taxi
5.68
6.35
2.67a1477
3.21aom
4.55
3.54ums
3.99mm
5.360m0
Table 5: Mean region size with the base predictor MQF?
9
8 7
6
PCP DR-CP
C-PCP
PCP
L-CP
PCP
C-PCP
L-CP C-HDR M-CP
DR-CP
C-PCP
DR-CP
C-HDR
STDQR
STDQR HD-PCP CopulaCPTS
C-HDR M-CP
STDQR
CP
HD-PCP
CopulaCPTS
CopulaCPTS
M.CP
Rank of CEC-X (x100)
Rank of CEC-Z (x100)
Rank of WSC
8
6
5 4 3 2
1
5 3 2
C-PCP CopulaCPTS M-CP C-HDR
DR.CP HD-PCP
PCP
C-HDR DR-CP
STDQR
STDQR
CopulaCPTS
C-PCP
PCP
M-CP
L-CP
L-CP
HD-PCP
Rank of Mean Region Size
Rank of Median Region Size
Figure 13: CD diagrams with the base predictor MQF² with 10 runs per dataset and method.
G.2. Distributional Random Forests
Figure 14 presents additional results for the base predictor Distributional Random Forests. Since this model does not rely on latent space, results for STDQR and L-CP are not included.
In terms of conditional coverage, the results align with those of MQF2, with C-PCP and C-HDR outperforming DR-CP PCP, and HD - PCP. Notably, M- CP achieves competitive conditional coverage, suggesting it pairs well with DDF-KDE. Similar to MQF2, all methods except for CopulaCPTS attain precise marginal coverage.
The median region size is normalized to a [0,1] range for each dataset to facilitate comparison. We observe that C-HDR generally achieves the smallest median region size, followed by DR- CP. The test time is the lowest for and CopulaCPTS
30

---

Multi-Output Conformal Regression
while C-PCP and HDR obtain the highest computation times.
CopUlaCPTS
C-HDR
HD-PCP
0.8
CEC-X (x110)
10"
(EC-Z((x10
10a
WSC
10-1
10-1
0.6
0.5
0.90
102
Marginal coverage
Median Rogi Sizee
0.8
Time (s)
107
0.80
0.4
0.75
10")
0.70
households scm20d
Ift
scmid
meps 21
meps 19
maps 20
Bio house
blog. calcofi
taxi
households scm2od
ff
m scmid mops 21 meps 19
meps house
bio blog data calcoft taxi
households scmzod
ff2
r
maps 19
meps 20 house
bio blog. data calcofi taxi
som maps
Figure 14: Conditional coverage metrics with the base predictor Distributional Random Forests across datasets sorted by size.
Figure 15 shows CD diagrams obtained with Distributional Random Forests as the base predictor. The results are consistent with Figure [4
5 4
3
2
DR-CP PCP HD-PCP CopulaCPTS
M-CP C-PCP C-HDR
DR-CP
C-PCP M-CP
DR-CP PCP
M-CP
PCP
C-PCP
HD-PCP
CopulaCPTS
HD-PCP C-HDR
CopulaCPTS
C-HDR
Rank of CEC-X (x100)
Rank of CEC-Z (x100)
Rank of WSC
7
6
5
4
3
2
7
6
5
4
3
2
1
5
4
3
2
1
CopulaCPTS M-CP C-PCP
DR-CP C-HDR HD-PCP
PCP C-PCP HD-PCP M-CP
C-HDR DR-CP CopulaCPTS
C-HDR C-PCP
CopulaCPTS DR-CP
HD-PCP
PCP
Rank Mean Region Size
Rank of Median Region Size
Rank of Time (s)
Figure 15: CD diagrams with the base predictor Distributional Random Forests based 10 runs per dataset and method.
G.3. Multivariate Gaussian Mixture Model
[6 additional results for the base predictor Multivariate Gaussian Mixture Model. Similarly to Distributional Random Forests, this model does not rely latent space and thus results for STDQR and L- are not included.
The conditional coverage also aligns with MQF², C-PCP and C-HDR outperforming DR RR-P. PCP, and HD PCP. M- HP and CopulaCPTS achieving intermediate conditional coverage. As expected, marginal coverage is precise for all methods except CopulaCPTS.
C-HDR often obtains the smallest median region size, while DR- consistently attains the best mean region size.
31

---

Multi-Output Conformal Regression
M-CP
CopulaCPTS
DR-CP
C-HDR
-
PCP
HD-PCP
C-PCP
0.9
10a
CEG-Z (x100)
10"
WSC
0.7
CEC-
10
10-1
0.6
10
0.90
Margina coverage
Median Region Size
0.85
0.6
Time (s)
0.80
10-)
0.75
0.70
0.0
households sem2od
riz
scmia
maps 21
meps 19
meips 20
house
bio blog data cak of
taxt
households som20d
semza maps 21
meps 19 meps house
bio blog cak
taxi
households scrzzod
i
meps 21 maps maps
house
bio blog data calcofi taxt
Figure 16: Conditional coverage metrics with the multivariate Gaussian mixture model base predictor across datasets sorted by size.
CD diagrams in Figure 17 are consistent with Figure 16
2
7
6
C-PCP C-HDR M-CP
C-PCP C-HDR M-CP
PCP
C-PCP C-HDR
DR-CP
DR-CP HD-PCP CopulaCPTS
DR-CP
HD-PCP
HD-PCP M-CP
CopulaCPTS
CopulaCPTS
Rank CEC-X (x100)
Rank CEC-Z (x100)
Rank of WSC
7
5
7
6
4
3
3
2
1
C-PCP CopulaCPTS
DR-CP
PCP
C-HDR DR-CP
C-PCP C-HDR HD-PCP
M-CP
HD-PCP
CopulaCPTS
PCP
DR-CP
C-HDR
C-PCP
CopulaCPTS
Rank of Mean Region Size
Rank of Median Region Size
Rank of Time (s)
Figure 17: CD diagrams based on Multivariate Gaussian Mixture Model parameterized by hypernetwork with M 10 and runs per dataset and method.
G.4. Impact of the number of samples
Figures 18 and 19 illustrate how conditional coverage, marginal coverage and region size change as a function ofK on all datasets. For a better comparison among datasets, the metrics CEC-X, CEC-Z, the median region size and the mean region size are normalized between 0 and with results averaged over 10 runs. Furthermore, the red line indicates linear regression fit, allowing to see the trend.
Conditional coverage metrics decreasing with K indicate that conditional coverage tends to improve with increasing number of samples. This is expected since an increasing number of Monte-Carlo samples allows a better estimation of the CDF of the scores in (13). Marginal validity is obtained with any K. However, small sizes K will lead to more duplicated conformity scores and thus a possibility of overcoverage. Median region sizes and mean region sizes also tend to decrease with K as the CDF approximation improves.
32

---

Multi-Oulput Conformal Regression
househalds
r11
scmld
meps_ 19
house
calcofi
taxi
Scm20d
f2
meps_21
meps_20 20
bio
blog_data
Relative CEC-X (x100)
2.0
Relative CEC-Z (x100)
2.0
0.8
0.6
0.5
10
30
100
300
10
30
100
300
10
30
100
300
K
Relative Median Region Size
K
Relative Mean Region Size
K
1.0
1.0
1.0
Marginal coverage
0.8
0.8
0.6
0.4
0.5
02
100
100
300
K
K
K
Figure 18: Evolution of conditional coverage, marginal coverage and region sizes of C-PCP as a function of the number of samples using the base predictor MQF?. The metrics CEC-X, and CEC-Z should be minimized, while the marginal coverage and WSC should approach a (indicated by the dashed black line). The red line, obtained by linear regression, indicates the general trend.
households
rf1
scm1d
meps_19
house
calcofi
scm20d
r22
meps 21
meps 20
blog data
Relative CEC-X (x100)
1.0
Relative CEC-Z (x100)
1.0
1.0
0.8
0.8
0.7
0.4
0.2
02
300
Relative Median Region Size
Relative Mean Region Size
0.9
0.8
0.6
0.7
Marginal
02
02
0.4
0.2
Figure 19: Reproduction of Figure 18 for C-HDR.
H. Comparison with Bonferroni correction
To better understand the prediction regions produced by FWER control methods, we provide a qualitative and quantitative comparison with Bonferroni correction. We consider Bonferroni correction applied to the scores CQR (see (16). similarly to M-CP. Figure 20 provides an illustrative example, and Table 6 provides results on this same dataset. This shows that Bonferroni computationally fast produces larger regions due to the rectangular shape.
33

---

Multi-Output Conformal Regression
Oracle
Bonferroni
M-CP
2
2
x
2
2
2220
2-
Y5
X
Y5
Y1
DR-CP
C-HDR
L-CP
2
x
-2
3 2 2
3 2-100-
2-22
X
Y1
X
Y1
Y1
Figure 20: Prediction regions for a bivariate unimodal dataset, conditional on a unidimensional input. The black, green, and yellow contours represent regions with nominal coverage levels of 20%, 40%, and 80%, respectively. The figure similar Figure 2 in the main text, with Bonferroni added as comparison. Both Bonferroni and M-CP are based on Conformal Quantile Regression (CQR) applied separately for each dimension.
Table 6: Detailed metrics for the unimodal heteroscedastic process from Figure 20 with 1 fixed to 0.8.
MC
Median Size
CEC-X
CEC-Z (x100)
WSC
Test time
Method
(x100)
Bonferroni
0.81 1813pom 08.006
9.07a.15
0.0241
0.0249
0.815
0.00339; 5.9c-05
M-CP
9.0037
8.62a.074
0.0240g.0031
0.01570.0031
0.7966.012
0,0959h.osa
DR-CP
0.7969.0019
6.83p.042
0.432p.019
0.403u.015
0.697 6.970.009
0.0557 7.00075
C-HDR
0.809
6.970.099
0.0129h.005s
0.0155g.0037
0.8 p..000
14.2n.
L-CP
0.798p.0024
8.06g.035
0.00586p,0009
0.00549
0.7940.0039
0.05840.0012
Comparison between C-PCP and CP2-PCP
In this section, we compare our proposed method, C-PCP, with the CP2-PCP method recently proposed by Plassier et al., 2025 More generally, we also compare the methods from the Cp² framework of Plassier et al., 2025 with our class of CDF- based conformity scores (Section 3.1 in the main text). Appendix [.1 we present the more general CP² framework using our own notation for clarity, with of CP2 and show the asymptotic equivalence with CDF-based methods. In Appendix 1.3 discuss the relationship between CDF-based and ca2-b methods.
particular case of CP² In Appendix I.2 we discuss the asymptotic properties
I.I. The CP2 framework
Let us define family of non-decreasing nested regions (R(;; HER such that NTER R(x; t) = 0 and ULER R(x;t) V. and Neses R(x; t') tx; Without loss of generality, these nested regions are expressed in terms of conformity score sw (x, 8) IR as follows:
R(x; t) {v E y sw (x. s t).
(134)
where sw (z, y) is continuous in y-
As the next step, we introduce family of transformation functions f=(A): : - IR parameterized by E R. It assumed
34

---

Multi-Output Conformal Regression
that for any 7, the function À (A) is increasing and bijective. Let ş E R be constant (e.g.. 4p 1). We also define the function (p(T) (p) and assume that - (T( is increasing and bijective.
As first step towards defining CP², we construct prediction region assuming knowledge of the conditional distribution Fyjx. For given imput E X, the prediction region is defined as:
Repa (z) -
(135)
where
= inf : ((X X = x) 21-00
(136)
implies that Repa(1) (x)) guarantees conditional coverage given Furthermore, using 134 and defining the random variable
W = sw sw,, (X,Y), Y), we can equivalently express (366 as
= inf {T : IP
S (+(() x = 2 - a} inf {T PP (o (s S
(137)
X = 2 aa)
(138)
Q (1
(139)
9. (QwIx=z a1))
(140)
where we used that De is increasing and bijective, with 96 (f+(z)) = In other words, Tx is the 1 - a quantile of 96 - (W).
However, in practice, Tュ cannot be computed directly because the true conditional distribution Fvia is unknown. Instead, can be estimated using sample {r(*), k E [K], drawn from the estimated conditional distribution Fyjx=a- If @fQ a) the 1 a quantile of the empirical distribution X EkE(K] 6aw (woy we can compute
f. = 9 . (1-
(141)
should be noted that this estimated prediction region loses the exact conditional and marginal coverage properties due the reliance the estimated conditional distribution. The following shows how conformal prediction restore some coverage properties.
From
134
using
135
we can write
Repa (x) = ly EV: (w S (p)
(142)
= {uv 5- (sw(z, v)) s )
(143)
where we used the invertibility of for any E IR.
Based on (43). Plassier et al., definee defined the following conformity score:
(x y) = (sw (s,((x)))) (x)
(144)
for which the corresponding prediction region Ray is given by
Rep= = ty E ) Scp= (I, w) S a).
(145)
where we used 2 from the main text.
As an example, taking (A()) = TÀ and & = 1, the conformity score becomes:
Scpz (Ix y) = sw (a, f)/Fa>
(146)
where 7z is defined in 141). Finally, we obtain CP2-PCH simply by replacing sw with spcr in (166).
35

---

Multi-Output Conformal Regression
1.2. Asymptotic properties
1.2.1. ASYMPTOTIC EQUIVALENCE PREDICTION REGIONS
In the following. we prove that the prediction regions generated by CP² (for any fr and () and CDF-based methods are identical in the oracle setting, asymptotically, |Deal - Specifically, for any X, both methods select the same threshold l Qw
X x) for the prediction region R(x; l1-a), which ensures coverage level of 1 - Proposition 6. Provided that the assumptions in Appendix I.1 hold, for any x E X, the prediction regions (x)) (for any choice of Jr and sp) and (C)F(x) are equivalent.
Proof. Using the fact that 9p (f+())) = T for any E IR and that By is increasing and bijective, we can write:
(x) = {v E y sw sw v) S (p))
(147)
(u E > 06 (sw (r. S Tx)
(148)
{y E V 94 (sw (I, v)) S 9 (Qw1x=((--a))
(149)
= [yy E y sw (a, y) S Qwix=x(1-o)
(150)
Let RcDF(x) denote the prediction region obtained using the conformity score SCDF [Deal] + As shown in Section SCDr (Y) ~1(0, 1), which implies q = 1 a. Therefore:
3.1
ReDF a)) = {v y (CDF(E,)) 1) S 11-a)
(151)
{y y Fw x=r(Sw (I, 4)) a}
(152)
{y V (w, Qwix= -a)
(153)
This shows that Repa (x) = (2) and that the threshold 1-aa = Qwix= (1 - a) is identical for both methods.
1.2.2. ASYMPTOTIC CONDITIONAL COVERAGE
Proposition 7. Provided that the assumptions in Section 5.2 the main text hold, specifically that Fylx=z = Fyl for all X. and |Deal oo, CP² achieves ACC K + oo.
Proof. Under these assumptions, have QWix=r QwIX=z which implies + To for all E X. Hence, the prediction region for CP² is given by:
(x) = tu E:: 8 (2, S sp).
Since this prediction region provides conditional coverage, it also ensures marginal coverage:
P(Y E Rxpx(X)) = -P(scpa (X, Y) p)
(154)
= =Ex [P(sp (X,Y) s x
(155)
Ex [1
(156)
=1aa. a.
(157)
Since q is the 1 - quantile of scp² ( Y). and as (Deall + oo, we have ş = y by definition. Therefore, since Repa (x)
achieves conditional coverage (see 136 the region Repa (x) also achieves ACC:
P(Y Re (X) X = x) P's X, X)Y)= < q X x) P( (X,Y) $ X = x))
(158)
(159)
> a
(160)
36

---

Multi-Output Conformal Regression
I.3. Relationship between CDF-based and CP2-based methods
natural question is whether there exists {fr)re and yp IR (with the assumptions introduced in Appendix I.1 such that CDF-based and CP2-based methods produce the same regions. In the simple case where the distribution of the base conformity score is in a location family, Proposition and Φ. However, the proof is not easily generalizable to a location-scale family. Further development of existing classes of conformal methods with ACC and their intersections is a promising avenue for future research. Interestingly, we discuss below that answering this question would also draw links between established univariate conformal methods.
g shows that the two methods equivalent for simple choice of fr
Analogy to univariate conformal prediction. To further clarify the distinction between CDF- and CP9-based methods, we draw analogy to the established univariate methods Dist-split (DS, Izbicki et al., 2020) and Conformalized Quantile Regression (CQR, Romano et al., 2019). Since CDF- and CP2-based methods calibrate one quantile intead of interval, we only consider the right-tail version of DS and CQR:
SECDF is analogous to DS but operates in the space of conformity instead of the output space V. DS uses the estimated conditional CDF of the output variable, Sps $ps(1,))) y) (y) transforming y based on its rank.
Scp= with difference adjustment analogous to CQR, and also operates in the space of conformity instead of the output space y. Note that Cp² with difference adjustment can be simplified to scpe(x, y) sw (a, v) Qw (1-a). Similarly,
CQR uses a score based on the difference from single estimated conditional quantile, SCQR (z, y) = y 1r aa).
jk== where (k) ~ PyJx== The difference lies in the way they transform sw (x. y) to obtain ACC. Recall that the conformity scores SECDF and scpa are given by
Both CDF-based and Cp2-based methods rely on sample (F(())
(ECDE v) = Σ (ssw(z; ((M)) I (w(2,)) v)) = (o)x
1
(161)
K
kE[K]
Scp (x, v) = f.. (sw (2) where F (Qw]x=((110)))
(162)
It is known that two conformal methods produce equal regions if and only if their conformity scores are equal after applying a strictly increasing function Φ R + i.e.:
SECDr y) = (x, Vor E X, N.yE E.
(163)
Given x E X, when K is finite, the conformity score SECDF(F, :) is discontinuous and is thus necessarily different from the conformity score Sxp=(x, :), which continuous. A more interesting setting is the case where K + 3o and SECDF( )) becomes continuous. We define the random variable W = sw (X, Y), with Y ~ Frix. Let Fwr and Q denote the conditional CDF and QF of W given X = r. The conformity scores are defined as follows:
BECDF (x. v) = (w (1,
(164)
ccpa (z. v) = ff_ (z, where î₂ = 966 (0
(165)
Thus, we require that
1= (ow n)) = + -FF (ow(a))))))) Vx E x
(166)
or equivalently
f: = ( o(Fw)=(u)
X
(167)
In Proposition we show that, in the particular where the conditional distributions { ( Izex belong to location family, there exists a simple choice of (f-)reR, ip E IR and strietly increasing & IR R such that the two methods are equivalent.
Proposition 8, Consider a scenario where all conditional distributions {Fw boex belong to a location family, i.c.,
Fi F(w - (.)) and (a) = F- (ax) + Pz:
(168)
for some continuous and strictly increasing base CDF F and location parameter ftg. The conformity scores SECDF and &cp" lead to the same prediction regions.
37

---

Multi-Output Conformal Regression
Proof. We will show that there a family of transformations (GTTER: y and strictly increasing Φ : IR R with the assumptions above such that, for any X and It,
(w))) = OlFi (w x))
(169)
Define the transformation function f, as:
==A)=
(170)
where t > 0, and define = a and ((A) =
The inverse transformations are:
- = -F(A-T),
(171)
and
94 (w) = w F-
(172)
Now, for X, compute
f = F- 1 a) + Ar - FF1((()) = Az.
(173)
Finally, we obtain the required equality
F f7. (w) = F(w - -x) =F(w (x)) Fwi((w) (w).
(174)
1.4. Empirical comparison
We perform a direct empirical comparison between CDF-based methods (C-PCP, C-HDR) and the corresponding Cp2 methods (CP2.PCP, CP2-HPD using both linear (-L) and difference (-D) adjustments from Plassier et al., 2025). Figure 21 shows that:
C-PCP performs comparably to CP3PCCPLL (best CP2 variant for PCP).
C-HDD performs comparably to CP-HPD-D (best CP2 variant for HPD).
Other CP2 variants (CPLPCP-D, CP-HPD--L) are generally outperformed by their CDF-based version.
38

---

Multi-Output Conformal Regression
C-PCP
CP2-PCP-D
CP1-PCP-L
C-HDR
.. CP*-HPD-D
CPI-HPD-L
0.84
Marginal coverage
0.82
CEC-X (x100)
100
0.80
10-1
0.78
0.76
10-2
0.9
0.8
CEC-2 (x100)
10
0.7
WS
10-1
0.6
0.5
10-2
households scm20d
rf2
rf1
scmld
meps 21
meps 19
meps 20
house
bio blog. data calcofi taxi
households scm20d
rf2
rfl
scmld
meps 21
meps 19
meps 20
house
bio blog data calcofi taxi
Figure 21: Comparison of CDF-base d methods and CP2 based methods.
The worsened conditional coverage of CP².PCP-D is an interesting observation that was not observed in the smaller scale study of Plassier et al., 2025 In the case of CP2-HPD-L, the poor conditional coverage is due to an incompatibility of the linear adjustment function with the (log-scaled) conformity score $DR-CP L y) = - -lo f(y (z), which can present negative values and thus a decreasing (instead of increasing) adjustment function (A) =λ.
This shows that our simpler SECDE formulation achieves the same practical benefits as scpa without the sensibility of choosing an adjustment function Sr.
J. Results on an image dataset
To better understand the behavior of prediction regions in high-dimensional spaces, we apply conformal methods to the CIFAR-10 dataset (Krizhevsky et al., 2014), which consists 32x32 RGB images, each labeled with one possible classes. We train a generative model conditioned on the image label, where y = [0, 1]3x32x32
3072) represents the image space, and X' (O..... ... , 9) (p- 1) represents the labels. The training, calibration, and test datasets contain 50,000, 1, 1,500, and 1 1,500 images, respectively. As noted in Angelopoulos and Bates, 2023 this calibration dataset size is sufficient ensure good marginal coverage.
Our generative model is a conditional Glow model (Kingma and Dhariwal, 2018 based on the implementation from Stimper et al., 2022 using 3-level multi-scale architecture with 32 blocks per level. Like MQF² (Appendix (22). this generative model is a normalizing flow and directly compatible with methods presented, except M- CP.P. For a direct comparison with M-CP, we compute quantiles based on samples from the generative model as in Appendix F2
The latent space of the conditional Glow model, due its multi-scale architecture, consists of three subspaces: 2 21 Zg x Zz, where 2, 4 Z2 122x8x8 and Z 16x16x16 As the distance function dz in the latent space, we use the maximum norm across the three spaces penalize high norms in any of them: dz(z) = max [1] l|, || =0]|2)>
where 2 = 21 X 22 X
Table 7 presents the metrics introduced in Appendix F44 All methods achieve marginal coverage despite the high
dimensionality y, which expected the marginal coverage distribution conditional on the calibration dataset
39

---

Multi-Output Conformal Regression
independent of d (Appendix E.1 The median of the logarithm of the region size is the smallest for C- and DR- CP. which matches results tabular datasets. The mean region size is not reported because it becomes infinity with machine precision. Instead, we report the mean of the logarithm of the region size, with similar conclusions to the median region size.
Regarding conditional coverage, as in other experiments, L-CP, CDHDR, C-PCP, and M- CP exhibit the smallest CEC-X and CEC-Z values, indicating superior conditional coverage. The WSC metric supports similar conclusions, with DR-CE CP and PCP being the least calibrated.
Table 7: Results obtained with a conditional Glow model on CIFAR-10 with 1 - = 0.9.
MC
Median Log Size
Mean Log Size
CEC-X (x (x100)
CEC-Z (x (x100)
WSC
Time
Dataset
Method
cifar10
M-CP
0.900
-7.100+033..
-7.05e+03
0.11
0.201
0.855n.012
0,465
DR-CP
0.9033.0ma
-8.30c+03 14e001
-8.33c+03 4e+00)
0.1522.030
0.325oo33
0.861
47.01,7e+01
C-HDR
0.902a.004)
-8.33e+03
-8.40c++0119e-00
0.0533
0.0629%..00
0.903u.000
4538,7e+01
PCP
0.899%.0035
-7.11c+113
1-44311
-7.06c=403
0.342o.070
O. 1.155.0000
0.825p.0098
20333.5e+01
HD-PCP
0,899
7.112++0348
-7.06e+031.2emm
0.359
0. 198
0.81 0.8199.0088
4066.9e+01
STDQR
0.898a.0046
-7.11c+0346
-7.06c+031,aem
0.357
0.205oom3
0.828p.01
204; 5e+01
C-PCP
0.900g.0036
-7.08c+0349
-7.04c+031. te++01
0.1188.021
0.0877a.023
0.880p.0067
408 9e+01
L-CP
0.900%.c033
-7.19c+13700
-7.15c+03
0.06688.0086
0. 1.19
0.877
47. 4..41,8c001
K. Full results
Tables and show the full results obtained with the setup described Section 6 Each metric is the mean over 10 independent runs. The standard emror of the mean is indicated as an index. For each dataset and metric, bold values indicate results statistically similar to the best performer (a = 0.05) according to a Z-test.
40

---

Multi-Output Conformal Regression
Table 8: Full results obtained with the setup described in Section
(Part 1),
MC
Median Size
CRC-X (x 100)
CEC-2 (x 100)
WSC
Test time
Dataset
Method
households
0.801
14,2
0.340uus
0.364uuz
0.779ummo
5.6%
CopulaCPTS
0.782o.0s4
12.3us7
0.524kn07
0.745ums
8.86
DR-CP
0.8024.0045
13.2
0.987a.10
1.88g.44
0.656mmis
0.225ux002
C-HDRR
0.807
10.60.3
0.209
0. 0.149ua20
0.795
1c.50
PCP
0.798,,aoss
20.5u
1.07omes
2. 2.35u15
0.632ams
5.48o.t
HD-PCP
0.800
15.6n
0.776une
1.38a.10
0.707am:4
5.76
STDQR
0.804
17.8u.41
0.91 8ourr
1. 1.970.98
0.677amu
8,45c.79
C-PCP
0.8033,0006
15.5u
0.179hu4s
0.120kun6
0.800umss1
11.2g.95
L-CP
0.800
18.6
0.204musn
0.118km
0.788
0.10
scm20d
M-CP
0.8006.0039
67. 67.6g.5
0.0682aa11
0.91
0.777a0000
8.08c.17
CopulaCPTS
0.8330:00s
1.12e+022.trr1
0.22
0.878uu4
0.802amm13
10.9621
DR-CP
0.799u.0oss
2.33e+02h
0.429pus4
2.72a.6
0.69 0.69Iums
0.560pots
C-HDR
0.806
42.0%.
0.159
0.102a017
0.796u006
9.42o.u8
PCP
0.798a,
1.05e+02
1.58
5 28
0.621nms6
6.276.19
HD-PCP
0.799..004
94.4, 40.4
0.504p.47
4. 4.7822
0.67
7.1 lc.42
STDQR
0.801
99.4110e01
0.540kumz
4. 4.86.17
0.620um
8.15
C-PCP
0.809
26.0s.
0.105mmm
0.0896mum4
0.789
14.3a.4s
L-CP
0.7966.0035
72.0 nc=00
0.1 0.166mm33
0.0873ams
0.786
0.0987ame9
f2
M-CP
0.797a.004
0.005470:007
0.2022n00
0. 0.96Mu.15
0.667am1
20.72.4
CopalaCPTS
0.785o.10
0.00555omms
0.299uus
1.18
0.635mm9
29.6u
DR-CP
0.799
0.949uu21
5.420.20
0.549nms
0,344o015
C-HDR
0.801
0.000690,.0a0
0.11
0.230hu47
0.732auus
21.52.
PCP
0.801
0.00700.ords
0.863
5. 5.95u.s
0.538
17.8
HD-PCP
0.800
0.0061 7oo000
0.776
5,58
0.563pme
18.30.4
STDQR
0.800
0.006246,001
0.788g.9
5. 5.67a.50
0.566mm25
25.74.1
C-PCP
0.8024
0.00262o.002
0.0925muts
0.169
0.732
38.64s
L-CP
0.800
0.00104g.0wsx
0.107
0.236
0.730mmm
0.0960huRR1
if
M-CP
0.797%.00s
0.005470xe7
0.202onz0
0.968u
0.667mms
20.91.
CopulaCPTS
0.785x.00
0.00555o.0r3
0.299
1.18.19
0.635pm9
29.841
DR-CP
0,799..02s
0.0021 5ec0000
0.949p.21
5,42a.30
0,5499..s
0.335o016
C-HDR
0.801
0.11
0.230
0.732
21.72.5
PCP
0.801e.a
0.00700..oss
0.863
5.95p
0.538pmm
17.9t4
HD-PCP
0.800
0.0061
70.0000
0.776n.19
5.580.40
0,563am29
18.42.4
STDQR
0.800
0.00624g.001
0.78%g11
5.670.50
0.566muz5
25.74.1
C-PCP
0.8024.0e11
0.00262p.0nz
0.0925mm
0.169hum7
0.732
38.849
L-CP
0.800
0
4.00104aa0mu
0.107am22
0.236m042
0.730
0.0976um67
semid
M-CP
0.796
0.528
L.2
2. 2.42nms
0.636mmr
85.6
CopulaCPTS
0.732o.011
0.323
1.720.00
3,49027
0.582amr7
47.82
DR-CP
0.7933.003
0.867aor7
1.500p07
5. 5.170.20
0.55%
0.584g0025
C-HDR
0.239
0.452ama
0.114mms
0.761ume
87.0g.ev01
PCP
0.795a.0054
0.698
1.770.1
8.1
0.516mm
5,530.44
HD-PCP
0.795
0.6844g002
1.75o.u
7 7.96222
0.530mm17
6.410.
STDQR
0.795
0.67
1.78a13
8.07
0.502
6.876.15
C-PCP
0.8033.0033
0.216
0.456mmm
0.154mmx
0,751um4
L-CP
0.799
0.197
0,463us59
0. 0.108
0.731
0,1022.005
me
M-CP
0.800
0.185mo13
0.926une
0.775mme
0.701amo
1.35e+02 nosat
CopulaCPTS
0.778o.004
0.171a014
0.957a.1
0. 0.693a099
0.684ao11
1.6le+02).9eii
CR-CP
0.803
0.227o013
3.75g.6
4.38
0.531umz
0.228oo11
C-HDR
0.807a.0045
0.132am24
0.437anms
0.260h
0.745amr
1.35e+02,
nosat
PCP
0.801a.031
0.359
3.170.13
3.750.44
0.550pm078
78.4
HDD-CP
0.8024.0024
0.2466.015
2.096.1
2. 2.18
0.601
78.74.3
STDQR
0.802a.022
0.283
2.60n12
2.97a
0.582amr
96.0ga
C-PCP
0.805
0 0.220
0.165
0.085
0.775mums
3e14
LCP
0.801
0.244.002
0.770u13
0.422
0.685uuu
0.1
meps_19
M-CP
0.803
0,214
0.702
0.622nus
0.70 0.709uuss
1.4e+02g
CopulaCPTS
0.804e.a22
0.595p4z
1.13026
0.926u29
0.72
1.77c+0225eitt
DR-CP
0.795omms
0.175mm1
3.91n1k
3.98071
0.501am
0.224mm11
C-HDR
0.807a.029
0.119
0,380mms
0.245pm9
0.753umm
1.4c+022
PCP
0.7946.003
0.3960.09
2.956.23
3.51a.3
0.542
99.9011
HD-PCP
0. 0.796
0,266mm
1.98g.14
2.05u33
0.583uuss0
1.00c+02, Tostt
STDQR
0.79 1o.9m2
0.307am4
2.63nax
2.95g.49
0.557mm:4
1.17++02, Roant
C-PCP
0.810
0 2.2
0.128
0.0757
0.797
2.4444023.emmt
L-CP
0.803
0.232
0.679
0.41.5m1
0.702muu
0.123
41

---

Multi-Output Conformal Regression
Table 9: Full results obtained with the setup described in Section
(Part 2).
Median Size
CEC-X (x0)0)
CEC-2 (x 100)
WSC
time
Dataset
Method
meps_20
M-CP
0.806
0.37 1oom
0.868
0.455
0.702ama
1. 1.93e+02;
1Ree+022.eem
lev01
CopulaCPTS
0.794k.mm
0.362o.0w
0,963
1,497
0.692ams
2.30c+02; 40401
DR-CI
0.8055mmm
0,223
3.52
2.79um6
0.530
0.227aow
C-HDR
0.805x.044
0.1 0.114g,82
0.439
0.122ams
0.745uni0
1.94e+022 Ic401
PCP
0.801amm
0,5356.04
2.83
2.42px7
0.544un0ss
1.20e+002 A0401
HD-PCP
0.804
0.4366
1.89
1.34
0.622amna
1.20e+02, 4e-07
STDQR
0 0.803a.048
0.47
2.45
1884,5
0.575pmm
1.40c+02)
now01
C-PCP
0. 0.806
0.34 Loow
0.186
0.0484
0.792auss
3. 1e130002 1o401
L-CP
0.799
0, 0,280p.aza
0.662au073
0.259hu81
0. 0.703pm14
0.127a.002
house
M-CP
0.801
1.170023
0.254g.023
0.1 0.190hu19
0.730kace
56.11
CopulaCPTS
0.8 22.0002
1.22
0,316
0.276mu27
0.750
60.73.6e00
DR-CP
0.801
0.664
0.895
1.200urs
0.627amm
0.283
C-HDR
0.807
0.651
0,388o.0s
0.1 0.114 1140013
0. 0.709
56.63.8e:01
PCP
0.801a.00
0.882
0.753
1. 1.14ms
0.643uamm
17.6a.
HD-PCP
0 0,003
0.680
0.694
0.789uus
0.649uucee
18.0g.9
STDQR
0.801am02
0.799
0.670
0. 0.788mms
0.649mm0m7
19.5ams
C-PCP
0. o000
0.858
0,275am0
0.0831am1
0.7 729
73.7
L-CP
0.802a.035
1.190003
0.174
0.0542au09
0.75 0.756
0.146
bio
M-CP
0.809
0,30sss
0.137a.003
0,253a013
0.764.00s5
1.27c+0260
CopulaCPTS
0.800a.0045
0.296aam2
0.1370.003
0.260no15
0.75
1.45c+0271
DR-CP
g.805
0.257
0.507
1.14p034
0.646uu066
0.5 1.5110000
C-HDR
808
0.218
0.0372a.0m3
0.0360a.055
0.794um854
1.29e+022n
PCP
0 0.802
0.3430007
0.567ax09
1.32am
0.628ruo22
1.27e+026.1
HD-PCP
0 0.800
0.259e
0,3522.00
0 0,8030020
0.673
1 1.27e+026.1
STDQR
0 w.123
0.269omus
0.389
0.912
0. 0.667ames
86.6
C-PCP
10
0.302mm0
0.0369
0.0404gm09
0.798umsz
2.54e+02, 22101
L-CP
0. 0.805
0 0,2676n0s1
0.02033.004
0.0198amm1
0.789
0.251aa11
blog_date data
M-CP
0.802a.040
0.170
0.292a.081
0.153am7
0.736nn12
3.06c++37neee
CopulaCPTS
0.8133007
0,0948
0.313
0,23
0.742an10
5. 13c4303 te402
DR-CP
0. 0.800.00114
0.00740.0006
1.06
1.50p,3
0.644ga229
0,515a00s
C-HDR
0. 0.809.0000
0.0155
0.237a06
0.0611a019
0.751am13
5.06e+003
PCP
0.801e.0033
0.14le
0.938
1.52a.s
0.643uxo
5. 5.74c+0279e
HD-PCP
0 0.8003.000
0. 1250.23
0.794gus
0.945a22
0. 0.660uucco
5.75c+022 5.7
STDQR
0.8 180007772
0.16 131633066
0.805
0,88
0.678pn12
5.84c+02a
C-PCP
0 0.004.000
0. 106
0.1633,049
0.1 113ass
0.7644m12
3.63e+013
L-CP
e.803
0.06760017
0.327
0.0624a.23
0.722am12
0.258
calcofi
M-CP
0.023
2.1 2.130004
0,433
0.446u0t6
0.734pp0
26.4
CopulaCPTS
1.
0.480aasa
0.492a04
0. 0.746
29.6
DR-CP
0.805
1.67
1,44
1666m
0.654pmot
0.529
C-HDR
0.8 0.805emm18
1.99cm0s
0,0294..12
0.0187
0.794p.001
27.7
PCP
0. 0.802a.016
2.33
1.64
1.79p41
0.638pnos4
26.5
HD-PCP
0.802a.033
1,89
0,980ga33
1.05mm
0. 0.6833.000
27.3
STDQR
0.799
1.97000
1.13am
1.23ams
0.67 1.676unoo
26.4a.
C-PCP
0. 8099.0000
2.8 10.022
0.03322,0093
O. 0.0253
0.806
52.923
L-CP
0,800
2.70
0.03324,199
0.0179A.m40
0. 0.792nn055
0.264a.12
taxi
M-CP
0.8022.0032
4.260000
0,0585p0004
0.0421
0.784
60.6g.6
CopulaCPTS
0.822acmm
4.72x111
0.1 0.114
0.0989
0.799
68.4g
DR-CP
0.805
2.624.099
0,383
0.451
0.707
0.539
C-HDR
0.809
2.62a.00
0.0388
1.044
0.793
61.94.
PCP
0. 0.804amm
4,03
0,341
0.399
0.715
60.2g.
HD-PCP
0. 0.805
3.1880000
0.1 0.194a.022
0.219um:2
0.750pnoss
61. 61.1gs
STDQR
0.805
3.63
0,203
0.2 0.224um:
0.748pnos
32.4g,1
C-PCP
0.807
4.02a004
0.0307a.0E3
0.03388
0.802auasa
1.21c402
L-CP
0,805
4.94x12
0.02644,0000
0.0196g.005
0.796uu846
0.243
42


Reviews:
[{"review_id": "rxdwItRRKJ", "reviewer": "Reviewer_zeNu", "summary": "This paper reviews latest developments of conformal prediction methods in multi-output regression tasks.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive evaluation and constructive feedback. We address your questions below.\n\n**Comparison of CP under different conditional density estimation (CDE) variants:**\n    We agree that evaluating CP methods across different CDEs is important for robustness. We performed extensive experiments (detailed in Appendix G) using a diverse set of CDE models: two normalizing flows (MQF² [33], conditional Glow [54]), Distributional Random Forests [12] (with KDE), and a GMM hypernetwork [23, 8]. Our key findings regarding the relative performance and properties of the compared conformal methods remained stable across these different base predictors, demonstrating the robustness of our conclusions.\n\n**Comparison between our normalizing flow approach and the VAE approach of STDQR**\n\nOur primary goal was a fair comparison among conformal methods. To achieve this, we standardized the base predictor architecture where appropriate. For STDQR [20], we replaced the original CVAE with a conditional normalizing flow (specifically, MQF² [33]). This ensures that comparisons between STDQR, L-CP, and other methods are not confounded by differences in the underlying generative model architecture (VAE vs. Flow). This adaptation was also motivated by the exact invertibility and direct density evaluation offered by flows, eliminating the noise associated with VAE sampling and the need for directional quantile regression in the latent space, as discussed in Appendix F.3 and recommended by Feldman et al. [20]. While a direct empirical comparison between the two models is outside the scope of this comparative study, our approach ensures a more direct comparison of the conformalization strategies themselves.\n\n**Pairing between CP methods and conditional density estimation methods**\n\nThis is a valid point. Ideally, each CP method could be paired with its optimal CDE. However, for a unified comparative study, this would introduce confounding factors, making it difficult to isolate the performance differences attributable to the CP methods themselves versus the CDE pairings. Our approach was to select a set of strong, representative CDEs (as listed above) and apply them consistently across all applicable CP methods. This ensures a fair comparison of the conformal methods within our unified framework. Appendix G shows that the relative rankings and conclusions are largely consistent across these CDEs, suggesting that our findings are not overly sensitive to a specific CDE choice within this representative set.", "claims_and_evidence": "Yes. The paper did a comprehensive overview and detailed analysis of various methods, categorized them into different variants, and compared the results both numerically and visually.", "methods_and_evaluation_criteria": "Yes. The paper used popular datasets and common evaluation criteria.", "theoretical_claims": "The theoretical claims are well supported by the proof in appendices.", "experimental_designs_or_analyses": "The experiments are well designed.", "relation_to_broader_scientific_literature": "The work is important and provides an overview of the multi-output conformal regression methods. It can lead to future developments in this field.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Strength:\n1. The paper is a comprehensive review of an important topic in conformal prediction. \n\nWeakness:\n1. Many CP methods rely on conditional density estimation techniques and can integrate with various approaches. Were the experiments conducted to compare the performance of CP methods under different conditional density estimation variants to ensure a comprehensive evaluation?\n2. Specifically, while the original ST-DQR paper used conditional VAE, this paper employs conditional normalizing flows. Did the authors perform a comparison between these two approaches?\n3. More broadly, for a fair comparison, CP methods should be evaluated based on their optimal performance when paired with different conditional density estimation methods. Were such evaluations conducted to ensure fairness and completeness?", "other_comments_or_suggestions": "No further suggestions.", "questions_for_authors": "NA", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethical_review_concerns": "NA"}}, {"review_id": "bLSQlcXw6z", "reviewer": "Reviewer_mLNW", "summary": "The paper considers conformal prediction for high-dimensional regression. While one can extend uni-dimensional regression algorithms to multi-dimensional ones, other algorithms that explicitly work in $\\geq 1$ dimensions also exist.\n\nThis paper introduces two conformity scores, C-PCP and L-CP, exploring their connection to existing methods. Further, the paper experimentally compares different algorithms on various synthetic and real-world datasets.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback.\n\n**The paper mentions that CP²-PCP is similar to C-PCP. What is the reason for not including it?**\n\nCP²-PCP was proposed concurrently with our submission period and published very recently (ICLR 2025). Hence, we were unable to include it in our empirical comparisons. However, recognizing its relevance, we have added an in-depth discussion in Appendix H comparing CP²-PCP with our CDF-based scores (including C-PCP) and highlighting the similarities and differences in their approaches to achieving conditional validity.\n\n**The novelty seems limited to the proposed conformity scores C-PCP and L-CP**\n\nIn addition to the proposed conformity scores C-PCP and L-CP, this work is the first (to our knowledge) to categorize systematically, implement (in a unified codebase for fairness), and empirically compare a broad range of multi-output conformal methods (marginal, density, sample, latent-based) within a single framework. This synthesis reveals crucial trade-offs (Table 1), motivates our proposed methods (C-PCP/L-CP each address specific gaps), and establishes theoretical connections, which were previously unexplored.\n\n**Unclear if/when proposed scores perform better; C-HDR seems best**\n\nC-HDR indeed performs very well, especially for region size when density estimation is accurate and feasible. However, C-PCP and L-CP offer distinct practical advantages in specific scenarios:\n*   **C-PCP:** Requires **no explicit density estimation or likelihood evaluation**, only samples from a generative model. This makes it applicable to models like diffusion models or GANs where densities might be intractable or unavailable, situations where C-HDR cannot be used.\n*   **L-CP:** Achieves asymptotic conditional coverage (ACC) with **significantly lower computational cost** (Figure 6, often >100x faster) compared to C-HDR and C-PCP, as it avoids per-instance sampling or complex density calculations.\nBoth C-PCP and L-CP achieve asymptotic conditional coverage (ACC), unlike DR-CP/PCP/HD-PCP/STDQR/M-CP/CopulaCPTS. Their region sizes are competitive, generally only surpassed by C-HDR (if applicable) and DR-CP (which lacks ACC guarantees). Table 1 and the Conclusion summarize these trade-offs.\n\n**Suggestions**\n\nThank you, these are excellent suggestions for improving clarity. We have implemented them as described in [this linked PDF](https://pdfhost.io/v/VUnLttqXdJ_rebuttal) and will incorporate other comments in the final version.\n\n**Questions:**\n1.  *Challenges in high-dim CP:* Primarily, the lack of a natural ordering (unlike 1D) makes simple extensions of univariate methods (like CQR directly) difficult. This motivates methods based on joint density (DR-CP, C-HDR), samples (PCP, C-PCP), or latent spaces (L-CP, STDQR) to capture complex dependencies and define multivariate regions.\n2.  *DR-CP empty contours (Fig 2):* This occurs when, for a given $x$ and a desired coverage level (e.g., $x=1$ and $1-\\alpha$ = 0.2), no region in the output space achieves a density value $\\hat{f}(y|x)$ high enough to be less than or equal to the constant threshold $-\\hat{q}$. This visually demonstrates DR-CP's lack of conditional coverage. Methods with asymptotic conditional coverage (C-HDR, C-PCP, L-CP) adapt their thresholds or regions based on x and avoid this issue asymptotically.\n3. *Claim of new classes:* We introduce generalizations of existing univariate concepts to the multivariate setting, forming broader classes:\n- CDF-based: Generalizes the univariate HPD-split/C-HDR score [31] to work with any base multivariate conformity score $s_W$ (Eq. 11), leading specifically to C-PCP when $s_W=s_\\text{PCP}$.\n- Latent-based: Generalizes univariate Distributional CP [13] to multivariate outputs using any invertible conditional generative model Q and any latent distance function $d_\\mathcal{Z}$ (Eq. 14), leading to L-CP. Feldman et al. [20] also uses a latent space but performs CP on grid samples mapped to the *output* space, not directly in the latent space like L-CP.\n4. *Practical takeaway from Prop 1:* An interesting practical takeaway follows from the fact that DR-CP and C-HDR are linked in the same way as PCP and C-PCP. Since DR-CP can be shown to asymptotically have the smallest average region size while C-HDR empirically has a smaller median region size, similar observations are expected for PCP and C-PCP. This is verified empirically: PCP has a smaller average region size across all base predictors, while C-PCP has a smaller median region size. We have added this insight to Section 5.3.\n5. *Distance for C-PCP/L-CP:* In our experiments, PCP (and thus C-PCP) uses Euclidean distance, and L-CP uses the Euclidean norm in the latent space. Both could be generalized to other metrics.\n6.  *i.i.d. vs exchangeable data:* Thank you for pointing this out, the analysis remains valid. We will relax this assumption.\n7. *Acronyms:* C-PCP = CDF-based Probabilistic Conformal Prediction. L-CP = Latent-based Conformal Prediction.", "claims_and_evidence": "This paper is a comparative study.", "methods_and_evaluation_criteria": "The paper uses a variety of synthetic and real-world datasets. Furthermore, the paper compares multiple algorithms for its study. The paper mentions that CP$^{2}$-PCP is similar to C-PCP (lines 195-197, column 1). What is the reason for not including it?", "theoretical_claims": "I checked the correctness of Proposition 1.", "experimental_designs_or_analyses": "The paper evaluates the algorithms on marginal coverage, conditional coverage, set sizes, and computational time. A few questions/suggestions:\n\n1.\tReplacing Fig. 5 with a table containing mean and median set sizes would help. The actual values convey more information than relative rankings.\n\n2.\tIs there a quantitative analysis for Section 5.1?\n\n3.\tFigs. 7 and 8 are difficult to compare.\n\n4.\tTables 3, 4, and 5 should bold the statistically significant values (accounting for the standard deviations).", "relation_to_broader_scientific_literature": "The paper compares different conformal prediction methods for high-dimensional regression. It also introduces 2 new conformity scores.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Strengths:\n\n1.\tThe paper considers conformal prediction for high-dimensional regression, which is of practical importance.\n\n2.\tThe paper compares different conformity scores on various datasets. This comparison includes empirical results and some theoretical links.\n\nWeaknesses:\n\n1.\tThe novelty seems limited to the proposed conformity scores C-PCP and L-CP.\n\n2.\tIt is unclear if and when the proposed conformity scores perform better than the existing ones. C-HDR seems to perform the best.", "other_comments_or_suggestions": "1.\tThe existing methods that C-PCP and L-CP build on should be discussed in detail in the main paper. For example, details of directional quantile regression [Feldman et al., 2023].\n\n2.\tThe non-conformity score need not be deterministic (lines 83-84, column 2). As long as the training data is independent of the calibration and test data, marginal coverage is guaranteed.\n\n3.\tThe conformity score by Sadinle et al. [2016] is for classification. Even though one can use it for regression, one should explicitly mention that.\n\n4.\tThe citations need to be corrected. For example, Sadinle et al. [2016] was in the Journal of the American Statistical Association in 2019.\n\nTypo:\n\n1.\t\"...$D_{cal} \\rightarrow \\infty$...\" $\\rightarrow$ \"...$|D_{cal}| \\rightarrow \\infty$...\" (line 162, column 2)", "questions_for_authors": "1.\tWhat are the challenges with extending conformal prediction to high-dimensional regression (Section 4)?\n\n2.\tWhat causes the absence of contours in Fig. 2 (lines 260-264, column 1)? The method should still construct a set.\n\n3.\tI am confused about the claim of the introduction of new classes of conformity scores (lines 124-126, column 2). Don't these classes already exist? The paper also provides existing methods for both classes.\n\n4.\tIs there a practical takeaway from Proposition 1?\n\n5.\tDo C-PCP and L-CP use Euclidean distance or a different distance?\n\n6.\tDoes the analysis change if the setup switches to exchangeable data (Section 2)?\n\n7.\tWhat do the acronyms C-PCP and L-CP stand for?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ijsbMHfkkp", "reviewer": "Reviewer_ReUg", "summary": "This paper performs a unified comparative study of existing conformal methods with different multivariate base models for constructing multivariate prediction regions. It generalizes two classes of conformity scores from the univariate to the multivariate case. Moreover, it conducts large-scale experiments comparing the different multi-output conformal methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and recognizing the value of our work.\n\n**Proposed scores seem not to have finite-sample conditional guarantees:**\n\nYou are correct. Our focus is on developing flexible conformal methods applicable to complex, modern generative models (density, sample, or latent-based) with minimal assumptions on the underlying data distribution or the base predictor's accuracy. Under such weak, distribution-free assumptions, achieving exact finite-sample conditional coverage is generally impossible, as proven by Barber et al. (2019) [5]. Instead, we provide asymptotic conditional coverage (ACC) guarantees (Appendix E.2) for C-PCP, L-CP, and C-HDR. This means that as the base model converges to the true distribution and the calibration set size increases, these methods provably achieve the desired conditional coverage level. Our experiments (Fig 3) further empirically validate that these methods also achieve superior empirical conditional coverage in finite samples compared to methods lacking ACC guarantees.", "claims_and_evidence": "Yes, the multivariate conformity scores introduced ensure asymptotic conditional coverage while maintaining exact finite-sample marginal coverage. The proofs are provided.", "methods_and_evaluation_criteria": "Yes, the comparative study and the large-scale experiments are evaluated on a wide range of conformal methods and datasets.", "theoretical_claims": "I did not go through the proofs for the theoretical claims.", "experimental_designs_or_analyses": "I read the descriptions of the experiments in the paper and did not see major issues.", "relation_to_broader_scientific_literature": "This paper generalizes two commonly used classes of conformity scores to the multivariate case. This is an important contribution as the classes of multivariate conformity scores are underexplored in the literature. Moreover, this paper provides an extensive comparison of existing conformal methods.", "essential_references_not_discussed": "Not that I know of.", "other_strengths_and_weaknesses": "Strengths:\n- The paper is well-written\n- This work studies multivariate conformal prediction, which is relatively underexplored. This can provide a valuable direction for future research.\n- Two multivariate conformity scores, which generalizes existing univariate conformity scores, are introduced. Theoretical analyses are provided. \n- This work gives a detailed and unified comparative study for existing multivariate conformal models.\n- The experiments cover a wide range of methods\n- The methods are implemented using a unified code base, ensuring fairness\n\nWeaknesses:\n- The proposed scores seem not to have finite-sample conditional guarantees", "other_comments_or_suggestions": "This work performs an extensive study on multivariate conformal methods and introduces two multivariate conformity scores. This paves the way for a research direction which is relatively underexplored.", "questions_for_authors": "No", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NnPkxAFZmy", "reviewer": "Reviewer_g7mK", "summary": "This paper provides an overview of multi-output conformal regression methods, putting them in a unified setting, and propose two new approaches based on scores CDF, that generalizes some previous methods, as well as a latent-based approach, that generalizes other families of approaches.\n\n## update after rebuttal\n\nI thank the authors for their replies. Overall the paper was quite enjoyable to me, and while as raised by other reviewers there are certainly ways to get further results, I do believe it makes a nice contribution with respect to the current literature on mutli-output conformal regression.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thorough, positive review and your support for acceptance. We appreciate you finding the paper enjoyable and the methods sound. We will correct the reference issues you kindly pointed out in the final version.", "claims_and_evidence": "Yes, the claims made about the new method, notably in terms of conditional coverage, are supported both by theoretical proofs (appendices E.2., that seem correct as far as I could tell with the time I could spend on it), as well as by extensive experiments (appendix F, notably).", "methods_and_evaluation_criteria": "The methods and evaluations criteria are classical metrics and arguments from conformal prediction papers dealing with multi-output regression, as well as some new unbiased metrics regarding the region size (discussed in appendix F). Overall the used methods appear to be sound.", "theoretical_claims": "I checked the claims as best as I could, and they appear correct to me. However, given the average size of ICML submissions (say, 30 pages with supplementaries included) and the number of reviews required by the confrence, my checks were not as thorough as they would have been for, say, a journal paper.", "experimental_designs_or_analyses": "I check the experimental designs, and all appeared sound to me. Graphs and results appear to be in-line with expectations.", "relation_to_broader_scientific_literature": "Yes, in the sense that providing reliable multi-output regression regions is a generally interesting problem.", "essential_references_not_discussed": "Most references I could think of were present in the paper, that appears to be very well educated in the field.", "other_strengths_and_weaknesses": "Overall, this is a quite enjoyable paper that wraps up many method of the literature and provides in addition some interesting new ones. I see no major weaknesses in the paper. My only regret is that this probably should have been a journal paper rather than a 6 pages paper with 30 pages of supplementaries, but I guess this is how things are done nowadays.", "other_comments_or_suggestions": "No specific comment, except that authors should double check references (some do not compile as expected, e.g., [37] and [61])", "questions_for_authors": "No specific questions to the authors, this is a nice work overall.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nwYeSyTOWN", "reviewer": "Reviewer_jr27", "summary": "The paper addresses the problem of applying conformal prediction to multiple continuous output prediction systems. The challenge is obtaining a small average region with a reasonable computation complexity. The paper presents a unified view of the current methods and proposes two novel approaches.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and suggestions for improvement.\n\n**Missing relation to family-wise error rate (FWER) control methods:**\n\nThank you for raising this point. We do briefly mention multiplicity control approaches in Appendix A (citing Timans et al., 2024 [59]). FWER methods (like Bonferroni or Holm corrections on univariate CPs) are indeed a related direction for constructing multivariate regions, typically by taking Cartesian products of univariate regions $\\hat{R}\\_i$. Their strength lies in leveraging well-understood univariate CP. However, our work focuses on methods that explicitly model the joint distribution $F\\_{Y|X}$ (using e.g., multivariate NFs, DRFs, GMMs, or generative sampling) to construct potentially non-rectangular regions that can better capture output dependencies and often result in smaller region sizes (as seen empirically). While FWER methods are valuable, especially when only marginal predictions are available, the methods we survey and propose aim for potentially tighter regions by directly modeling the multivariate structure. We have revised Appendix A to better clarify this distinction.\n\n**Conclusion unclear/short:**\n\nWe agree that the conclusion was too brief and lacked practical guidance. We have expanded it to provide clear, actionable takeaways based on our unified study and experiments, summarizing the trade-offs:\n*   **If only univariate margin models are available/desired:** Use M-CP or CopulaCPTS (fast, but hyperrectangular; may poorly capture dependencies).\n*   **To minimize *average* region size (potentially sacrificing conditional coverage):** Use DR-CP (if density $\\hat{f}$ available), else STDQR (if invertible latent model Q available), else PCP (if only samples $\\hat{Y}$ available).\n*   **To achieve good *conditional coverage* and small *median* size:** Use C-HDR (if density $\\hat{f}$ available), else L-CP (if invertible Q available), else C-PCP (if only samples $\\hat{Y}$ available).\n*   **For good conditional coverage with lowest *computational cost* among ACC methods:** Use L-CP (if invertible Q available).\n\n**Questions**\n1.  *Relation to FWER:* Please see our response above regarding the distinction (modeling joint vs. combining univariate CPs, region shapes, and dependency handling). To better understand the prediction regions produced by FWER control methods, we also provide a qualitative and quantitative comparison with Bonferroni correction on our dataset from Figure 1 in [this linked PDF](https://pdfhost.io/v/gqCAunzaRU_rebuttal-jr27). The method is fast computationally but produces larger regions due to the rectangular shape, similarly to M-CP.\n\n2.  *Need high-dimensional/unstructured data experiments:* We agree on the importance of high-dimensional settings. To this end, we included an experiment on CIFAR-10 image generation (Appendix I), where the output is a 3x32x32 image ($d$=3072). We used a conditional Glow model [34, 54]. The results (Table 3) largely confirm our findings from tabular data: C-PCP, L-CP, and C-HDR demonstrate superior conditional coverage (WSC, CEC) compared to other methods even in this high-dimensional setting. We agree this experiment deserves more visibility and will ensure it is also mentioned in the introduction.\n\n3. *How is the density $f(y|x)$ known?*\nCrucially, we do not assume the true density $f_{Y|X}$ is known. All density-based methods (DR-CP, C-HDR, HD-PCP) rely on an estimate $\\hat{f}(y|x)$ learned from data. In our experiments, $\\hat{f}$ is provided by models like normalizing flows (MQF², Glow), DRF+KDE, or GMMs (detailed in Appendix F.2). Conformal prediction takes this potentially imperfect estimate $\\hat{f}$ and calibrates the resulting prediction regions (using the calibration set) to provide valid finite-sample marginal coverage guarantees, regardless of how accurate $\\hat{f}$ is. Methods like C-HDR, L-CP, and C-PCP additionally provide asymptotic conditional coverage guarantees as long as the estimator $\\hat{f}$ (or sampler $\\hat{Y}$, or latent map $\\hat{Q}$) converges to the true data generating process.", "claims_and_evidence": "yes", "methods_and_evaluation_criteria": "yes", "theoretical_claims": "yes", "experimental_designs_or_analyses": "Yes", "relation_to_broader_scientific_literature": "The main contribution is a unified treatment of previously proposed methods.", "essential_references_not_discussed": "There is a closely related research direction of multiple testing in parallel in conformal prediction and family-wise error rate control. Popular methods are Bonferroni and Bonferroni-Holm  corrections of the rectangular\nshaped region.", "other_strengths_and_weaknesses": "The topic is important, and this paper is well-written and provides an up-to-date comparative review of current methods. It can be a good review paper for people interested in the topic.\n\nAfter comparing many methods it is not clear what the conclusion is.  The conclusion section is short and non-informative. \n\nThe presentation could be improved. The notation is sometimes too dense, making it challenging to follow.\n\nThe introduction could provide a clearer motivation for the need for multi-output conformal methods.", "other_comments_or_suggestions": "no", "questions_for_authors": "How is it related to family-wise error rate control methods? Popular methods are Bonferroni and Bonferroni-Holm corrections of the rectangular. I was surprised that a review paper didn't even mention this line of research.\n\nWhile the tabular datasets are informative, additional experiments on real-world high-dimensional or unstructured data (e.g., images, text) would better demonstrate the practical impact of the proposed methods.\n\nFor the methods DR-CP, C-HDR, and others, It is not clear how you know the density f(x|y) from the data.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper provides a highly valuable and thorough contribution to the field of multi-output conformal prediction by systematically categorizing and unifying existing methods. Its primary strength lies not in introducing highly novel techniques but in offering a clear, well-structured framework that synthesizes prior work, making it easier for researchers to understand the landscape of this important area. The authors go beyond a simple survey by conducting detailed empirical evaluations across diverse datasets, supported by rigorous theoretical analysis. This combination of breadth and depth ensures the work serves as both a useful reference and a practical guide for future research.

The writing is clear and polished, making the paper accessible even to those less familiar with multi-output conformal prediction. The inclusion of a unified codebase further enhances the quality and impact of the work, promoting reproducibility and providing a practical tool for researchers and practitioners alike. The authors’ efforts to standardize comparisons across methods and conduct experiments under consistent conditions demonstrate a commitment to scientific rigor and fairness.

While the proposed methods, C-PCP and L-CP, are not highly novel, they are thoughtfully integrated into the broader framework, and the authors provide an honest and balanced discussion of their trade-offs relative to existing approaches, such as C-HDR. Reviewer concerns around novelty, finite-sample guarantees, and other comparisons were addressed comprehensively in the rebuttal, with the authors clarifying key points and improving the conclusion to include actionable insights and practical takeaways.

Overall, this paper exemplifies high standards of academic writing, organization, and execution. It fills an important gap in the literature by systematically analyzing and categorizing multi-output conformal prediction methods while providing a valuable resource for both new and experienced researchers in the field.
Decision: Accept (poster)