CONTRA: Conformal Prediction Region via Normalizing Flow Transformation
Authors: Zhenhan Fang
Aixin Tan
Jian Huang
Venue: ICLR 2025 Poster (2025)
Keywords: Uncertainty Quantification
Conditional Density Estimation
Conformalized Quantile Regression
Coverage Probability
Multi-dimensional Prediction
Primary area: generative models
TL;DR: Utilizing normalizing flows to construct conformal prediction regions aligned with the estimated conditional density.

Abstract: Density estimation and reliable prediction regions for outputs are crucial in supervised and unsupervised learning. While conformal prediction effectively generates coverage-guaranteed regions, it struggles with multi-dimensional outputs due to reliance on one-dimensional nonconformity scores. To address this, we introduce CONTRA: CONformal prediction region via normalizing flow TRAnsformation. CONTRA utilizes the latent spaces of normalizing flows to define nonconformity scores based on distances from the center. This allows for the mapping of high-density regions in latent space to sharp prediction regions in the output space, surpassing traditional hyperrectangular or elliptical conformal regions. Further, for scenarios where other predictive models are favored over flow-based models, we extend CONTRA to enhance any such model with a reliable prediction region by training a simple normalizing flow on the residuals. We demonstrate that both CONTRA and its extension maintain guaranteed coverage probability and outperform existing methods in generating accurate prediction regions across various datasets. We conclude that CONTRA is an effective tool for (conditional) density estimation, addressing the under-explored challenge of delivering multi-dimensional prediction regions.

Published as a conference paper at ICLR 2025
CONTRA: CONFORMAL PREDICTION REGION VIA NORMALIZING FLOW TRANSFORMATION
Zhenhan Fang, Aixin Tan
Jian Huang
Department of Statistics and Actuarial Science
Department of Applied Mathematics
The University of lowa lowa City, IA 52242, USA (zhenhan-fare fang, aixin-tame mixin-tanjeuiowa..
The Hong Kong Polytechnic University Hongkong, China
j. huang@polyu. edu hk
ABSTRACT
Density estimation and reliable prediction regions for outputs are erucial super- vised and unsupervised learning. While conformal prediction effectively generates goverage-guaranteed regions, it struggles with multi-dimensional outputs due to reliance one-dimensional nonconformity scores. To address this, introduce CONTRA: CONformal prediction region via normalizing flow TRAnsformation. CONTRA utilizes the latent spaces of normalizing flows to define nonconformity scores based on distances from the center. This allows for the mapping of high- density regions in latent space sharp prediction regions in the output space, surpassing traditional hyperrectangular or elliptical conformal regions. Further, for scenarios where other predictive models are favored over flow-based models, we extend CONTRA to enhance any such model with a reliable prediction re- gion by training simple normalizing flow on the residuals. We demonstrate that both CONTRA and its extension maintain guaranteed coverage probability and outperform existing methods in generating accurate prediction regions across var- ious datasets. We conclude that CONTRA is an effective tool for (conditional) density estimation, addressing the under-explored challenge of delivering multi- dimensional prediction regions.
INTRODUCTION
In unsupervised and supervised learning, density and conditional density estimation are frequently used to communicate the inherent variability the potential outcome and assess the reliability of estimations and predictions (Hall & Yao 2005 Guhaniyogi al. 2014; Dalmasso al. 2020). These methods can report regions that aim to capture the true outcome with a specified probability. such 90%, based the estimated density. However, the actual coverage rate depends the underlying model assumptions and does not come with any inherent guarantees.
Conformal prediction is way produce regions with guaranteed coverage rate of future outcomes non-asymptotically and free of model assumptions (Papadopoulos et al., 2002] Vovk et al., 2005 Lei
al.) 2018, Romano et al. 2019 Lei al. 2015 Izbicki al.) 2020 Chernozhukov et 2021
Briefly, the most popular split conformal approach partitions the whole dataset into a proper training set and a calibration set. predictive model is trained using the former, and non-conformity scores are calculated on the latter that measure the deviation of the predictions from the actual outputs. Quantiles of the non-conforyy scores are used threshold values, producing prediction regions that capture future outcome with desired probabilities. More accurate predictive models usually yield tighter conformal prediction regions.
The conformal prediction literature has primarily focused on one-dimensional output settings, partly because the conformal idea depends on quantiles. Among the relatively few multi-targeting meth- ods in the literature, most restrict region shapes to be boxes ellipsoids. For complex conditional distributions of the output, such as those with multi-modes or unequal-tails boxes or ellipsoids encompass low-density areas to ensure validity. and become inflated. Two recent approaches that allow flexible shape for calti-dimensional conformal regions are the spherically transformed directional quantile regression (ST-DQR) (Feldman al. 2023) and the probabilistic
different directions.

---

Published as a conference paper at ICLR 2025
conformal prediction (PCP) (Wang al., 2022). But their prediction regions composed of nu- merous balls, which often lead to highly irregular boundaries and more disconnected regions than desirable, hindering interpretability.
We propose CONTRA, a new method that produce prediction regions for multi-dimennional outputs. CONTRA stands for CONformal prediction region via normalizing flow TRAnsformation. Here, normalizing flows (NF) (Dinh et al. 2016 Huang et al. 2018 Chen et al. [2019] Kingma & Dhariwal
2018 Papamakarios al. 2021) and conditional normalizing flows (CNF) (Winkler al. 2019) generative methods that provide samples, hence density estimates, for the output.
An example is shown in Figure I for the prediction region of dropoff locations of Taxi (y) given arbitrary pickup point (x, shown as blue pin) in New York city. This is a situation where the con- ditional density of the outcome is multimodal. Ten different conformal predictions were displayed. The five methods with restricted shapes in (e) and (f) yield large and unnatural prediction regions. The CONTRA region in (b), a variation of CONTRA in (c), and the PCP and ST-DQR methods based on either NF diffusion models in (d), (g) and (h) produced comparable areas with similar sizes. But both CONTRA prediction regions in (b) and (c) are much more connected than that of PCP and ST-DQR, providing the most effective and interpretable results that maintained theoretical rigor.
KDE (unsafe)
(b) CONTRA
(c) ResCONTRA
(d) NF + PCP
(c) Ellipses
(f) Rectangles
(g) DM PCP
(h) DM + ST-DQR
Figure I: NYC Taxi data. Prediction regions of drop-off' location given a pickup coordinate (blue pin). (a) shows a 90% high-density region estimated via kernel density estimates (KDE) based on dropoff locations from the 200 nearest pickups the blue pin. While there guarantee of coverage, the shape of this region offers an informal visual reference to help users assess conformal regions. (b)-(f) show various conformal predictions based on NE, trained on 3600 samples and calibrated on 1200. The methods include our proposed CONTRA and DesCONTRA: PCP and ST-DQR (the latter result is not shown due to its similarity to the PCP); NLE and RCP that always lead to elliptical regions; and MCQR, Dist-Slpit, and CQRppn that always lead to rectangular regions; (g) and (h) show the PCP and ST-DQR conformal regions based diffusion model, in particular, Denoising Diffusion Probabilistic Model 300 epochs.
(Ho al. 2020) trained with 800 timesteps, learning rate 0.001 and
We briefly explain here how CONTRA works. Given an input, CNF learns bijection that maps latent variable from simple base distribution, like standard multivariate Gaussian, to the outcome variable (details in Section 2 After training a CNF on the proper training set, CONTRA maps high- density regions (HDRs) of the base distribution to the outcome, while the split-conformal approach adjusts these regions in the latent space by comparing latent representations of the calibration set to the expected base distribution. Consequently, CONTRA's prediction region is a bijection of single high-density area, naturally aligning with the output's conditional distribution, circumventing the need union many regions, seen in PCP and ST-DQR.
While How-based models excel in many contexts, there are times when other predictive approaches preferable. For these situations invent variation of CONTRA called ResCONTRA, which enhance any predictive model with a reliable prediction region by training simple NF on the resid- uals. The ResCONTRA approach builds on Colombo (002)), where one-dimensional residuals were
2

---

Published as a conference paper at ICLR 2025
transformed follow base distributions like Uniform Gaussian using relatively simple bijec- tions. ResCONTRA extends this idea to nulti-dimensional outputs and introduces the concept of symmetric base distribution in Rº en route to define ideal non-conformity score. In describ- ing ResCONTRA in section 3.4, we also highlight the need for a three-step calibration procedure to ensure exchangeability that leads to coverage guarantee. Compared to CONTRA, which trains single model, ResCONTRA seems less efficient as it trains two models on smaller data sets. But ResCONTRA has the potential to excel in situations, say, when [E(y)x) (x) is highly intricate and can be better learned with techniques like XGBoost rather than an underlying bijection. Some empirical comparisons of the two proposed methods can found Section and Appendix
We provide details of CONTRA and ResCONTRA in Section 3 review existing multi-output confor- mal predictions Section 4 We demonstrate that CONTRA and ResCONTRA achieve the desired coverage probability (e.g., 90%) with high accuracy in Section They outperform shape-restricted methods by providing smaller regions that better reflect the true density. And they have smoother boundaries and improved interpretability compared the flexible regions of PCP and ST-DQR.
2 CONDITIONAL NORMALIZING FLOW
Let (X, Y) denote a random vector with unknown probability density (or mass) function Pxy y). supported on S C ((++)), For each x E X' C I,, the support of Y is denoted by V/ C R Suppose data set D {+), )= has been drawn from pxy (x, y). The objective of conditional density estimation method is to find an estimate i' x(-Ixx ) of the true conditional density (-)x(-Ix) for all x E X' based on D. In a parametric approach, where candidate models are = {pY|x,o; 0 E }].
the maximum likelihood principle leads to pvpx (-1x) pxix,a where Ô is given by:
n
i = arg max 0
II (y)1)) = arg max Σ i=1
log pY)x,,(y/(xx.)
(1)
The CNF approach specifies P to be such that, each P (-Ix) is a density of a q-dimensional random vector Y that be transformed via differentiable bijection to a simple random vector Z, such as a g-dimensional Gaussiar]] We denote this transformation by
y = (z)
such that for any x X, ((.x)) is differentiable bijection from 2 to Vxx with inverse function (-,xx. Let det denote the determinant function. The
variable technique implies:
Pr(x,(ly|x)) p2((tt (y, x)) det 0l (y.
θy
can be extremely challenging to directly constructing an expressive enough collection of transfor- mations, (to, 0 E 0). capture the true conditional density, Py]x (y|x). Fortunately, the compo- sitional nature of bijective functions allows complex transformations to be built from simpler ones: Lo(:; x)
bm,0 0 6m-1,0 ... (,, Accordingly,
Σ log det Itg.
log p log Pz (z) -
(1-1-
1=1
where Zo z, Zi t1,0 (Z) for / m. and Zm to x) y. And equation | can be rewritten as:
i arg max
[og pz pz(z.) Em log det
0t, az." (n
From cquation it's essential to specify transformations t, that are easy to evaluate Jacobian determi- nants for. Several structured approaches, including autoregressive flows (Huang et al.) (018). linear flows (Kingma & Dhariwal 2018), and residual flows (Chen et al., 2019). facilitate such tractability. In this paper, employ affine transformation model known realNVP (Dinh et al. (016)., which leverages coupling layers (Dinh et al. 2014) to improve computational efficiency. Further details on
reaINVP are included in Appendix C
'The general CNF framework allows the base distribution of to depend on and A. We intentionally used free of that the latent z different values of are comparable and can be pooled for calibration.
3

---

Published as a conference paper at ICLR 2025
3 CONTRA: CONFORMAL REGION VIA NORMALIZING FLOW TRANSFORMATION
Having trained CNF model, 1g. it is tempting to perform naive conditional density estimation as follows. Given any x of interest, define the naive (1 - a) prediction region to be
(x+4() {y : y 1g(2, xa+1), 2 Bı-a),
where Bi-a is the q-dimensional ball that the highest probability region of size (1 (1-a) a) in the latent space for the standard Gaussian. However, the actual coverage rate of this region can vary greatly depending on the quality of y = tg(z, (.. Xn+1), Z ~ Ng I) as an approximate sampler for PYIxxx..1 Approximation emrors can be high for X values not well represented in training.
Given x we would like methods to construct ((x,+1) C IR that satisfies the marginal coverage guarantee:
P'yn+1 E ((x++1)) 2 a.
3.1 CONFORMAL REGIONS: FROM THE LATENT SPACE TO THE OUTPUT SPACE
Two general ideas to achieve equation 3 are full conformal prediction (Vovk et al. 2005) and split conformal prediction (Papadopoulos al. 2002 Lei al. 2018 2015 Angelopoulos & Bates 2021]. Full conformal prediction requires recalculating non non-conformitt scores across the entire dataset for each new instance. Split conformal prediction partitions the whole dataset into two dis- joint sets: the proper training set, D2 y.) E } of size nı; and the calibration set, D2 = y.) E I2] of size n2 = n n1. We choose the split approach for computational efficiency in developing CONTRA.
Recall (=( +] S + RR stands for the CNF model with a standard Gaussian base trained from D1. Given any (x, y). there is a latent representation of the output, which we denote by = 1 x).
(y,
Denote the collection of latent representations for D2 by
Eeal {2. IR 2; 1 (y i E T2)
(4)
Let r1-a denote the [(1 - ( + 1)] smallest member of [1|2:112, E T2], where || is the Euclidean norm. We define the conformal ball size (1 - a) to be
A = (z TR. || S ri-a).
(5)
Then E contains at least - a 100% of the points in Zeal- By the inflation of quantiles lemma (Romano et al., 2019), assuming the points in D2 and the point to predict, (X4t1,n411) are ex- changeable, we have
P(2n+1 E É) = P( +1 S ri-n)) 2 a.
The CONTRA prediction region is defined to be the mapping of E in the output space:
C(xa+1) = (B xx+1).
The algorithm for producing CONTRA prediction region is summarized Algorithm Its coverage guarantee is stated in Proposition in Appendix A.1
Algorithm 1 Conformal Region via Normalizing Flow Transformation (CONTRA)
Input:
I: Data (s), E RR x Ro.
2: Miscoverage level a E [0, 1).
3: A CNF algorithm with a standard Gaussian base distribution.
4: A point that needs prediction region for its output, yn+1-
Procedure:
1: Randomly split {(x,,
y:), into two disjoint sets D1 and Dz-
2: Fit fg by A(Dr).
3: Obtain Zeal in equation 4
4: Compute ri-o and define E in equation
5: Compute ((xx+1) lg( (F, xn+1).
Output
A prediction region for y is given by ((x++1).
4

---

Published as a conference paper at ICLR 2025
3.2 PRACTICAL ASPECTS OF PRESENTING CONTRA OUTPUTS
There is flexibility in implementing steps and 5 of Algorithm and in displaying the prediction region. One way is to get samples inside E and map them to the output space. Samples can be random or deterministic, using methods like Monte Carlo, grid points, or Quasi Monte Carlo. Another less costly way to only sample from the boundary of É, which will map of Ò in the output space. The validity of this approach hinges on Proposition Iin Appendix [A. which states that boundaries of sets are preserved under homomorphisms, including the CNF transformation.
3.3 CONNECTEDNESS AND VOLUME OF CONTRA
There is universal criteria for comparing different prediction regions, provided they have guaran- teed coverage. However, smooth boundaries and smaller volumes are generally considered desirable characteristics.
Connectedness, Fewer disconnected sets and smoother boundaries for a prediction region means predictions close each other are more likely to be classified the same way, either inside or outside of the prediction region. These properties contribute to robust and interpretable inferences in practice.
The prediction region Ò of CONTRA is indeed closed and connected due to the following.
Proposition 1. (James 2000 Chap.3) Suppose E C Z is closed and connected, and is a homeo- morphism. Then {(E) y is also closed and connected.
Volume calculation. Assessing, calibrating, and comparing different prediction regions necessitates tool to calculate the volume of the regions, which can be challenging for irregular multi-dimensional shapes. We propose an approximation method for the volume of C(x.+1). using by-products from the training process of CNF, Note that
Vol (C(x)+1)) = (ce((J,(z))) da Vol(E) / [de((Jr,(z)) (ol(E)/(dee(urre))))
dz.
Vol(E)
The above integral can be approximated using a Monte Carlo estimator based on a random sample, {zu / drawn uniformly from É, with density 1/Vol
B
Vol (C(xn+z)) Vol(E)
1
Σ
det((J,, (2+))
B
3.4 EXTENSION TO WORK WITH OTHER PREDICTION MODELS: RESCONTRA
The CONTRA method proposed so far restricts the fitted model be NF. We now present vari- ation, ReSCONTRA, to enable conformal prediction for any user-chose prediction methods. As explained in the introduction, this is inspired by an idea in Colombo (2024). We split the whole dataset into D1, D2 and Dz. Here, D, is used to train the user-chosen point estimator, f. Residuals, r = y. - fxx,). for points in the latter two sets are calculated. Then the standard CONTRA is applied to Dē the calibration set respectively in the split conformal framework. Denote the NF that transforms the residuals in D to latent z by t- and the conformal ball Ex is calibrated with Dj. Finally, the prediction region for new data point is given by
and D = [(x,,r)) r.)) which serve the proper training and
Ca (xx+1) = f(xx+1) + (E",
It's straightforward to see that marginal coverage guaranteed, as
P(yn+1 E (* =P(()) PP E ("(E",xati) 21 a.
Remark: If the distribution of z were perfectly independent of x, Corollary 2.6 of Colombo 2024)) suggests that conditional coverage probability of yn+1 given any can be achieved. However, since perfect independence is not achievable practice with finite data, they also derived in Theorem 2.7 theoretical bound for the potential reduction in conditional coverage probabilities. This bound
5

---

Published as a conference paper at ICLR 2025
depends the deviation between the learned distribution of and the ideal base distribution that is independent x. Since both CONTRA and ResCONTRA aim to transform z to have symmetrical distributions free of x, both methods are expected approximately achieve the desired conditional coverage similar to those shown in Coloobo (2002)) is an on-going work to analyze the conditional coverage probabilities of CONTRA and ResCONTRA both theoretically and empirically, with the challenging but important goal of deriving practically useful bounds.
4 RELATED WORK
In contrast CONTRA, traditional conformal prediction methods construct prediction regions di- rectly in the output space. We review some of them below.
4.1 CONFORMAL APPROACHES THAT TARGET MULTI DIMENSIONLL OUTPUTS
Limited conformal methods in the literature target munti-dimensional outputs directly. An example is the robust conformal prediction (RCP) (Johnstone Cox (021), which employs the global covari- matrix produce ellipsoid region. The normalized locally ellipsoid (NLE) (Messoudi al.. 2022) extends RCP by incorporating local covariance matrices, making the region adaptive to x.
Two recent methods, PCP and ST-DQR, are closer to CONTRA as they don't restrict shape of the prediction regions. First, the PCP method was proposed by Wang et al. (2022). Given a generative model for estimating the conditional density, PCP defines the non-conformity score of data point (x.y y:) as
S;
min ISRCK
ly: -/-l,
where (y.).k=z is sample of the output generated from the estimated conditional density, Then the prediction region at X,+1 is obtained by generating a sample {5 === and form
K
K
T (x 1)= = Uty lly - shtalla S s1
k=1
where S1-0 is the [[] [1- a + 1)7-th smallest member of {st, E 12). Note that each PCP prediction region is the union of K balls. They have flexible, but often irregular, disconnected shapes, and are sensitive to the choice of K and a. This brings challenge to interpreting the regions.
Next, the ST-DQR method was proposed by Feldman et al. (2023). It learns an r-dimensional latent representation of the output, e.g., using the conditional variational auto-encoder (CVAE). The latent variable is encouraged to follow unimodal distribution, so that methods like the directional quantile regression (DQR) are applicable to form convex probability regions for it. Samples are generated in the output space corresponding to points in the latent region. Calibration and the final prediction region are created similarly to the PCP.
ST-DQR and our CONTRA are similar in that they both depend on latent representations, but differ substantially in the latent representation. CONTRA uses bijection and are often able to learn the latent variable follow the Gaussian reference distribution rather closely. In contrast, the latent variable in ST-DQR is typically of lower dimension than the output, and its distribution is only coarsely similar to reference. And this is why additional steps like DQR are needed to form latent probability regions in ST-DQR, but our CONTRA can directly use HDR of the reference Gaussian. When calibrating the regions, ST-DQR (and PCP) had introduce yet another step: union balls around generated samples, and calibrate the common radius of the balls. Whereas CONTRA directly calibrate radius of the Gaussian HDR. Afterall, these methods have their strengths and weaknesses, but each is valuable addition to the user's toolkit.
4.2
ONE- DIMENSIONAL CONFORMAL APPROACHES AND UPGRADING THEM TO CAPTURE MULTI- DILLII-DIMESIONNL OUTPUTS
The literature conformal predictions for one-dimensional output is rather rich and we list a few state-of-the-ar methods, Locally adaptive split conformal prediction (Lei et al. 2018) enhances

---

Published as a conference paper at ICLR 2025
reliability by adapting to local data variability, while Dist-split (Izbicki et al. 2020 constructs in- tervals based on the estimated cumulative density function. However, both methods tend to provide broader intervals when the underlying conditional density is not symmetric. Conformalized quantile regression (CQR) (Romano et al. 2019) offers asymmetric intervals using two conditional quantile estimators, but it remains inefficient for multi- distributions. This limitation be addressed by CD-split (Izbicki et al. 2022). which can generate discontinuous prediction intervals, but may encounter stability issues.
Although these methods are designed for handling one-dimensional output, they can be building blocks for valid prediction regions for multi-dimensionn outputs. Bonforonni method (Bonferron) Dunn (961) is a simple but conservative way to do so. Alternatively, we provide in section new method called multi-target conformalized quantile regression (MCQR) that lifts tradi- tional one-dimensional CQR to work for multi-dimensional outputs.
1936
4.2.1
4.2.1 MCQR: MULTI-TARGET CONFORMALIZER QUANTILE REGRESSION
As in CQR, we train on D1 a pair of lower and upper quantile estimators, Qj and Qÿ, for each dimension of the output = (y...... Vd). Given a weight vector w = (W11, W12.--. Wo1 Wg2), the choice of which to be discussed later, we define non-conformity of point (x as:
{ WSI (Q)(()) - -y)), -yylys-------- - -3(x))
8 max 100nn
Then we set the prediction region be
MMCQ (MCQR(Xn+1)> = y (x(+1)
81-a
S yj S Q) (xn+1) + 81-a . = 1, 1.....
Wy1
Wy2
For any choice of the weight vector w = (WIL, W12, Wgl, W22), the MCQR satisfies the marginal which we prove Appendix A.3
coverage guarantee defined in equation 3
Finally, since ÖMCOR
(x,+1) is box, its volume is simply
Vol( (MMOR
(xn+1))= = II (x(x))+))))) - 22 (xn+1) + j=1
S1-a
$1
Wiz
Wil
Coming back to the problem of how to specifyy the weight vector w, solution is to minimize the average volume of prediction regions in the calibration set, that is
Σ Vol (MMCOR(xx))
arg min w
712 iEf2
5 EXPERIMENTS
In this section, we systematically compare the performance of the proposed CONTRA and ResCON- TRA to that of other conformal prediction methods reviewed in section including PCP, NLE, RCP, MCQR, Dist-split and CQR. The last two were designed for one-dimensional outputs. We employ the Bonferroni approach to produce valid multi-dimensional regions, labeled Dist-split and CQR respectively. Throughout the experiments, the miscoverage rate is set at = 0.1 for each prediction region, hence nominal coverage rate of 90%.
Experiments are conducted on four synthetic and six real datasets. For each real dataset, input vari- ables x were standardized before model training. There are many ways to implement CNF, all of which would be compatible with our CONTRA method. Here, we applied RealNVP due to its sim- plicity and computational efficiency. Specifically, we used 6 to 10 coupling layers for each CNE, Each coupling layer (details in Appendix C) involves two neural networks, each including 2 hidden layers with 512 hidden units and ReLU activation function. Optimization is done with the Adam gradient descent method (Kingma & Baj 2014) with learning rate of 10- for most of cases. Training epochs are mostly set to be 200. To implement ResCONTRA, we selected support vector regression as the predictive model.
7

---

Published as a conference paper at ICLR 2025
We also studied the impact of underfitting and overfitting CONTRA regions, well as the impact of data sizes. Practical guidelines are provided for interesting readers in Appendices and
All models were trained varied from 1 to minutes, depending on the datasets and RealNVP structures used.
A30 GPU with 32 GB of memory. The training time for cach CNF
5.1 SYNTHETIC DATA ANALYSIS
We experimented with four setups where YJX are mixture- -aauur--aaasan,, spiral, moon and ring- shaped, respectively. We present results for the first two cases below while leave the latter two in the Appendix [ due to limited space. Data generation details can be found in Appendix The sample sizes of the training, calibration and testing sets of the mixture Gaussian and spiral examples are 3375, 1125, and 500, respectively. We further split the training set for ResCONTRA into 60% for a support vector regression model and 40% for the first calibration to train a normalizing flow, For additional simulation studies, see Appendix D
Table shows the empirical coverage probabilities of each method, averaged over 20 replications, with standard errors the brackets. Each replication is different split between the proper training set and the calibration set, while the test set is fixed. All conformal methods achieved the nominal level of 0.9. In both setups, CONTRA and PCP are based on the same CNF and generated the small- est two prediction regions. ResCONTRA delivered comparable results, with discrepancies largely because the "user-chosen". The other methods produced significantly larger regions.
Table 1: Coverage and volume for 2-dimensional 90% conformal regions in two synthetic datasets. Each table entry is the average of results 20 random splits, with standard error in the parentheses. The method that achieved the smallest volume in boldface.
Metric
CONTRA
ResCONTRA
PCP
NLE
RCP
MCQR
spli.
CQR...
Mixt.
Covenge
Volume
64.92(1.015)
(.910))
90.94(1.727)
0.900.003)
21.32(0.845)
0.9100.005)
129.891 10.299)
0.89(0.005)
131.12(12.378)
(.90/0.0022)
109.51(0.412)
0.89/0.003)
112.3/00.555
0.90y0.002)
112.300..63)
0.900.0022)
Spiral
Coverage
Volume
23.70(1.622)
0.91(0.003)
41.73(2.275)
0.9000.003)
23.22(1.154)
0.91(0.002)
68.97(0.917)
0.91(0.003)
68.16/0.943)
0.9110.002)
64.69 (0.199)
0.91/0.003)
65 (5.821/.157)
0.91/0.003)
65.42/0.202)
0.1(0.0002)
26-
5
25
15
10
110
(a) CONTRA
(b) ResCONTRA
(c) PCP
(d) Ellipses
(c) Rectangles
Figure 2: Prediction regions of two-twweeeee dimensional outcome given specific value from the test set. Colored lines show the boundary various prediction regions. case of PCP, 3000 disks are shown. Orange points show a random sample size 2000 from true conditional distribution given x.
Beyond numerical metrics, we visualize and compare these conformal regions. Figure shows the prediction regions at arbitrarily chosen x value from the test set. CONTRA, ResCONTRA and PCP construct prediction regions that align well with the true HDRs. The remaining four methods are limited generating elliptical rectangular regions, hence contain extensive areas of lower density. Among the top three performers, PCP constructs regions from a union of circles. While increas- ing K improves the approximation the true underlying distribution, also raises computational
8

---

Published as a conference paper at ICLR 2025
costs and creates more complex boundaries. After all, the synthetic data experiments suggest that the proposed CONTRA and ResCONTRA methods outperforms others. They consistently achieve nom- inal coverage, capture true conditional density with flexible shapes, and maintain smooth, connected regions for robust inference.
5.2 REAL DATA ANALYSIS
We next applied CONTRA, ResCONTRA and competing methods to form prediction regions for six datasets from the public-domain. Each dataset was partitioned into 60% training, 20% calibration, and 20% testing. For ResCONTRA, the 60% training portion was further split into 60% for training and 40% for the first calibration. Summaries of their dimensions and sizes are given in Table 2 The first example has a one-dimensional output, and we check if methods for multi-dimensional output work effectively in this special case. Then, we took on examples with two-dimensional output with various input dimensions and training/calibration set sizes. Finally, we tackled a challenging task of density estimation for a dimensional output, a significant hurdle for traditional approaches.
Table 2: Summary of dataset structures in real data experiments.
Dataset
dim(x)
dim(y)
my (training)
n2 (calibration)
ni (lest)
Bio
9
1
3000
1000
1000
Taxi
2
2
3600
1200
1200
Energy
2
460
154
154
2D RF
20
2
5403
1801
1801
SCM20D
61
2
3000
1000
1000
4D RF
20
4
5403
1801
1801
Briefly, Bio focuses on the physicochemical properties of protein tertiary structures, derived from CASP 5-9 experiment (Rana, 2013). It includes nine predictors that provide information the structural and geometric properties of molecules, with the aim of predicting the size of the residue. Taxi contains longitude and latitude details for pick-up and drop-off location in New York for the year 2016. The goal is to predict the (conditional) distribution of the drop-off location given any pick-up location. This is also the example featured in Figure Energy is used train a predictive model for heating and cooling loads given eight muilding-relateed predictors like relative compactness, roof area, surface area and so on (Tsanas Xifara 2012). 2D RF and 4D RF are based on the same river flow dataset, which comprises more than one year of hourly flow observations from eight sites within the Mississippi River network (Spyromitros- Xioufis et al., 2016). Twenty observations of river network flows at various sites and past time points are used to predict the flows 48 hours into the future at 2 and 4 sites, respectively. SCM2OD contains 5000 records from the 2010 Trading Agent
Competition in Supply Chain Management (Spyromitros-iioute Xiouuf et al. 2016). The goal is predict the mean price 20-days into the future.
Table 3: Coverage and volume for multi-dimensione dimensional 90% conformal regions in seven real datasets. Each table
entry is the average of results over 20 random splits, with standard error in the parentheses. The method that achieved the smallest volume is in boldface.
Metric
CONTRA
ResCONTRA
PCP
NLE
RCP
MOQR
Dist-split... um
CQR,..
Bio
Cowerge
Volume
13.23(0.385)
0.9000.002)
12.64(0.178)
0.90)0.000)
12.54/0.072)
0.90(0/003)
13,49/0 132)
(189(0.00)
12.55(0.07)
0.900.002)
Taxi
Coverage
0.89.0.802)
0.8900.000
0.8900.022
0.9000.001)
0.90/0.002)
0.89(0.002)
0.90/0.002)
0.110.002)
Volumet'x 100
8.1(0.290)
9,06/0.251)
8.95(0.132)
10.7700.188)
12.21(0.222)
10.86/0.283)
14.65(0.378)
1244(0.266)
Energy
Coverage
Volume
18.24(1.269 269)
0.8700.006)
22.76(4 238)
0.87:0.000)
16.4000.901)
0.5800.006)
19.14(1.995 (42)1962)
0.8660.006)
26.04(2.888)
0.5500000
25.39(2.235 238)
0.84(0 007)
32,12(2.35 383)
0.84(0.010)
27.73(2.343)
0.8(0.008)
2D RF
Coverage Volume
5.29:0.074)
0.91(0.012)
10.3410.207)
0.9000 002)
0.910.0022
7.19001144
(5.50(0.307)
0.9100-002)
23.85(0.644)
0.900.000
12.14(0.187) 1855
0.91 001)
0.92/00001)
12.30/0.19))
(.22)
SCM20D
Valumet x 10*)
Coverage
0.89:0.001)
6.30x0.120)
0.89/0.000)
8.91(002007)
0.89000022)
6.97(.114)
10.37(0.599)
0.90V0.002)
0.89/0.002)
6.52(0.153)
0.8910.002
(.1000115)
0.90/0.003)
8.82(0.231)
0.9(000002)
8. 8.50002206)
4D RF
Valumet xx 10=)
Coverage
0.90 (0.003)
0.59:0.043)
0.89 (0.002)
1.54(0.098)
0.89 (0.003)
1.13 (0.057)
0.89 (0.002)
1044022)))
52.93(11403)
0.90 (0.002)
10.92(0.597)
0.89(0.002
20.40 (1.219)
0.91(0.001)
12.6800.665)
0.92/0.002)
Table 3 reports the empirical coverage probability and volumn of various prediction regions. The coverage rates generally very accurate around the nominal level 0.9. Slightly lower coverage

---

Published as a conference paper at ICLR 2025
rates that range from 0.84 to 0.88 are observed across the different methods for the energy dataset. This is not surprising given the marginal coverage guarantee established for these conformal methods meant to deliver the nominal coverage when averaged over all possible training and calibration sets. Hence, smaller training and calibration sizes led to a bigger chance of seeing large deviation from the nominal level and relatively low empirical coverage rate out of 20 replications. In addition, for Dist-split pon a reduction in coverage is observed when the dimension of x is high. For Bio, see that CONTRA, ResCONTRA and PCP are as good as state-of-the-art methods designed for one-dimensional problems. In multi-dimense ditensiinnnii experiments, CONTRA, ResCONTRA and PCP constructed considerably more compact prediction regions than their competitors, with CONTRA keeping the smallest volume with most cases. Overall, the CONTRA-type methods excel in these numerical metrics and offer smoother regions than other top performers, providing a significant ad- vantage.
6 SUMMARY AND DISCUSSION
The main contribution of our work is the proposal of a new multi-dimennional conformal prediction method, CONTRA, that allows reliable conditional density estimation. We demonstrated with var- ious examples that CONTRA surpasses other methods in offering compact prediction regions with flexible shapes and smooth boundaries for casy interpretation.
Alongside CONTRA, introduce two new methods, ResCONTRA and MCQR, catering users with their own preferred methods. ResCONTRA extends any user-selected point predictors with valid prediction regions. MCQR extends CQR, a popular one-dimensional conformal prediction method, to create prediction boxes that aim for optimal compactness among box-shaped conformal regions.
Note that most conformal prediction methods are intrinsically one-dimensional: points in the cal- ibration set are projected into line via properly defined one-dimensional conformal score. An empirical quantile of these scores serves as the threshold on the line, and is projected back to the output space to form prediction regions. Any choice of the projection will lead marginal cover- age guarantee under the exchangeability condition. Quality projections are key to generating desir- able prediction regions, which exhibit properties like small volume, minimal disconnected sets, and smooth boundaries that allow robust and interpretable conclusions in weal-world applications. CON- TRA is an intuitive and effective approach for defining such projections. that possess latent representations of the output, and project points to a line based on the density of their latent representation with respect to the standard Gaussian. This reduces to using the distance of the latent representation to the origin in the q-dimensional space as the non-conformity score.
leverages NF and CNF
One challenge of CONTRA is it requires training a bijection that transforms q-dimensional output to its latent representation. Better trained bijection leads to more compact CONTRA regions. In theory, the ground truth bijection always exists (Bogachev et al., 2005). But in practice, more complex distributions with higher q require more data points and increasingly complicated structures in an NF (or CNF) to learn the true bijection well (Durkan et al. 2019). In the case where data size is small and NF (or CNF) is hard to train, one possible solution is to inspect the latent representations of the training and the calibration points and check their distribution against the standard multivariate Gaussian. Patterns of the deviation such bias in certain directions, correlation among different directions, asymmetry/skewness kurtosis and so on could potentially be modeled and utilized to adjust the shape of Ê to maintain 90% empirical coverage and be transformed back to obtain adjusted prediction regions in the output space, all while retaining the coverage guarantee.
The current version of CONTRA uses Gaussian base distribution for the latent representation. This limits its application to model continuous outcomes. Given the wide application of supervised learn- ing with multi-dimensionaa discrete or mixed-type outcomes, there is an urgent need to perform un- certainty quantification in these settings by developing effective conformal prediction methods. Since multi-dimensional classification models as well as one- me--imensioonl
there is wealth of literature conformal methods for discrete outcomes, their confluence has the potential to drive rapid advance- ments in this direction.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Anastasios N Angelopoulos and Stephen Bates, A gentle introduction to conformal prediction and distribution-fiee uncertainty quantification. arXiv preprint arXiv:2107.0751, 2021.
Vladimir Igorevich Bogachev, Aleksandr Viktorovich Kolesnikov, and Kirill Vladimirovich Medvedev. Triangular transformations of measures. Sbornik: Mathematics, 196(3):309, 2005,
Carlo Bonferroni. Teoria statistica delle classi calcolo delle probabilita. Pubblicazioni del Istituto Superiore di Scienze Economiche e Commericiali di Firenze, 8:3-62, 1936.
Ricky TQ Chen, Jens Behrmann, David K Duvenaud, and invertible generative modeling. Advances in Neural Information Processing Systems, 32, 2019.
Jacobsen. Residual flows for
Victor Chernozuukov, Kaspar Wüthrich, and Yinchu Zhu. Distributional conformal prediction. Pro- ceedings of the National Academy of Sciences, 18(48))2100774118, 2021.
Nicolo Colombo, Normalizing flows for conformal regression, 2024. URL hitps://arxiyy
org/abs/ 200p:///400033465533
N. Dalmasso, Pospisil, A.B. Lec, R. Izbicki, P.E. Freeman, and A.I. Malz. Conditional density estimation tools in python and with applications to photometric redshifts and likelihood- cos-
mological inference. Astronomy and Computing, 30: 100362, January 2020. ISSN 2213-1337. doi: 2sccmm.019.10033..
URL http::/dkx doi orw.orggg//...o 1016/1..11111 a.a.aassom....... 2019.
100362
Laurent Dinh, David Krueger, and Yoshua Bengio. Nice: Non-linear independent components esti-
mation. arXiv preprint arXiv:1410.8516, 2014.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. arXiv preprint arXiv: 2016.
Olive Jean Dunn. Multiple comparisons among means. Journal of the American statistical associa- tion, 56(293):52-64, 1961.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Ad- vances in neural information processing systems, 32, 2019.
Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated multiple-output quantile regression with representation learning. Journal of Machine Learning Research, 24(24):1-48, 2023.
Rajarshi Guhaniyogi, Shaan Qamar, and David Dunson. Bayesian conditional density filtering.
arXiv preprint 1401:344 2014.
Peter Hall and Qiwei Yao. Approximating conditional distribution functions using dimension reduction. The Annals of Statistics, 33(3), June 2005. ISSN 0090-5364. doi: 10.1214/
009053604000001282 URL URLpttp:///xx
doi wwt.:xx.oooooo////.o/o http:////...0///0000000000000000
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851, 2020.
Chin-Wei Huang. David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregressive
flows. In International Conference on Machine Learning, pp. 2078-2087. PMLR, 2018
Rafael Izbicki, Gilson Shimizu, and Rafacl Stern. Flexible distribution-fii conditional predictive bands using density estimators. In Silvia Chiappa and Roberto Calandra (cds.), Proceedings of
the Twenty Third International Conference Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 3068-3077, PMLR, 26-28 Aug 2020. URL
wtteee/s/w110//22ock220000 html
Rafael Izbicki, Gilson Shimizu, and Rafael B Stern. Cd-split and hpd-split: Efficient conformal regions in high dimensions. Journal of Machine Learning Research, 23(8)):1-32, 2022.
Munkres James. Topology. Prentic Hall of India Private Limited, New delhi, 2000.
11

---

Published as a conference paper at ICLR 2025
Chancellor Johnstone and Bruce Cox. Conformal uncertainty sets for robust optimization. In Con- formal and Probabilistic Prediction and Applications, pp. 72-90. PMLR, 2021.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization.
CoRR, abs/1442.6980,, 2014. URL https:: sew.sissccccnn...rrr....mm htt.://cccccooa...ooo/oooppppaap 6628106
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1xl convolutions. Advances in neural information processing systems, 31, 2018
Jing Lei, Alessandro Rinaldo, and Larry Wasserman. A conformal prediction approach to explore functional data. Annals of Mathematics and Artificial Intelligence, 74:29-4: 2015.
Jing Lei, Max G'lll,, Alessandro Rinaldo, Ryan Tibshirani, and Larry Wasserman. Distribution- free predictive inference for regression. Journal of the American Statistical Association, 113(523): 1094-1111, 2018.
Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Ellipsoidal conformal inference for multi-target regression. In Conformal and Probabilistic Prediction with Applications, pp. 294- 306. PMLR, 2022.
Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19-23, 2002 Proceedings 13, pp. 345-356. Springer, 2002.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lak- shminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Ma- chine Learming Research, 22(57):1-64, 2021. URL http: // fmmr..r.. wwwpgppppp///////2
1028. html
Prashant Rana. Physicochemical Properties of Protein Tertiary Structure. UCI Machine Learning Repository, 2013. DOI: https///ol.org/113333///000
Yaniv Romano, Evan Patterson, and Emmanuel Candes, Conformalized quantile regression. Ad- vances in neural information processing systems, 32, 2019.
Eleftherios Spyromitros-Xinn Xiouffs, Grigorios Tsoumakas, William Groves, and Ioannis Vlahavas. Multi-target regression via input space expansion: treating targets inputs. Machine Learn- ing, 104(1):55-44 February 2016. ISSN 1573-0565. doi: 10.10077119994--1-55556 URL http: //tx dx doi www.oogg//......... 1001 1ttp:/111111001111111111011 016-5546--
Athanasios Tsanas and Angeliki Xifara. Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools. Energy and Buildings, 49:560- 567, 2012. ISSN 0378-7788. doi: ott.a//c.................o..o//a //www.ccieneedi.eccc
URL https:
Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in a random world, volume 29, Springer, 2005,
Zhendong Wang, Ruijiang Gao, Mingzhang Yin, Mingyuan Zhou, and David M Blei. Probabilistic conformal prediction using conditional random samples. arXiv preprint arXiv:2206.0655 065544 2022.
Christina Winkler, Daniel E. Worrall, Emiel Hoogeboom, and Max Welling. Learning likelihoods with conditional normalizing flows. ArXiv, abs/1912.00042. 2019. URL https: //api.
semant was.annttscchllllo htt.:co........ooooooooooooooo/.
12

---

Published as a conference paper at ICLR 2025
A PROOF
A.1 MARGINAL COVERAGE GUARANTEE OF CONTRA
Proposition (A simple Corollary the marginal coverage theorem of split conformal methods) Suppose the sample points in D2 and (Xu+1, are exchangeable. Then for any CNF model tär- the corresponding conformal ball É based on D2 satisfies
TP (yn-1 tg( (x+1)) 2 21-a.
Further, if 2. E Zoal and bn41 are almost surely distinct, the above probability is bounded above by a + #
Proof.
Since la is trained by CNP, differentiable bijection. Hence
2n+1 E E (g(2n+1, ( xn+1)
- Yu+1 E tgg xn+1).
Therefore
Plyn+1 E tg(E,xnn1)) = P(an+1 EJ 2 1-a
A.2 SET BOUNDARIES REMAIN BOUNDARIES UNDER НOMОMORPHISM
The following result is well-know, but we provide a proof to self-contained. Many background materials can be found in (James) 2000). Let U and V be two topological spaces. (They will be played by and Y respectively CNF.)
Definition 1. functionn
V is called homomorphism is bijective and continuous, and
t- is continuous.
Proposition 3. Ift: :4 + V is homomorphism and E C U, then t(0E) = Ot(E), where OE the boundary of E.
Proof. The boundary of F in 21 can be expressed by OE = E. n E where E. denotes the closure of E and E denotes the complement E.
Since t is homomorphism, we have
((E9) (t(E))
Therefore,
t i
Applying & to the boundary of E. obtain
((EE) - [ B9) (E (8) =((E))
Hence, t(0E) (t(E).
A.3 GUARANTEED COVERAGE PROOF OF MCQR
Proposition Suppose (x;, I...., + 1, are exchangeable, then the prediction region ČMMOOR (x.+1) constructed by the MCOR algorithm satisfies
P yntz CMCCRR (x (xx)))((x,
21-a.
(6)
If the non-conformity scores, s5.==1.... n + 1, are almost surely distinct, then the probability in the left hand side of equation is also bounded above by 1 m2+I.
13

---

Published as a conference paper at ICLR 2025
Proof. First, show that y CMCQR (x) is equivalent to $ S1-a-
$ $1-a
max {1011 (x)(x)) - Ub) u2/2 (1) - (7 s 81-a
j=llnng
(0) - vs) -
$1-a ^ U5 - (x) S
$1-0 for j = 1.....9 w2
Wij
('((x)
$1-a
Sn & (x( + $1-a for j = 1, 9
Wji
Since non-conformity scores 8;, i 12 and Sn+1 exchangeable. Applying Lemma] |
P( I 81 s1-a| : -
calibration residuals &; and sn+1 are almost surely distinct,
P(sn+1 S s1-al(x,, :EE 11) 1
n2 +
Taking expectation over D1, the marginal coverage is guaranteed.
Lemma (Inflation of quantiles, Romano al. 2019) Suppose Zi...., Zut1 are exchangeable random variables. For a E 1),
P S (((atn+1))a) 2 a,
Where [[a((n+]]].)) is [a(n + 1)1 smallest value in Zi...., Lns or a-th empeical quantile of Zi Ze Moreover, if the random variables Zi,..., Zn+1 are almost surely distinct, then also
P(Zn+1 S Z ([o(++)]] s sa n
B SYNTHETIC DATA STRUCTURE
A model with mixture Gaussian error term.
Yi = 3xXX, 5x3 4X, - 6X2 + +7+
Y22 X,X, X3 3X, 8
where = [Y, Y2 and & [y
[X, X21 ~ N), # = [-2.0, -1.5"" 0.3N (10.0)" 0.I + J2)) 0.4N (15, 5)", 1.5(I2 J2)) + 0.3N (110,07 0j 12). J2 is a 2 by 2 matrix with all elements equal to one.
2. A model with spiral curve error term
Yi 2x 3X3 +5, X1X = X2X, 4X3 + 3 +7
where X has the same structure as in Model 1; E1 ((0c (0 sin( 0.12), where 0 (0,2 2T).
3. A model with a moon curve error term. The specification of Y and X are the same as the previous setup, And the error term follows a moon-shaped distribution, Ei ~ N (c)s(()) 0.13), 82 N (sin(0), 0.1²), 0 (T.
4. A model with a ring error term. The specification of Y and X are the same as the previous setup, ε cos( 82 (0), where 22 U( 2 O U (0,,2T))
14

---

Published as a conference paper at ICLR 2025
C COUPLING LAYER
For vector E 189, partition into two subspaces: (V1,) E R x R and a bijection function g(-) 129-91 91 Define
211 = Yh
212 g(y12 (y,,, x))
In particular, let g(yi;; m(,,, Y1, O (yp, x)) e(y,, x) (Winkler et 2019), then the determinant of Jacobian is exp( 91 j= :1y,,,))))
The transformation é composed by multiple coupling layers. To prevent the composition of two consecutive coupling layers from reducing to the identity function, we can switch the roles of the two subsets (Dinh et al., (2166). For instance, consider scenario where aim to transform y into z via w. In this case, we have:
Z1, =g ((0(w,,m(w)g,)) 212 W/,
! Wh yn
WI, = g (yyimmyy,
D ADDITIONAL SIMULATION STUDIES
Table Coverage volume 2-dimensional 90% conformal regions two synthetic datasets. Fach table entry is the average of results over 20 random splits, with standard error in the parentheses. The method that achieved the smallest volume is in boldface.
Metrie
CONTRA
ResCONTRA
PCP
NLB
RCP
MCQR
Dist-split
CQR
Coverage
0.91(0.002)
0.9000.003)
0.91(0.003)
0.91(0.003)
0.91(0.003)
0.90/0.003)
0.91(0.003)
0,91(0.033)
Volume
1.56(0.020)
2.65(0.062)
1.66/0.015)
3.32(0.040)
3.28(0.037)
2.50(0.014)
2.63(0.012))
2.62/0.013)
Ring
Coverage
0.90(0.004)
0.91(0.002)
0.90(0.002)
0.91(0.008)
0.91(0.003)
0,89(0.004)
0.89c0.003)
0.89(0.003)
(x10P)
4.57(0.405)
1.97000300)
6.14(1.537)
5.17(0.009)
5.19(0.010)
5.19/0.009)
300
25
-7
10
-10
-20
-10
-20
-20
-2B
-20
-10
-BC
-22
-00
-10
(a) CONTRA
(b) ResCONTRA
(c) PCP
(d) Ellipse
(c) Rectangle
Figure 3: Prediction regions of a two-dimensional outcome given a specific x value from the test set. Colored lines show the boundary of various prediction regions. In the case of PCP, K 3000 disks are shown. Orange points show random sample of size 2000 from the true conditional distribution of y given x.
E THE IMPACT UNDERFITTING AND OVERFITTING ON CONTRA
Using Normalizing Flow model that is either too complex too simple result in conformal regions that are overly sensitive or excessively large. A straightforward guideline to mitigate these overfitting and underfitting issues is to evaluate how closely the latent variables from the calibration
15

---

Published as a conference paper at ICLR 2025
(a) Over-dispersed
(h) Under-dispersed
(c) Nearly Gaussian
(d) Underfitting
(e) Overfitting
(f) Regular fitting
Figure 4: CONTRA prediction regions and latent for the calibration set under three different NF models. for drop-off location given specific pickup coordinate (blue pin) for the NYC taxi data. Data size equals 4800, with 75%-25% training-calibration split. and (d): underfitting NF with 2 coupling layers and hidden units per layer, trained for 50 epochs 4 seconds. (b) and (c): An overfitting NF with 16 coupling layers and 1024 hidden units, trained for 500 epochs in 205 seconds. (c) and (f): A regular fitting NF with coupling layers and 256 hidden units, trained for 200 epochs in 17 seconds.
set set resemble a random sample from the standard Gaussian. The top row of Figure provides a visual check for bias, overdispersion, corresponding to the rightmost picture had the best out-of-sample performance among the three, and would be the best choice for constructing CONTRA regions. This confirmed in the bottom row of Figure 4 In addition to the visual check, classical metrics can be used to quantify the deviation between the sample and the standard Gaussian, aiding model tuning.
underdispersion from which users tell that the trained model
F IMPACT OF DATA SIZE ON CONTRA AND ITS MAIN COMPETITOR
(a) (x11 m2) (225, 75)
(b) (vi. m2) (900,300)
(c) (m1, v2)) (3600, 1200)
Figure 5: The impact of data size prediction regions of drop-off off location given specific pickup coordinate (blue pin) for the NYC taxi data. The top row shows results of our proposed CONTRA method; the bottom row shows results ST-DQR based on a diffusion model. Data sizes used for each column are 300, 1200, and 4800, respectively, with a 75%-25% caliing-calibration split.
We examined the performance of CONTRA and the ST-DQR based on the diffusion model across various sample sizes, with a particular focus on smaller sample sizes that pose challenges for both NF and diffusion model training. The resulting prediction regions for both methods were consistent with
16

---

Published as a conference paper at ICLR 2025
expectations and appeared reasonable. As the sample size increased from very small to moderate, uncertainty decreased, leading to smaller prediction regions. However, as the sample size continued to grow, the prediction region stabilized, remaining approximately the same size, which correctly reflects the inherent uncertainty in predicting outcomes for new subjects.
G A COMPARISON BETWEEN CONTRA AND RESCONTRA
This section compares the two proposed methods, CONTRA and ResCONTRA. We use empirical examples to confirm the following heuristics: When NF is effective learning y/x, the one-step CONTRA method generally outperforms the two-step ResCONTRA due to its larger training set. Whereas when the relationship between y and x, such as E(ylx) f(x), is highly complex but the residual y f(x follows a relatively simple distribution, ResCONTRA tends to provide better prediction regions. This is because ResCONTRA divides the original training data into two subsets. One subset is used to train an estimator for f(x), using any preferred method that specializes for this task of point estimation. The second subset then used model the distribution of the residual using a relatively simple NF.
Above, we provided a general guideline for choosing between CONTRA and ResCONTRA, In prac- tice, users don't need to decide in advance. They can try both methods and use the method that lead to the prediction region that better suits their needs.
3.87
226
(a) CONTRA
(b) Res- SVR-10
(c) SVR-
(d) Res -Gs-X
(e) Res-XGB-6
Figure 6: Comparing CONTRA to ResCONTRA in terms of the latent variable, z, of the calibration set on two different datasets. Plots of that closely resemble the standard bivariate Gaussian distribution result in smaller calibrated radii, r.9, and better conformal prediction regions. Five conformal prediction methods are implemented: (a) CONTRA with a 10-layer NF; (b) and (c) ResCONTRA with SVR followed by 10- and 6- layer NFs, respectively; (d) and (e) ResCONTRA with XGBoost followed by 10- and 6-layer NFs, respectively. For CONTRA, the sizes for training, calibration, and testing are (3375, 1125, 500). For ResCONTRA, the training set further split into 60% to train point estimator and 40% to train NE.
CONTRA
SVR-10
SVR-6
XGB-10
XGB-6
Mixture-Gaussian (%)
0.71
2.67
1.69
2.67
3.38
Multiplicative Gaussian (%)
1.24
0.27
0.09
0.00
0.00
Table 5: The percentage of points out of the [-7.5, 7.5] x [-7.5, 7.5) range.
For illustration, we apply CONTRA and four different implmentations of ResCONTRA to two exam- ples. The first example the same Kexture-Gaussiar example in Section 5.1 And the second exam- ple specifies more complex relationship between y and x, described in equation 7 For each method that learns the conditional distribution of y given x, we examine values of the latent variable, z, of the calibration set. Resemblance of the distribution of z to the bivariate standard normal indicates good sample learning of the distribution of y xx. Figure 6 displays digupayss the values of z for the two exam-
ples in the top and bottom rows, respectively. All 10 plots show the [-7.5, 7.5] 1-7.5, 7.5] region of the latent space, while 5 the proportions of z falling outside the shown range. Over- laying on each plot is the theoretical 90% highest density region of the bivariate standard Gaussin
17

---

Published as a conference paper at ICLR 2025
(in green) and the smallest circle centered at the origin that contains 90% of the of the calibration set, with radius denoted by 7. From the first row of Figure 6 and Table we can see that CONTRA (a) is the best performer for the example with relatively simple model and more complicated error; and from the second row, we can see that ResCONTRA methods, especially (c) and (e) that combine tools specialized for point estimation and relatively simple NF for the residuals, the best perform- ers to capture the complex relationship between y and x while providing small r that leads to small conformal regions.
Details of the relatively complex dataset.
This is the dataset that corresponds the results in the bottom row of Figure 6 and Table 5
Yi eeee 62 3X2 + 0.5X X4X - 15xg + 0.7x, 0.3 + sin (11) + 5
{ 1/22 - x1 + 4x3 XaX4e2 + 2Xc.Xre1e + 0.6Xs 1.2x.ef + cos(X10) +
(7)
where Y ,Y, Y2],, X~ N(,, μ is 10-dimensional mean vector with each component independently drawn from the Uniform[- 10, 10] distribution, and e [21, ~ N (0, I2).
H AGREEMENT BETWEEN CONTRA PREDICTION REGIONS AND HPDS
(a) Mixt. Gaussian
(b) Moon
(c) Gamma
(d) Ring
(e) Spiral
Figure 7: CONTRA conformal regions of the output variable y given some fixed x value for five setups with coverage levels 50%, 70% and 90%, respectively. The orange points are random samples size 2000 from each of the true conditional distribution of y given x. We can see the CONTRA regions of various levels properly capture the high-density regions in each case.
18


Reviews:
[{"review_id": "0ruxkbxzfI", "reviewer": "Reviewer_vtsj", "summary": "The authors propose a new conformal prediction method for multi-dimensional regression. It is based on fitting a (conditional) normalizing flow in the space of outputs. Authors perform numerical experiments on public datasets that illustrate certain benefits of the proposed method.", "questions": "- Why don't you study other density estimation methods? \n\n- Performance of normalizing flows heavily depends on the size of the dataset. How will it influence the results of the proposed method? I think that this should be studied.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The proposed method makes a lot of sense as flexible estimator of conditional label distribution is obviously required to achieve optimal set size, especially in multidimensional case. \n- The results are adequately explained.\n- Experimental results show certain benefits of the proposed approach in terms of prediction set size.", "weaknesses": "- The proposed approach boils down to application of the normalizing flows to the split conformal prediction framework. There are no technical difficulties on this way. Any other conditional density estimator can be used here instead.\n\n- All the theorems in the paper are standard (not sure, why even proofs are provided). There seems to be no theoretical contribution for the present paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "FwDkQC4yIO", "reviewer": "Reviewer_e6SY", "summary": "This work addresses the limitation of most of the current conformal prediction algorithms which struggle with multi-dimensional outputs. The proposed method CONTRA, uses conditonal normalizing flows (CNFs) to learn a bijection, mapping a latent variable form a simple base distribution to the output variable, and then mapping the high density regions of the base distribution to the output. Using Conformal Prediction ( particularly the split conformal prediction algorithm) , CONTRA then calibrates the results to provide marignal coverage gaurantees of CP on including the true prediction region.", "questions": "1. with respect to point 2 in the weaknesses, how do competing methods, particularly PCP, are expected to perform in cases that CONTRA is not able to perform well. A discussion ~ potentially a brief experiment on synthetic data ~ where we understand the trend of performance deterioration in CONTRA, and the effectiveness of the solutions proposed on page 10,  could illicit a more clear choice to follow in whether one should try to apply the proposed solutions or go with another method.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The writing is clear, and to the point. The contributions are presented honestly and placed well with respect to the prior work. \n2. The experiments are elaborate, and show the main metrics in CP: efficiency in terms of the size ( volume ) of the prediction sets and coverage. Both synthetic and real world datasets are considered, and the results are consistent across the two settings ( synthetic vs real-world). \n3. the output regions of CONTRA do seem to visually support the authors' claims in that they are interpretable and have smoother boundaries compared to the related work. In most settings/datasets CONTRA seems to be outperforming most other baselines ( except PCP ). \n4. They address an important and practical problem (multi-modal outputs ) which extends the applicability of Conformal Prediction. \n5. Their methods seems simple to follow and implement.\n6. They introduce variation of the method, to adopt to different choices of the fitted model for the conditional density estimation.", "weaknesses": "1. PCP ~ ( a competing method) is able to outperform both variations of CONTRA and ResCONTRA in some cases in terms of providing smaller prediction regions while maintaining coverage guarantees.\n2. the authors mention instances in which CONTRA could deteriorate in performance. Particularly when its hard to learn bijection well. but no results are provided to illicit the deterioration trends or to show the remedy solutions discussed by the authors can help.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "KVixGLdIRP", "reviewer": "Reviewer_rPVJ", "summary": "CONTRA is a novel method for multi-dimensional output prediction that utilizes normalizing flows to define conformal prediction regions. By mapping high-density regions in latent space to output space, it generates more precise and connected prediction regions, avoiding the oversized or irregular regions produced by traditional methods. The paper also introduces ResCONTRA, an extension that can enhance any predictive model with reliable prediction regions. Experiments demonstrate that CONTRA not only guarantees the desired coverage but also produces prediction regions with smaller volumes and smoother boundaries, outperforming existing methods on both synthetic and real-world datasets.", "questions": "None", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The manuscript is generally easy to follow and includes a sufficiently comprehensive discussion of relevant prior works. Experimental results are presented well.  The proposed solution has solid theoretical foundations.", "weaknesses": "- In line 186, how to construct the prediction region? \\hat{E} is a q-dim ball, which contains infinitely many vectors, making its implementation unclear.\n- More discussion about CNF will be valuable.  \n        1. Providing details about choosing q_1 in Appendix D will be better. Equal division? Moreover, for Bio dataset, the dim-y is only 1, that can not be split.   \n         2. A possible limitation is that normalizing flows is notoriously difficult to train in high-dimensional spaces. Extending CONTRA to handle high-dimensional outputs (e.g., image generation) would significantly enhance its practicability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "fmBxQsVbPm", "reviewer": "Reviewer_bEdm", "summary": "The authors propose a technique to compute valid prediction regions in multiple-output regression tasks. The approach is based on learning a Normalizing Flow that transforms the output conditional probability into a 'simpler' base distribution. The invertibility of the flow allows mapping valid convex intervals of the feature space, e.g. rectangles, to non-convex smooth regions in the original space.", "questions": "- What is the difference between the proposed approach and methods that learn a feature space map before applying CP, e.g. STDQR?\n- What is the difference between using a normalizing flow as a prediction model and a model for the residuals? From Figure 1, the latter seems to produce smoother regions. Is this always the case?\n- Can the method be extended to tasks with many more outputs, e.g. images?\n- In practice, requiring the map to be invertible may make Normalizing Flow hard to train or less flexible than unconstrained models. Have you tried Normalizing Flow with different structures?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Using Normalizing Flow to improve Conformal Prediction is intriguing and flexible. \n- Normalizing Flows and their likes are increasingly popular for training generative models. The proposed method could be further improved by leveraging advances in that domain.\n- Multiple-output Conformal Prediction is an important but often overlooked problem.", "weaknesses": "- The use of Normalizing Flow is not theoretically justified. The introduction should explain the differences compared to learning more standard feature maps.\n- The authors should describe in detail the overlaps between their work and Colombo 2024, which is not mentioned with other related work but seems closely related to the proposed approach.\n- It is unclear why the proposed approach is superior to existing multiple-output Conformal Prediction approaches.\n- The authors should give more details about the invertible flow, $t_\\theta$, and discuss possible overfitting or underfitting problems. Especially in comparison with other non-invertible probability maps (e.g.diffusion models).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)