  Quantile regression (QR) is a powerful tool for estimating one or more
conditional quantiles of a target variable $\mathrm{Y}$ given explanatory
features $\boldsymbol{\mathrm{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 $\boldsymbol{\mathrm{Y}}$ given the
features $\boldsymbol{\mathrm{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.
