Fast Nonlinear Vector Quantile Regression
Authors: Aviv A. Rosenberg
Sanketh Vedula
Yaniv Romano
Alexander Bronstein
Venue: ICLR 2023 poster (2023)
Keywords: optimal transport
quantile regression
vector quantiles
uncertainty quantification
multi-output regression
conformal prediction
software

Abstract: $$
\newcommand{\rvar}[1]{\mathrm {#1}}
\newcommand{\rvec}[1]{\boldsymbol{\mathrm{#1}}}
$$
Quantile regression (QR) is a powerful tool for estimating one or more conditional quantiles of a target variable $\rvar{Y}$ given explanatory features $\rvec{X}$.
A limitation of QR is that it is only defined for scalar target variables, due to the formulation of its objective function, and since the notion of quantiles has no standard definition for multivariate distributions.
Recently, vector quantile regression (VQR) was proposed as an extension of QR for vector-valued target variables, thanks to a meaningful generalization of the notion of quantiles to multivariate distributions via optimal transport.
Despite its elegance, VQR is arguably not applicable in practice due to several limitations:
(i) it assumes a linear model for the quantiles of the target $\rvec{Y}$ given the features $\rvec{X}$;
(ii) its exact formulation is intractable even for modestly-sized problems in terms of target dimensions, number of regressed quantile levels, or number of features, and its relaxed dual formulation may violate the monotonicity of the estimated quantiles;
(iii) no fast or scalable solvers for VQR currently exist.

In this work we fully address these limitations, namely:
(i) We extend VQR to the non-linear case, showing substantial improvement over linear VQR;
(ii) We propose {vector monotone rearrangement}, a method which ensures the quantile functions estimated by VQR are monotone functions;
(iii) We provide fast, GPU-accelerated solvers for linear and nonlinear VQR which maintain a fixed memory footprint, and demonstrate that they scale to millions of samples and thousands of quantile levels;
(iv) We release an optimized python package of our solvers as to widespread the use of VQR in real-world applications.

Published as a conference paper at ICLR 2023
FAST NONLINEAR VECTOR QUANTILE REGRESSION
Aviv A. Rosenberg! 13t, Sanketh Vedula Yaniv Romano 12, and Alex M. Bronstein 1.3
'Deparmmen of Computer Science, Technion
e Department of Electrical and Computer Engineering, Technion
Sibylla, UK
(avivr, sanketh, yr romano, bron (cs. technien. Equal Contribution
ABSTRACT
Quantile regression (QR) is powerful tool for estimating one or more conditional quantiles of a target variable Y given explanatory features X. A limitation of QR is that it is only defined for scalar target variables, due to the formulation its objective function, and since the notion of quantiles has standard definition for multivariate distributions. Recently, vector quantile regression (VQR) was proposed
extension of QR for vector-valued target variables, thanks to a meaningful generalization the notion of quantiles to multivariate distributions via optimal transport. Despite its elegance, VQR is arguably not applicable in practice due to several limitations: (i) it assumes a linear model for the quantiles of the target Y given the features X; (ii) its exact formulation is intractable even for modestly-sized problems in terms of target dimensions, number of regressed quantile levels. or number of features, and its relaxed dual formulation may violate the monotonicity of the estimated quantiles; (iii) no fast or scalable solvers for VQR currently exist. this work we fully address these limitations, namely: (i) We extend VQR to the non-linear case, showing substantial improvement over linear VQR; (ii) We propose vector monotone rearrangement, a method which ensures the quantile functions estimated by VQR are monotone functions; (iii) We provide fast, GPU-accelerated solvers for linear and nonlinear VQR which maintain fixed memory footprint, and demonstrate that they scale to millions of samples and thousands of quantile levels; (iv) We release an optimized python package of our solvers to widespread the use of VQR in real-world applications.
INTRODUCTION
Quantile regression (QR) (Koenker Bassett 1978) is a well-known method which estimates a conditional quantile of target variable Y. given covariates X. A major limitation of QR is that deals with a scalar-valued target variable, while many important applications require estimation of rector-valued responses. A trivial approach is estimate conditional quantiles separately for each component of the vector-valued target. However this assumes statistical independence between targets, a very strong assumption rarely held in practice. Extending QR to high dimensional responses is not straightforward because (i) the notion of quantiles is not trivial to define for high dimensional variables, and in fact multiple definitions of multivariate quantiles exist (Carlier et al. 2016); (ii) quantile regression is performed by minimizing the pinball loss, which is not defined for high dimensional responses.
Seminal works of Carlier et al. (2016) and Chernozhukov et al. (017) introduced notion of quantiles for vector-values random variables, termed vector quantiles. Key their approach is extending the notions of monotonicity and strong representation of scalar quantile functions to high dimensions, i.e.
Co-monotoniciry: (Qy(u) - Qx)))) T (u u') 2 0, Vu, u' [0. 1]a Strong representation: Y (U),
(1)
U'O,1]" 1]dd
(2)
where Y is d-dimensional variable, and Qy [0, 11dh + Rd is its vector quantile function (VQF).
Moreover, Carlier et al. (016)) extended QR to veluor targets, which leads to vector quantile regression (VQR). VQR estimates the conditional vector quantile function (CVQF) Q from

---

Published as a conference paper at ICLR 2023
Wh
(b)
Figure I: (a) Visualization of the vector quantile function (VQF) and its a-contours, high- dimensional generalization of a-confidence intervals. Data was drawn from 2d star-shaped distribution. Vector quantiles (colored dots) overlaid on the data (middle). Different colors correspond to a-contours, each containing 100 (1 2a percent of the data. The VQF Qv(u)= (Q)(u), Q2(u)] is co-monotonic with u (u, ug); Q1. . Q2 are depicted as surfaces (left, right) with the corresponding vector quantiles overlaid. On Q1, increasing i1 for a fixed 22 produces a monotonically increasing curve, and vice versa for Q2- (b) Visualization of conditional vector quantile functions (CVQFs) via a-contours. Data was drawn from joint distribution of (X, Y) where Y X æ has star tar-shaped distribution rotated by z degrees. The true CVQF Qvix changes non- oon-linearly with the covariates X, while E Y [ remains the same. This demonstrates the challenge of estimating CVQFs from samples of the joint distribution. Appendix [ppendix[[providdee provvides further intuitions regarding VQFs and CVQFs, and details how the a-contours are constructed from them.
samples drawn from P(X,Y) Y), where Y is a d-dimensional target variable and X are k-dimensional covariates. They show that function Q which obeys co-monotonicity ( and strong represen- tation 2 exists and is unique, as a consequence of Brenier's polar factorization theorem (Brenier)
1991). Figure provides a visualization of these notions for a two-dimensional target variable.
Assuming linear specification Qyix(u; a) B(u) T x + a(u). VQR can be formulated as an optimal transport problem between the measures of YJX and U, with the additional mean- independence constraint E [U] = E [X] The primal formulation of this problem is a large-scale linear program and is thus intractable for modestly-sizen problems. A relaxed dual formulation which is amenable to gradient- based solvers exists but leads to co-monotonicity violations.
The first goal of our work is to address the following limitations of Carlier et al. (2016) 2020): (i) the linear specification assumption on the CVQF, and (ii) the violation of co-monotonicity when solving the inexact formulation of the VQR problem. The second goal of this work is to make VQR an accessible tool for off-the-shelf usage on large-scale high-dimensions there are available software packages to estimate VQFs and CVQFs that can scale beyond toy problems. We aim to provide accurate, fast and distribution-free estimation of these fundamental statistical quantities. This is relevant for innumerable applications requiring statistical inference, such as distribution-fre. uncertainty estimation for vector-valued variables (Feldman et al. 2021) 1 hypothesis testing with vector-valued test statistic (Shi et al. 2022), causal inference with multiple interventions (Williams & Crespi 2020), outlier detection (Zhao et al. 2019) and others. Below we list our contributions.
datasets. Currently
Scalable VQR. We introduce a highly-scalable solver for VQR in Section 3 Our approach, inspired by Genevay et al. (2016) and Carlier et al. (2020), relies on solving new relaxed dual formulation of solvers which maintain a constant memory footprint regardless of problem size. We demonstrate that our approach scales to millions of samples and thousands of quantile levels and allows for GPU-acceleration.
the VQR problem. We propose custom stochastic- wwwnstti-graiient--aa
Nonlinear VQR.
To address the limitation of linear specification, in Section 4 we propose nonlinear vector quantile regression (NL-VQR). The key idea is fit nonlinear embedding function of the input features jointly with the regression coefficients. This is made possible by leveraging the relaxed dual formulation and solver introduced in Section 3 We demonstrate, through synthetic and real-data experiments, that nonlinear VQR can model complex conditional quantile functions substantially better than linear VQR and separable QR approaches.
Vector monotone rearrangement (VMR). In Section 5 we propose VMR, which resolves the co-monotonicity violations in estimated CVQFs. We solve optimal transport problem to rearrange
2

---

Published as a conference paper at ICLR 2023
the vector quantiles such that they satisfy co-monotonicity. We show that VMR strictly improves the estimation quality of CVQFs while ensuring zero co co-monotonieity violations.
Open-source software package. We release feature-rich, well-tested python package vqu" implementing estimation of vector quantiles, vector ranks, vector quantile contours, linear and nonlinear VQR, and VMR. To the best of our knowledge, this would be the first publicly available tool for estimating conditional vector quantile functions at scale.
2 BACKGROUND
Below we introduce quantile regression, its vector-valued targets. We aim to provide brief and self-contained introduction to the key building blocks of VQR. For simplicity, we present the discretized versions of the problems.
based formulation and the extension to
Notation. Throughout, Y, X denote random variables and vectors, respectively; deterministic scalars, vectors and matrices are denoted as y. æ, and X. P(x,Y) denotes the joint distribution of the X and Y. 1N denotes an N-dimensional vector of ones, O denotes the elementwise product, and II f.} is an indicator. We denote by N the number of samples, d the dimension of the target variable, k the dimension the covariates, and the number of vector quantile levels per target dimension. Qy|x (u; æ)) is the CVQF of the variable Y|X evaluated at the vector quantile level u for X xt
Quantile Regression. The goal of QR is to estimate the quantile of a variable YJX. Assum- ing a linear model of the latter, and given u (0, 1), QR amounts to solving (bu, a.)) arg ming E [Pm(Y x x)]. where b, J On is the u-th quantile of Y|X æ. and (z(2), known the pinball loss, is given by pu(2) max max(0, =} + (u 1)z. Solving this problem produces an estimate of Qyix for a single quantile level u. In order to estimate the full conditional quantile function (CQF) Qyix(u), the problem must be solved at all levels of u with additional monotonicity constraints, since the quantile function is non-decreasing in The CQF discretized at quantile levels can be estimated from N samples (æt, 3i P(x,Y) by solving
N
Σ
Σ 1=1
Pulys - by T x; au)
min B,a
s.t. Vi,
b x.
Tx: + Ou b.. + aus
where B and a aggregate all the b, and Gu, respectively. We refer to eq. (3) as simultaneous linear quantile regression (SIQR). This problem is undefined for vector-valued Y, due the inherently 1d formulation of the monotonicity constraints and the pinball loss (z))..
Optimal Transport Formulation. Carlier et (2016) showed that SLQR 33 can be equivalently written as an optimal transport (OT) problem between the target variable and the quantile levels, with an additional constraint of mean independence. Given N data samples arranged as y E RN X and T quantile levels denoted by = [:: T ......I T we can write.
Ig N xk
max 1120
Hy
s.t. II 1
(4)
II = 11
IIX = X
[3]
where II the transport plan between quantile levels u and samples (x, y). with marginal constraints
1 IN. / 1 Ir and mean-independence constraint X = + 17 1 IN T X. The dual variables &p = T a and B D-TT-- B, where D T is first-order finite differences matrix, and R. B E Tx contain the regression coefficients for all quantile levels, Refer to appendices A.1 and A.2
for full derivation of the connection between SLQR 3 and OT
Vector Quantile Regression. Although the OT formulation for SLQR is specified between Id measures, this formulation is immediately extensible to higher dimensions. Given vector-valued
targets E R arranged in Y E N their vector quantiles are also in Rd. The vector quantile
can he installed with pip install vqr; source available at https://gt http//gtttb.cc.o com vi vistalab-tedt betcchnio//qq
3

---

Published as a conference paper at ICLR 2023
levels are sampled uniform grid on [0, 1]d with T evenly spaced points in each dimension, resulting in T d-dimensional vector quantile levels, arranged as U E IRT xd The OT objective can be written as E - = II @ S, where S E R TxN,, and thus can be naturally extended to d > 1 by defining the pairwise inner- inner-product matrix = UY E R XN The result is a d-dimensional discrete linear estimate of the CVQF, Q xx; which is co-monotonic (I) with for each æ. Appendix detalis how the estimated CVQF obtained from the dual variables in the high-diment -ige--imensional case.
3 SCALABLE VQR
To motivate our proposed scalable approach to VQR, we first present the exact formulations of VQR and discuss their limitations. The OT-based primal and dual problems are presented below.
T
ui T
y,
min & T + :p_ + (BT X)
TT20
v..p.
(5) s.t. Vi.j
s.t. II 1T
[gb]
(6)
III = μ
[p]
P: + B, I)j + 2/
2 us T
yj
II- X
[g]]
Here v RN, 4 ITT and B E IRT xk are the dual variables, μ # v 1 the marginals and X -
X is the mean-independence constraint. The solution of these linear programs results in the convex potentials (b y) and ((u) x + <p(u), where the former is discretized over data samples and the latter over quantile levels. The estimated CVQF is then the transport map between the measures of the quantile levels and the data samples, and is given by @ x) = V. (B(u) T x + } (see also Appendix A.6). Notably, co-monotonicity of the CVQF arises due to it being the gradient of convex function.
The primal VQR problem (5) has T N variables and T (k + N constraints. Solving the above linear programs with general-purpose solvers has cubic complexity in the number of variables and constraints and is thus intractable even for modestly-sized VQR problems (see fig. A9). Numerous works proposed fast solvers for variants of OT problems. A common approach for scaling OT to large samples is to introduce an entropic regularization term - Ei, II, log II., to the primal objective. As shown by Cuturi (2013). this regularized formulation allows using the Sinkhorn algorithm which provides sigmihcant speedup. Other solvers have been proposed for OT variants with application-speciff regularization terms on II (Ferradans al. 2013 Solomon et al. 2013 Benamou al. 2015). However, the addition of the mean- mnean-independence constraint on II renders these solvers inapplicable for VQR. For example, the Sinkhorn algorithm would require three projections per iteration: one for each marginal, and one for the mean-indeppnddence constraint. Crucially, for the latter constraint the projection onto the feasible set has no closed-fore -fom solution and thus solving VQR via Sinkhorn is very slow in practice.
Limitations of existing methods.
Relaxed dual formulation. Another approach for scaling OT is by solving its dual formulation. Genevay et al. (016)) showed that the relaxed dual version of the standard OT problem is amenable to stochastic optimization and thus scale to large samples. The relaxed dual formulation of VQR is obtained from eq. (6) (see appendix ().4). and can we written as
T
Σ exp (1 (u:Tu y, -3
min T (T X) + Σ pr log 4.8
BTT,, xx (3,5)))))))))))
where 5 controls the exactness of the the objective; as decreases, the relaxed dual more closely approximates the exact dual. Of note are some important properties of this formulation: (i) it encodes the mean- mean-indess -naeeeddeeendeece constraint in the objective; (ii) it is equivalent version of the VQR OT primal (5) (see appendix A.5)): (iii) it is an unconstrained convex objective amenable to stochastic- -rradient based optimization. The key drawbacks of this approach the linear specification of the resulting CVQF in x, and the potential violation of co- co-monotonicity due to relaxation of the constraints. We address the former limitation in Section 4 and the latter in Section S
entropic-reguularized
SGD-based solver. The relaxed dual formulation of VQR 7 involves only Tdd le + N optimization variables and the objective is amenable to GPU-based acceleration, as it involves only dense-matrix
4

---

Published as a conference paper at ICLR 2023
300
Time (mins.)
Time (mins.)
0.8
fyxill
0.28
vxx kxxll:
0.6
200
0.4
100
0.24
0.22
0.00
0.25
0.50 0.75 1.00 1.25 1.50 1.75
2.00
10
20
30 40
50
60
70
80
90
100
Number of samples
le6
Quantile levels (per dimension)
Figure 2: Proposed VQR solver scales to large N and T with improving accuracy. Computed on MVN dataset with d 2 and l 10. Blue curve shows KDE-L1 distance, defined in section Left: Sweeping N; T 50 and Bx 50k. Right: Sweeping T; N 100k, and By 2500.
multiplications and pointwise operations. However, calculating this objective requires materializing two T X N inner-proet proouct matrices, causing the memory footprint to grow exponentially with d. This problem exacerbated when is high-dimensional, as it must also be kept memory. These memory requirements severely limit the problem size which can be feasibly solved on general- purpose GPUs. Thus, naive stochastic gradient descent (SGD) with insufficient for obtaining a scalable solver.
of samples 7 is
To attain constant memory footprint w.r.t. T and N, sample data points from {(x), 32) 3,)
together with vector quantile levels from {u;} TMd and evaluate eq. 7 on these samples. Contrary to standard SGD, when sampling quantile levels, select the corresponding entries from / and μı; likewise, when sampling data points, we select the corresponding entries from sb. Thus, by setting fixed batch size for both data points and levels, solve VQR with a constant memory footprint irrespective of problem size (fig. 2 Moreover, we observe that in practice, using SGD adds only negligible optimization error the solution, with smaller batches producing more error as expected (Appendix fig. [Ala).. Refer to Appendix B for convergence analysis of this approach and to Appendix G for implementation details.
NONLINEAR VQR
In this section we propose an extension to VQR which allows the estimated model to be non-linear in the covariates. Carlier et al. (2016) proved the existence and uniqueness of CVQF which satisfies strong representation, i.e.. ( x) Qy (U; and is co-monotonic w.r.t. U. In order to estimate this function from finite samples, they further assumed a linear specification in both and æ, given by the model Q (u; x>) = (u)) T a + a(u). This results in the OT problem with the mean- independence constraint (5). However, in practice and with real-world datasets, there is no reason to assume that such specification is valid. In cases where the true CVQF is a non- non-linear function of æ this model is iis-ifecdf and the estimated CVQF does not satisfy strong representation.
Extension to nonlinear specification. We address the aforementioned limitation by modelling the CVQF
non-linear function of the covariates parametrized by O, i.e., x (u) T (x) + auu). We fit 0 jointly with the regression coefficients B, The rationale is to parametrize an embedding go (a)) such that the model is Q better specified than O+ the sense of strong representation. Encoding this nonlinear CVQF model into the exact formulations of VQR (eqs. 5) and (6)) will no longer result in linear program. However, the proposed relaxed dual formulation of VQR 7 can naturally incorporate the aforementioned non-linear transformation of the covariates, as follows
T
Σ exp (w: yj B, T. go())) w3y) ) (8)
min &6 v uir (B Go(X)) +e Σ us log $33,0
where Gx I go (X) the empirical mean after applying go to each sample.
We optimize the above objective with the modified SGD approach described in section 3 The
estimated non-linear CVQF can then be recovered by QNL (u:z) V. {B(u) ge (x
p(u)).
crucial advantage of our nonlinear CVQF model that ge may embed x E R into a different dimension, K o k. This means that VQR can be performed, e.g., lower-dimensional projection
5

---

Published as a conference paper at ICLR 2023
Ground truth
VQR
CVAE
NL-VQR
Figure 3: NL-VQR quantitatively and qualitatively outperforms other methods in conditional distribution estimation. Comparison of kernel- density estimates of samples drawn from VQR, CVAE, and NL-VQR, on the conditional- els were trained with N 20k sam- ples; for VQR bers depict the KDE-LI metric. Lower is better.
V
V
V
1-277 X=2.1
dataset. Mod-
102 0.17 0.30
was used. Num-
V
V
0.32
of x or on a "lifted" higher-dimensional representation. For low-dimensional X which has highly- nonlinear relationship with the target Y, lifting it into a higher dimension could result in the linear regression model (with coefficients B and a) be well-specified. Conversely, for high-dimensional X which has intrinsic low-dimensionality (e.g,, an image), ge may encode projection operator which has inductive bias, and thus exploits the intrinsic structure of the feature space. This allows one to choose Go suitable for the problem at hand, for example, convolutional neural network in case of an image or a graph neural network if x: lives on graph. all our experiments, we indeed find that learned non-linear feature transformation substantially improves the representation power of the estimated CVQF (see fig. 3 Appendix figs. A13 and A14).
Simultancous Quantile Regression (SQR). SQR is the task of simultaneously resolving multiple conditional quantiles of a scalar variable, i.e. a possibly non-linear version of eq. (3) Nonlinear SQR has been well-studied, and multiple approaches have been proposed (Brando et al. 2022; Tagasovska
& Lopez-Pazz 2019). We highlight that SQR can be performed as special case of NL- VQR (8) when d = 1. See Appendix fig. A12 for an example. Other SQR approaches enforce monotonicity of the CQF via different ad-hoe techniques, while with NL-VQR the monotonicity naturally emerges from the OT formulation explained in section 3 thus arguably a preferable approach. To our knowledge, OT-based SQR has not been demonstrated before.
Prior attempts nonlinear high-dimensional quantiles. Feldman et al. (2021) proposed approach for constructing high-dimensional confidence regions, based on a conditional variational autoencoder (CVAE). crucial limitation of this work is that it does not infer the CVQF; their
resulting confidence regions are therefore not guaranteed to be quantile contours (see fig. Tab. because
they do not satisfy the defining properties of vector quantiles (eqs. (1) and (2)
VECTOR MONOTONE REARRANGEMENT
Solving the relaxed 7 may lead to violation of co-monotonicity in QY since the exact constraints in eq. (6) are not enforced. This is analogous to the quantile crossing problem in the sealar case, which can also manifest when QR performed separately for each quantile level (Chernozhukov et al. 2010) (Appendix fig. A12d). In what follows, we propose way to overcome this limitation.
Consider the case of scalar quantiles, i.e., d 1. Denote @ (u; æ)) as an estimated CQF, which may be non-monotonic u due to estimation error and thus may not valid CQF. One may convert Qyix (u; into a monotonic Qxix (u; æ)) through rearrangement as follows. Consider a random variable defined Vix Q (U; where U U(0, 1). Its CDF and inverse CDF given by (): = I QQi S du. and Qxy (u;z:) = inf {v: PXix(U: x)) 2 m Q (u; x) is the true CQF of YIX and is thus necessarily monotonie. It can be shown that Q (u; x) is no worse estimator for the true CQF Qxx (w; than Q (u; æ) in the nn-norm sense (Chernozhukov et al. 2010). In practice, this rearrangement is performed in the scalar case by sorting the discrete estimated CQF. Rearrangement has no effect if Q (u; x)) is already monotonic.
Here, we extend the notion of rearrangement to the case of d > 1. As before, define Ý == Q az) where U 1[0, 1]4, and Q x (1; ar) is the estimated CVQF. If it is not co-monotonic,

---

Published as a conference paper at ICLR 2023
10"1
6 10-2
MV
10-2
No VMR
No VMR
3 10-2
VMR
VMR
10-1
10
10-3
10
10
10-2
10-3
10
VMR improves strong representation (QFD; left), and completely eliminates tonicity violations (MV; right). VMR allows for smaller values, thus improving accuracy (lower QFD) without compromising monotonicity (zero MV). Without VMR, MV increases (right, red) as
Figure 4:
the problem becomes more exact (smaller E). MVN dataset; N 20k, T = 50, d = 2, k 1.
as defined by eq. I we can compute a co-monotonic Q x (u; æ) by calculating the vector quantiles
of Y = z, separately for each We emphasize that this amounts to solving the simpler vector quantile estimation problem for specific æ, as opposed to the vector quantile regression problem.
Let Q = [bl, E T xd be the estimated CVQF xxpx (u; æ) sampled at 7 levels U = [[,,]; E
RT xd by solving VQR. To obtain Q& x (w; x) sampled at U we solve
max Tay 20
Σ 1
Ti. U; 9i
s.t. II
7
and then compute Q = 7 UQQ We define this procedure as vector monotone rearrangement
(VMR). The existence and uniqueness of co-monotonic function, mapping between the measures of
Mrxx(U U: x) is due to the Brenier's theorem as explained in section Brenier 1991 McCann 1995; Villani 2021). VMR can be interpreted as "vector sorting" of the discrete CVQF estimated
U and x
with VQR, since it effectively permutes the entries of Q such that the resulting Q is co-monotonic
with U. Note that eq. (9) is an OT problem with the inner-produc matrix as the ground cost, and
simple constraints on the marginals of the transport plan II. Thus, VMR ( significantly simpler to solve exactly than VQR (5). We leverage fast, off off-the-shelf OT solvers from the POT library (Flamary
al. 2021 Bonneel et al. 2011) and apply VMR post- post-processing step, performed after the estimated CVQF Q (u; x) is evaluated for a specific x.
In the 1d case, quantile crossings before and after VMR correction can be readily visualized (Appendix figs, Al2d and Al2eh. For > 1, monotonicity violations manifest as (u,-u, ( (u,,)) < 0 (Appendix lig. [AIT. Figure demonstrates that co-monotonicity violations are completely eliminated and that strong representation strictly improves by applying VMR.
6 EXPERIMENTS
We use four synthetic and four real datasets which are detailed Appendices D.I and D.2 Except for the MVN dataset, which was used for the scale and optimization experiments, the remaining three synthetic datasets were carefully selected to he challenging since they exhibit complex nonlinear relationships between X and Y (see e.g. fig. Tbp. We evaluate using the following metrics (detailed in Appendix E}: (i) KDE-LI, an estimate of distance between distributions: (ii) QFD, a distance measured between an estimated CVQF and its ground truth; (iii) Inverse CVQF entropy: (iv) Monotonicity violations; (v) Marginal coverage; (vi) Size of a-confidence set. The first three metrics serve as measure of strong representation (2). and can only used for the synthetic experiments since they require access to the data generating process. The last two metrics are used for real data experiments, as they only require samples from the joint distribution. Implementation details for all experiments can be found in Appendix G
6.1 SCALE AND OPTIMIZATION
To demonstrate the scalability of our VQR solver and its applicability large problems, we solved VQR under multiple increasing values N and T, while keeping all other data and training settings fixed, and measured both wall time and KDE-LI. We used up to N 2 10y data points, d dimensions, and T = 100 levels per dimension (I = 1002 levels in total), while sampling both data points and quantile levels stochastically, thus keeping the memory requirement fixed throughout. This
7

---

Published as a conference paper at ICLR 2023
enabled to run the experiment on a commodity 2080Ti GPU with 11GB of memory. Optimization experiments, showing the effects of e and batch sizes on convergence, are presented in Appendix FF]
Figure 2 presents these results. Runtime increases linearly with N and quadratically with T, can be expected for d = 2. KDE-LI consistently improves when increasing N and T. showcasing improved accuracy in distribution estimation, especially more quantile levels estimated. To the best of our knowledge, this is the first time that large-scale VQR has been demonstrated.
6.2 SYNTHETIC DATA EXPERIMENTS
Here our goal is to evaluate the estimation error (W.r.t the ground-truth CVQF) and sampling quality (when sampling from the CVQF) of nonlinear VQR. We use the conditional-banane rotating stars and synthetic glasses datasets, where the assumption of a linear CVQF is violated.
Baselines. We use both linear VQR and conditional variational autoencoders (CVAE) (Feldman et al. 2021> as strong baselines for estimating the conditional distribution of Y|X. We emphasize that CVAE only allows sampling from the estimated conditional distribution; it does not estimate quantiles, while VQR allows both. Thus, we could compare VQR with CVAE only on the KDE-L1 metric, which is computed on samples. To the best of our knowledge, besides VQR there is other generative model capable of estimating CVQFs.
The results of the conditional- experiment, comparing VQR, NL-VQR and CVAE, are presented in fig. 3 Appendix table [A2 and fig. [A13] Two-dimensional KDEs of conditional distributions for three values of x are depicted per method (fig. B. Qualitatively, sampling from either NL-VQR or CVAE produces accurate estimations of the ground truth distribution, when compared linear VQR. This indicates that the linear VQR model is mis-specified for this data, which is further corroborated by the entropy metric (table A2 The entropy of the inverse CVQF is lower for linear VQR, indicating non-uniformity and therefore mis-specification. Quantitatively in terms of the KDE-LI metric, NL-VQR outperforms CVAE by 23%, and in terms of QFD, nonlinear VQR results in 70% lower error than linear VQR.
Conditional Banana.
Rotating Stars. The results of the rotating stars experiment with the same baselines as above, are presented in Appendix fig. A14 and table A3 As illustrated in fig. Tb this dataset is highly challenging due to the non-trivial dependence between X and the tails of Y Y.. We observe that NL-VQR qualitatively and quantitatively outperforms both CVAE and linear VQR by a large margin. It is therefore the only evaluated method that was able to faithfully model the data distribution.
These results highlight that NL-VQR is significantly better at estimating CVQFs compared to VQR. Another point of note is that the CVAE model required 35 minutes train, while VQR and NL-VQR were trained within less than 4 minutes. NL-VQR therefore outperforms CVAE with almost an order of magnitude speedup in training time, and outperforms VQR while requiring similar runtime. Additional synthetic experiments, showcasing OT-based SQR and VMR are presented in Appendix F
6.3 REAL DATA EXPERIMENTS
VQR has numerous potential applications, as detailed in section I Here we showcase one immediate and highly useful application, namely distribution-fiee uncertainty estimation for vector- vector-valued targets. Given P(x.y) and a confidence level a E (0, 1). the goal is to construct a conditional a- confidence set C. (x) such that it has marginal coverage of - at, defined as IP [Y E C,, = 1 a. A key requirement from an uncertainty estimation method is to produce a small C. Cxxx)) which satisfies marginal coverage, without any distributional assumptions (Romano et al.) 2019).
Baselines. We compare nonlinear VQR against (i) Separable linear QR (Sep-QR): (ii) Separable nonlinear QR (Sep-NLQR); (iii) linear VQR. For Sep-QR and Sep-NLQR the estimated CVQFis is
Soo (u; x)
QQQR (uz;
QQR (uai x))
(10)
where @ ()::)) is obtained via 1d linear or nonlinear quantile regression of the variable Y,. These separable baselines represent the basic approaches for distribution-fres uncertainty estimation with vector- valued variables. They work well in practice and the size of the a-confidence sets they produce serves an upper bound, due to their inherent independence assumption.
Evaluation procedure. The key idea for comparing distribution free uncertainty estimation ap- proaches is as follows. First, nominal coverage rate a* is chosen. An a-confidence set is then
8

---

Published as a conference paper at ICLR 2023
constructed for each estimation method, such that it obtains the nominal coverage (in expectation). This is similar to calibration in conformal prediction (Sesia & Romano 2021) since a is calibrated to control the size of the a-confidence set. Finally, the size of the a-confidence set is measured as the evaluation criterion. Appendix G contains experimental details and calibrated a values.
Results. Across the four datasets, NL-VQR acheives 34-60% smaller a-confidence set size compared to the second-best performing method (section and fig. [A15
New the confidtecce region (sover buittae)
Requealed Coverage
table A4). The reason for the superior performance of NL- V QR is its ability to accurately capture the shape of
70
60
the conditional distribution, leading to small confidence
sets with the same coverage (section 6.3 and fig. A16).
N.VVR
Sep-
Sep-QR
SOFTWARE
2.5
NL-VOR e1 .000
We provide our VQR solvers as part of a robust, well-tested python package, vqz (available in the supplementary ma- terials). Our package implements fast solvers with support for GPU-acceleration and has a familiar sklea rn-style API. supports any number of dimensions for both in- put features (k) and target variables (d), and allows for arbitrary neural networks to be used the learned non- linear feature transformations, 9o. The package provides tools for estimating vector quantiles, vector ranks, vector quantile contours, and performing VMR refinement step after fitting VQR. To the best of our knowledge, this would be the first publicly available tool for estimating conditional vector quantile functions at scale. See Ap- pendix H for further details.
NL Sep QR (2.19)
VOR (4.43)
2.
yjx
15
0.0
-2.0
-15
-1.0
-0.5
0.0
0.5
1.0
1.5
Figure NL-VQR produces sub- stantially smaller confidence sets than other methods at a given confidence level. Top: Mean C. size; Bottom: Ex- ample C. (zx for specific x. Calculated on house prices test-set. Full details in figs. A15 and A16
CONCLUSION
In this work, we proposed NL-VQR, and a scalable ap- proach for performing VQR in general. NL-VQR over- comes key limitation of VQR, namely the assumption that the CVQF linear in X. Our approach allows mod- elling conditional quantiles by embedding X into space where VQR is better specified and further to exploit the models like CNNs). We proposed new relaxed
structure of the domain of X (via domain-specinn dual formulation and custom solver for the VQR problem and demonstrated that we can perform VQR with millions of samples. As far as we know, large scale estimation of CVQFs even VQFs has not been previously shown. Moreover, we resolved the issue of nigh-dimensiooal quantile crossings by proposing VMR, a refinement step for estimated CVQFs. We demonstrated, through exhaustive synthetic and real data experiments, that NL-VQR with VMR is by far the most accurate way to model CVQFs. Finally, based on the real data results, argue that NL-VQR should be the primary approach for distribution-fre uncertainty estimation, instead of separable approaches which assume independence.
Limitations. As with any multivariate density estimation task which does not make distributional assumptions, our approach suffers from the curse of dimensionality, especially in the target dimension. Overcoming this limitation requires future work, and might entail, e.g., exploiting the structure of the domain of Y. This could be achieved by leveraging recent advances in high-dimensional neural OT (Korotin et al. 2021). Another potential limitation is that the nonlinear transformation (a)) shared across quantile levels (it not a function of u), though evaluating whether this is truly limiting assumption in practice requires further investigation.
In conclusion, although quantile regression a very popular tool, vector quantile regression is arguably far less known, accessible, and usable in practice due to lack of adequate tools. This is despite that fact that is a natural extension of QR, which can be used for general statistical inference tasks. We present the community with an off-the-shelf method for performing VQR in the real world. We believe this will contribute to many existing applications, and inspire a wide-range of new applications, for which it is currently prohibitive or impossible to estimate CVQFs.

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGEMENTS
Y.R. was supported by the Israel Science Foundation (grant No. 729/21). Y.R. thanks Shai Feldman and Stephen Bates for discussions about vector quantile regression, and the Career Advancement Fellowship, Technion, for providing research support. A.A.R., S.V., and A.M.B. were partially supported by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 863839), by the Council For Higher Education Planning Budgeting Committee, and by the Israeli Smart Transportation Research Center (ISTRC).
REFERENCES
C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa. The quickhull algorithm for convex hulls. ACM Trans. Math. Softw., 22(4):469483, dec 1996. ISSN 0098-3500. doi: 10.1 1051135/235815 235821. hLL.:tt6e/.oo.cccr0000.1442/255. 235821
org.
Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré. Itera- tive bregman projections for regularized transportation problems. SIAM Journal on Scientific Computing, 37())A1111AA11. 2015.
Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich. Displacement inter- polation using Lagrangian mass transport. ACM Transactions on Graphics, 30(6): n 2011. doi: 10. 145/2070781 2.241992 URL httpt:/pp:/ha inria..i............ fr ft/afff-0ha622222
Axel Brando, Joan Gimeno, Jose Rodríguez-Serrano, and Jordi Vitrià. Deep non-crossing quantiles through the partial derivative. arXiv preprint arXiv:2201.12845 2022.
Yann Brenier. Polar factorization and monotone rearrangement of vector- valued functions. Commu- nications pure and applied mathematics, 44(4):375-417, 1991.
Guillaume Carlier, Victor Chernozhukov, and Alfred Galichon. Vector quantile regression: An optimal transport approach. Annaly of Statistics, 4463):1 2016. ISSN 00905364. doi: 10.1214/15-ASS1401.
Guillaume Carlier, Victor Chernozhukov, and Alfred Galichon. Vector quantile regression beyond the specified case. Journal of Multivariate Analysis, 161:96-102, 2017. ISSN 10957243, doi: 10.1016/j jmva.
Guillaume Carlier, Victor Cheroozukov, Gwendoline De Bie, and Alfred Galichon. Vector quantile regression and optimal transport, from theory to numerics. Empirical Economics, 2020. ISSN 14358921. doi: 10. (0.10077.0008810200009999
Benjamin Charlier, Jean Feydy, Joan Alexis Glaunès, François-Davia Collin, and Ghislain Durif. Kernel operations the gpu, with autodiff, without memory overflows. Journal of Machine Learning Research, 22(74):1-6, 2021. URL
html
Victor Chernozhukov, Iván Fernández- Val, and Alfred Galichon. Quantile and probability curves without crossing. Economelrica, 1893) 2010.
Victor Chernozhukov, Alfred Galichon, Marc Hallin, and Mare Henry. anngeeanto depth, quantiles, ranks and signs. The Annals of Statistics, 45(1):223 256, 2017. doi: 10.1214/ 16-A0$1450. http:htt:///da wrg/10...000.. 1214/16 2
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Steven Diamond and Stephen Boyd. CVXPY: Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1-5, 2016.
Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated multiple-output quantile regression with representation learning. arXiv preprint arXiv:21 2021.
10

---

Published as a conference paper at ICLR 2023
Sira Ferradans, Nicolas Papadakis, Julien Rabin, Gabriel Peyré, and Jean- SannFFaançoss Aujol. Regularized discrete optimal transport. In International Conference on Scale Space and Variational Methods in Computer Vision. pp. 428-439. Springer, 2013.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Lactitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, levgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer. Pot: Python optimal transport. Journal of Machine Learning Research, 22(78):1-8, 2021. URL httthttt:://ttt http://mmloorrr/reeeeeeeeee/e www.orgpppppp///pppoo 20 451, html
Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach. Stochastic optimization for large- scale optimal transport. Advances in neural information processing systems, 29, 2016.
Saeed Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming, SIAM Journal on Optimization, 23(4):2341-236 2013.
Elad Hazan. Lecture notes: Optimization for machine learning. arXiv preprint arXiv:1909.03550, 2019.
Roger Koenker and Gilbert Bassett, Regression Quantiles. Econometrica, (633) 1978. ISSN 00129682. doi: 10.2307/1913643.
Alexander Korotin, Lingxiao Li, Aude Genevay, Justin M Solomon, Alexander Filippov, and Evgeny Burnaev. Do neural optimal transport solvers work? continuous wasserstein-2 benehmark. Advances in Neural Information Processing Systems, 14:11499-14600 2021.
Robert J McCann.
Existence and uniqueness of monotone measure-preserring maps. Duke Mathe- matical Journal, 80(2)309-3323, 1995.
Gabriel Peyré and Marco Cuturi. Computational optimal transport. Foundations and Trends in Machine Learning, 116-66)355-607, 2019.
Tyrrell Rockafellar. Conjugate duality and optimization. SIAM, 1974.
Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized Quantile Regression. In H Wallach, H Larochelle, A Beygelzimer, F tetquuuoteeee Alché-Buc, E Fox, and R Garnett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Matteo Sesia and Yaniv Romano. Conformal prediction using conditional histograms. Advances in Neural Information Processing Systems, 34:6304 6614,, 2021.
Hongjian Shi, Mathias Drton, and Fang Han. Distribution-free consistent independence tests via center-outward ranks and signs. Journal of the American Statistical Association, 117(537):
395-410, 2022. doi: 1000000001121449.22220 1782223. 01621459.2020.1799
Justin Solomon, Fernando De Goes, Gabriel Peyré, Mareo Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. Convolutional wasserstein distances: Efficient optimal transportation on geometric domains. ACM Transactions on Graphics (ToG), 34(4):1-11, 2015.
Natasa Tagasovska and David Lopez-Paz. Single-model uncertainties for deep learning. Advances in Neural Information Processing Systems, 32, 2019.
Cédric Villani. Topics in optimal transportation, volume 58. American Mathematical Soc., 2021.
Pauli Virtanen, Ralf Gommers, Travis Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, J Carey, ilhan Polat, Yu Feng. Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientifie Computing in Python. Nature Methods, 17:261-272 2020. doi: 10.1038/s4 192-019-0666-22
11

---

Published as a conference paper at ICLR 2023
Justin R. Williams and Catherine M. Crespi. Causal inference for multiple continuous exposures via the multivariate generalized propensity score, 2020. URL httpttsss/aaa..v..rr///////// 13767
Yue Zhao, Zain Nasrullah, and Zheng Li. Pyod: A python toolbox for scalable outlier detection. Journal of Machine Learning Research, 20(96))1-7, 2019, http:://mmm..rrg/ v20/19-00 011. html
12

---

Published as a conference paper at ICLR 2023
A DERIVATIONS
In this section, we present the following derivations:
Formulating one-dimensional QR as correlation- -axiaiiiiimmaiiiiiatiiinnn problem (Appendix A.1
Rephrasing correlation maximization one-dimensional optimal transport (Appendix A.2)
Extension of the OT-based formulation for QR to the multi-dimennional targets (Ap- pendix A.3)
Relaxing the dual formulation of the OT-based VQR problem (Appendix A.44)
The equivalence between the entropic-regullrizze version of the OT primal and the relaxed dual formulation (Appendix A..3.
Calculating the conditional (vector) quantile functions from the dual variables (Ap- pendix A.6
The derivations in this section not rigorous mathematical proofs; instead they are meant show an casy-to-follow way to obtain the high-dimensional VQR relaxed dual objective (which we eventually solve), by starting from the well-known one-dimensional QR based the pinball loss.
A.1 QUANTILE REGRESSION AS CORRELATION MAXIMIZATION
The pinball loss defined above be written as,
u2, >> 0
Pu(z)
(u 1)2,
8 8
u12 - + 2, 2>0
(u 1)2, 250
zt (u 1)z
where 2+ A m0) Note that also define - ف max -2) which later. Given the
continuous joint distribution Px,y, and assuming linear model relates the u-th quantile Y to X, then performing quantile regression involves minimizing
min Ex (V - by Tx aw) b,.
+ - (1- (V bw TX aw)
(11)
Define Z(au; bu) S b X aus the above problem be written
min Ex,Y [anu b.)+ (1 - u) ((a,, b.)] bu
(12)
Define Pu A Z(u, b and N₂ A ((a b.) as the positive and negative deviations from the true u-th quantile, then notice that (i) P. 2 O and N. 2 0; (ii) P. N, ( bw). Introducing Pa and N,, as slack variables, we can rewrite the above optimization problem as
min Exy [Pu - (1 (a,, b.)l ag bu, Nu
S.t.
Pv 2 0, N, 20
P. N., Z(an, bu), with multiplier [V.].
Solving the above problem is equivalent to solving the Lagrangian formulation
min max Ex,y [Pu. u) ((a. bw.) V. (Pu. -o. Z( bw))l
b..
s.t., P. 2 O,No 20.
13

---

Published as a conference paper at ICLR 2023
Substituting Zu (au bu) Y - b. X a, get
min Pu,Nu,,Ru..b V.
max
Ex, Pu - 261-111 u) (Y - b. _ X- 15-aw) an) - -Va (Pu - No -Y+ x + au)
S.L,,
N,
In the first term, Y can be omitted because it is independent of the optimization variables, and in the second term, Y can be separated, this yields
max Ex [V.Y] Vu
min Ru. bu
Ex,y Pa (u-1) () T an) Va (Pa N₂ + ba X + am)]
s.l., P. 2 0, No 20
Rewriting the terms by separating au, bu, Pu, Nu yields
max
Ex [Va]+
min b
Ex,y [a (1 V.)) N.V. aulV -01 (1- (1-4) u)) by T (V,X (1- -1-n)xx]
Vu
s.l., Pa 2 0, No 2.
By treating Pu, Nu, Ou, by as Lagrange multipliers we can obtain the dual formulation. Since Pa, Na, 2 they must participate incquality constraints, and since aus by are unconstrained, they participate in equality constraints. Thus,
Σx, [V.Yu] V.
s.t. Vu 20 0
[Nu]
V. S 1
[Pu]
Σx [Va] = (1-u)
[au]
Ex.Y [V,
u)E X
[bu]-
Note that (i) complementary slackness dictates that we always have Pu, 2 0, and that (ii) the last two constraints can be seen together implying the mean independence condition of V E [X].
Solving for multiple quantiles simultaneously can be achieved by minimizing the sum of above optimization problem Yu e 10, 1). In addition we demand monotonicity of the estimated quantiles, such that u' - b., X a, 2 b.. X + a.' According to the KKT conditions, the solution to the above problem must satisfy complementary slackness, which in this case means
(1 - Vu =0
V.NN 0.
Following the definitions of Pu, No. we have the following relation between the estimate quantile b,. X au and Vu:
Y > by T TX + ag
=
Pa >
V =
Y < b T X + 0g
= Na >
V, = 0
Y by X au
P. No =
0 <
Since Px,Y is continuous distribution, then for any aus bu. we have P Y b. TX am] 0. Thus in practice we can ignore this case and write
I(Y 2 by xx
(13)
Therefore the monotonicity constraint translates to Vn as
n' b. T au 2 b,d X Oud
V, < Var
14

---

Published as a conference paper at ICLR 2023
Adding the monotonicity constraint Vu to the gives rise dual problem:
max V.
Ex.Y
[V.Yu] du
s.t.
Vu 20 0
[N
Vu S
[P.]]
Ex [Va] = (1 - u)
[au]
E [VaX] (1 u) [X]
[bu]
2 u' Va S Vie
Let us now consider dataset ((a,, y:) of ii.d. samples from Px,y where E R and E R. We are interested in calculating quantiles for levels u1 = 0 < U2 <... < ur I 1. Furthermore, denote the sample mean ã E RR where T2 = E L.a/N. By discretizing the above problem, we arrive
max Vr.
Σ Σ
Vr
s.t. V 2 0
[N
V 1
[Pr.]]
N
Σ
Vr,1 = (1 - ur)
[a,]
N im1
N Vr, 2x)
Σ
[b,
VT, S VT-1i: VI,
Now we can vectorize the above formulation by defining the matrix V E RRNN containing elements Vr, the first-order finite differences matrix
1
0
1
1
D
0
-1
E IR. TRT
0
-1 1
and the vector of desired quantile levels, u [ui.... Wr] T. The monotomicity constraint therefore becomes v, D > 0 Vi = I,...,N . where v, is the i-th column of V. We also denote the covariates E R Nxk and response vector y E RVV. Thus, the problem vectorizes as,
matrix X
max 1T T Vv V
s.t. V7,, 20 0
[N,,
V7 1
[Pul
V1y = (17 u)
[a]
N
1
VX (IT u.
[B]
N
V T
D20.
D
where a E R. and E Txk are the dual variables which contain the regression coefficients per quantile level. We can observe the following for the above problem:
15

---

Published as a conference paper at ICLR 2023
1. For any random variable, its zeroth quantile is smaller than equal to any values the variable takes. In particular, for r = 0 we have Vo, = 1 Vi because from eq. (13) we know that Vo, is the indicator that ys is greater that the zeroth quantile of Y.
2. The last constraint V D'D 2 O, which enforces non-increasing monotonicity along the quantile level, i.c., V I V VT 2 T', also enforces that Vr, 2 0.
Therefore, we can observe that the first (V+,; 2 0), second (V+,i S 1) and last (V T D 2 0) constraints are partially-rednnde and can be condensed into only two vectorized constraints: (1) V D > 0 which ensures monotonicity and non-negativity of all elements in V; (2) V DIr S Ix which enforces that Vo,; < vi. Notice that condensing the constraints in this manner comes with the inability to interpret the meaning the Lagrange multipliers Pr..: N₁ However, the advantage that we have less constraints in total and the interpretability of the multipliers a, B is maintained. Thus, we arrive at the following vectorized problem,
max 1T Vy
s.t.
V1N (17 u)
[a]
VX = =(1r u1)
[
N
V 2 0
V D1T SIN-
A.2 CORRELATION MAXIMIZATION AS OPTIMAL TRANSPORT
Following Carlier et 2016)) the above problem can be ce-oomuulted as an Optimal Transport problem, i.e. the problem of finding a mapping between two probability distributions.
Assume we are now interested in estimating the quantiles of Y at T uniformly-sampii levels. T In the above problem, one can decompose the objective as T
[u,.... 4r] T T
2
1r VV (1r D DD (D 1T) T V)y.
we then denote
II D VVERRTX
N
1 T
1 2
T T TERT. T
D
T'
then we can write the objective NT Hy. In addition, denote
1 D (IT u) = ERT V 1 N
17
INERN
which represent (respectively) the empirical probability measure of the quantile levels u and the data points (X, y) (we choose both measures be uniform), Now, by using the decomposed objective, and by multiplying the first two constraints by on either side, we obtain the following equivalent problem:
T fl
max u
s.t. II N = 11 lT
[D-Ta]
(14)
IIX uv X= 1
D-TBB
14 カ=v= N 1N
Note that (i) we ignore the normalization constants where they do not affect the solution; (ii) the Lagrange multipliers are scaled by D since the constraints scaled by D T
16

---

Published as a conference paper at ICLR 2023
The interpretation of this formulation is that seek transport plan II between the measures of the quantile levels and the target Y[X, subject to constraints on the plan which ensure that its marginals the empirical measures μ and и, and that mean independence E [ = E [X] holds. Each individual entry i in this discrete plan is the probability mass attached to (ui, U5) in this optimal joint distribution.
A.3 EXTENDING THE OPTIMAL TRANSPORT FORMULATION TO VECTOR QUANTILES
We now wish to deal with the case where the target variable can be vector. Observe that for the scalar case, we can write the OT objective
Hl ΣΣ
H.y. 11 3 S,
i=1==
where S RRTN is matrix of pairwise products, i.e. S u.Ylj. and @ denotes the Hadamard (elementwise) product.
Now let assume that y; TRd for any d 2 1, and thus target data will now be arranged Wxd The quantile levels must now also d-dimensional, since we have a quantile level dimension for each data dimension. We will choose a uniform grid on [0, 1] on which to compute our vector quantiles. Along each dimension of this grid we sample equally-spacee points, giving us in total 7 points in dimensions.
To keep the formulation of the optimization problem two dimensional, we arrange the coordinates of these points as the matrix U E RT xd Thus, we can naturally extend the pairwise product matrix S to the multi-dimensional case using d-dimensional inner anner-product between each point the quantile
level grid U and each target point in Y. This yields the simple form S UYT,, where S E TT which can be plugged in to the above formulation (14) and solved directly. Thus we obtain eq. (5)
A.4 RELAXING THE EXACT OPTIMAL TRANSPORT DUAL FORMULATION
Recall the exact dual formulation of the OT primal VQR problem (5)
min TT + 40 1 + ir (B X)
$.
s.t. Vi,j:
ii B xz + 2/ 2 uu yj
The above problem has unique solution CCrllerr al. 2016). Thus, first-order optimality conditions for each yps yield
P: = max y, - 3, x, b/y
Substituting the optimal y into the dual formulation results an unconstrained but exact min-max problem:
min 2b v + tr (8 X) $
Σ
Me max {u: - 3i B;xT,
13 -
(15)
==1
We relax this problem by using smooth approximation for the max operator, given by
J max(#:) elog
Σ exp
Plugging the smooth approximation into eq.
15 yields the relaxed dual in eq.
A.5 EQUIVALENCE BETWEEN REGULARIZED PRIMAL AND RELAXED DUAL
Adding an entropic regularization term the ba-dd primal formulation of the VQR problem 5 and converting into a minimization problem yields,
min II
III, S) + + (II, log II
s.t.
II 1ps
+
(16)
II P
| [
IIX = X
[-B]
17

---

Published as a conference paper at ICLR 2023
where RN, / RTT X E R are defined in section Note that the entropic regularization term can be interpreted as minimization of the KL-divergence between the transport plan II and the product of marginals, i.e., DRL.
III pM
In order to show the equivalence between this regularized problem and the relaxed dual 77. the key idea is to use the
duality theorem (Rockafellar 1974). This allows us write the dual of an optimization problem in terms of the convex conjugate functions of its objective. To apply this approach, we reformulate eq. (16) into an unconstrained problem of the following form:
(AA"W) + 9" (W) max WEW
min
-f(-V) g(AV),
VEV
where we define pair of operators A + W and A W + V adjoint to each other; V + IR, I V + I and W + I, g" W + I are pairs of convex conjugate
functions. In our problem W W and V V are the vector spaces W W 1 X and
== V* = RTT x RW x
We define the operator
A'II II 11d, IN IIX)
and an indicator function,
[0, a,
ia (2)
00, z + a.
We can now represent eq. (66) as an unconstrained problem,
min II
II, - S) + (TI, log IT) + X) (A
(17)
To derive A, the adjoint operator of A", we can write
((bb, ((b),,,)).AAIII,
(wb, II lrw) + IIN) + (B, IIx)
tr T II 114) + to (4 III,) + tr (3 IIX)
= (llTab T .I) + (GlN ,,II + (3XT, ,I)
(A(b. (.8), I. w.
Thus, A(v,,4,,9) L T + +XXT We then define the functions
(AIIII) = x). (A'II) g" (III = (II, S) + (II, log II),
and their corresponding convex conjugates are therefore given by.
f1v.4.. 3) ((1), (14.,,,300 3), (v, μr X)), (zb, v) + (40 (B,X)
Wij Sij -1
g(W) = exp
E
Using fr and g", we can write 17 as minw (A*WW (W)) then by the Fenchel- + Rockafeller duality theorem, we get the equivalent dual form maxv (-s(-v)- ((VV) Substi- tuting V = 1--,--p. -)) (the dual variables of eq, (16)), converting to minimization problem, and omitting constant factors in the objective, we get
T
min T + 1 (3 P.
ΣΣ exp
(S E
1, 0: 13, 23)) (18)
18

---

Published as a conference paper at ICLR 2023
Now write yp in terms of $, 3 by using first-order optimality condition of the above problem.
Taking the derivative w.r.t. (p and setting to zero, we have:
0 = - exp
Pi
N
exp
1
(Sig - 1/ BB xx,),))
)=1
10
H4
exp
Ej cxp (: (Sig - 1/ - -i xa)))
(19)
5
ii = log
Σ exp
(Sij w/i - Bi az))
(20)
14
c
Finally, substituting eqs.
19 and
20
into eq.
18
and omitting constant terms yields,
min vb To + tr (3 x) + Ms log Σ exp
33, (Sa
43,)
(21)
4.
where S.y u y, Thus, we obtain that eq. 21 is equal to eq.
In summary, we have shown that the relaxed dual formulation of the VQR problem that we solve 7 is equivalent to an entropic-regularized version of the OT-based primal formulation eq. (5).
A.6 EXTRACTING THE VECTOR QUANTILE REGRESSION COEFFICIENTS
The dual variables obtained from the OT formulation (eq. (4) are sp E RT and 3 E RT In the case of scalar quantiles (d 1) we can obtain the conditional quantile function from the dual variables by applying the T x T finite differences matrix D defined in appendix A.1
@ x) [D (Bx +)),
This is effectively taking the u-th component the first-order discrete derivative where is one of the discrete quantile levels.
In the vector case, have T discrete d-dimensional quantile levels. Equivalently, the relation between the dual variables and the quantile function is then
Y (u': x)) V (ax "plll'
Here Bz + is in Tdd and its gradient with respect evaluate its gradient at o ne of the discrete levels u' to obtain the d-dimensional quantile. As explained in section the expression Bæ + 60 is convex in u, and thus the co-monotonicity of the estimated CVQF is obtained by virtue of being the gradient convex function.
Vu (3x + +p) is in RT We then
19

---

Published as a conference paper at ICLR 2023
B CONVERGENCE OF VQR
Below we mention a few comments regarding the optimization and approximation error of VQR.
Linear VQR. The relaxed dual formulation of VQR, presented in eq. (6). is unconstrained smooth convex minimization problem. Thus, this case, gradient descent is guaranteed to converge to the optimal solution up to any desired precision. Moreover, given fixed quantile levels, objective can be written as the minimization of an expectation under the data distribution P(x,v). Under these criteria, performing SGD by sampling i.id. from the data distribution, is known converge the global minimum (Hazan 2019). Chapter 3). We refer to Section 5.4 in (Peyré & Cuturi] 2019) (and references therein) for further analysis and details regarding the convergence rates of SGD and other stochastic optimization methods applied to this problem. Specifically, stochastic averaged gradient (SAG) was shown to have improved convergence rates compared SGD Genevay et (2016).
However, in practice we find that SGD converges well for this problem (see Figure fig. AT0a. and we use it here for simplicity.
Appendix
Nonlinear VQR. In the case of non-lineat VQR, as formulated in eq. S). the optimization is over non-convex objective, and thus no convergence to global minimum is guaranteed. Convergence analyses for non convex objectives only provide guarantees for weak forms of convergence. For e.g.. the analysis methodology introduced by [Ghadimi Lan (2013) can used to show that under the assumption of uniformly bounded gradient estimates, the norm of the gradient of the loss function decreases on average O(t-1/2), when t - co, where is the iteration. This can be viewed as a weak form of convergence that does not guarantee convergence to a fixed point, even to local minimum.
Approximation error as E - O. The nonlinear VQR formulation (8) produces an estimate T
By decreasing E, one can make the problem more exact, as we have shown in practice (fig. AT0b). Following the approach of Proposition A.1 in Genevay et al. (2016), it can be shown that this estimate approaches the true CVQF as E - 0.
20

---

Published as a conference paper at ICLR 2023
C CONDITIONAL VECTOR QUANTILE FUNCTIONS
C.1 INTUITIONS ABOUT VECTOR QUANTILES
To interpret the meaning of CVQFs, us consider the 2-dimensional case, where we assume Y = (Y,, Y2)". Given specific covariates vector æ (#1...., ..). and level u = u2)
may write the components of the conditional vector quantile function as,
Q (u; x) = Q (u; x))
= 1 (my; Y2 (ug; Q (ug: Jx (u1; x)
Q2 u: x)
where QY,IY,, (u; y, x) denotes the scalar quantile function of the random variable Y, at level u, given Y, = y, X 2. Thus, for example, the first component Q1 (u; a)) is a 2D surface where moving along un for fixed u2 yields scalar, non-decreasing function representing the quantiles of Y1 when Y₂ is at a value corresponding to level u2 (see Figure A6b). In addition, the vector quantile function is co-monotonic with the sense defined by eq. (I. For higher dimensions, becomes more involved as the conditioning in each component is on the vector quantile of all the remaining components of the target Y (Figure/A6c)
C.2 CONDITIONAL QUANTILE CONTOURS
useful property of CVQFs is that they allow natural extension a-confidence intervals high dimensional distributions, which we denote a-contours.
Vector quantile contours.
Formally, we define the conditional contour of Y X - at confidence
level 2 x (z), given by
2 = {Qyix(u;=))
E Ua
(22)
d
U₂ = U
i=1
u₁ = {uu.....,
u, E [a, 1- a), u fa,1-a3]
where u-; denotes any component of u (ui...., ua) except for u4 Note that the contour can be calculated not only using the true CVQF. (x;x)x;;))) as above, but also using an estimated CVQF, Oxix(u;;))
For simplicity consider the case. By fixing e.g. i, be one of {a, a} and then sweeping over U2 (and vice versa), we obtain a set of vector quantile levels U' corresponding to the confidence level (see fig. la left and right). Mapping the set Ua back to the domain of Y by using the values the CVQF Q uu; æ) along it, we obtain a contour of arbitrary shape, which contains 100 (1 - 2a) percent of the d-dimensional distribution (see fig. middle and fig. A7 right).
Separable quantile contours. In contrast to vector quantile contours obtained by VQR, which can accurately model distributions with arbitrary shapes, using separable quantiles produces trivial box-shaped contours. Consider a CVQF estimated with separable quantiles presented in eq. (10).
Due to the dependence of each component of Q (u; æ) only on the corresponding component of u, the shape of the resulting quantile contour will always be box-shaped (see fig. A7 middle).
Such trivial contours result in inferior performance for applications of uncertainty estimation, where a confidence region with a the smallest possible area for a given confidence level is desired.
21

---

Published as a conference paper at ICLR 2023
(wn)
1.4
1.2
1.0
0.8
0.
0.4
0.2
0.0
-02
02
00
14
20
28
32
38
.4
9.50
58
.62
68
74
80
.88
92
39
(a)
Q; Up]
Oz(0, Uz!
uy=2.02
ky=0.02
11
0,10.1
0.-0.34
0,-0.34
8,50.5
uy50.50
0,6066
0100.6
0,80.82
Mj =,
0,90.9
0ym0.96
0.0
0.2
0.4
0.5
0.5
1.0
0.0
0.2
0.5
0.8
U2
20
11.5
11.0
25
11.0
10.5
0.
10.5
100
22
(b)
(c)
Figure A6: Visualization (a) 1d, (b) 2d, and (c) 3d vector quantile functions estimated on the
MVN dataset. In each plot, (u(u) the ith component of the vector quantile Q plotted over all quantile levels u. The number of quantile levels calculated was 50, 25 and 10 for 1d, 2d and 3d quantiles respectively. For 2d quantiles (b), the top plot shows multiple monotonic quantile curves of one variable while keeping the other at a fixed level.
22

---

Published as a conference paper at ICLR 2023
Figure A7: Visualization of a-quantile contours constructed using separable quantiles and vector quantiles. All the contours presented contain ~ 90% of the points (i.e., a = 0.025). Left (top bottom): samples drawn from bivariate normal, heart-shaped, and star-shaped densities. Middle & Right: black triangles overlaid on the density constitute the a-qunntile contour estimated using separable quantiles and vector quantiles, respectively. The bivariate normal density has zero mean, unit variance and correlation of -0.7. Vector quantile contours accurately capture the shape of distribution whereas the separable quantile contours are always box-shaped due the assumption of independence.
23

---

Published as a conference paper at ICLR 2023
D DATASETS
D.1 SYNTHETIC DATASETS
Below we describe the four synthetic datasets which were used in the experiments.
MVN. Data generated from linear model y Ax + 7, where Ui0 11s A E Raxk is random projection matrix and is multivariate Gaussian noise with a random covariance matrix.
This dataset was introduced by Feldman et al. (2021). and inspired by similar dataset in Carlier et al. (0017). The target variable Y E R² has banana-shaped conditional distribution, and its shape changes non-trivially when conditioned X E R (see the left-most column of the panel presented in Figure A13). The data generating process defined as follows
Conditional-
continuous-valuue covariate
X U,0.8,3.2) ~~ m], & 2m], U[-0.1,0.1]
3 3
UiO, 1
11
7 BX
Yo - cos (Z) + 1) +r sin (ω) + sin (X), Y1 + cos (Ø),
cooss
and then Y [Yo. Y.
Synthetic glasses. This dataset was introduced by Brando et al. (2022). The target variable Y E RR has bimodal conditional distribution. The mode locations shift periodically when conditioned continuous-valued covariate ~ U[0, 1] (fig. Al2ah. The data generating process defined as
21 3TX
Z2 T( + 3X)
~ Beta (a = 0.5,0 5 1)
~ Categorical 1)
Y: 5 (21) 2.5 +
Y2 sin 2.5
Y = 7Y, + +11++(1100 (1-7)
Rotating star. In this dataset, the target variable of Y E R² has a star-shaped conditional distribution, and its shape is rotated by X E R degrees when conditioned on a discrete-valuee valued X, taking values in [0, 10, 20, ... 60). Data is generated based
600 600 binary image of star. See the first column Figure A144 to visualize conditional distributions as a function of X. Since conditional distributions differ only by rotation, E [Y[X] remains the same for all X. However, the shape of the distribution changes substantially with X, especially in tails. Thus, this dataset is challenging candidate for estimating vector quantiles which must also represent these tails in order to properly
recover the shape of the conditional distribution.
D.2 REAL DATASETS
We perform real data experiments on the b blog_data, bio, house, meps__0 datasets obtained from (Feldman et al. 2021). The original real datasets contained one-dimensional targets. [Feldman et al.] (2021) constructed an additional target variable by selecting a feature that has high correlation with first target variable and small correlation to the other input features, so that hard to predict. Summary of these datasets is presented in table A1
Dataset Name
N
k
Additional Response
blog_data
52397
279
Time between the blog post publication and base-time
bio
45730
8
2
F7 Euclidean distance
house
21613
Latitude a house
meps_20
17541
138
Overall rating of feelings
Table Al: Information about the real data sets.
24

---

Published as a conference paper at ICLR 2023
E EVALUATION METRICS
E. METRICS FOR SYNTHETIC DATA
The quality of an estimated CVQF can be measured by how well it upholds strong representation and co-monotonicity (eqs. ① and (2). However, measuring strong representation requires knowledge of the ground-trune truth joint distribution P(x.Y) or its conditional quantile function. With synthetic data, adherence strong representation can be measured via several proxy metries. The second key property, violations of co-monotonicity, can be measured directly. Below we describe the metrics we use evaluate these properties on estimated conditional quantile functions.
Entropy of inverse CVQE. strong representation holds, the inverted CVQF, x XY must result in a uniform distribution, when evaluated on samples drawn from the true conditional distribution. As measure of uniformity, calculate normalized entropy of the inverted CVQE The inversion procedure done as follows.
1. Sample L evaluation covariate vectors, {am) values random from the ground truth distribution Px.
2. For each 2t,
24- (b) For each ym.,, find the level of the closest vector quantile w.r.t the Euclidean distance, i.e.
m=1 (a) Sample M points (ym.)) } m=1 from Y
um arg min w
ym - @xi
(c) For E {I....,TT) TT] denote &; as the number of times that the quantile level is found in {um,) M m= and calculate p; c//M.
(d) Calculate the entropy hy (exp (hr) 1) 1).
E; pi logp; and the normalized entropy by
3. Report the CVQF inverse entropy F, Σ, hu-
Note that,
1. The entropy metric normalized such that its values are in [0, 1] where corresponds to uniform distribution, and 0 corresponds to a delta distribution.
2. Some non-uniformity arises due to the quantization of the quantile level grid into discrete levels per dimension. We report reference entropy, calculated uniform sample of M quantile-level grid points.
3. We used L = 20 and M 4000.
Distribution distance (KDE-LI). Since the CVQF fully represents the conditional distribution P(xxx) can serve as a generative model for it. We used inverse-transform sampling to generate data from the estimated conditional distribution though the fitted VQR model, Q (u; x), follows.
1. Sample L evaluation covariate vectors, { values random from the ground truth distribution Px.
2. For each æt.
(a) Sample M quantile levels { m=1 uniformly.
(b) Generate M samples from the estimated distribution of Y using the estimated
CVQF: { (úm.) @Qxxx(ummt 2i1)
m=1" M
(c) Sample additional M points ymm 1 m=1
from the ground truth conditional distribu-
tion Py P
(d) Calculate Kernel Density Estimate (KDE) fyix=a; from the VQR samples (ūm,))-
25

---

Published as a conference paper at ICLR 2023
(c) Calculate the KDE fy Xyp== from the ground truth samples { (mm..)
3. Calculate the KDE-L1 metric as
L
Σ L
Frxx-m - BXX==m 11.
The KDEs are calculated with 100 bins per dimension. An isotropic Gauassian kernel was used, with = 0.1 for the conditional-banany dataset and 0.035 for the star dataset. We used pykeops (Charlier et al., 2021) for fast implementation of high-dimensional KDEs. We used L = 20 and M 4000.
Quantile function distance (QFD). This metric measures the distance between true CVQF,
Oxix and an estimate for it obtained by VQR, Qyxx.
1. Sample L evaluation covariate vectors, {ιi)/= values at random from the ground truth distribution Px.
2. For each æt.
(a) Sample M points {ym.))mm)
M from the ground truth conditional distribution P
(b) Estimate an unconditional vector quantile function on fymm m=11 quantile estimation, not regression. Denote the estimated unconditional vector quantile function as x This serves as proxy for the ground-truth conditional quantile function.
) i.e. perform vector
(c) Denote the estimated conditional quantile function evaluated at æ;; Q Q
(d) Compute the normalized difference between them elementwise over each the Td discrete quantile levels, i...
d = 0 @X x= 1
3. Calculate the QFD metric as Σt= di.
We used L 20 and M 4000.
Percentage of co-monotonicity violations (MV). This value can be measured directly. Given an estimated vector quantile function (uu), with T levels per dimension, there are in total 22dd quantile level pairs. Thus, we measure
I (m: (u: ( - @(1,)) < <0)
T
E.2 METRICS FOR REAL DATA
With real data, only finite samples from the ground truth distribution are available. In particular, since real datasets have continuous X, this means that never have more than sample from each conditional distribution Y X æ. Strong representation can therefore not be quantified directly with the above metrics. Instead, opt for metrics which can be evaluated finite samples and are specifically suitable for our chosen application of distribution- uncertainty estimation (section 6.3)
Size of a conditional a-confidence set (AC). This metric approximates [Ca (æ)] i.c. the size of an a-confidence set constructed for a specific covariate æ. Lower values of this metric better, as they indicate that the confidence set is a better fit for the shape of the data distribution because, intuitively. the implication is that the same proportion of the data distribution can be represented by smaller region. The metric is computed as follows.
1. Estimate CVQF (u: a:), R e.g. via VQR (7). NL-VQR 8 separable quantiles 10)
26

---

Published as a conference paper at ICLR 2023
2. Choose test covariate, denoted by x².
3. Construct the corresponding conditional a-contour, x (x"), as defined by eq- (22).
using the estimate
4. Construct the corresponding conditional a-confidence set C, from the contour
xx=") by calculating convex hull of the points within the contour.
5. Calculate the value the metric as the volume of the resulting convex hull (area for 2).
We note that using convex hull is only linear approximation of the true (a)) defined by the points in the contour. For reasonable values T, we find it to be good approximation, and it is an upper-bound on the area/volume of the true C, (x). In practice we scipy with qhull (Virtanen) Barber et al.) 1996) to construct d-dimensional convex hulls and measure their volume.
et al. 2020
Marginal Coverage (MC). This metric measures the proportion of unseen data points that are contained within the the conditional a-confidence sets, C. (z), that was obtained from given estimated CVQE. computed as follows.
1. Estimate a CVQF Qvpx (u; x) using one of the aforementioned methods on training set sampled from PyIx-
2. Denote disjoint held-out test
{ { 1イ TT p
Prix.
3. Measure the marginal coverage
Nr NT Eu E Co()))
where Ca(x]] is constructed as convex hull of the conditional a-contour 2x (x7 constructed using the estimated x a)) (eq. (2
27

---

Published as a conference paper at ICLR 2023
F ADDITIONAL EXPERIMENTS
F.1 SCALE
Scale with respect to d and k. Figure A8 demonstrates the runtime of our solver with respect the number target dimensions d and covariate dimensions
Time (mins.)
(mins.)
10
20
0
10
10*
larget dimensions
Covariate dimensions
Figure A8: Effect of d and on runtime. VQR solver runtime shown as a function number of target dimensions d (left) and covariate dimensions k (right). Calculated on the MVN dataset with N 10k; for experiment 7 10; for k experiment T and exponentially with as expected, due to having T quantile levels. For k. we can see that runtime remains relatively constant even for hundreds dimensions, then increses due to memory constraints.
Runtime scales
Comparison to ceneral-purpase solver. To showcase the scalability of our solver, we compare its runtime against an off-the-shelf linear program solver, ECOS, available within the CVXPY package (Diamond & Boyd, 2016). These results are shown in Figure A9 As expected, our custom VQR solver approaches the performance of the general purpose solver which achieves slightly more accurate solutions due to solving exact problem. However, the results demonstrate that with linear program solver, the runtime quickly becomes prohibitive even for small problems. Crucially, is important note that linear program solvers only allow solving the linear VQR problem, while our solver supports nonlinear VQR.
250
VQR
CVX
VOR
CVX
0.60
70
VOR
CVX
VQR
CVX
1.8
1.6
60
200
0.55
50
1.4
Time (min)
Timo (min)
0.50
122
KDE
30
1.02
1.00
20
0.8
50
0.40
10
0.6
0.35
0.4
2000
4000
6000
8000
10000
10
20
30
40
50
N (samples)
7 (quantile levels per dimension)
Figure A9: Our proposed solver is orders of magnitude faster than general-purpose linear program solver, and approaches the exact solution the problem size increases. VQR and CVX (using ECOS LP solver) solver runtimes (solid lines) are shown function of number of samples N (left) and quantile levels per dimension T (right). Dashed lines indicate KDE-LI as measure of solution quality. CVX obtains the ideal solution due to solving exact problem. Calculated on the MVN dataset 2, 10. Left: T 30; Right: N = 2000.
F.2 OPTIMIZATION
To further evaluate our proposed VQR solver, we experimented with various values & and batch sizes of samples and quantile levels, denoted Bx and BT respectively. Figure 10a shows the effect of the mini-batch size in both and T on the accuracy of the solution. We measured the QFD metric, averaged over 100 evaluation values. These results presented in Figure|A100 As expected we observe improved accuracy when increasing batch sizes (both By and Br) and when decreasing 8.
28

---

Published as a conference paper at ICLR 2023
6o0500
Bv=625
104
6g-1000
109
Bv-1250
6,=2500
6,-2500
04-3000
Rp-10000
10-1
10
Iterations
Iterations
(a) Effect of batch sizes Bx and By on optimization. Computed on MVN dataset with d and
k 1. The vertical axis is the QFD metric defined in section Left: Sweeping over Bx: T = 25.
Right: Sweep over Br; T = 50 produces 502 total levels.
1-1000
--0 0022
c=0,00
c-0
100001
0-2
02
Iterations
Iterations
(b) Decreasing in the relaxed dual improves strong representation. Y-axis: QFD metric. Computed on MVN dataset with N 20k, k 1. Right: T 50, d 2: Left: T 100, 1.
Figure A10: Optimization experiments showing the effect of and both batch sizes on convergence.
F.3
VMR
Figure 4 shows the effectiveness of our VMR procedure. We can see that when applying VMR, there are no monotonicity violations, even for low values of E, and, moreover, QFD metric improves. Figures A12 and A12c demonstrate how VMR can remove monotonicity violations, even when the distribution cannot be estimated correctly (in this case due to linear VQR). Finally, fig. All visualizes the co-monotonicity violations in 2d.
0
20
20
40
60
60
80
25 50
75
25 50 75
Quantile level index pairs
Quantile level index pairs
Figure All: Visualization of monotonicity violations in 2d. Left: co-monotonicity matrix, showing the value of (u; - (Q(u.) - Quu,)) for all i.j, after fitting VQR without VMR with d = 2, T = 10 on the MVN dataset. Right: quantile-level pairs i,j where co-monotonicity is violated. VMR resolves all these violations.
F.4 NONLINEAR VQR
Here we include additional results comparing linear and nonlinear VQR for various datasets. Fig- ure AT2c compares linear and nonlinear VQR the synthetic glasses dataset. Figures A and table A3 present the results of linear and nonlinear VQR on the conditional-banana and rotating star datasets for many values of æ.
29

---

Published as a conference paper at ICLR 2023
10
10
10
5
-5
0.2
0.4
0.6
0.8
0.4
0.6
0.8
0.2
0.4
0.8
X
X
(a) Ground truth joint distribution
(b) VQR
(c) VQR + VMR
10
10
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
X
X
NL-VQR
NL-VQR VMR
Figure A12: VQR and NL-VQR can be used for Simultaneous Quantile Regression (SQR), and VMR eliminates quantile crossings. SQR on the synthetic glasses dataset (a; gray points show ground-truth distribution). Linear VQR fails to correctly model the conditional distribution (b, c; gray points sampled from linear VQR), while nonlinear VQR reconstructs the shape exactly (d, e; gray points sampled from nonlinear VQR). In both cases case, VMR successfully eliminates any quantile crossings (c, e).
Table A2: Quantitative evaluation of linear VQR, CVAE and nonlinear VQR the conditional- banana dataset. Evaluated on N = 4000 samples with = 50 levels per dimension. The KDE-LI, QFD and Entropy metrics defined in section Arrows t/4 indicate when higher/lower better; Numbers are mean + std calculated over the 20 values of æ. Reference entropy is obtained by uniformly sampling the 2d quantile grid.
CVAE
VQR
NL-VQR
KDE-L1
0.175 + 0.021
0.873 + 0.175
0.135
0.024
Quantile Function Distance (QFD) Inverse CVQF Entropy (Ref: 0.70)
0.179 + 0.034
0.055 + 0.023
0.309 + 0.077
0.560 + 0.014
Train time (min)
35
3.5
4
30

---

Published as a conference paper at ICLR 2023
Ground truth
VQR
C-VAE
NL-VQR
x=1.0
1.26
0.21
0.19
x=1.4
0.74
0.15
0.15
V
V V
x=1,8
1.01
0.17
V
VVV
V
V
X=2.2
0.97
0.17
0.12
X=2.6
V
V
V
V
0.58
0.17
0.13
x=3.0
V
V
V
0.96
0.19
0.11
Figure A13:
Qualitative evaluation of VQR, CVAE and NL-VQR the conditional- dataset. Depicted are the kernel-density estimates of samples drawn from each of these methods. Models were trained with N 20k with 50 quantile levels per dimension. Numbers depict the KDE-LI metric. Better viewed as GIF provided in the supplementary material.
Table A3: Quantitative evaluation of linear VQR, CVAE, and nonlinear VQR on the rotating stars dataset. Evaluated N = 4000 samples with T 50 levels per dimension. The KDE-LI, QFD and Entropy metrics are defined in section Arrows f/4 indicate when higher/lower better; Numbers mean + std calculated uniformly sampling the 2d quantile grid.
the 20 values of æ. Reference entropy is obtained by
CVAE
VQR
NL-VQR
KDE-L1
0.95 + 0.04
0.35 + 0.06
0.20 + 0.01
Quantile Function Distance (QFD) Inverse CVQF Entropy (Ref: 0.69)
0.18
0.04
0.08
0.01
0.51 = 0.02
0.59 + 0.02
Train time (min)
35
3.5
31

---

Published as a conference paper at ICLR 2023
Ground truth
VQR
C-VAE
NL-VOR
0.41
0.20
x=10
X=
0.35
x=30
0.42
0.93
0=40
0.31
0.99
X=50
0.28
0.95
0.22
x=60
0.38
1.03
0.19
Figure A14: Qualitative evaluation of linear VQR, CVAE and nonlinear VQR the rotating stars dataset. Depicted are the kernel-density estimates samples drawn from each of these methods. Models were trained with N 20k with T = 50 quantile levels per dimension. Numbers depict the KDE-LI metric. Better viewed as GIF provided in the supplementary material.
32

---

Published as a conference paper at ICLR 2023
Dataset: blo
Dataset: blog
conficeecc region
Requested Coverage
Requested Coverage
Area of the confirence bester)
80
Marginal Coveragee
80
Margina Coveragi
70
70
50
NL-VOR
NL-VOR VOR
Sep-NIOR
40
Datasel: House Prices
Datasel: MEPS-20
Requested Coverage
Requested Coverage
Area confidencc
(lower is better))
80
Marginal Coverage
tf confidience (lower better)
90
Marginal coverage
70
70
50
NL-VOR
Sep-NLOR
Sep-QR
40
NL.VOR
VQR
Sep-NLOR
Sep
40
Figure A15: Nonlinear VQR produces substantially smaller a-confidence sets in real data experiments. Comparison nonlinear VQR (NL-VQR), separable nonlinear QR (Sep-NLQR). linear VQR (VQR), and linear separable QR approaches (Sep-QR) on the four real datasets presented in table AT The confidence sets are constructed for each method in such a way to produce same desired level of marginal coverage (blue dashed line). Error bars generated from 10 trials, each with a different randomly-sampler train/test split.
MLAGE-O (190)
SE-GR-1300 +3,00 M.Sts--n
VOK (021)
1.0
-0.
-1.
-2.
-2.
MLACk (L90
N_Sip-on 10.
2.0
814409 (-40)
8L-ca-on (L
MIVIO VOR Sep-OW
2.
-35
-25 -.0
22
Figure A16: Nonlinear VQR produces a-quantile contours which model the shape of the conditional distribution and have significantly smaller area compared to other methods. Each column depicts a-quantile contours for two test covariates (top, bottom) sampled from bi (left), house pri ces (middle), and blog_da og_data (right) datasets. Separable QR approaches produce box-shaped confidence regions due to the assumption of statistical independence, whereas VQR and nonlinear VQR produce contours with non-trivia shapes that adapt the test covariate hand (each column, top vs bottom). The areas of the a-confidence sets are reported relative to NL-VQR in the legend.
33

---

Published as a conference paper at ICLR 2023
G IMPLEMENTATION DETAILS
Solver implementation details. The log-sum-exp term in eq. E decreases, which can be mitigated by implementing it as
becomes numerically unstable as
logsumexp, = logsumexp,
max + m
Scale and Optimization experiment. For both scale and optimization experiments, run VQR for 10k iterations and use learning rate scheduler that decays the learning rate by a factor of 0.9 every 500 iterations if the error does not drop by 0.5%.
Synthetic glasses experiment. We set N = 10k, T = 100, and = 0.001. We optimized both VQR and NL-VQR for 40k iterations and use learning rate scheduler that decays the learning rate by a factor of 0.9 every 500 iterations if the error does not drop by 0.5%. In NL-VQQR, go, we used a 3-layer fully-connected network with each hidden layer of size 1000. We used skip-connections, batch-norm, and ReLU nonlinearities between the hidden layers. For NL-VQR and VQR, set the initial learning rate to be 0.4 and 1, respectively.
Conditional Banana and Rotating Star experiments. We draw N 20k samples from P P(X.Y) and fit T = 50 quantile levels per dimension (d = 2) for linear and nonlinear VQR (NL-VQR).. Evaluation is performed by measuring all aforementioned metrics on 4000 evaluation samples from 20 true conditional distributions, conditioned on = [1.1, 1.2,
3.0]. For NL-VQR, as ge we use a small MLP with three hidden layers of size (2, 10, 20) and a ReLU non-linearity. Thus go lifts X into 20 dimensions on which VQR is performed. We set ε = 0.005 and optimized both VQR and NL-VQR for 20k iterations. We used the same learning rate and schedulers as the synthetic glasses experiment.
Real data experiments. In all the real data experiments, we randomly split the data into 80% training set and 20% hold-out test set. Fitting performed on the train split, and evaluation metrics, marginal coverage and quantile contour area (calculated as reported in appendix E.2). are measured on the test split. We repeat this procedure 10 times with different random splits. The baselines we evaluate are NL-VQR, VQR, separable nonlinear QR (SeR-NLQR), and separable linear QR (Sep-QR). In the separable baselines, two QR models are fit separately, one for each target variable. For both nonlinear separable QR and NL-VQR baselines, we choose to be an MLP with three hidden layers. In the case of NI.- VQQ, hidden layers sizes are set to (100, 60, 20) and for nonlinear separable QR, they are set to (50, 30, 10). All methods were run for 40k iterations, with learning rate set to 0.3 and и = 0.01. We set = 50 for NL-VQR and VQR baselines, and T = 100 for separable linear and nonlinear QR baselines. The parameter that determines the a-quantile contour construction calibrated separately for each method and dataset. The goal of the calibration procedure is to achieve consistent coverage across baselines in order to allow for the comparison of the a-quantile contour areas. Table A4 presents the calibrated a's that were used.
Dataset
NL. VQR
VQR
Sep-NLQQ
Sep-QR
a
Cov.
Area
a
Cov.
Area
a
Cov.
Area
a
Cov.
Area
bio
0.05
85.00
1.57
0.03
82.87
2.49
0.01
85.13
3.18
0.02
85.76
3.59
aata
0.05
84.17
1.19
0.05
81.50
3.41
0.025
84.60
3.02
0.025
71.10
4.73
house
0.05
82,29
1.69
0.06
80.75
3.39
0.025
82.51
2.53
0.03
82.48
4.76
meps_ 20
0.02
74.13
1.81
0.06
74.82
4.04
0.01
74.73
2.78
0.01
75.77
9.70
Table A4: Choice of a for different real datasets. Nominal (requested) coverage for bi and -at_ was 85%. For house and meps_ 20 datasets, it was 82% and 75%, respectively.
Machine configuration. All experiments were run on a machine with an Intel Xeon E5 CPU, 256GB of RAM and an Nvidia Titan 2080Ti GPU with 11GB dedicated graphics memory.
34

---

Published as a conference paper at ICLR 2023
H SOFTWARE
We provide our VQR solvers as part of a comprehensive open-source python package vęr which installed with pip install vqr. The source code is available at /https://github.c. com/ vistalab Sestalab-techhion///S. twchnion/vgrreeee and it contains several jupyter notebooks which demonstrate the of our package.
import numpy as np from ver Import VectorquantileRgrreess from w.w.sollves..rgulaaiiee...
import Regulaarrzeeiiaal@@@oo...eer
# Create the VQR solver and regressor.
vgr_solver ReguaarieeDaalyBSS0 verbose=True, T=50, epsilcn=le-2, eum_pcocs=1 1r=0.9
vgr VectorQuant (erelaantleeeeeeesseeee TeRegressor ver
Pit the model some training data.
vgr. oqr.fit(ee (X_rrai, Y_train)
Sample from the fitted conditicnal distribution, given specific x. Y sampled = vgr, sample (n=100, x=X_test[0]
calculatc coverage on the samples. cov_sampled vqr. coverage (X_sampled, x=X_test [0), alpha=alpha)
Listing 1: Minimal example of using the vqz library. Demonstrates instantiation of solver and regressor, fitting VQR, sampling from the fitted CVQF and coverage calculation.
We note that although reference Matlab implementation of linear VQR by the original authors exists) this implementation relies on solving the exact formulation (5) using general-purpose linear program solver. We found it be prohibitively slow to use general-purpose linear program solvers for VQR (fig. A9p. For example, we noted runtime of 4.5 hours for N 7k, d = 2, k = 10, T = 30, and for N 10k the solver did converge even after hours. Thus, this approach is unsuitable even for modestly-sized problems. Moreover, when using such solvers, only linear VQR can be performed. For comparison, our proposed solver converges to a solution of the same quality on problems of these sizes in less than one minute and supports nonlinear VQR.
35


Reviews:
[{"review_id": "s_jsmIrNA-", "reviewer": "Reviewer_fWcC", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper introduces an efficient, nonlinear model for vector quantile regression (VQR). VQR is the problem of learning the quantiles of a multivariate response variable Y, conditioned on set of covariates x. The approach extends prior work by allowing for a nonlinear function of x and implementing vector rearrangement to guarantee monotonicity of simultaneous quantile estimates.\n\nI am not familiar with Optimal Transport, so unfortunately, a fraction of the paper's key ideas were unintelligible to me.", "summary_of_the_review": "I really like this topic and feel that it deserves more attention. Unfortunately, the audience for this paper would be relatively small. However, I believe it extends the state of the art quite a bit for this problem, and as such it probably deserves to be published.", "strength_and_weaknesses": "Strengths:\n* If the authors' literature review is taken to be comprehensive, then they have advanced the state of the art for vector quantile regression (VQR) by quite a large amount: they introduce the first (partially) nonlinear model for VQR, they introduce a vector variation of monotonic rearrangement, and they provide an efficient and scalable implementation.\n\nWeaknesses:\n* If I understand Section 4 correctly, the model does not allow for interactions between x and u. This seems to imply that the model can only model situations where each conditional distribution over Y (for a given x value) is exactly the same as the others (for all other x values), up to an affine scale-and-shift transformation. This seems like an overly-restrictive model, but I guess prior work, using a fully linear model in both u and x, is even more restrictive.\n* Vector quantile regression seems like quite a niche topic. In my experience in corporate research, even scalar-valued quantile regression is little-known and little-used.\n* There doesn't seem to be a comparison to, or a mention of, Bayesian methods or generative modeling, which seems like a very natural fit for this problem.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "rbfoDkTTF1", "reviewer": "Reviewer_2xhV", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this work, motivated by the limitation that there are no fast or scalable solvers for vector quantile regression (VQR), the authors provide a highly-scalable solver for VQR that relies on solving its relaxed dual formulation. In addition, the authors propose vector monotone rearrangement (VMR), which serves as a refinement step and resolves the co-monotonicity violations in estimated conditional vector quantile functions (CVQFs). Moreover, motivated by the limitation of VQR that it assumes a linear model for the quantiles of the target given the features, the authors propose nonlinear vector quantile regression (NL-VQR). Based on extensive experiments on both synthetic and real data, the authors demonstrate the superiority of VMR and NL-VQR in modeling CVQFs accurately. ", "summary_of_the_review": "Overall, I think that this is a good work with sufficiently novel methods and interesting numerical results.", "strength_and_weaknesses": "Strength:\n\nThe studied problem is of sufficient interest to the community, and the authors' work seems to be a significant contribution to the field of VQR. \n\nWeaknesses:\n\nAlthough the experimental results are sufficiently promising, I hope that similar to some closely relevant works such as Carlier et al. (2016) and Chernozhukov et al. (2017), the authors can provide some theoretical justifications for their proposed methods.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "9MX-jqbfNQy", "reviewer": "Reviewer_u5QX", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a scalable solver for non-linear multivariate quantile regression. While vector quantiles have been defined and studied in previous work, the authors address several limitations: developing a scalable solver, extending the linear regressor model to allow non-linear models and making sure monotonicity is ensured.", "summary_of_the_review": "I believe this is an interesting contribution even though writing and structure could be improved significantly. Also, I am not too familiar with the background on vector quantiles, but appreciated the conformal prediction results.", "strength_and_weaknesses": "Strengths:\n- Clean abstract in terms of the limitations addressed and the included contributions.\n- Good visualization/illustration of vector quantiles in Figure 1.\n- Brief summary of previous work in introduction and discussion of related work throughout the paper.\n- Even though I am not too familiar with this line of work, vector quantiles seem to be relevant in several niche applications and this paper provides a more scalable solver. With the increased interest in distribution-free uncertainty estimation, I can imagine this becoming and interesting tool.\n- The extension to non-linear models seems to be a particularly important contribution given the experiments and I can imagine this to be impactful in this line of work.\n- I appreciate the experiments on conformal prediction that show the benefit of this method.\n- Providing an open-source implementation.\n\nWeaknesses:\n- The paper is generally more difficult to follow, mainly due to two reasons: the content being very technical and dense (requiring a lot of familiarity with the topic or time to go through equations by trial and error) and missing structure. For example, in section 3, there is a high-level outline/overview missing. The primal discussion, as example, seems less relevant and not be used (more like a comment on related work), while the dual formulation and SGD approach is actually used in the proposed method.\n- The discussion of SQR and prior attempts in section 4, while interesting, seems less important in the main paper and distracted me. I’d prefer this space to be spent on more structure and a bit more guidance through the technical content.\n- When less familiar with this line of work, the experimental setup could also use more details – e.g., describe at least one metric properly and give a high-level idea of the datasets.\n- The details of the contribution over Genevay 2016 in terms of the solver is a bit unclear to me. As of my understanding the SGD optimization setup of that paper is used with the “trick” of obtaining constant memory footprint. Is that correct? This could be made clearer throughout the paper. In the abstract/intro it sounds like a novel solver is proposed instead.\n- I couldn’t find the exact non-linear model that is used in the main paper and would be interested in a discussion of the constraints here. Are generally any differentiable non-linear models supported?\n- I am also missing a baseline in terms of runtime. In Figure 2, how would previous solvers perform? When would they stop working?\n- While I appreciate the conformal prediction experiments, this part is difficult to follow and judge from the main paper as the results are exclusively found in the appendix.\n- Generally, this paper and the experiments are difficult to read without the appendix and I think this hurts the paper. I would appreciate if the authors could prioritize experiments and discuss a subset where the results are actually in the main paper, instead.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster