Multivariate Conformal Selection
Authors: Tian Bai
Yue Zhao
Xiang Yu
Archer Y. Yang
Venue: ICML 2025 poster (2025)
Keywords: Uncertainty Quantification
Machine Learning
Model-free Selective Inference
Conformal Inference
Multivariate Responses
Primary area: general_machine_learning->supervised_learning
TL;DR: We propose Multivariate Conformal Selection (mCS), a framework extending Conformal Selection to multivariate settings, achieving finite-sample False Discovery Rate control and enhanced selection power with applications in drug discovery and beyond.

Abstract: Selecting high-quality candidates from large datasets is critical in applications such as drug discovery, precision medicine, and alignment of large language models (LLMs). While Conformal Selection (CS) provides rigorous uncertainty quantification, it is limited to univariate responses and scalar criteria. To address this, we propose Multivariate Conformal Selection (mCS), a generalization of CS designed for multivariate response settings. Our method introduces regional monotonicity and employs multivariate nonconformity scores to construct conformal $p$-values, enabling finite-sample False Discovery Rate (FDR) control. We present two variants: $\texttt{mCS-dist}$, using distance-based scores, and $\texttt{mCS-learn}$, which learns optimal scores via differentiable optimization. Experiments on simulated and real-world datasets demonstrate that mCS significantly improves selection power while maintaining FDR control, establishing it as a robust framework for multivariate selection tasks.

Multivariate Conformal Selection
Tian Bai
Yue Zhao²
Xiang Yu? Archer Y. Yang
Abstract
In these settings, true test responses (e.g., binding affinity or alignment score) are often unavailable, requiring selec- tion to rely machine learning model predictions. Since selected targets drive downstream decisions, quantifying selection uncertainty is essential for maintaining efficiency. Conformal selection (Jin & Candès, 2023) provides a model- agnostic framework for selection with uncertainty quantifica- tion by extending conformal prediction ( (Vovk et al., 2005) to multiple hypothesis testing, using conformal p-values (Bates et al., 2023) and multiple testing corrections. It has shown promise in real-world drug discovery (Bai et al., 2024) and LLM alignment (Gui et al., 2024).
Selecting high-quality candidates from large datasets is critical in applications such as drug discovery, precision medicine, and alignment of large language models (LLMs). While Confor- mal Selection (CS) provides rigorous uncertainty quantification, it is limited to univariate responses and scalar criteria. To address this issue, we pro- pose Multivariate Conformal Selection (mCS), generalization of CS designed for multivari- ate response settings. Our method introduces regional monotonicity and employs multivariate nonconformity scores to construct conformal p- values, enabling finite-sample False Discovery Rate (FDR) control. We present two variants: mCS-dist, using distance-based scores, and mCS learn which learns optimal scores via differentiable optimization. Experiments on sim- ulated and real-world datasets demonstrate that mCS significantly improves selection power while maintaining FDR control, establishing it as a ro- bust framework for multivariate selection tasks.
However, existing conformal selection methods are limited to univariate responses and selection targets of the form y > c, where e a user-defined threshold. Many real-world applications require selection based multiple interdepen- dent criteria. For instance, LLM outputs must simultane- ously satisfy alignment requirements such faimesss, safety, and correctness (Bai et al., 2022), which are better repre- sented as multivariate alignment scores. This highlights the need for a principled selection method with uncertainty quantification multivariate settings.
In this paper, we extend the conformal selection framework to multivariate response settings. To ensure sampe false discovery rate (FDR) control in our procedure, we gen- eralize the concept of monotonicity (Jin & Candès, 2023) to regional monotonicity for the nonconformity function. We propose two types of nonconformity scores that satisfy this property: (i) distance-based nonconformity scores for regular-shaped and convex target regions, and (ii) learning- based method for optimizing nonconformity scores. The latter approach leverages a loss function that penalizes ei- ther the smooth selection size or the conformal p-value to learn an optimal score. This method is particularly effec- tive when the dimension of responses is large or when the target region is irregular or nonconvex. Through experi- ments on simulated and real-world datasets, both variants of mCS demonstrates enhanced selection power over baseline methods while ensuring finite-sample FDR control.
1. Introduction
Selecting subset of promising candidates from a large pool is crucial across various scientific and real-world applica- tions. In drug discovery, researchers search vast chemical spaces to identify compounds with strong effects, such as high binding affinity to specific target (Szymański al., 2011; Scannell al., 2022; Sheridan et al., 2015; Zhang al., 2025). Similarly, precision medicine aims to identify positive individual treatment effects (Lei & Candès, 2021), and post-hoc certification of large language model (LLM) outputs seeks to retain only trustworthy generations that meet user-defined criteria (Gui et al., 2024).
'Department of Mathematics and Statistics, McGill Univer- sity. Montreal, Canada 'Department of Mathematics. University of York, York, UK MRRL, Merck & Co., Inc., Rahway, NJ, USA "Mila- Quebee Al Institute, Montreal, Quebec, Canada. Corre- spondence to: Archer Y. Yang <archer megill.ca>.
2. Background and Related Work
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
Problem Setup We let x E RP represent the p- dimensional features, and let y E Rd denote the d-

---

Multivariate Conformal Selection
dimensional multivariate response variables. We consider training dataset urain { yi yi and test dataset 1 Diest {In+i)j=1 where the corresponding test responses 1--+ are unobserved. We further assume that the combined set of training and test samples (æi, y=li=1 1a+m drawn i..d. from an unknown, arbitrary distribution Dxxy.
with selection tasks. Specifically, the potentially complex shapes of the multivariate CP sets C1-a may be incompati- ble with the pre-defined target region R.
Even in the simpler cases, such as when the response is binary, using CP for selection introduces a multiplicity issue (Jin & Candès, 2023). In this context, CP only controls the per-comparison error rate (PCER), which differs from the FDR. PCER also measures the Type-I error, and is defined as (1) with the denominator replaced by m, the size of the test dataset. By definition, the PCER is always smaller than the FDR, making it a less stringent error control criterion. As a result, procedures controlling PCER may fail to meet the stricter requirements of FDR control.
We formulate the selection problem as follows: Given a predefined d-dimensional closed region R S RR, our goal to identify a subset of indices _ { m) from Drest such that as many test observations j S satisfy Dn+j E R possible. while controlling the FDR (Benjamini & Hochberg, 1995) below a user-specified level q. The FDR is defined as the expected proportion of false discoveries (j E S but Un+s $ R) among all selected observations
Conformal Selection Conformal selection (CS) (Jin & Candès, 2023) is model-aggnstic selection framework that guarantees finite-sample FDR control. However, CS only considers the univariate response case (d = 1) and assumes that the selection region takes the form [c, +oo), where c is predefined threshold. In this setting, CS formulates one hypothesis test per candidate:
|S n Hol [S]
FDR = E
(1)
where Ho [j Un+j & R). with the convention that 0/0 in the fraction above. This criterion measures the overall Type-I error rate of the selection procedure.
The overall Type-II error of selection can be quantified by the power, defined as the expected proportion of desirable observations (Un+) E R)) that are correctly selected,
Hoj : yn+s Se vs. Hig : yo+s > c.
Rejecting the null hypothesis Haj indicates that the j-th test sample is selected, as its response is deemed to exceed the threshold c.
[S I H4l 1711
Power =
(2)
CS uses nonconformity scores to guide its selection process. A nonconformity measure quantifies how atypical (or non- conforming) an observation is, based on the relationship between inputs and responses. For calibration samples, the nonconformity scores V; V( y:). and for test sam- ples, they are Vn+j (a c), where e replaces the un- observed Un+j. These scores are then used compute con- formal p-values through a rank-based comparison of Va+s against the calibration scores Vi,..., Vn- A lower rank of Vn+, relative to the calibration scores provides stronger evidence for rejecting Hoj.
where H4 = {j Untj E R3. The Type-II error of selec- tion is therefore (1 Power). An ideal selection procedure should aim to maximize the power while keeping the FDR below the specified nominal level.
Conformal Prediction Conformal prediction (CP) (Vovk et al., 2005) is a popular framework for uncertainty quantifi- cation that constructs prediction intervals on a per-sample basis. Assuming exchangeable calibration and test data, CP provides prediction sets C1-a(2)) (x) with sample
coverage guarantees for a E (0, 1):
To determine the final selected subset S, CS applies the Benjamini-Hochberg (BH) procedure (Benjamini & Hochberg, 1995; 1997), widely used method for control- ling FDR in multiple testing setting, to the set of conformal p-values. The use of the BH procedure ensures that the overall Type-I error rate kept below the specified level.
P2 E (z)) 2 a.
Although CP was originally designed for univariate re- sponses, numerous studies have proposed multivariate gen- eralizations (Kuleshov et al., 2018; Bates et al., 2021; Mes- soudi et al., 2022: Johnstone & Cox, 2021; Feldman et al., 2023; Park et al., 2024; Klein et al., 2025). However, these multivariate CP methods are not directly applicable to our selection problem. The primary objective of CP construct- ing confidence sets for predictions- - does not naturally align
3. Multivariate Conformal Selection
In this section, we introduce the key concepts, procedures, and theoretical foundations of multivariate Conformal Se- lection (mCS).
Laterr, we will relax the i.i.d. assumption exchangeability conditions.
First, for each d-dimensional multivariate response Untjz
2

---

Multivariate Conformal Selection
mCS performs the following hypothesis test:
(practical) conformal p-values:
Pj py= Di= 1( Putst + U Σ 1(V PutaB) })
Hoj: Dn+s E R vs. Hiz Bn+s E R
n ++
(4)
where R represents arbitrary pre-defined closed target region in R.. Multivariate responses yn++ can represent ei- ther regression classification outcomes. In this paper, we focus on the more challenging regression setting. For dis- cussion of the classification case, please see Appendix B.1. The mCS process consists of three main steps:
Standard results from conformal inference ensure that the oracle conformal p-values pž follow the Unif(0, 1) distribu- tion, in particular implying that they are conservative in the sense that pj as random variable has super-uniform dis- tribution on [0, 1]. satisfying P(rj S a) S a (Bates et al., 2023). To guarantee that the practical conformal p-values Pj also maintain this conservativeness, the nonconformity function V must satisfy a property that we introduce as regional monotonicity. With this property the p 's in the resulting mCS procedure will ensure the control of the FDR. We formally define the regional monotonicity as follows: Definition 3.1 (Regional Monotonicity). A nonconformity score V : X x y/ - R satisfy the regional monotone property if V(x, y²) < V(x, y) for any x E EX, y E R and y E R.
1. Training: Construct a multivariate predictive model j for y. This model can be obtained using any suitable machine learning algorithm.
2. Calibration: Build regionally monotone multivari- ate nonconformity function based on i., and evaluate this function on the calibration dataset and test dataset. Subsequently, we compute the conformal p-values for each test sample.
3. Thresholding: Apply the Benjamini Benjanini-Hochberg (BH) procedure as in the original CS procedure to the set of conformal p-values, yielding the final selection set S.
Definition 3.1 leads to the conservativeness of Pj- The fol- lowing proposition formalizes this result, with proof avail- able in Appendix A.1.
Both the training step and the calibration step rely on the labeled dataset Durin- In the case where model j is already available, mCS can be directly applied using all labeled data for calibration. Otherwise, the training data is divided into two subsets: one for model training (the proper training dataset) and the other for calibration. For simplicity, in the following discussion,
Proposition 3.2. Given that the calibration data { yi) together with the j-th data test data point (In+j Un+j) are exchangeable for j E (1....,m}. re- gionally monotone nonconformity scores V ensures that the conformal p-value pj defined in (4) is conservative in the following sense,
assume that j is already available and all training data {a;;, J are used for calibration so Dtrain-
P(pj S a (rand) j E Ho) S a, for all E E0, 1).
(5)
that Deal=
The conformal p-values (Bates et al., 2023) are used to per- form the hypothesis tests. the true responses {ym++) 1 were observed, the oracle conformal p-value would be de-
Remark 3.3 (Clarification on conservativeness) The con- servativeness described in (5) differs from the conventional notion of statistical conservativeness, which is not condi- tional on the event j Ho. Due the inherent randomness in リж+s within our hypothesis tests, an unknown depen- dency exists between the p-values and the event j E Hoo Consequently, the standard form of conservativeness does not hold in this context.
fined as
P E <
Vnng
+U,(+ 1(V,
Vo+g))
n ++1 1
(3)
Remark 3. 3.4 (Univariate monotonicity as special case). Re- gional monotonicity generalizes the univariate monotonicity concept introduced in CS (Jin & Candès, 2023), which was originally defined for nonconformity scores increasing in their second argument. The original definition restricted to the univariate case, where the target region R. is spec- ified as (c, +20). For a univariate nonconformity score V X' x ) R. monotonicity implies regional monotonic- ity. Specifically, for any y E R = (-oc, c] and y' E Rt it holds that V(x, v) S V( 3/) for all x E X. More- over, even in the univariate case, the original definition of monotonicity across the entire domain X' is unnecessary and monotonicity across the regions R and R suffices.
where V; = V(x,, v,) for i 1,..., + and V is multivariate nonconformity function based on A. and Uy- ~ (n, 1) is an independent random variable for the tie- breaking of the nonconformity scores. We defer the specific choice of to later sections.
The evaluation of the oracle p-value p is infeasible, be- cause that in the above definition (3), the computation of Vn+j V (Zn+jj )n+t)) requires knowledge of the unob- served response Yntj- To address this issue, we replace Vn+, with Vnts = (I Tn+s), where rnts is an ar- bitrarily chosen point in the target region R, yielding the
3

---

Multivariate Conformal Selection
Once a regionally monotone multivariate nonconformity score is defined, conformal p-values can be computed us- ing (4), Leveraging these conformal p-values, mCS again applies the Henjamini-HHochee procedure (Benjamini & Hochberg, 1995) to construct a selection set S. The com- plete approach is outlined in Algorithm I.
nonconformity score are follows:
(regull)r) Di(z, = D2(z, R²) = inf
=| - 8 sllp.
sER
(7)
(cclipped)
Di(2, R²) = M 1 E R OR),
D2(2, R²)= inf aERe
|| - 8 pp:
(8)
Algorithm 1 mCS: Multivariate Conformal Selection
where M is a large constant that serves as a relaxation of infinity, approximately satisfying M 2D2(((x), R*). We discuss the role of M in the following paragraphs. These scores generalize the signed error score and the clipped score (Jin & Candès, 2023), respectively.
Input: Calibration data Deat = {( H.) test data Dieu = {I+++) 1m jj=1- target target region R. FDR nom- inal level q E (0, 1), regionally monotone nonconfor-
mity score V X x y/ + R.
Output: Selection set S.
1: Compute Vú = (z,, v,) for i = 1..... n. and Vats = V( r5) for j = 1,..., m with rn+j E R. 2: Construct conformal p- v-luuu } as in (4).
This formulation enables the decomposition of V( y) into two terms. The first term D1 inherently ensures the regional monotonicity of the score V (Definition 3.1), as Dj is distance function satisfying 0 = Di(y', R*) & (y. R)) for any y' E Re and E R. Meanwhile, the second term D2 measures the distance between the predicted responses of the points and R this term is designed to increase as ((æ) moves away from RE, ensuring that the test points with the predicted responses having large distance from R will yield smaller test scores Vn+j. According to (4), these points will have smaller conformal p-values and are more likely to be rejected (selected) by the BH procedure. Note that when the predictive model f outputs an estimated condi- tional distribution P(y) Ø)- - as in classification, conditional density estimation, or Bayesian models - the second term D2 can be replaced by the predicted probability of being in the target region: y/ E R. i.e. f E (t)dP(y)a)) (z). This serves the same purpose: points with high predicted probability of satisfying the selection criterion will receive lower scores and are more likely to be selected.
3: (BH procedure) Compute the BH selection threshold
k* max(k E Z200 ΣØ= 1 (p, S qk/m) 2 k).
and return the selection set as S = [j : Pj S qk" /m).
The following theorem shows that Algorithm controls FDR in finite sample. The proof can be found in Ap- pendix A.2.
Theorem 3.5. Suppose V is a regionally monotone non- conformity score, and for any j E (1....,m), m), the random variables Vi....,Vo, Vn, Vn+j are exchangeable conditioned on (Vu+elezj. Then, for any E (0, 1), the output S of mCS satisfies FDR < g.
4. Choices of Nonconformity Score
While the previous sections demonstrate that FDR- controlled selection can be performed using any regionally monotone nonconformity score, the selection power heav- ily depends on the quality of the chosen score. Although related studies have explored this issue in the context of CP (Romano et al., 2019; Kivaranovic al., 2020; Sesia & Candès, 2020), there has been limited focus on the choice of scores specific for CS.
We note that computing Vn+s = (( rwty) in (4) requires choosing a point Tn+j E R to ensure FDR con- trol. Although any point in R works, selecting rn+s on the boundary ÖR is optimal for maximizing selection power for both regular and clipped scores. In this case, for any fixed En+jj, the test score Vn+j achieves its minimum, ensuring it to be smaller than a larger proportion of calibration scores Vi = V( V:). For example, with the clipped score, if Tn-j E OR, then the first term Di becomes 0 in Vn+j. while for any calibration samples with y, E R, we have Dj = MM (except when y, lies exactly on ØR, which occurs with zero probability.) This yields smaller p- p-values p, and enables more samples to be selected. We require M to be large for the same reason.
In this section, we introduce two types of nonconformity score that satisfies the regional monotonicity. As a result of Theorem 3.5, applying Algorithm with the purposed scores would guarantee FDR control.
4.1. mCs Distance- based Scores
For multivariate selection, we propose distance- listance-based non- conformity scores of the following form:
The following asymptotic analysis provides theoretical ba- sis for preferring the second, clipped score (8) over the first score (7) when choosing a distance-based nonconformity score. This result extends the original CS framework (Jin & Candès, 2023) to the multivariate setting, as formalized in the following theorem:
V(a, y) (y R (2(())) R
(6)
where ğJ is a trained predictive model of y. and D1 and D2 are distance functions. Two examples of distance-based
4

---

Multivariate Conformal Selection
Theorem 4.1. Let V be any fixed regionally monolone nonconformity score, and suppose ((r,, De) p= are ex- changeable from distribution xxxY. Let (xt, y) denote random pair also drawn from D DxxY, and define F(V, u) P(V(2, y) < <) + u. P(V y) = for any U E I and u E [0. 1]. Assuming the choice of Tutj = r E R is fixed, define
based on maximizing the realized FDR, is provided in Ap- pendix B.2, along with further discussions Theorem 4.1. For empirical comparison of the performance of the two scores (7) and (8), refer to Appendix C.2.1.
4.2, mCS-learn: Learning-based Nonconformity Scores
t- sup { [0. 1]: P(F(V(Z, T), S $q). (9)
The two distance- based nonconformity scores introduced in the previous section offer straightforward and practical solutions for many scenarios. However, their effectiveness, particularly in the design of the second distance term D2(-), is limited in some cases. For example, our numerical simu- lations indicate that mCS- st. would only achieve subop- timal power when R is nonconvex set; see Appendix C.2.3 for further details. Furthermore, when the target region R is irregular, constructing a closed-form distance function can be challenging, leading to higher computational costs and potential inaccuracies.
Suppose that, for any sufficiently small e >> 0, there exists some E (L" -e, ()) such that P(P(V(rr))0)3) S q. Then the output S of Algorithm 1 from input (();, v,)), U (En+isye satisfies
lim FDR (()( R*) and
lim Power = (x) & yy, (10) 8,
P(y R)
To address these challenges, we propose an alternative method mCS-learn, which leverages loss function that penalizes either the smooth selection size or conformal p- value to learn an optimal nonconformity score within the following family:
We omit the proof of Theorem 4.1, it closely mirrors Proposition 7 in Jin & Candès (2023). An intuitive planation can be found in Appendix B.2. Leveraging a characterization of the BH procedure (Storey et al., 2004), the theorem establishes that the asymptotics of FDR and power can be precisely achieved by replacing each term in (1) and (2) with its population counterpart. Notably, in this context, l" represents the population version of the BH rejection threshold for the p-values.
(x, y) = M -fy E R - fol.,,, (1I)
where M is large constant and fo x' - R is flexible function parametrized by 0, that can be chosen from a specific machine learning model class, such as kernel ma- chines, gradient boosting models or neural networks, etc. The first term, an indicator function identical to D.(.) in (8), ensures regional monotonicity in Definition 3.1 and boosts selection power, suggested in Section 4.1. The second term generalizes the distance term D2(-) from (6), offer- ing a more expressive framework for constructing optimal nonconformity scores.
Theorem 4.1 indicates that the second, clipped score (8) is preferable to the first score (7) for achieving higher power. According to (10), since the value V(x, r) = in || - 8| is identical for both scores assuming r E aR, it suffices to compare their asymptotic BH thresh- olds t". The score with the larger 1" achieves higher asymp- totic power and is therefore more effective. In the definition of tt in (9), the fraction can be rewritten as
The following result demonstrates the expressiveness of the family (11) by showing that it can include the optimal nonconformity score for any selection task. The proof is provided in Appendix A.3.
P(F(V U) t)
Gv(()= = (2,
P(F(V(a: r),U) 5)
((V(a, y) <)
Let {(x, } be sampled i.i.d. from a distribution, and assume a fixed choice of Tn+j = T. Un- der Algorithm 1, for any nominal FDR level q and target region R, there exists a function f such that the score con- structed using f in (11) maximizes the number of selected samples (and thus the power) among all scores with FDR control.
Proposition 4.2.
((V(x, s)
when s is the inverse of F(:U) U) at f, assuming it exists. The first equality follows from the fact that ((V(x,)) U) ~ Un, 1), while the second arises from the monotonicity of F with respect to its first argument v. To maximize t in (9) (equivalently, the inverse s") ), an effective score should yield larger V(æ, y) relative to V(æ, r), thereby reducing (t for a fixed f. This, in turn, results in larger /" when computing sup(t: (t) < q} in (9).
Remark 4.3 (Subfamilies of the score class). A notable subfamily of (11) is
M 1yy E R (QR))----)))) - fe((()) R),
(12)
This criterion is precisely satisfied by the clipped score. An alternative justification for favoring the clipped score,
where fo depends solely on the prediction ((x). This sub- family includes the clipped distance-based score from Sec-
5

---

Multivariate Conformal Selection
tion 4.1 as a special case. Here, regional monotonicity is guaranteed for any constant M, but the score also be generalized further by allowing f to depend æ as well:
v in (11). Here, D and Dj D-iiinn serve the cali- bration dataset and test dataset respectively, to obtain the smoothened conformal p-value for training.
Next, to smooth the BH procedure, we first apply the soft- sorting operation to the smooth p-values i7 to obtain their corresponding ranks uf.
M (y E RR aR] - (1 ((ar);; R).
In contrast, the broader family defined in (11) offers greater flexibility by incorporating y in the second term. However, this added expressiveness requires a sufficiently large M to preserve regional monotonicity. Specifically, M must satisfy M > fo (z. y; R) In practice, M is chosen to be sufficiently large to ensure that this inequality holds across the entire dataset.
af soft-ran (
(14)
and then compute the softened selection size (()) as
5()) ( log Σ E000 where of (1) a df//m' - D9
Remark 4.4 (Incorporating pretrained models). This could be achieved in several ways, such as using the predictions of f as inputs to fo. or train fo as a prediction on top of when both models are implemented neural networks. While such practice does not increase the expressiveness of the score family, it often facilitates the training of fa. as p(x) estimates y and is very informative for selection. Since fo can directly learn the data and the selection task, lea-n can still perform well when f is poorly fitted; see Appendix C.2.4.
T
In the above equation, the sigmoid function with tempera- ture coefficient T in (ii) serves as a smooth approximation of the indicator function 1 (D8 < qaf./m") while the log- sum-exp function in (i) approximates the element-wiss max function (ax (afss......... an s ). This formulation ensures that (0) closely approximates the BH selection size, as the selection size is defined by the largest rank with a p-value below the threshold.
To identify an optimal function within the family (11), we introduce a differentiable loss function that mimics the inher- ently non-differentiable mCS procedure. The "hard" sorting and ranking operations the mCS workflow are replaced with their smooth, differentiable counterparts (Blondel et al., 2020; Cuturi et al., 2019). We adopt the implementation in- troduced in Blondel et al. (2020), with lg regularization and regularization strength set to 0.1. The resulting loss function is then used for a chosen machine learning method to train fo. Specifically, we partition the calibration data into three batches Deal= = DJ J-urin U Dy-val U D' where D y-urain and Dasval are used for training and validating fo, respectively. Upon completion of training and validation, Algorithm can be applied with D' the calibration dataset for 8 to generate the final selection set S.
To maximize the selection size, the loss function for learning the mCS-learn score can be defined as
L (0) S()).
(15)
While above formulation of the loss I22 is intuitive in its attempt to approximate the final selection size using differen- tiable functions, the inclusion of two soft sorting steps may reduce numerical stability and impede the training process.
In the mCS procedure, the BH procedure is solely intended for multiplicity correction for FDR control, and thus not necessarily required in the formulation of the loss function, whose primary objective
learn the function fe. A sim- pler yet effective alternative loss function directly penalizes the smooth p-values pt through the following loss function:
Training Step. The training loss function is defined based on the smoothed selection size. We denote the softened rank of an element E A within the set A by soft-rank(44 (a; A). We randomly partition D into two subsets tmii and Dymmin., and we assume D D-mmin = {( y=) and D Di-arinn = "( 1 })
m
(2
0 D3 [I >
E R) 1- E ER)] (16)
where p is minimized when the j-th sample is deemed desirable, as indicated by the first term. The second term, scaled by a balancing coefficient 7, ensures that the p-value is not uniformly small but becomes relatively larger for less favorable samples. This approach eliminates the need to es- timate the BH threshold via a secondary soft-ordering (14), leading to improved numerical stability and enhanced over- all performance. A comparison of the two loss functions can be found in Appendix C.2.5. After computing the loss in each epoch, we can follow the standard backpropagation procedure to train 0.
for notational sim- plicity. We define the smooth conformal p-value F for j == 1,...,,m m' as
soft-rank (V
(g (V U (D n' 11
p3
(13)
where Vi0 : (0 v,) are computed on D, Darrain and De
#+ = V (2n Fn(+z) are computed on Dy iminn, with
6

---

Multivariate Conformal Selection
Validation Step. To avoid overfitting Sa. we perform an additional model selection procedure using hold-out dataset y-val Specificallyy, for each epoch & I,...,T of the backpropagation procedure, we apply K random partitions on y-vvll to obtain D f-vall and f-vat2 for k I..... K. For each k we then apply Algorithm I with the setting Deal = J-vall and Diest" == D f-val2, and record the validation power Pio (t) We then compute the average validation power Power(t)= = E Pt(t)/K.
adapted from CS address the multivariate case. Later, we employ these adapted methods as baselines and compare their performance against our proposed method.
In the scenario when the target region R C R is rectangular, the overall selection criterion can be decomposed, allowing each dimension to be evaluated independently. By applying CS separately to each dimension, we obtain d selection sets S Sa. where each Sa contains observations satisfy- ing the k-th corresponding marginal criterion. The final selection set S is then given by the intersection of these individual sets, S = nd S SR. We refer to this approach as CS iint..
In the end, we select from {1..... T to be the epoch with the highest average validation selection power Power(t), and deploy the associated model fo; for the final selection.
It can be shown that Cs_int, when each CS subroutine is conducted a nominal level q, fails to control the FDR at or below q. This issue is analogous to the intersection hypothesis testing (IHT) problem in statistics. A common approach to address this is the Bonferroni correction (Dunn, 1961), which adjusts the nominal level of each subroutine to a lower threshold g/d. Then obtain the intersection of individual sets. However, this method is widely recognized as being overly conservative (Pereger, 1998; Westfall & Young, 1993). We refer to the Bonferroni-adjusted CS_int as CS ib. An alternative to the Bonferroni correction is to adaptively account for intersection hypothesis testing, Rather than predefining the nominal levels for subroutines, we determine suitable values by validating on a hold-out datasct. This approach necessitates additional data splits to construct the hold-out set, and we refer to it as CS_is.
Finally, Algorithm 2 details the complete learning procedure for MCCS-eearn scores. The code for reproduction can be found at
com Ti -inn-Baaaa /mas.
Algorithm 2 MCCS- Learning Procedure
Input: Training data D-trin validation data y-va target region R. FDR level y E (o, 1), other hyperparameters. Output: Trained nonconformity function fe.
1: Initialize parameters 0 = 0g-
2: for epoch ( t= 1,..,, T do
Training Step
3: Randomly partition y-uain into two disjoint sub- sets Da-rriin and D f-rain2--
4: Use the current fe to obtain Vi from D-rtrint and
Di + from D --min..
5:
Compute the smooth conformal p-values Dy (13).
Beyond considering each dimension separately, another nat- ural adaptation of CS involves transforming the response vector before applying CS. Specifically, each response y: is converted to a binary indicator reflecting whether it meets the selection criterion, defined as ÿ = 1 [3, E R). Under this transformation, the new selection threshold can be set to c 0, as y; > is equivalent to y; E R. We refer to this approach as bi.
6: Compute the loss function using (15) or (16).
7:
Back- pookkgropeette to update the parameters = O₂.
Validation Step
8: Apply K random partitions on D f-val to obtain DiA) and {A) J-vat2 for each & = 1, .,...,,.. K. For each k, apply Algorithm 1 for score func- tion v with the setting Deal == D Dras := Dia) f-vai2 and compute validation power Pk Compute Power(t) = Σk Pt(t)/K.
9:
5.2. Numerical Results
We compare the performance of mCS dt, wws-mmmsaaaaae (abbreviated mCS-d and mCS-1 respectively) against the baseline methods outlined in Section 5.1.
10:
11: end for
12: Determine i = arg max, Power(t) and return for-
In our data generation processes, covariates æ are sampled uniformly from Unif(-1, 1)" where is the covariate di- mension, and the responses y are generated as y = f((x) +c, where j denotes the regression function and represents noise drawn from either a multivariate Gaussian or multi- variate t-distribution. By varying the regression function, the size of response dimensions, and the choice of Gaus- sian or heavy-tailed noise, we create a range of selection problems with differing levels of difficulty.
5. Simulation Studies
5.1. Baseline Methods
While the standard CS approach is originally designed for univariate settings and cannot be directly applied to multi- variate selections, appropriate adaptations can be devised. In this section, we introduce several naive methods directly
We consider two selection tasks where the target region R.
7

---

Multivariate Conformal Selection
is defined as:
ing various factors of the selection problem and the mCS algorithm are provided in Appendix C.2.
Task 1. The (shifted) first orthant, R = {y: : yk 2 ck vk] Task 2. A sphere centered at c. R E 1y c||2 r).
Table I:
Observed FDR for Task (shifted first orthant R).
The nominal FDR level is g = 0.3.
These two specific tasks are particularly relevant in applica- tions, as they simulate scenarios where (1) d criteria must be simultaneously satisfied or (2) an instance must be suffi- ciently close to specific point e to be deemed acceptable.
Setting C$_int CS.ib CS.is bi
mCS-0
mCS
0.773
0.000
0,280 0.240
0.277
0.251
0.801
0.000
0.133 0.300
0.315
0.266
0.724
0.000
0.204 0.295
0.264
0.278
In our simulation, the coefficients с and r for each selec- tion problem are chosen to ensure that approximately 15% to 35% of the data points lie within R across all six data- generating processes. Detailed descriptions of the data gen- eration process, model specifications, and the specific values of the coefficients are provided in Appendix C.1. We first train a support vector regression model ja using 1000 data points, and use an additional labeled dataset of 1000 samples to construct selection sets for different methods in compar- ison. We evaluate the selection power and FDR using test dataset of size 100. We adopt the clipped score (8) for mCS-dist and adopt the loss function in (16) with balanc- ing coefficient y = 0.5 for -earn. For mCS-mmasrarre the calibration data is split to D y-train D-val and Deai with ratio 111:1, and the model fe is formulated as a two-layer MLP with batch normalization. The response dimension is set to be d 30, and nominal FDR level is set at q 0.3. Number of iterations for validation is set to K- 100. The selection process is repeated across 100 iterations, with new dataset generated independently for each iteration.
4
0,811
0.000
0.255 0.309
0.277
0.315
0,810
0.000
0.300
0.266
0.308
0.239
0.778
0.000
0.374
0.245
0.287
0.258
Table 2: Observed power for Task 1 (shifted first orthant R).
Setting cs int
CS ib
CS is
bi
mCS-d mCCS-
1.000
0.000
0.406
0.324
0.555
0.325
1.000
0.000
0.012
0.069
0.104
0.108
3
1.000
0.000
0.039
0.059
0.068
0.102
4
1.000
0.000
0.124
0.194
0.324
0.198
1.000
0.000
0.019
0.035
0.060
0.042
1.000
0.000
0.101
0.027
0.046
0.034
Table 3: Observed FDR and power for Task 2 (spherical R).
The nominal FDR level is q = 0.3.
FDR
Power
Setting
bi
mCS-
mCS-1 1
bi
MCS-d
MCS-
0.255
0.265
0.279
0.636
0.760
0.534
0.260
0.263
0.273
0.343
0.405
0.421
0.216
0.223
0.263
0.115
0.134
0.179
0.316
0,286
0,254
0.192
0.333
0.180
Table and Table 2 summarize the experimental results for the first selection task. As shown in Table I, Cs_int substantially violates FDR control. cs_is s provides only approximate FDR control, and in scenarios such as Setting 6, the FDR control may be compromised. In each setting, the red numbers in Table 2 indicate the highest achieved power among all methods that properly control the FDR. Among the four remaining methods that always maintain valid FDR control, our two proposed methods consistently achieve the best and second-best power, outperforming the baseline methods under all settings.
5
0.307
0.291
0.273
0.137
0.170
0.189
6
0.292
0.283
0.207
0.055
0.063
0.061
6. Real Data Application
We apply the mCS framework to drug discovery, selecting drug candidates with desirable biological properties. Each candidate corresponds to a chemical compound, where the feature vector encodes structural chemical characteris- tics, and the multivariate response y represents biological properties. The multidimensional nature of y reflects the need to evaluate multiple biological criteria before advanc- ing a compound. Ensuring FDR control improves down- stream processes, such as wet-lab validation.
Table 3 summarizes similar results for the second selec- tion task. Baseline methods Cs_int and CS_ib are not included as they are not applicable to non-rectangular target regions. Figure shows the realized FDR and power curves across varying nominal FDR levels, ranging from 0.05 to 0.5 in increments of 0.05. Results are shown exclusively for Setting 3 due to space constraints. Among the methods compared, mCS-dist iss, mCS-learn, and bi demonstrate consistent FDR control, as their respective curves remain below the black dashed line indicating the nominal FDR threshoid. Notably, higher power across all nominal levels.
We employ an imputed public ADMET dataset compiled from multiple sources (Wenzel et al., 2019; Iwata al., 2022; Kim et al., 2023; Watanabe et al., 2018; Falcón-Cano et al., 2022; Esposito et al., 2020; Braga et al., 2015; Alia- gas et al., 2022; Perryman et al., 2020; Meng et al., 2022; Vermeire et al., 2022), comprising n = 22805 compounds with d 15 biological assay responses. We focus on three different selection tasks with the following target regions:
-earn also achieves consistently
Additional simulation results and ablation studies explor-
Task 1. The (shifted) first orthant, R = { y/k 2 Ck Vk),
8

---

Multivariate Conformal Selection
Task FDR
Task 2 FDR
Task Power
Task
0.5
0.7
1.0
0.5
0.6
0.4
0.5
0.3
0.4
0.3
0.2
0.2
0.1
0.1
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.5
0.1
0.2
0.3
0.4
0.5
Nominal FDR level
Nominal FDR level
Nominal FDR level
Nominal FDR level
bi
mCS-d
mCS-I
CS CS_int
CS_ib
CS_is
Figure 1: Observed FDR and power across varying nominal levels for Task 1 (shifted first orthant It) and 2 (spherical R).
Task 2. A sphere centered at c, R ty: cll2 S r). Task 3. The complement of sphere centered at R {y: |y c|22 r.).
Table 5: Observed power of methods with real data.
Task
q
CS cs CS_ib cs_is
bi
mCS-d d mCS-1
0.3
0.993
0.000 0.019 0.000
0.006
0.010
2
0.3
0.000
0.278
0.086
We included the two tasks introduced the numerical sim- ulations (Section 5), and we also designed the third task to evaluate the performance of various methods under a non- convex target region with real data. Each of the selection task is designed so that approximately 15%-30% of the compounds qualify as acceptable. Further details on the dataset and problem setup, including response descriptions and cutoffs, can be found in Appendix D.
0.3
0.000
0.410
0.431
0.5
1.000
0.003
0.225
0.000
0.433
0.193
0.5
0.000
0.759
0.515
0.5
0.000
0.449
0.589
mCS-leart exhibits superior performance under noncon- vex target regions, corroborating the results presented in Appendix C.2.3. We note that although the baseline method bi performed well in the simulation settings, it barely se- lected any compound in the current task. This outcome may be attributed to the suboptimal performance of the underly- ing binary classification model, which achieved an Fl score of only 0.31.
In this experiment, the underlying predictor ji is specified as drug property prediction model from the DeepPurpose. Python package (Huang et al., 2020) with Mo rgan drug encoding. We train the model using Yraain 12000 samples, provide neat = 8000 samples for calibration and reserves the remaining data of size Thest = 2805 as test data, We keep the configuration and hyperparameters of the methods unchanged as in Section 5. Two nominal levels q 0.3 and q = 0.5 are considered, and the selection processes for each method are repeated across 500 iterations.
7. Conclusion
We propose multivariate conformal selection, an extension of conformal selection to multivariate response settings. Our experiments demonstrate that mCS significantly improves selection power while maintaining rigorous FDR control, outperforming existing baselines across simulated and real- world datasets. The flexibility of mCS makes it a valuable tool for selective tasks involved in diverse fields including drug discovery. Looking forward, we anticipate that the mCS framework can be further extended to handle addi- tional practical challenges, including settings with hierarchi- cal or structured responses. By addressing these challenges, mCS has the potential further enhance its applicability and impact across diverse scientific and industrial domains.
Table 4: Observed FDR of different methods with real data.
Task
q
cs_int CS.ib cs.is bi mCS-d mCS-1
1
0.3
0.760
0.000 0.303 0.038 0.304
0.275
2
0.3
0.000
0.300
0.293
3
0.3
0.084
0.301
0.296
I
0.5
0.782
0.393
0.496
0.040
0.499
0.488
2
0.5
0.000
0.499
0.498
3
0.5
0.084 0.501
0.497
Tables 4 and Table 5 summarize the FDR and power of the competing methods respectively. Consistent with our numerical simulation, the methods bi, mCs-dist, and all demonstrate valid FDR control. Among the methods guaranteed to control FDR, mCS-dist and mCS-learn consistently achieve the best and the second- best power across all tasks and nominal levels. Notably,
Impact Statement
This work does not present any foresecable societal conse- quence.

---

Multivariate Conformal Selection
Acknowledgments
C. H. Pred-herg: A novel web-accessible computational tool for predicting cardiac toxicity. Molecular Informat- ics, 34(10))698-701, 2015.
This work was supported by Natural Sciences and Engineer- ing Research Council (NSERC) Discovery Grant (RGPIN- 2024-06780) and FRQNT Team Research Project Grant (FRQ-NT 327788).
Cuturi, M., Teboul, O., and Vert, J.-P. Differentiable ranking and sorting using optimal transport. Advances Neural Information Processing Systems, 32, 2019.
References
Dunn, O. Multiple comparisons among means. Journal of the American Statistical Association, 56(293):52-666 1961.
Aliagas, I,, Gobbi, A., Lee, M.-L., and Sellers, B. D. Com- parison of logp and logd correction models trained with public and proprietary data sets. Journal of Computer- Aided Molecular Design, 36(3):253-26,, 2022.
Esposito, C., Wang, S., Lange, U. E., Oellien, F., and Riniker, S. Combining machine learning and molecular
Bai, T., Tang, P. Xu, Y., Svetnik, V., Khalili, A., Yu, X., and Yang, A. Conformal selection for ef- ficient and accurate compound screening in drug discovery. ChemRxiv, 2024. doi: 10.26434/ Dhemxxix-2024-p33ph-22 URL URLhttps://chemx/t hhee//engggeeccmmmxx///ttttttiie 67cffbbie8812211122224ffff
dynamics to predict p-glycoprotein substrates. Journal of Chemical Information and Modeling, 60(10))4730-4749 2020.
Etemadi, N. and Kaminski, M. Strong law of large num- bers for 2-exchangeable random variables. Statistics probability letters, 2843):24--200 1996.
Falcón-Cano, G., Molina, C., and Cabrera-Pérez, M. Á. Reliable prediction of caco-2 permeability by supervised recursive machine learning approaches. Pharmaceutics, 14(10):1998, 2022.
Bai, Y., Jones, A., Ndousse, K., Askell, A., Chen, A., Das- Sarma, N., Drain, D., Fort, S., Ganguli, D., Henighan, T., et al. Training helpful and harmless assistant with rein- forcement learning from human feedback. arXiv preprint arXiv:2204.05866. 2022.
Feldman, S., Bates, S., and Romano, Y. Calibrated multiple- output quantile regression with representation learning. Journal of Machine Learning Research, 24(24): 2023.
Bates, S., Angelopoulos, A., Lei, L., Malik, J., and Jor- dan, M. Distribution-fie free, risk-controlling prediction sets. Journal of the ACM (JACM), 18(6):1-34, 2021.
Gui, Y., Jin, Y., and Ren, Z. Conformal alignment: Knowing when to trust foundation models with guarantees. Ad- vances in Neural Information Processing Systems, 2024.
Bates, S., Candès, E., Lei, L., Romano, Y., and Sesia, M. Testing for outliers with conformal p-values. The Annals of Statistics, 51(1), February 2023. ISSN 0090-5364. doi: 20.1214//2--os2222 URL http: //tp:///xx./ doi. orr// 10. 12..121/220 / AOS2244.
Heid, E.,, Greenman, K. P., Chung, Y., Li, S.-C., Graff, D. E., Vermeire, F. H., Wu, H., Green, W. H., and McGill, C. J. Chemprop: machine learning package for chemical property prediction. Journal of Chemical Information and Modeling, 64(1):9-17, 2022
Benjamini, Y. and Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical So- ciety. Series (Methodological), 57(1):289-300, 1995. ISSN 00359246. URL httt:://wwtt httpp:/www.jssor..co jstor.org/ stable/ / 2246101..
Huang, K., Fu, T., Glass, L. M., Zitnik, M., Xiao, C., and Sun, J. Deeppurpose: deep learning library for drug- target interaction prediction. Bioinformatics, 36(22-23): 5545-5547, 2020.
Benjamini, Y. and Hochberg, Y. Multiple hypotheses testing with weights. Scandinavian Journal of Statistics, 24(3):
Iwata, H., Matsuo, T., Mamada, H., Motomura, T., Mat- sushita, M., Fujiwara, T., Maeda, K., and Handa, K. Predicting total drug clearance and volumes of distri- bution using the machine learning- learning-mediated mediated multimodal method through the imputation of various nonelinical data. Journal of Chemical Information and Modeling, 62 (17):4057-4065, 2022.
407-418, 1997. doi:
00072. URL htt.ss::/onniiielliiaaaaaats wiley. do! com/doi abs 10. 111.1.1111111111
com
00072.
Blondel, M., Teboul, O., Berthet, Q., and Djolonga, Fast differentiable sorting and ranking. In International Con- ference on Machine Learning. pp. 950-959, PMLR, 2020.
Jin, Y. and Candès, E. J. Selection by prediction with con- formal v-luuus. Journal of Machine Learning Research, 24(244):1-41 2023.
Braga, R. C., Alves, V. M., Silva, M. F., Muratov, E., Fourches, D., Lião, L. M., Tropsha, A., and Andrade,
10

---

Multivariate Conformal Selection
Johnstone, C. and Cox, B. Conformal uncertainty sets for robust optimization. In Conformal and PndPPre Probabilistic Prediction and Applications, pp. 72-90. PMLR, 2021.
Sesia, M. and Candès, E. J. A comparison of some confor- mal quantile regression methods. Stat, 9(1):226,, 2020.
Sheridan, R. P., McMasters, D. R.,, Voigt, J. H., and Wildey, M. ecounterscreening: using qsar predictions to priori- tize testing for off-target activities and setting the balance between benefit and risk. Journal of Chemical Informa- tion and Modeling, 55(2):231 2015.
Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., Yu, B., al. Pubchem 2023 update. Nucleic Acids Research, 51(D1): D1373-D1380, 2023.
Kivaranovic, D., Johnson, K. D., and Leeb, H. Adaptive, distribution-free prediction intervals for deep networks. In International Conference on Artificial Intelligence and Statistics, pp. 4346 4356. PMLR, 2020.
Storey, J. D., Taylor, J. E.,, and Siegmund, D. Strong con- trol, conservative point estimation and simultaneous con- servative consistency of false discovery rates: unified approach. Journal of the Royal Statistical Society Series B: Statistical Methodology, 66(1):187-205) 2004.
Klein, M., Bethune, L., Ndiaye, E., and Cuturi, M. Mul- tivariate conformal prediction using optimal transport. arXiv preprint ar arXiv:2502.03609 2025.
Szymański, P., Markowicz, M., and Mikiciuk Hikiciuk-Olasik, E. Adaptation of high-throughput screening in drug discov- ery- toxiiooicologaaal screening tests. International Journal of Molecular Sciences, 13(1):427-452 2011.
Kuleshov, A., Bernstein, A., and Burnaev, E. Conformal prediction in manifold learning. In Conformal and Proba- bilistic Prediction and Applications, pp. 234-253. PMLR, 2018.
Vermeire, F. H., Chung, Y., and Green, W. H. Predict- ing solubility limits of organic solutes for wide range of solvents and temperatures. Journal of the American Chemical Society, 144(24):10885-1777) 2022.
Lei, L. and Candès, J. Conformal inference of counter- factuals and individual treatment effects. Journal of the Royal Statistical Society Series B: Statistical Methodol- ogy, 83(5):911-933 2021.
Vovk, V., Gammerman, A., and Shafer. G. Algo- rithmic Learning in a Random World. Springer, 2005. URL https: htttps:/aap..smmmtiiiiiee www.nticscoollal.. org/CorpusID: 1111088325
Meng, J., Chen, P., Wahib, M., Yang, M., Zheng, I.., Wei, Y., Feng, S., and Liu, W. Boosting the predictive perfor- mance with aqueous solubility dataset curation. Scientific Data, 9():71,, 2022.
Watanabe, R., Esaki, T., Kawashima, H., Hatsume-Kitatani. Y., Nagao, C.. Ohashi, R., and Mizuguchi, K. Predict- ing fraction unbound in human plasma from chemical structure: improved accuracy in the low value ranges. Molecular Pharmaceutics, ((11)):502-5333 2018.
Messoudi, S., Destereke, S., and Rousseau, S. Ellipsoidal conformal inference for multi-target regression. In Con- formal and Probabilistic Prediction with Applications, pp. 294-306. PMLR, 2022.
Wenzel, J., Matter, H., and Schmidt, F. Predictive multi- task deep neural network models for adme-tox properties: learning from large data sets. Journal of Chemical Infor- mation and Modeling, 59(3):1253-1) 1263, 2019.
Park, J. W., Tibshirani, R., and Cho, K. Semiparametric conformal prediction. arXiv preprint arXiv:2411.02114 2024.
Perneger, T. V. What's wrong with bonferroni adjustments. BMJ, 31 12311199 1998.
Westfall, P. H. and Young, S. S. Resampling Based Multiple Testing: Examples and Methods for p-value Adjustment, volume 279. John Wiley & Sons, 1993.
Perryman, A. I.., Inoyama, D., Patel, J. S., Ekins, S., and Freundlich, J. S. Pruned machine learning models to predict aqueous solubility. ACS Omega, 5(27):16562- 16567, 2020.
Yang, K., Swanson, K., Jin, W., Coley, C., Eiden, P., Gao, H., Guzman-Perez, A., Hopper, T., Kelley, B., Mathea, M., et al. Analyzing learned molecular representations for property prediction. Journal of Chemical Information and Modeling, 59(8) 33888, 2019.
Romano, Y., Patterson, E., and Candes, E. Conformalized quantile regression. Advances in Neural Information Processing Systems, 32, 2019.
Zhang, K., Yang, X., Wang, Y., Yu, Y., Huang, N., Li, G., Li, X., Wu, C., and Yang, Artificial intelligence in drug development. Nature Medicine, pp. 1-15, 2025.
Scannell, J. W., Bosley, J., Hickman, J. A., Dawson, G. R., Truebel, H., Ferreira, G. S., Richards, D., and Treherne, J. M. Predictive validity in drug discovery: what is, why it matters and how to improve it. Nature Reviews Drug Discovery, 21(12):915 931, 2022.

---

Multivariate Conformal Selection
A. Technical Proofs
A.1. Proof of Proposition 3.2
Proof of Proposition 3.2. Since V is fixed, the nonconformity scores Vi...., result from conformal inference Vovk et (2005, Proposition 2.4), the oracle p-value pj defined as in (3) is uniform distributed with value ranging in (0, 1), and P(y) S <a) a. This gives
,, Vn and Vn+s are exchangeable. By a standard
P Σ a and j E Ho) S a.
When the null hypothesis Hoj is true, Un+i E R Since ra+j E R, by the regional monotone property, we have Vn+5 (a Un+j) S Vi (nts) = Vn+j. We then have Pj S Pj by definition, and
P(P S a and j E Ho) - P(P) S a and j E Ho) S a.
A.2. Proof of Theorem 3.5
Proof of Theorem 3.5. We adapt the proof of Theorem 6 in (Jin & Candès, 2023). In the proof that follows, we fix index j E (1...,,m) m.. For notational simplicity, only in this proof we deal with the deterministic conformal p-values,
1 Pj = n
n
Σ
1( S Pn+sll
We note that the deterministic conformal p-values are only valid when the scores (Vi) (Va+))7= have no ties almost surely, and therefore we also make this assumption. We highlight that the validity of the random conformal p-value does rely on this statement. We also define the corresponding deterministic oracle p-values,
pj
1+
HV, S Vn+t3]]
n
This version p-values are also conservative by standard results in conformal inference (Bates et al., 2023; Jin Candès, 2023). For l,..., m, define set of slightly modified p-values,
p = 1 Σ 1( Va+T + 1(v++s S P.+ n n 1+
Also define S(a,...., am S (1....,ml m) as the rejection index set obtained by the Benjamini-Hochberg procedure, from p-values taking values ai,...y am- Then, the output mCS is S(pi, PPm.)
In the sequel, will compare S(p Pm) to
S(p (1)
11) P1 P 1411 .P)))
(1)
for the test sample j that is falsely rejected, i.c.
E S. Yutj E R.
First, on this event, we have Vn+j V(zntj U Bn4g) V(2n43, Vnt, and p) S Py. For the remaining valuuss (1) 1 hez we consider two cases:
If V P Vatg. then p 2 Py. In this case, we also have V.n 2 Vatj as Vntj 2 Vnng. This implies
Σ 1(v P.4t + 11Vnts $ Vatil] m+1
1[, < Vati + +1] = Pt-
1
12

---

Multivariate Conformal Selection
(ii) If If Patg. then p
S Pj. We also have
a Vn+1 +
1[
1[, & Putst +1] = Pj.
n
n
n +
2=1
Since j E S, by the definition of Senjamini-Hochber procedure ( E S pt has smaller rank when ordering the p-values.
To summarize, if we replace Pm.) by (p) ......P--11111 : Pj. Pi...-....... Pin ) on the event {j E SnY RR), ), such replacement does not modify any of those p-values pt if they satisfied pt Z Pj- Also, for all p-values pt with pt S Pj including j itself (1 = j). their replaced values pt are still no greater than Pj- Since all p-values are no larger than their original values after the replacements, the size of rejection set must not decrease. On the other hand, since j E S and no p-values larger than ア are modified, no new hypotheses can be rejected by the new set of p- values. We conclude that
,,,""
())
1 Sj = S ()) ................................ ...........P......... Pjt1 . S(p...... .mm. = S
on the event (Ya+s E j E S). By decomposing the FDR we have
Σ m 1(Y,t) E R E $)
FDR = E
(8)
m
S E.je SH
=
m
qk
m )
-[11s]]-= = gh
E[11S)]1= STH
The last equality is again by the property of the Benjamini-Hochber, procedure. Also, by its step-up nature, sending p1 to zero does not change the rejection set if the corresponding hypothesis of p3 is rejected, i.. on the event {j E SS}. We have
(2) = .........PPPPPP................. ... p( s(( 1,
9)
p()) 87
Thus,
1( 41571/m} Σg E
1 1if S (/5)mm//mm
FDR <
o[ E S}]] = ==
E
(1,
max 1S)
By definition, (0( Jizz is invariant after permuting (1i) U (Vrtz)- Since (Vi) U (Vn+s} are exchangeable conditioned on (V~+theg (izj, the distribution of ( Jizy is independent from the ordering of {Vi)== U ((n+z) conditioned on the (unordered) set [V,.... Vn, Vn+ U Vin4shees. Also, conditioned on (Vatt1es, S only depends on {pt () Jızj which is in turn only dependent on the unordered set [Vi,..., Vn, Vn+ and Pž only depends on the ordering of (v] U
{Vn+s). ). This implies that -o is independent on p) conditioned on the (unordered) set [Vi...., Vn, Vn+zl and (V+t1445. Therefore, by the conservative property of conformal p-values and conditional independence,
So0 m
PP Ppp
[Vi,..., Vn, Vatal U (Pote)ers)
m
13

---

Multivariate Conformal Selection
By the law of total expectation,
Lip S {05)
LPp S 4/5j-+0l/m)
E
E
E
[V. Vn Vn+zl u (uttiets))
max(1,15)
max
S1S}-tol
E 1115 # #0}
[V2...., Vny Vn+sl U (Vatit+s
ml m|
q
m
Since when S7-u 0, 1(p 2
(/0j-=o///) ///)m 0. Now, the proof is concluded by summing over every j ....., m.
A.3. Proof of Proposition 4.2
Proof of Proposition 4.2. For any specific dataset and selection problem, let Vop denote the optimal nonconformity score that controls the FDR and maximize the number of selected samples. Define V, (V y₁) for calibration samples and propi yopl (Ta+jj Tn+ÿ) for test samples. Consider the following score W within the family (11):
n+j
W(r, y) = M 1 E R OR) + (x,
Picking rn+j = r E ÖR, the test nonconformity scores Wn+; satisfy:
Wn+s = WIEn+js = Vop (an Tn+s) = porm
and for each calibration score,
W; ((x,)))) 2 VO y,) = V,
Therefore, by replacing Vou with W, the calibration nonconformity scores Wi,... Wn, may increase, while the test nonconformity scores n+ remain unchanged. According (4), this leads decrease in each conformal p-value pj- which in turn increase the value k* in Algorithm 1. This means that would select more samples, by the definition of S. Consequently, we conclude that there must exist score within the family (11) that achieves the optimal selection size. A similar argument applies to the maximization selection power, thereby concluding the proof.
B. Deferred Discussions
B.1. Diseussions on Classification Responses
In our paper, we choose to focus primarily on the regression response case for multivariate conformal selection. In fact, the selection problem for classification responses (univariate or multivariate) can be directly reduced the univariate conformal selection framework introduced by Jin & Candès (2023):
For the univariate classification setting, suppose the response space is composed of classes y = _ C2 with target region R = - Ck (with < K). Then, by defining binary response: h = 1s E R). the original selection problem directly translates into a univariate conformal selection problem, where we select samples with ÿ; =
For multivariate classification case, e.g. suppose that bivariate responses y; are drawn from a joint class space y (p(())(2)) (C( c( and the target region R = U(e.)ez(((() c(2)). Again, we define binary indicator 9: = 11y E RJ, converting the original multivariate selection problem into a standard univariate selection task:
Hoj Dnts < 0.5 versus H1 Dnts 2 0.5.
Since P(j: 1) P(y, R), there is direct correspondence between the multivariate and univariate nonconformity scores:
V( = M 1 e R] - P(y E R|x)
and
V( M 10 2 0.5) (x)))
14

---

Multivariate Conformal Selection
where j((x)= Pi l.). Moreover, regional monotonicity is simply the usual monotone condition of univariate conformal selection: V( 1) V -=
Thus, every classification- based selection task can be naturally and effectively solved using existing univariate conformal selection methods.
B.2. Additional Discussions on Theorem 4.1 and Advantages of The Clipped Nonconformity Score
In this section, we first provide interpretations about Theorem 4.1 and offer an alternative perspective on the advantages of the clipped nonconformity score.
We first note that the (practical) conformal p-values defined in (4) can rewritten as
Pj
11
< Vats} n
U,
3 1[V, Puts + + 1
1 U,
2m1
By the Seivenko-Caantllly theorem under exchangeability (Etemadi & Kaminski, 1996), as n + 00,
1
sup 1RR
LI
<tt+
LIV.
t m+1* U, (V((2, y) < - U, - (.V =t) 3 0.
71
n+1
1=1
1=1
Replacing & the above by Vn+j = V (z r) yields
Pj 3 (x, < ((m4)).)) r)) + U, (VV(v, v) (((nss,)) = (V(
1 F(((x, r), U).
Therefore, F(V (x, r), U) is also conservative random variable.
The quantity tn' in Theorem 4.1 can also be viewed as the asymptotic counterpart of a corresponding finite-sample quantity. characterization of the BH procedure (Storey al., 2004) states that the rejection set = {j P; S ST"}. where the BH rejection threshold 7" defined as
T" sup { EE []
1 Σm (0, s saf
By the fact that J Σラ Hp < P(pj 1) 2 0 a m - x (which follows due to the the weak law of large numbers for triangular arrays), and our earlier characterization of F(V (x. r), Uy) which implies Pp,, I =t) - r() as n - oo, the fraction the definition of 7* satisfies (in the limit n, m - 0o)
-n2 EE Hp 5t
r)
r), U) St)
This establishes 1" as the asymptotic limit of the BH rejection threshold 7", which further implies
1 (ps ST (F) 1*) and
(p, - and Un+j E R) B U) R).
Theorem 4.1 then gives the asymptotic version of FDR and power:
S L ynts E R$
1 P, ST yntz E RCH
FDR = E
[S]
: S TTH
1
1(p < T" and Uatj E R
E
1 E 1p, ST*
U)) R
U))
15

---

Multivariate Conformal Selection
[8n (i Dn+s E RH
Ki Pj:Ps S ST (i Yn+s E R}]
Power = E
Híj Untj E RH]
Kí : ynts E R
-==1 1 (p) T and Un+j E R]
E
ンー= 1 1yn R}
((F(V r). U) St . E R)
U)
P(y E R)
Theorem 4.1 suggests that the clipped score should be preferred from the perspective of maximizing the realized (asymptotic) FDR (which is still controlled below the nominal level). The asymptotic FDR is bounded through the following chain of inequalities:
F)) R
() (()()
U))
() (()((
(17)
P(P(V(E)))
By the definition of t", the final term in (17) is at most the nominal level g. the bounds in this chain of inequalities could be further tightened, the realized FDR would more closely align with the nominal level q, effectively allowing for greater selection power as more of the FDR budget could be utilized.
Using the clipped score tightens the inequality (+) in (17). Assuming the score function V is clipped, we observe that
(( (*) ())
((F(V a U) < 1" E P(F(V( v), U) S7, S R) P(F(V(Z,))
())
is approximately zero, since for y E R. the clipped score satisfies V( y) = M. implying that P(F(V y), U) S y R) & 0, because that F(u, u) is monotone with respect to its first argument v. Notably, the large constant M is originally introduced conformal selection as a relaxation of infinity, and these approximations indeed hold. This phenomenon is verified through our extra simulations in Appendix C.2.1.
C. Additional Details for Numerical Simulations (Section 5)
C.1. Data Generating Processes and Configuration of the Selection Tasks
The configuration of the true regression function y and noise term e Rd can be found in Table 6.
For the noise term, the degree of freedom for the t-distribution scenario is set to v = 3, and the scale matrix Σ is specified as a matrix with diagonal entries of 0.5 and off-diagonal entries of 0.05. The second column denotes the k-th output entry of the true regression function μ4, which relates to the k-th coordinate yk of the response y and the k-th coordinate e of the noise e as yk = [p(=)] + Gg. In all of our simulations, we take p = 10 (recall that æ Unif(-1, 1)" but y E Rd)). If, when generating y the value of de (the b-th coordinate of a:) for l p needed, we take 3eg = x(((-1) med p)+1-
Table 7 summarizes the values of coefficient vectors (c and r) that define the selection tasks for each response dimension. Within each vector, all entries share the same value. For instance, if is listed as 1, it indicates that = 11, 1).
The six settings we consider differ in two key aspects that influence the difficulty of selection. First, the regression function exhibits varying degrees of nonlinearity across settings. Specifically, Settings and 4 are linear, while Settings 2 and have weak nonlinearity, and Settings and exhibit strong nonlinearity. The degree of nonlinearity affects the predictive accuracy of the estimated regression function á. which in turn impacts the selection difficulty. Second, the distribution of the noise term differs across settings, leading variations in the conditional variance. In Settings 2, and 3, follows multivariate Gaussian distribution, whereas in Settings 4, 5, and 6, it follows multivariate (-distribution. The conditional variance cov(E) = Var(y 2:) in the latter three settings is /ー E = 3E, which is three times larger than that of the former three settings, thereby increasing the difficulty of selection.
To better illustrate the data distributions, Figure 2 presents example scatter plots of the response vector y for the six settings, with the response dimension set for visualization purposes. Higher-dimensional cases are not displayed as they are
16

---

Multivariate Conformal Selection
Table 6: True regression functions and noise distributions
Setting
[n(-)l
f0x+1 Ik+2 +
NO,
2k
2
+ x xx42 +
20)
14 > 07- 1 x > 0.5) (0.25 + xx+2)
3
E)
1(24Ik11 50} 1 v 0.5} (Øb+2 - 0.25) + 0.75
Same Setting
(0,
Same as Setting 2
(2, E)
Same as Setting 3
(0(0 Σ)
Table Coefficients for Selection Task 1 and Ck. (Task 1)
Response Dimension
Cp. (Task 2)
Tk (Task 2)
2
1.5
5
0.2
2.6
10
-0.2
2
4.1
30
0.6
2
7.5
more challenging interpret. Settings within the same column share the same conditional expectation Ely æ), while those in the lower exhibit greater noise and more dispersed scatter patterns due the multivariate I-distributed noise e. The target regions for the two tasks, defined based the coefficients in Table 7, are highlighted red and yellow shaded areas in the plots.
C.2. Extra Simulated Experiments
C.2.1. COMPARING Two DISTANCE DISTANCE-BASED BASED SCORES
In Tables 8 and 9, compare the performance of the two distance-based scores (7) and (8) in dist.. All other configurations are kept the same as in Section 5.
Table 8: Observed FDR of mCS dist with score (7) and (8) for Task and
Task
Task 2
Setting
Score (7)
Clipped Score
Score (7)
Clipped Score (8)
1
0.154
0.277
0.102
0.265
2
0.269
0.314
0.127
0.263
3
0.241
0.265
0.150
0.223
0.202
0.278
0,226
0.286
5
0.282
0.308
0.240
0.292
6
0.265
0.287
0.220
0.283
Across wide range of tasks and experimental settings, the clipped score (8) consistently demonstrates superior performance.
17

---

Multivariate Conformal Selection
Setting 1,
Setting nonlinearty, Gaussian
Setting
linear, Gaussian
strong nonlinearit, Gaussian
8
2
0-
0-
0-
-4
-2
Setting 4, Einear, (-distributed noise
Setting 5,
Setting 6,
weak nonlinearity, d-distributed noise
strong nonlinearity, E-distributed noise
6
0
-4
2
0
2
Figure 2: Scatter plots of 1,000 i.i.d. samples from our data generating processes. Each point represents a sample, with the r- and y-axes corresponding its responses entries y and y2: respectively. The red and the yellow shaded areas represents the two target regions for the two selection tasks.
Table 9: Observed power of mCS dist with score (7) and (8) for Task I and 2.
Task
Task
Setting
Score (7)
Clipped Score (8)
Score (7)
Clipped Score (8)
0.305
0.555
0.440
0.760
2
0.069
0.104
0.208
0.405
3
0.046
0.068
0.075
0.134
4
0.191
0.324
0.222
0.333
5
0.052
0.060
0.093
0.170
0.038
0.046
0.048
0.063
Although both scoring methods guarantee valid FDR control, the clipped score tends to exhibit higher observed FDR. This finding is consistent with the asymptotic analysis in Appendix B.2, which suggests that the clipped score can leverage a larger portion of the available FDR budget.
C.2.2. VARYING RESPONSE DIMENSION
The main article included results with response dimension d 30. In this section, we access the performance of MCCSSSA and mCs-learn in lower dimensional response settings {2, 10), examining whether mCS-lean continues to outperform competing methods for smaller values of d. To simplify the analysis, we focus on Setting 3 for these experiments. All other configurations are kept unchanged.
Tables and 11 show that even for smaller d, mCS-leard and mCS-dist continue to outperform the other baseline methods, achieving the best and the second-ber best power respectively (while controlling FDR). Notably, the performance gap of mcs--ee a over competing methods becomes more pronounced the response dimension increases, suggesting that
18

---

Multivariate Conformal Selection
Table 10: Observed FDR of different methods for lower response dimensions.
Task
Task
CS-int
CS_ib
CS_is
bi
ICSS-d
MCS-1
bi
MCS-d d
mCS
0.244
0.000
0.142
0.303
0.268
0.281
0.270
0.299
0.302
0.706
0.000
0.255
0.290
0.324
0.294
0.246
0.290
0.302
10
0.735
0.000
0.342
0.311
0.276
0.310
0.257
0.265
0.311
30
0.724
0.000
0.204
0.295
0.264
0.278
0.216
0.223
0.263
Table 11: Observed power of different methods for lower response dimensions.
Task
Task
c$_int
c$_ib
CS_is
bf
mCS-d
MCS-1
mCS-d
mCS-1
2
0.029
0.000
0.022
0.032
0.047
0.063
0.050
0.068
0.069
1.000
0.000
0.067
0.041
0.049
0.083
0.050
0.057
0.072
10
1.000
0.000
0.055
0.044
0.058
0.069
0.062
0.069
0.106
30
1.000
0.000
0.039
0.059
0.068
0.102
0.115
0.134
0.179
is particularly advantageous for high-dimensional settings.
C.2.3. NONCONVEX TARGET REGION
To examine the suitability of our methods for tasks involving more irregular target regions, we conduct experiments on two additional tasks which the target region nonconvex:
3. The complement of the (shifted) orthant, { yk 2 C yk]
{y : yk < ck for some k},
4. The complement of sphere centered at c, R {y 11 cl12 r).
To ensure that a reasonable proportion (15%-35%) of responses y fall within the selection region, the coefficients cj and defined differently from their counterparts in Task and 2. Table 12 presents the specific values, where each scalar, following the convention in Table 7, represents vector whose entries are all equal to that scalar. All other setups, configurations, and model hyperparameters are the same as in Section 5.2.
Table 12: Coefficients for Selection Task and
Response Dimension
Ck (Task 3)
Cj (Task 3)
Tk (Task 4)
2
0.5
2
5
0.8
10
1.1
30
1.6
9.5
Table 13 and Table summarize the observed power and FDR for the two additional tasks respectively, where the dimension 30 and nominal FDR level is q 0.3. For the nonconvex tasks, mCS-di (iG-distt itst demonstrates inferior performance compared to the baseline bi. In contrast, -ccs-lea rn performs comparably to bi in relatively easy settings and surpasses bi in more challenging scenarios, such as Settings 3 and 6. Consequently, mCS -learn emerges as the preferred choice when R is nonconvex or otherwise irregular.
C.2.4. OTHER PRETRAINED MODELS fi
We also test the performance of various methods when the underlying model fi is less predictive. To simulate this effect, we train i as a linear model instead of support vector machines in Section All other setups remain unchanged.
19

---

Multivariate Conformal Selection
Table 13: Observed FDR of different methods for Task and 4.
Task
Task
Setting
bi
MCS C
MCS
bi
mCS
MCS-1
0.321
0.267
0.280
0.327
0,298
0.304
2
0.312
0.228
0.332
0.381
0.271
0.286
3
0.334
0.257
0.276
0.385
0.276
0.275
4
0.314
0.328
0.277
0.298
0,298
0.284
0.273
0.266
0.279
0.305
0.288
0.313
0.237
0.255
0.303
0.335
0.263
0.265
Table 14: Observed power of different methods for Task 3 and 4.
Task
Task
Sctting
b1
mCS-d
mC3-1
bi
mCs-d
MCS-1
0.305
0.052
0.255
0.763
0.311
0.555
2
0,008
0.011
0.017
0.098
0.015
0.089
3
0.018
0.005
0.017
0.043
0.012
0.044
0.307
0.050
0.241
0.589
0.233
0.446
0.024
0.001
0.034
0.114
0.017
0.117
0.023
0.001
0.040
0.036
0.016
0.062
Table 15 and Table 16 summarize the FDR and power of different methods when j is fitted linear model. Although overall performance declines for most procedures, ICS-diss st and mCS-learm remain comparatively stronger than the baselines. In particular, mCS-les ea rm appears less sensitive to the imperfect linear thanks its data-adaptive nonconformity score yo, which is learned based on the model μ and the observed data.
Table 15: Observed FDR of different methods for lower response dimensions.
Task
Task
Setting
CS int
ib
bi
mCS-d
mCS-1
bi
mCS-d
mCS-7
1
0.776
0.101
0.266
0.274
0.264
0.267
0.257
0.273
0.276
2
0.798
0.000
0.230
0.269
0.299
0.288
0.245
0.222
0.275
3
0.742
0.000
0,210
0.312
0.235
0.206
0,232
0.223
0.287
4
0.803
0.000
0.267
0.318
0.290
0.282
0.307
0.286
0.260
5
0.813
0.000
0.174
0.306
0.279
0.272
0.273
0.298
0.285
6
0.779
0.000
0.268
0.259
0.306
0.264
0.292
0.267
0.279
Table 16: Observed power of different methods for lower response dimensions.
Task
Task
Setting
CS int
CS ib
CS is
bi
mCS-d
mCS -
bi
mCS-d
mCS-1
1.000
0.204
0.467
0.195
0.550
0.379
0.063
0.842
0.583
2
1.000
0.000
0.024
0.077
0.071
0.089
0.386
0.337
0.419
3
1.000
0.000
0.026
0.066
0.074
0.079
0.129
0.138
0.148
4
1.000
0.000
0.203
0.120
0.309
0.186
0.040
0.408
0.188
5
1.000
0.000
0,016
0.036
0.048
0.051
0.176
0.135
0.162
6
1.000
0.000
0.028
0.042
0.039
0.048
0.056
0.056
0.073
20

---

Multivariate Conformal Selection
C.2.5. COMPARISON OF LOSs FUNCTIONS mCS-learn
Here compare the performance of mCS -learn with loss L1 (15) and loss L2 (16). The response dimension is d 30, and the nominal FDR level is q 0.3.
Table 17: Observed FDR of mCs- with loss I and L2-
Task
Task
Setting
L1
L2
L1
L2
0.255
0.251
0.276
0.279
0.279
0.267
0.278
0.273
0.299
0.278
0.331
0.263
0.274
0.315
0,280
0,254
0.210
0.239
0.278
0.273
6
0.296
0.258
0.361
0.207
Table 18: Observed power of mCS-learn with loss L and Lz.
Task
Task
Setting
L
L
L1
L2
0.040
0.325
0.051
0.534
2
0.024
0.109
0.111
0.421
0.022
0.102
0.041
0.179
0.034
0.199
0.019
0.180
0.009
0.042
0.041
0.189
0.012
0.034
0.032
0.061
As noted in the main article, employing the L2 loss obviates the need for an additional round of smooth ranking, thereby enhancing both numerical stability and training efficacy. This explains the superior power observed in Table 18 when mCS-leard utilizes the L2 loss.
C.2.6. COMPARISON OF SCORE FAMILIES IN mCS learn
When learning fo in mCS -larn, various forms of data may be utilized as inputs, giving rise to different score families. The simplest approach draws exclusively on the covariates æ as features, making it applicable even in scenarios where the model fi is not available. Alternatively, incorporating the response y expands the family (11). When j is available, its predictions (ax)) can also be included as input. Although this does not increase the overall expressiveness of the family, can accelerate the training process. In this section, we evaluate the performance of the following score families with varying complexity:
1. Covariate only: M 1(y R) 0R} fo(æ; R).
2. Prediction only: M 1(y R aR) fo(i((ar); R).
3. Covariate and Prediction: M 1g R (R)R) - fo ji(m);; R).
4. Full Family (12) M 1 E R 0R) - -fe(0, y: R).
5. All available information: M 1 E (R) fo (ee.f(().))))) y: R).
As discussed in the main article, families 4 and incorporate as an input, which compromises the regional monotonicity of the score unless M is sufficiently large. Table 19 summarizes the FDR and power of mCS-leard with different families, with Setting 3. For both Task and Task 2, families 4 and violated FDR control when M = 103. Among the three subfamilies that maintained valid FDR control, family exhibited the best performance.
21

---

Multivariate Conformal Selection
Table 19: Observed FDR and power of mCS learn with different score families.
FDR
Power
Family
Task 1
Task
Task
Task 2
1
0.298
0.237
0.108
0.165
2
0.260
0.291
0.096
0.175
3
0.278
0.263
0.102
0.179
4
0.711
0.594
0.972
0.807
5
0.667
0.494
0.782
0.536
D. Additional Details for Real Data Application (Section 6)
D.1. Overview of Drug Discovery Data and Contiguration of the Selection Tasks
The drug discovery dataset we used in Section 6 is compiled from various public sources (Wenzel et al., 2019; Iwata et al., 2022; Kim al., 2023; Watanabe et al., 2018; Falcón-Cano et al., 2022; Esposito et al., 2020; Braga al., 2015; Aliagas al., 2022; Perryman al., 2020; Meng et al., 2022; Vermeire et al., 2022). Because the integrated data contained missing values, we employed Chemprop (Yang et al., 2019; Heid et al., 2023) impute these entries. The resulting imputed dataset was then used in all subsequent experiments. The processed dataset contains n 22805 data points, and can be found at https: / / thub. com/ Taanm/iinn-aaaa an Ba mcs.
We list the names, units and cutoffs of the responses (only relevant the first selection task) in the imputed dataset in Table 20, and provide detailed descriptions of their biological significance and drug discovery relevance in Figure 4. Recall that in the first task consider target regions of the shape R {y y Z c vk), where the cutoffs defining the selection problem. Figure 3 shows the distribution of these 15 responses, with vertical red lines indicating their corresponding cutoffs. Approximately 21 1% of the test dataset compounds exceed all 15 thresholds, thereby qualifying for selection.
the values cj
For the second task, the target region is defined as sphere {y 1ay - c|l2 S r). For convenience, we take the center the sphere the same as the cutoffs с₁ task 1, and let 2.4. Under this definition, approximately 24% of the compounds qualify for selection. Similarly, for the third task where the target region is the complement of a sphere {y |y 1yy-ch cll2 S r'}, adopt the same center and set r' 3.4. This choice ensures that 18% of the compounds qualify for selection.
Table 20: List responses in the drug discovery dataset.
Name
Unit
Cutoff ce
CL microsome human
log (mL/min/kg)
CL m_mmcrosomeemouss
log 10 (mL/min/kg)
4
[L.microsome.ral
log (mL/min/kg)
CL.total .og
log 10 (mL//in/kg)
0.5
CL total human
log 10 (mL/min/kg)
monkeey
log 10 (mL/min/kg)
0.5
rot
log (mL/min/kg)
CYP2C8 inhibition
log 10 (nMolar IC50)
3.5
CYP2C9 inhibition
log (nMolar IC50)
3.5
BYP2D6.inhibitie
log 10 (nMolar IC50)
3.5
CYP3A4 inhibition
log 10 (nMolar IC50)
3.5
Papp Caco2
log 10 (10-%cm/s)
0.8
Pgp human
log (efflux ratio)
-0.2
bERG binding
log 10 (nMolar IC50)
3.5
LogD pH 7.4
log 10 (M/M)
22

---

Multivariate Conformal Selection
Cl_ microsomee human
hmmm//oosomeemeee mouse
Cl_ www.miccoomeemm
Cl_totalddg
1500
1250
2300
1000
2000
1000
750
500
500
1000
250
250
0
CI total humaa
Cl_ Lattttall_aaat monkey
CI. total rat
CYP2C8 Inhibition
2500
2500
2000
1500
1500
1500
1500
1000
1200
1000
500
Inhibition
Papp.
2000
2000
2000
1500
1500
1000
500
500
human
hERG binding
LogD pH
1500
25.0
1500
1090
1000
500
Figure 3: Histograms of 15 responses the drug dataset. The vertical red lines denote the corresponding cutoffs (for the first task) for each response.
23

---

Multivariate Conformal Selection
Figure 4: Description and drug discovery relevance of the responses.
mimirooooome human:
Intrinsic metabolic clearance human liver microsomes. (Help predict rate liver metabolism in human through the vitro in vivo correlation.)
CL.mit microsomee mouse:
Intrinsic metabolic clearance in mouse liver mierosomes. (Help predict rate of liver metabolism in human through the in vitro in vivo correlation.)
CL microsome-
Intrinsic metabolic clearance in rat liver microsomes. (Help predict rate of liver metabolism in human through the in vitro in vivo correlation.)
CL total dog:
Total body clearance measured in vivo in dogs. (Drug exposure in this species, crucial for determine human dose regimens through translational modeling.)
Co.total human:
Total body clearance measured in vivo in humans (Drug exposure in this species, crucial for determine
human dose regimens through translational modeling.)
CL total monkey:
Total body clearance measured in vivo in monkeys. (Drug exposure in this species, crucial for determine human dose regimens through translational modeling.)
CL.total
Total body clearance measured in vivo in rats. (Drug exposure this species, crucial for determine human dose regimens through translational modeling.)
CYP2C8 inhibition:
Inhibition potential against human CYP2C8 enzyme. (Assessment of risk of drug-drug interactions (DDIs)
involving drugs metabolized by CYP2C8. Important for safety assessment.)
CYP2C9 inhibition:
Inhibition potential against human CYP2C9 enzyme. (Assessment of drug-drug interactions (DDIs)
involving drugs metabolized by CYP2C9 (e.g., warfarin). Important for safety assessment.)
CYP2D6 inhibition:
Inhibition potential against human CYP2D6 enzyme. (Assessment of drug-drug interactions (DDIs) involving drugs metabolized by CYP2D6 (e.g., some antidepressants, beta-blockers). Important for safety assessment.)
CYP3A4 inhibition:
Inhibition potential against human CYP3A4 enzyme. (Assessment of drug-drug interactions (DDIs) involving drugs metabolized by CYP3A4 (a very common pathway). Important for safety assessment.)
Papp_Caco2.
Apparent permeability across Caco-2 cell monolayer. (Potential for oral absorption by crossing the intestinal wall. High Papp suggests good absorption is likely.)
Pgp human:
Interaction with human P-glycoprotein (Pgp) efflux transporter. (Assess if the drug is pumped out of cells by Pgp, potentially limiting oral absorption and brain penetration.)
bERG .indding
Binding/inhibition the hERG potassium channel. (Assessment of risk of cardiac toxicity (QT prolonga- tion, arrhythmias). A critical safety screen; high binding is a major safety concern.)
LogD 7.4:
Logarithm the distribution coefficient at pH 7.4. (Predicts drug's lipophilicity (fat vs. water solubility) under physiological conditions. Influences absorption, distribution, membrane crossing, and CNS penetra- tion.)
24

---

Multivariate Conformal Selection
D.2. Extra Experiments Real Data
D.2.1. NUMERICAL STABILITY DIFFERENT METHOSS
In this section. we test the numerical stability of mCS-dist and experiments where sample data randomly partition data for each iteration, fix the pretrained model ji, the calibration data caal and the test data Dest This is to avoid the variability from data sampling pretraining, and only consider the inherent variance of the methods. Due to the high computation cost of retraining fo, we also keep fo the same across different iterations. We keep all other setups and configurations unchanged as in Section 6.
on real data. To do this, unlike previous
Table 21: Observed standard error of FDRs for different methods with real data.
Task
9
C8.int
CS_ib
CS.is
mCS-d
mCS-1
1
0.3
0.000
0.000
0.248
0.000
0.000
0.000
2
0.3
0.000
0.000
0.002
3
0.3
0.000
0.021
0.004
1
0.5
0.000
0.000
0.032
0.000
0.001
0.004
2
0.5
0,000
0,000
0.002
3
0.5
0.000
0.028
0.005
Table 22: Observed standard error of powers for different methods with real data.
Task
9
int
Cs
bi
mCS-d
mCS-1
0.3
0.000
0.000
0.022
2
0.3
0,000
0,000
0.001
0.3
0.5
0.000
0.056
2
0.5
0,000
0.000
0.004
0.5
0.000
0.011
0.005
For all methods, the standard errors for both FDR and power remain low, with the exception of CS is nominal level of q 0.3. This phenomenon likely arises from the unstable nature of the subroutine nominal level search at lower user- specified nominal levels. The failure of cs is to control the FDR at this low nominal level (Figure 1) further corroborates this interpretation.
25


Reviews:
[{"review_id": "D8VGXFL0a9", "reviewer": "Reviewer_QDhQ", "summary": "This paper proposes multivariate conformal selection (mCS), which extends conformal selection to multi-response settings by introducing the concept of regional monotonicity, generalizing univariate monotonicity, and defining multivariate non-conformity scores. mCS guarantees finite-sample false discovery rate (FDR) control. Two types of non-conformity scores are proposed. The first, mCS-dist, is based on predefined distances. The second, mCS-learn, involves a term that is learned via gradient-based optimization where the sorting has been replaced with soft (differentiable) sorting. Experiments on simulated and real-world datasets show that mCS enhances selection power while maintaining FDR control.\n\n## update after rebuttal\nThe new experiments regarding hyperparameters, numerical stability, and calib/test data sizes are convincing. I recommend the paper's acceptance.", "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": "> The method requires an extra split... in higher dim.\n\n**A1**: From our experience, the (multivariate) conformal selection procedure is insensitive to the size of calibration data - the calibration scores only affect the resolution of the p-values, which is typically sufficient when $|D\\_{cal}| \\geq 100$. \n\nTo support this point, we ran an additional experiment using $\\texttt{mCS-learn}$ with real-data Task 2 at  $q=0.5$. We set $|D\\_{f-train}|=6400, |D\\_{f-val}|=800$ and $|D\\_{test}|=200$, respectively, and we consider different $|D'\\_{cal}|$ shown below:\n\n|$n\\_{cal}$|FDR|power|\n|-:|-:|-:|\n|100|0.486|0.594|\n|500|0.501|0.603|\n|3000|0.499|0.598|\n\nAs shown, even with a small calibration set, the method maintains FDR control and strong power.\n\n> Why is a further split of the training set necessary...\n\n**A2**: In short, this is because the computation of smoothed p-values require both calibration and test data; therefore to train $f_\\theta$ based on the performance of a certain score, we would need to split the training set into two, which serve as calibration data and test data respectively. \n\n> What was the exact Algo used for soft ranking...typically very sensitive...explain\n\n**A3**: We adopted the implementation in Blondel et al. (2020), with $\\ell_2$ regularization and strength set to 0.1. \nPreliminary tests showed that within a reasonable range, different regularization strengths produce similar overall behavior in mCS. While cross-validation could optimize this parameter further, its exact value is not central to our primary contributions. To maintain clarity and simplicity, we chose not to include additional tuning steps for this hyperparameter.\n\n> The 2nd loss introduces another hyperparameter $\\gamma$...\n\n**A4**: This value was chosen based on a series of preliminary experiments. Similar to the previous question, while we could cross-validate on $\\gamma$, for the simplicity of our presentation, we fix $\\gamma$ to 0.5 in our main procedure.\n\n> In Algo 2, line 8, what is $k$...\n\n**A5**: To accurately estimate the power of the selection rule given by $V^{\\theta}$ (defined as an expectation in Eq. 2), we average empirical selection power over multiple random partitions of $D\\_{f-val}$. Each random partition produces different power values, even with the same $V^{\\theta}(x,y)$. We use $k=100$ partitions, which provides a stable and accurate approximation. We clarified this point in the updated version of our manuscript.\n\n> For the simulated study... bigger test set than 100.\n\n**A6**: We appreciate the suggestion. While a test size of 100 may seem small for general machine learning tasks' evaluation, conformal selection performance is generally insensitive to test set size (as noted in Sec 3.1 of Jin and Candes, 2023). Additionally, we repeated our experiments over 100 independently generated datasets, ensuring the results are stable. This stability is also supported by our simulation results provided in response to subsequent comments.\n\n> Please report the std. errors...\n\n**A7**: For both simulation and real-data experiments, each iteration involves sampling a new dataset, introducing variation in data and model training. Reporting std errors under these conditions would reflect primarily data variability rather than the inherent stability of the selection methods.  \n\nHowever, we agree numerical stability is important. Thus, in Appendix D, we provide an additional real-data experiment with a fixed dataset across iterations. Below are the observed std dev. of FDR and power (over 100 iterations) from this experiment: \n\nFDR/Power:\n\n|Task|$q$|CS\\\\_int|CS\\\\_ib|CS\\\\_is|bi|mCS-d|mCS-l|\n|-:|-:|-:|-:|-:|-:|-:|-:|\n|1|0.3|0.000/0.000|0.000/0.000|0.248/0.022|0.000/0.000|0.000/0.000|0.000/0.000|\n|2|0.3|-|-|-|0.000/0.000|0.000/0.000|0.002/0.001|\n|3|0.3|-|-|-|0.000/0.000|0.021/0.005|0.004/0.002|\n\nThese results confirm strong numerical stability for our proposed methods. Note that for $\\texttt{mCS-learn}$, the trained model $f_\\theta$ was fixed, so variability from model training is not considered here.\n\n> The choice of $r_{n+j}$...?\n\n**A8**: We have added a detailed discussion about the choice of $r\\_{n+j}$, please refer to our response under A5 for Reviewer 2 (PjuQ). \n\nFor both the generalized signed score and the clipped score, choosing $r_{n+j}$ on the boundary is already optimal. In this case, the first term $D_1$ is minimized (equal to 0), so no further optimization is needed to enhance selection power.\n\n> How ... to adapt when ... outputs a conditional distribution $\\hat{P}(Y|X)$?\n \n**A9**:  When the model $\\hat\\mu$ outputs an estimated conditional distribution $\\widehat{P}(y|x)$, the second term $D_2$ can be replaced by the predicted probability of being in the target region: $y \\in R$, i.e. $\\int \\boldsymbol{1}\\\\{\ny\\in R\\\\} d\\widehat{P}(y|x)$. This serves the same purpose: points with high predicted probability of satisfying the selection criterion will receive lower scores and are more likely to be selected.", "claims_and_evidence": "- The extension of conformal selection to the multi-response setting is a novel and significant contribution.\n- The framework guarantees finite-sample FDR control. The theoretical arguments seem correct.\n- The proposed non-conformity scores are well-motivated.\n- There were some ambiguities in the definition and implementation of the `mCS-learn` algorithm for me to assess it fully (please see comments in \"Methods And Evaluation Criteria\" section)\n- Overall, the clarity of the presentation could be significantly improved. I suggest moving up some of the experimental details (such as the definition of settings, tasks) to the main text so that the figures and tables can be interpreted from information given in the main text. \n- Both `mCS-dist` and `mCS-learn` both perform well against baselines in the simulated setting. All of the figures and tables in the main text are on the simulated setting, however, and not the real data application.\n- If the authors are planning to make the aggregated and imputed ADMET dataset public, that would be a significant contribution worth mentioning up front.", "methods_and_evaluation_criteria": "I had multiple questions and concerns about the form of the `mCS-learn` method as well as its presentation:\n- The method requires an extra split of the calibration set into training, validation, and proper calibration sets. I would like to see more discussion on the signal lost from a reduced calibration set size, which I'd presume would be more significant in higher dimensions (many targets).\n- Why is a further split of the training set necessary in line 3 of Algorithm 2?\n- What was the exact algorithm used for soft ranking? The paper cites both Blondel et al. 2020 and Cuturi et al. 2019, but does not indicate the exact algorithm nor the hyperparameters it would introduce. \n- Related to the above, soft sorting algorithms are typically very sensitive to the regularization parameter and small numerical errors may impede training. From Algorithm 2, it seems like only $\\theta$ was chosen via a held-out validation split and the sorting regularization wasn't. If it was predetermined, please explain the procedure. \n- The second loss function (Eq 16) introduces another hyperparameter $\\gamma$. How was this determined?\n- In Algorithm 2, line 8, what is $k$? Why does repeated application of mCS on $\\mathcal{D}_{f\\textrm{-val}}$ lead to different power values?\n- Section 5.2: what was $p$?\n- It is very difficult to interpret the tables of results. The \"settings\" are only defined in the Appendix -- it would also help to refer to them by their respective characteristics in the main text (for example, setting 1 = \"linear, Gaussian noise\").  Please also indicate the nominal level in the table caption.", "theoretical_claims": "The proof of Theorem 3.5 seems correct.", "experimental_designs_or_analyses": "- For the simulated study, as there is full control over the data generation, the authors can use a bigger test set than 100 for cleaner evaluation\n- Please also report the standard errors of the metrics across runs in the tables.\n- Please share the anonymized code so the details of the implementation can be checked", "relation_to_broader_scientific_literature": "- Extending CS to multi-response settings is a significant contribution. The paper introduces the concept of regional monotonicity, which generalizes univariate monotonicity, required to guarantee FDR control. \n- Even in the univariate setting, regional monotonicity generalizes the threshold-based framework of Jin and Candes 2023 to arbitrary sets of intervals.\n- The authors provide finite‐sample theoretical guarantees (Theorem 3.5) using similar arguments as Jin and Candes 2023.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "- Algorithm 1, line 1: $r_j$ is confusing as it reuses $j$ index --> $r_{n+j} \\in R$\n- Please explain the role of $M$ up front when it's introduced, right after Eq 8 and also Eq 11\n- In Algorithm 2, please include a line that assigns the value to $\\theta_t$ (currently introduced without definition)\n- In Algorithm 2, $t^*$ is overloaded with the previous discussion of Theorem 4.1\n- Fig 1: Please make the figure axis labels and legend font larger\n- Fig 1: Next to task number, please include a few-word description of what the task is designed for", "questions_for_authors": "- The choice of $r_{n+j}$ (choosing it on the boundary of $R$ that would be most informative for the given $x$) also seem to offer room to improve selection power, as long as the choice doesn't violate exchangeability. Is this correct? If so, did the authors consider simply optimizing for this within the `mCS-dist` variant?\n- Why is it necessary to further split $\\mathcal{D}_{f\\textrm{-train}}$ into two?\n- How does the method have to adapt when the underlying predictor outputs a conditional distribution $P(Y|X)$ rather than a point estimate?", "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": "gmmgZdrv9S", "reviewer": "Reviewer_PjuQ", "summary": "This paper addresses the important problem of multivariate conformal selection. While multivariate conformal prediction is relatively well studied, this appears to be the first work on multivariate selection tasks.\n\n\n## update after rebuttal: I did not change my score as the recommendation is already to accept.", "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": "> It would be good to see an evaluation on classification outcomes as well.\n\n**A1**: \nIn our paper, we chose to focus primarily on regression tasks because they represent a more challenging setting for conformal selection.\n\nIn fact, the selection problem for classification (univariate or multivariate) can be directly reduced to the univariate conformal selection framework introduced by Jin and Candes (2023). This reduction explains why we did not include separate evaluations specifically for classification scenarios, as such evaluations would essentially revisit methods already established in the literature.\n\nTo clarify this point, let us briefly demonstrate the reduction explicitly: \n\n* For the univariate classification setting, suppose the response space is composed of classes $\\mathcal{Y} = \\cup_{k=1}^K C_k$ with target region $R = \\cup_{k=1}^{s} C_k$ (with $s < K$). Then, by defining a binary response: $\\tilde{y}\\_{i} = \\boldsymbol{1}\\\\{y\\_{i} \\in R \\\\}$, the original selection problem directly translates into a univariate conformal selection problem, where we select samples with $\\tilde{y}\\_{i} = 1$.\n\n* For multivariate classification case, e.g.  suppose responses $\\boldsymbol{y}\\_{i}$ are drawn from a joint class space $\\mathcal{Y} = (\\mathcal{Y}^{(1)},\\mathcal{Y}^{(2)}) = \\cup_{k,\\ell} (C^{(1)}\\_k, C^{(2)}\\_\\ell)$, and the target region $R = \\cup_{(k, \\ell) \\in \\mathcal{I}}(C^{(1)}\\_k, C^{(2)}\\_\\ell)$. Again, we define a binary indicator $\\tilde{y}\\_{i} = \\boldsymbol{1}\\\\{\\boldsymbol{y}\\_{i} \\in R \\\\},$ \nconverting the original multivariate selection problem into a standard univariate selection task: \n$$H_{0j}: \\tilde{y}\\_{n+j} < 0.5 \\text{\\quad versus \\quad } H\\_{1j}: \\tilde{y}\\_{n+j} \\geq 0.5.$$  \n\nSince $P(\\tilde{y}\\_{i} = 1)=P(\\boldsymbol{y}\\_{i} \\in R)$, there is a direct correspondence between the multivariate and univariate nonconformity scores:\n$\nV(\\boldsymbol{x},\\boldsymbol{y}) = M\\cdot\\boldsymbol{1}\\\\{\\boldsymbol{y}\\in R\\\\} - \\widehat{P}(\\boldsymbol{y}\\in R|\\boldsymbol{x}) \n$\nand\n$\nV(\\boldsymbol{x},\\tilde{y}) = M\\cdot\\boldsymbol{1}\\\\{\\tilde{y}\\geq 0.5\\\\} - \\tilde{\\mu}(\\boldsymbol{x}) \n$\nwhere $\\tilde{\\mu}(\\boldsymbol{x})\\equiv \\widehat{P}(\\tilde{y}=1|\\boldsymbol{x})$. Moreover, regional monotonicity is simply the usual monotone condition of univariate conformal selection:\n$V(\\boldsymbol{x},\\tilde{y}=0)\\leq V(\\boldsymbol{x},\\tilde{y}=1).$\n\nThus, every classification-based selection task can be naturally and effectively solved using existing univariate conformal selection methods. We will clarify this point explicitly in our final manuscript.\n \n\n> The following paper seems also related \"Klein et al. (2025)\".\n\n**A2**: We have now cited this paper.\n\n> Figure 1 is too small to read properly. \n\n**A3**: We have now increased the font size and improved the clarity of Figure 1.\n\n> Could you expand more about tie breaking? ...\n\n**A4**: In our work, we follow the standard practice in conformal methods and break ties randomly.\nIf instead, one used a deterministic worst-case tie-breaking rule-always ranking the test score behind tied calibration scores-the conformal p-value becomes:\n$\np_j^{dtm} = \\frac{1}{n+1}(1+\\sum_{i=1}^n \\boldsymbol{1}\\\\{V_i \\leq \\widehat{V}_{n+j}\\\\}).\n$\nThis rule ensures that $p_j^{dtm} \\geq p_j$ for every test sample, making the test deterministic and conservative (though not uniformly distributed). Consequently, applying the BH procedure to these deterministic p-values yields a fully reproducible selection rule, but with a slight reduction in statistical power compared to random tie-breaking. We will clarify this trade-off explicitly in the final manuscript.\n\n> Footnote 2: The choice of $r_{n+j}$ needs to be exchangeable in which sense? ....\n\n**A5**: In the original manuscript, footnote 2 mistakenly suggested that the choice of $r_{n+j}$ required some form of exchangeability. After careful reconsideration, we see clearly now that no such condition is needed. The point $r_{n+j}$ can be chosen arbitrarily from the region $R$. \nFrom the viewpoint of our proof, the critical step is the regional monotonicity condition. This condition alone ensures the relationship $p_{j}^\\* \\leq p_j$ under the null hypothesis. The oracle p-values $p_{j}^{*}$ are uniform precisely because of the exchangeability of the calibration data. The choice of \n$r_{n+j}$, therefore, does not influence the oracle p-value and imposes no extra exchangeability constraints.\nWe have revised the manuscript accordingly to clearly state this simpler and accurate assumption.\n\n> How would one define regional monotonicity for a classification problem?\n\n**A6**: Please see Q1. \n\n> In Theorem 4.1 why do you need to assume i.i.d.?...\n\n**A7**: The proof of Theorem 4.1 indeed relies crucially on the strong law of large numbers. This result requires independence (or at least certain mixing conditions). Exchangeability alone is not sufficient here. Please see more detail in Appendix B.", "claims_and_evidence": "The claims are supported by proofs and simulation studies.", "methods_and_evaluation_criteria": "As this is the first paper on multivariate conformal selection there does not seem to be any benchmark data sets available. \n\nIt would be good to see an evaluation on classification / binary outcomes as well.", "theoretical_claims": "The proofs look plausible.", "experimental_designs_or_analyses": "The experimental setup is well explained in the supplementary material.", "relation_to_broader_scientific_literature": "The following paper seems also related\n\nKlein, Michal, et al. \"Multivariate Conformal Prediction using Optimal Transport.\" arXiv preprint arXiv:2502.03609 (2025).", "essential_references_not_discussed": "None came to mind", "other_strengths_and_weaknesses": "Figure 1 is too small to read properly. Perhaps move Task 2 to the supplementary?", "other_comments_or_suggestions": "See above", "questions_for_authors": "Could you expand more about tie breaking? What if one used a worst-case rule for ties instead of breaking them at random, so that one has a deterministic test procedure? \n\nFootnote 2: The choice of $r_[n+j}$ needs to be exchangeable in which sense? Regarding its coordinates, or ensuring exchangeability with the calibration data? \n\nHow would one define regional monotonicity for a classification problem? \n\nIn Theorem 4.1 why do you need to assume i.i.d.; would exchangeability suffice?", "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": "b3mtcYvC9u", "reviewer": "Reviewer_1RGn", "summary": "This paper introduces Multivariate Conformal Selection (mCS), a new approach for selecting candidates in settings with multivariate responses, such as drug discovery and large language model alignment. Unlike traditional Conformal Selection, which is limited to single-variable outputs, mCS extends the framework by leveraging multivariate nonconformity scores and ensuring False Discovery Rate (FDR) control. The authors propose two variants: mCS-dist, which uses distance-based scores, and mCS-learn, which optimizes selection criteria through differentiable learning.", "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": "> On page 12, just before \"The last equality is again by\", $p\\_{j}^\\* \\in \\mathcal S\\_{j}^\\*$ should be replaced with $j \\in \\mathcal S\\_{j}^*$... since the p-value corresponding to $j$ in $S^{\\*}\\_{j\\to0}$ is 0.\n\n**A1**: Thank you very much for pointing out the typo and the error in our proof. We have corrected the typo accordingly.\n\nLet us clarify the corrected reasoning clearly. In the original manuscript, we aimed to show:\n$$\n\\text{FDR} \\leq \\sum\\_{j=1}^m \\mathbb{E}\\bigg[ \\frac{ \\boldsymbol{1} \\\\{p\\_j^{\\*} \\leq q|\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}| / m \\\\}}{\\max(1, |\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}|)} \\bigg]\n$$\nwhich is correct, yet we used a wrong path to prove this inequality. We claimed:\n$$\n\\text{FDR} \\leq \\sum\\_{j=1}^{m} \\sum\\_{k=1}^{m} \\frac{1}{k} \\mathbb{E}\\big[ \\boldsymbol{1}\\\\{|\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}| = k\\\\} \\boldsymbol{1}\\\\{j \\in \\mathcal{S}\\_{j \\to 0}^{\\*} \\\\} \\big] \n= \\sum\\_{j=1}^m \\mathbb{E}\\bigg[ \\frac{ \\boldsymbol{1} \\\\{p\\_j^{\\*} \\leq q|\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}| / m \\\\}}{\\max(1, |\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}|)} \\bigg].\n$$\nwhich is incorrect, since as pointed out, $\\mathcal{S}\\_{j \\rightarrow 0}^{\\*} := \\mathcal{S}(p\\_{1}^{(j)}, \\dots, p\\_{j-1}^{(j)}, 0, p\\_{j+1}^{(j)}, \\dots, p\\_m^{(j)})$'s $j$-th argument is 0, not $p\\_{j}^{\\*}$.\nHowever, consider the following chain of inequalities:\n$$\n\\text{FDR} \\leq \\sum\\_{j=1}^{m} \\sum\\_{k=1}^{m} \\frac{1}{k} \\mathbb{E}\\big[ \\boldsymbol{1}\\\\{|\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}| = k\\\\} \\boldsymbol{1}\\\\{j \\in \\mathcal{S}\\_{j}^{\\*} \\\\} \\big] \n= \\sum\\_{j=1}^m \\mathbb{E}\\bigg[ \\frac{ \\boldsymbol{1} \\\\{p\\_j^{\\*} \\leq q|\\mathcal{S}\\_{j}^{\\*}| / m \\\\} }{\\max(1, |\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}|)} \\bigg] \\leq \\sum\\_{j=1}^m \\mathbb{E}\\bigg[ \\frac{ \\boldsymbol{1} \\\\{p\\_j^{\\*} \\leq q|\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}| / m \\\\}}{\\max(1, |\\mathcal{S}\\_{j \\rightarrow 0}^{\\*}|)} \\bigg].\n$$\nThis corrected argument allows us to properly continue the proof as originally intended. We have updated the manuscript to clearly reflect this correction.\n\n\n> In the proof of Theorem 3.5, could the authors comment on this sentence... but I was not able to understand it.\n\n**A2**: Thank you for pointing this out. In fact, the assumption that $V\\_1, \\dots, V\\_n, \\widehat{V}\\_{n+1}, \\dots, \\widehat{V}\\_{n+m}$ have no ties is unnecessary. We can simply handle the case where $\\widehat{V}\\_{n+l} = \\widehat{V}\\_{n+j}$ by combining it directly with the scenario originally labeled (i), where $\\widehat{V}\\_{n+l} > \\widehat{V}\\_{n+j}$. After this adjustment, the proof proceeds without difficulty.\n\nWe have clarified this point and updated the manuscript accordingly.\n\n> In Section A.1, it would be good to provide a reference when it is stated: \"By a standard result from conformal inference...\"\n\n**A3**: We added a reference (Vovk et al., 2005) to corroborate the claim.\n\n> Other Comments Or Suggestions\n\n**A4**: We have corrected all the typos and addressed all the minor suggestions. We highly appreciate your careful review, which has helped us improve the clarity and quality of our manuscript.", "claims_and_evidence": "Through both simulated and real-world experiments, the study demonstrates that mCS enhances selection accuracy while rigorously controlling FDR. These numerical results support the thereotical claim of the paper.\n\nI have some concerns on the proof of one of the main results of the paper (please refer to the section Theoretical Claims).", "methods_and_evaluation_criteria": "Proposed methods and/or evaluation criteria  make sense for the problem at hand.", "theoretical_claims": "I would be graeful is the authors could clarify the following aspects on the proof of Theorem 3.5:\n\n- On page 12, just before \"The last equality is again by\", $p^*_j \\in \\mathcal S^*_j$ should be $j \\in \\mathcal S^*_j$ (same later in the sentence and in the next equation). With these corrections, I seems to me that the equality before \"By definition, {p_l }_{l \\neq j} is invariant after...\" is not true anymore. Indeed, the event \"j \\in S^*_{j\\to0}\" is equal to \" 0 \\leq q|S^*_{j\\to0}/m\" since the p-value corresponding to j in S^*_{j\\to0} is 0.\n\n- In the proof of Theorem 3.5, could the authors comment on this sentence: \"Above two cases happen with probability 1 since there are no ties almost surely.\" ? I was not clear to me why there are no ties almost surely since V. It might be related to the fact that the proof is derived considering diterministic p-values (and thus the terms in Eq.(3) for the tie-breaking of the nonconformity scores are not there), but I was not able to understand it.", "experimental_designs_or_analyses": "The authors provide analysis of their method considering different dimensions of the output space, different properties for the target region and discussed results on both simulated and real data. The empirical evaluation is solid.\n\nI would have been interested to have a more concrete/detailed description of the biological meaning the output space for the application on real data.", "relation_to_broader_scientific_literature": "The paper extends the work from Jin & Candès \"Selection by prediction with conformal p-values\" to the multivariate response setting. Some key concepts originally introduced by  Jin & Candès are adapted to the multivariate setting.", "essential_references_not_discussed": "Relevant related works are cited as far as I am aware.\n\nThere is another line of work proposing conformal methods to control the FDR on the set of true edges in graph when considering a link prediction problem (such as \"Conformal link prediction for false discovery rate control\", Ariane Marandon). These works are not cited by the authors but - while being connected - they adress a different task than the one considered in the paper and thus are not essential.", "other_strengths_and_weaknesses": "-", "other_comments_or_suggestions": "Here are some typos:\n\n- \"but replaces the the denominator\" (sec. 2)\n\n- In algo 1, at line 1, r_j should be r_{n+j}.\n\n- on the left side of page 5, F (V (x, y), U ) ~ U(0,1) and not (-1,1)", "questions_for_authors": "I would be grateful if the authors could clarify my question on the proof of Theorem 3.5.", "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": "None"}}]

Author rebuttal:

Meta-review: This paper introduces Multivariate Conformal Selection (mCS), extending Conformal Selection to multivariate response settings. By leveraging regional monotonicity and multivariate nonconformity scores, mCS ensures finite-sample False Discovery Rate control. Two variants—distance-based and learned scores—demonstrate strong performance in simulations and real-world tasks, improving selection power while maintaining rigorous uncertainty guarantees. All reviewers provided positive feedback. The concerns raised were mostly addressed in the rebuttal.
Decision: Accept (poster)