domain: stat math

text:
Eustasio Del Barrio

Alberto González-Sanz
email: alberto.gonzalezsanz@math.univ-toulouse.fr

Marc Hallin
email: mhallin@ulb.ac.be

Nonparametric Multiple-Output Center-Outward Quantile Regression *

Keywords: Multiple-output regression, Center-outward quantiles, Optimal transport

), we are considering the problem of nonparametric multiple-output quantile regression. Our approach defines nested conditional center-outward quantile regression contours and regions with given conditional probability content irrespective of the underlying distribution; their graphs constitute nested center-outward quantile regression tubes. Empirical counterparts of these concepts are constructed, yielding interpretable empirical regions and contours which are shown to consistently reconstruct their population versions in the Pompeiu-Hausdorff topology. Our method is entirely non-parametric and performs well in simulations including heteroskedasticity and nonlinear trends; its power as a data-analytic tool is illustrated on some real datasets.

Introduction

Quantile regression, single-and multiple-output

Forty-five years after its introduction by Koenker and Bassett (1978), quantile regressionarguably the most powerful tool in the statistical study of the dependence of a variable of interest Y on covariates X = (X 1 , . . . , X m )-has become part of statistical daily practice, with countless applications in all domains of scientific research, from economics and social sciences to astronomy, biostatistics, and medicine. Unlike classical regression, which, somewhat narrowly, is focused on conditional means E[Y |X], quantile regression indeed is dealing with the complete conditional distributions P Y|X=x of Y conditional on X = x. Building on that pioneering contribution, a number of quantile regression methods, parametric, semiparametric, and nonparametric, have been developed for an extremely broad range of statistical topics, including time series, survival analysis, instrumental variables, measurement errors, and functional data-to quote only a few. Sometimes, a simple parametrized regression model allows for a parametric approach, yielding, for instance, linear quantile regression. In most situations, however, parametric models are too rigid and a more agnostic nonparametric approach is in order. We refer to Koenker (2005) for an introductory text and to Koenker et al. (2018) for a comprehensive survey.

In single-output models (univariate variable of interest Y ), this nonparametric approach is well understood and well studied, and the history of non-parametric estimation of conditional quantile functions goes back, at least, to the seminal paper by [START_REF] Stone | Consistent Nonparametric Regression[END_REF]. The results are much scarcer, however, in the ubiquitous multiple-output case (d-dimensional variable of interest Y, with d > 1), and the few existing ones are less satisfactory-the simple reason for this being the absence of a fully satisfactory concept of multivariate quantiles.

A major difficulty with quantiles in dimension d > 1, indeed, is the fact that R d , contrary to R, is not canonically ordered. A number of attempts have been made to overcome that issue, the most remarkable of which is the theory of statistical depth. That theory has generated an abundant literature which we cannot summarize here-we refer to [START_REF] Serfling | General notions of statistical depth function[END_REF] or [START_REF] Serfling | Quantile functions for multivariate analysis: approaches and applications[END_REF](Serfling, , 2019a,b) ,b) for general expositions and authoritative surveys.

Several depth concepts coexist. The most popular of them is Tukey's halfspace depth [START_REF] Tukey | Mathematics and the picturing of data[END_REF], but all depth concepts (including the most recent ones: see, e.g., Konen and Paindaveine (2022)) are sharing the same basic properties. Tukey's halfspace depth characterizes, for each distribution P over R d (for simplicity, assume P to be Lebesgueabsolutely continuous) depth regions D P (δ) (resp., depth contours D P (δ)) as collections of points with depth (relative to P) larger than or equal to δ (resp., equal to δ), δ ∈ (0, 1/2]. Depth regions are convex, closed, and nested as δ increases, and have been proposed as notions of quantile regions and contours-an interpretation that is supported by the L 1 nature of Tukey depth (Hallin, Paindaveine and Šiman, 2010).

Among the merits of this interpretation is that it has imposed the idea that quantiles, in dimension d ≥ 2, should rely on some center-outward ordering with central region of depth δ = 1/2 rather than a southwest-northeast extension of the classical univariate "leftto-right" linear ordering of the real line. Unfortunately, depth regions fail to satisfy the quintessential property of quantile regions: the P-probability content P D P (δ) of the quantile region D P (δ) indeed very much depends on P. This is not a minor weakness: the univariate median Y P 1/2 of an absolutely continuous distribution P, for instance, is characterized by the fact that P (-∞, Y P 1/2 ] = 1/2 irrespective of P-who would call median a quantity Y P med such that P 1 (-∞, Y P1 med ] = 0.4 while P 2 (-∞, Y P2 med ] = 0.6 ? None of the depth concepts in the literature is satisfying that essential property of quantiles, though, and depth regions, therefore, cannot be considered as bona fide quantile regions.

A review of the various existing multiple-output quantile regression models (linear and nonparametric, depth-based and others) can be found in Hallin and Šiman (2018). A nonparametric quantile regression model based on a directional form of Tukey depth is developed in Hallin et al. (2015) but suffers the same lack of control over the probability contents of the quantile regions involved as the depth-based quantile concept itself. So does also the directional concept of M-quantiles proposed by Merlo et al. (2022).

Recently, based on measure transportation ideas, new concepts of quantiles in dimension d > 1 have been introduced in [START_REF] Chernozhukov | Monge-Kantorovich depth, quantiles, ranks and signs[END_REF][START_REF] Hallin | On distribution and quantile functions, ranks and signs in R d : a measure transportation approach[END_REF][START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF] under the names of Monge-Kantorovich depth and center-outward quantile function. Center-outward quantile functions define nested closed regions C P (τ ) and continuous contours C P (τ ) indexed by τ ∈ ([0, 1) such that, for any absolutely continuous P, P C P (τ ) = τ irrespective of P. Unlike the previous concepts, thus, and unlike the depth-based ones (as proposed in Hallin, Paindaveine and Šiman (2010) or Kong and Mizera (2012)), these measuretransportation-based quantiles do satisfy the essential property that the P-probability contents of the resulting quantile regions do not depend on P. Moreover, the corresponding quantile regions are not necessarily convex and, as shown in Figure 1, they are able to cap- ture the "shape" of the underlying distribution. We refer to [START_REF] Hallin | Measure transportation and statistical decision theory[END_REF] for a survey of measure-transportation-based center-outward quantiles, the dual concepts of multivariate ranks and signs, and their many applications in inference problems [START_REF] Ghosal | Multivariate ranks and quantiles using optimal transportation and applications to goodness-of-fit testing[END_REF] Motivated by this long list of successful applications, we are proposing in this paper a novel and meaningful solution, based on the concept of center-outward quantiles, to the problem of nonparametric multiple-output quantile regression. Namely, for a pair of multidimensional random variables (X, Y) with values in R m × R d (Y the variable of interest, X the vector of covariates) and joint distribution1 P, we define (Section 2.2) the center-outward quantile map

Q ± of Y conditional on X = x as u ∈ S d → Q ± (u| x) ∈ R d (1.1) (S d the open unit ball in R d ),
with the essential property that, letting

C ± (τ |x) := Q ± (τ S d |x) τ ∈ (0, 1), x ∈ R m , (1.2)
we have

P Y ∈ C ± (τ | x) X = x = τ for all x ∈ R m , τ ∈ (0, 1), and P, (1.3) justifying the interpretation of x → C ± (τ | x)
as the value at x of a regression quantile region of order τ of Y with respect to X. For τ = 0,

C ± (0 | x) := τ ∈(0,1) C ± (τ | x) (1.4)
yields the value at X = x of the regression median x → C ± (0 | x) of Y with respect to X. The same conditional quantile map characterizes nested (no "quantile crossing" phenomenon) "regression quantile tubes of order τ " (in R m+d )

T ± (τ ) := x, Q ± (τ S d x) x ∈ R m , τ ∈ (0, 1) (1.5)
which are such that

P [(X, Y) ∈ T ± (τ )] = τ irrespective of P, τ ∈ (0, 1). (1.6)
For τ = 0, define

T ± (0) := (x, y)) x ∈ R m , y ∈ C ± (0 | x) = τ ∈(0,1) T ± (τ ) (the graph of x → C ± (τ | x))
; with a slight abuse of language, also call T ± (0) the regression median of Y with respect to X.

None of the earlier attempts to define multiple-output regression quantiles-neither the depth-based definitions in Hallin et al. (2015) or Paindaveine and Šiman (2011), the directional concepts of marginal M-quantiles [START_REF] Breckling | M-Quantiles[END_REF] considered by Merlo et al. (2022), nor the measure-transportation-based approach proposed by [START_REF] Carlier | Vector quantile regression: An optimal transport approach[END_REF]-is characterizing quantile regions that satisfy requirements (1.3) and (1.6).

Carlier, [START_REF] Carlier | Vector quantile regression: An optimal transport approach[END_REF] deserves special attention, though, as the first attempt to break with directional and depth-based approaches to multiple-output quantile regression by means of innovative measure transportation ideas. Focusing on linear quantile regression, their method is based on a concept of multivariate quantile functions defined over the unit cube [0, 1] d rather than the unit ball S d . While yielding (under the assumption of a linear regression) an asymptotic reconstruction of the distributions of Y conditional on X = x, however, their choice of the unit cube does not directly allow for the definition of quantile regions of given order similar to the center-outward regression quantile regions C ± (τ | x) ⊂ R d or the quantile regression tubes T ± (τ ) ⊂ R m+d .

As for the depth-based quantile regions, moreover, they are necessarily convex, even for distributions with highly non-convex shapes as in Figure 1. To circumvent this convexity problem, [START_REF] Feldman | Calibrated multiple-output quantile regression with representation learning[END_REF] propose a clever transformation of the data turning its distribution into a latent one with convex level sets by fitting a conditional variational auto-encoder [START_REF] Sohn | Learning Structured Output Representation using Deep Conditional Generative Models[END_REF]). The probability contents of the quantile regions resulting from this machine-learning type of "convexity reparation," however, remain out of control (they still depend on P).

Figure 2 provides, for m = 1 and d = 2, a visualization of the regression median and quantile tubes of orders τ = 0.2, 0.4, and 0.8, along with some of the corresponding conditional regression quantile contours

C ± (τ |x) := x, Q ± (τ S d-1 x) (S d-1 the unit sphere in R d ), for the model Y = Y 1 Y 2 = sin( 2π 3 X) + 0.575 e 1 cos( 2π 3 X) + X 2 + e 3 2 2.3 + 1 4 e 1 + 2.65 X 4 (1.7)
where X ∼ U [-1,1] , e = 1 + 3 2 sin(πX/2) 2 v, v ∼ N (0, Id), and X and v are mutually independent. Actually, since exact values cannot be obtained for these population concepts, a very large sample of n = 400, 040 observations was generated from (1.7), and the consistent estimation procedure described in Section 3 was performed to obtain the picture. Note the non-convexity of the conditional contours for τ = 0.8, the non-linearity of the regression median, and the marked heteroskedasticity of the regression.

Outline of the paper

The paper is organized as follows. Section 2 is dealing with the population concept of conditional center-outward quantile map and the resulting center-outward regression quantile contours and regions, Section 3 with their estimation. In Section 3.1, we show how to construct empirical quantile contours and regions, the consistency of which is established in Section 3. Numerical results are provided in Section 4. Monte Carlo experiments in Section 4.1 show the ability of our method to pick heteroskedasticity, nonlinear trends, and the shape of conditional distributions; comparisons also are made with the results of the depthbased method of Hallin et al. (2015). Real datasets are analyzed in Section 4.2, illustrating the power of our method as a data-analytic tool. Section 5 concludes with references to the recent literature on the numerical aspects of optimal transportation and perspectives for future developments.

Nonparametric center-outward quantile regression

Notation

For convenience, we are listing here the main notation to be used throughout the paper. Unless otherwise stated, we denote by (Ω, A, P) the triple defining the underlying probability space. Let d be the d-dimensional Lebesgue measure, B d the Borel σ-field, and P(R d ) the space of Borel probability measures on R d . The support of a probability P ∈ P(R d ) is denoted as supp(P); its closure as supp(P). Throughout, (X, Y) denotes an R m+d -valued random vector with probability distribution P = P X,Y ∈ P(R m+d ), m-dimensional Xmarginal P X and d-dimensional Y-marginal P Y . The distribution of Y conditional on X = x is denoted as P Y |X=x

2 . The open unit ball, the closed unit ball, and the unit hypersphere in R d are denoted by S d , S d , and S d-1 , respectively. We denote by U d the spherical uniform over S d -that is, the product of a uniform distribution over the unit hypersphere S d-1 (for the directions) and a uniform distribution over [0, 1] (for the distance to the center).

Conditional center-outward quantiles, regions, and contours.

Let us provide precise definitions for the concepts we briefly presented in the Introduction and properly introduce center-outward quantiles, regions, and contours.

For any P ∈ P(R d ), denote by Q ± = ∇ϕ and call center-outward quantile map the (Lebesgue-a.e.) unique gradient of a convex function ϕ : S d → R such that Q ± (U) ∼ P, for any U ∼ U d -in the measure transportation convenient terminology, Q ± is pushing P forward to U d , which we denote as Q ± #P = U d . This only defines Q ± at ϕ's points of differentiability (recall that convex functions are differentiable at almost every point in the interior of their domain: see Theorems 26.1 and 25.5 in [START_REF] Rockafellar | Convex Analysis[END_REF]). At ϕ's points of non-differentiability u, let us define Q ± (u) as the subdifferential ∂ϕ(u) of ϕ, namely,

Q ± (u) = ∂ϕ(u) := y ∈ R d | for all z ∈ R d , ϕ(z) -ϕ(u) ≥ y, z -u , u ∈ S d ;
then, Q ± is an everywhere-defined set-valued function. Slightly abusing the notation, we also write ∂ϕ for the set of all points (u, y) ∈ R m+d such that y ∈ ∂ϕ (u). With this notation, we can introduce the concepts of conditional center-outward quantiles, contours, and regions. Definition 2.1. Call conditional center-outward quantile map of Y given X = x the centeroutward quantile map u → Q ± (u|X = x), u ∈ S d of P Y |X=x , x ∈ R m . The corresponding conditional center-outward quantile regions and contours of order τ ∈ (0, 1) are the sets

C ± (τ x) := Q ± τ S d |X = x and C ± (τ x) := Q ± (τ S d-1 |X = x) ,
respectively. The conditional center-outward quantile maps also characterize (see Definitions (1.4), (1.5), and (1.6)) conditional medians C ± (0 | x) and regression quantile tubes T ± (τ ).

When no confusion is possible, we also write Q ± (u|x) for Q ± (u|X = x). The terminology center-outward regression quantile region, contour, and median is used for the mappings x → C ± (τ x), x → C ± (τ x), and

x → C ± (0 | x), x ∈ R m .
Recall, however, that, in the absence of any assumptions on the conditional probabilities P Y |X=x , the mappings u → Q ± (u |X = x) typically are set-valued, see [START_REF] Rockafellar | Variational Analysis[END_REF]. Whenever continuous, single-valued functions (typically, on the punctured unit ball S d \ {0}) are needed, we will make the following assumption.

Assumption (R)

For P X -a.e. x ∈ R m , the conditional distribution P Y |X=x admits, with respect to the Lebesgue measure, a density p Y |X=x with convex support supp(P Y |X=x ); moreover, for every R > 0, there exist constants

0 < λ x R ≤ Λ x R < ∞ such that λ x R ≤ p Y |X=x (y) ≤ Λ x R for all y ∈ supp(P Y |X=x ) ∩ R S d .
(2.1)

In the classical single-output case (d = 1), consistent estimation of conditional quantiles similarly requires the continuity of the conditional quantile maps (see [START_REF] Stone | Consistent Nonparametric Regression[END_REF]). In dimension d > 1, the continuity of center-outward quantile maps follows from assumptions similar to Assumption (R)-see [START_REF] Figalli | On the continuity of center-outward distribution and quantile functions[END_REF] and del Barrio, González-Sanz and Hallin (2020).

Empirical center-outward quantile regression

We now proceed with the construction of empirical versions of the conditional center-outward quantile concepts defined in Section 2.2 and their consistency properties.

Empirical conditional center-outward quantiles

Let (X, Y) (n) := (X 1 , Y 1 ), . . . , (X n , Y n ) be a sample of n i.i.d. copies of (X, Y) ∼= P XY . In this section, we develop an estimator of the conditional center-outward quantile maps u → Q ± (u|X = x), x ∈ R m . Our estimator is obtained in two steps: in Step 1, we construct an empirical distribution of Y conditional on X = x and, in Step 2, we compute the corresponding empirical center-outward quantile map.

Step 1. For each value of x ∈ R m , our estimation of the conditional distribution of Y conditional on X = x involves a sequence of weight functions w (n) : R m(n+1) → R n measurable with respect to x and the sample X (n)

:= (X 1 , . . . , X n ) of X observations, of the form (x, X (n) ) → w (n) x, X (n) := w 1 (x; X (n) ), . . . , w n (x; X (n) ) (3.1)
where w

(n) j : R m(n+1) → R, j = 1, . . . , n are such that w (n) j (x; X (n) ) ≥ 0 and n j=1 w (n) j (x; X (n) ) = 1 a.s. for all n. (3.2)
We refer to a function w (n) of the form (3.1) satisfying (3.2) as a probability weight function and define the empirical conditional distribution of Y given X = x as

P (n) w(x) := n j=1 w (n) j (x; X (n) )δ Yj , (3.3)
where δ Yj is the Dirac function computed at Y j . Following [START_REF] Stone | Consistent Nonparametric Regression[END_REF], we say that the sequence w

(n) is a consistent weight function if, whenever (X, Y ), (X 1 , Y 1 ), . . . , (X n , Y n ) are i.i.d., where Y is real-valued and such that E|Y | r < ∞ for r > 1, E n j=1 w (n) j (X; X (n) )Y j -E(Y |X) r -→ 0 as n → ∞. (3.4)
Step 2. To estimate the conditional quantiles, consider a regular grid

G (N ) of U d consisting of N gridpoints denoted as G (N ) 1 , . . . , G (N )
N . The number N here is arbitrarily chosen as factorizing into a product of integers of the form N = N R N S + N 0 with N 0 = 0 or 1. That regular grid is created as the intersection between -the rays generated by an

N S -tuple u 1 , . . . , u N S ∈ S d-1 of unit vectors such that N S -1 N S
j=1 δ uj converges weakly to the uniform over S d-1 as N S → ∞, and -the N R hyperspheres with center 0 and radii j/(N R + 1), j = 1, . . . , N R , along with the origin if N 0 = 1. Based on this grid, we define the sequence of discrete uniform measures

U (N ) d := 1 N N j=1 δ G (N ) j ∈ P(R d ), N ∈ N over G (N ) and require that N → ∞ with both N R → ∞ and N S → +∞. By construc- tion, U (N ) d
converges weakly to U d as N → ∞. Note that imposing N 0 = 0 or 1 is not a problem, since N , unlike n, is chosen by the practitioner. Having N 0 = 0 or 1 yields the fundamental advantage that all points of G (N ) have multiplicity one and that Corollary 3.1 in [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF], to be used below, applies.

Our estimation of the conditional center-outward quantile maps relies on the optimal transport pushing U

(N ) d forward to P (n)
w(x) -more precisely, adopting (since typically N = n) the Kantorovich formulation of the optimal transport problem, on the solution of the linear program (solvable using efficient numerical methods such as the auction or Hungarian algorithms-see [START_REF] Peyré | Computational Optimal Transport: With Applications to Data Science[END_REF] and references therein)

3 min π:={πi,j } N i=1 n j=1 1 2 |Y j -G i | 2 π i,j , s.t. n j=1 π i,j = 1 N , i ∈ {1, 2, . . . , N }, N i=1 π i,j = w (n) j (x; X (n) ), j ∈ {1, 2, . . . , n}, π i,j ≥ 0, i ∈ {1, 2, . . . , N }, j ∈ {1, 2, . . . , n}.
(3.5)

Here, any N n-tuple π := {π i,j |i = 1, . . . , N j = 1, . . . , n} satisfying the constraints in (3.5) represents a transport plan from U

(N ) d to P (n) w(x) -that is, a discrete distribution over R d ×R d with marginals U (N ) d and P (n) w(x) . Let π * (x) = {π * i,j (x)|i = 1, . . . , N j = 1, .
. . , n} be a solution of (3.5) (an optimal transport plan). Theorem 2.12(i) in [START_REF] Villani | Topics in Optimal Transportation[END_REF] implies that its support supp(π

* (x)) := {(G i , Y j )|π * i,j (x) > 0} is cyclically monotone, 4
hence is contained in the graph of the subdifferential of a convex function. Therefore, the idea is to construct a smooth interpolation of π * (x) that maintains this property.

Note that, for any gridpoint G i , i ∈ {1, . . . , N }, the constraints in (3.5) imply that there exists at least one j ∈ {1, . . . , n} such that (G i , Y j ) ∈ supp(π * (x)). Since more than one such j may exist, we choose the one which "gets the highest mass" from G i , and in case of ties we choose the "smallest" one: namely, let

T * (G i | x) := arg inf y : y ∈ conv {Y J : J ∈ arg max j π * i,j (x)} , (3.6)
where conv(A) denotes the convex hull of a set A. Since conv

{Y J : J ∈ arg max j π * i,j (x)} is closed and convex in R d , (3.6) defines a unique T * (G i | x). Due to the cyclical monotoni- city of supp(π * (x)), there exists a convex function ϕ * (• x) : R d → R with subdifferen- tial ∂ϕ * (• x) such that, for all 1 ≤ i ≤ N , ∅ = {Y j : (G i , Y j ) ∈ supp(π * (x))} ⊂ ∂ϕ * (G i x).
Since sub-differentials are convex sets, this entails

conv{Y j : (G i , Y j ) ∈ supp(π * (x))} ⊂ ∂ϕ * (G i x). Consequently, {(G i , T * (G i | x)) : i = 1, . . . , N } is
cyclically monotone and satisfies the assumptions of Corollary 3.1 in [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF]. This implies the existence, for all x, of a continuous cyclically monotone map u → Q

(n) w,± (u x), say, interpolating the N -tuple (G 1 , T * (G 1 | x)) , . . . , (G N , T * (G N | x)), i.e., such that Q (n) w,± (G i | x) = T * (G i | x) for i = 1, . . . , N .
In particular, we proceed as in [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF] by choosing the smooth cyclically monotone interpolation with largest Lipschitz constant. This continuous map u → Q (n) w,± (u x) from S d to R d will be called the empirical conditional center-outward quantile function of Y given X = x. It defines the empirical center-outward regression quantile regions and contours

C (n) w,± (τ x) := Q (n) w,± τ S d x and C (n) w,± (τ x) := Q (n) w,± τ S d-1 x , τ ∈ (0, 1) (3.7)
which we are proposing as estimators of C ± (τ x) and C ± (τ x), respectively. The intersection τ ∈(0,1) C (n) w,± (τ x) yields the empirical conditional center-outward regression median region. The definition of empirical regression quantile tubes

T (n) w,± (τ ) := x, Q (n) ± (τ S d x) x ∈ R m , τ ∈ (0, 1)
naturally follows.

Remark 3.1. Note that the results of this section and the next one still hold for any continuous map with cyclically monotone graph satisfying

(u i , Q (n) w,± (u i x)) ∈ conv ({Y j : (u i , Y j ) ∈ supp(π * (x))}
) for all i = 1, . . . , N. The reason for choosing the "smallest" y in (3.6) is to have a "universal criterion."

4 Recall from [START_REF] Rockafellar | Convex Analysis[END_REF]

that a set S ⊂ R d × R d is cyclically monotone if any finite sub- set {(x k 1 , y k 1 ), . . . , (x kν , y kν )} ⊂ S, ν ∈ N satisfies ν-1 =1 y k , x k +1 -x k + y kν , x k 1 -x kν ≤ 0, where •, • stands for the scalar product in R d .

Consistency of empirical conditional center-outward quantiles, regression quantile regions, and regression quantile contours

The objective of this section is to justify the definitions of Section 3.1 by showing the consistency of the empirical quantile regions and contours defined in (3.7). The asymptotic behavior of these regions and contours, quite naturally, depends on the regularity of the conditional distributions involved. In fact, as discussed before, when Assumption (R) does not hold, the population conditional quantile maps are not necessarily defined for every u ∈ S d , but only for a set of U d -probability one. Consistency results can be obtained despite this a.s. definition, provided that population quantile maps are extended into set-valued maps. The following theorem shows5 , under mild assumptions, that any possible limit of

Q (n) w,± ( • |x) asymptotically belongs to the set Q ± ( • |X = x). Theorem 3.2. Let (X, Y), (X 1 , Y 1 ), . . . , (X n , Y n ) be pairs of i.i.d. random vectors with values in R m × R d
and let w (n) be a consistent sequence of weight functions. Then, for every u ∈ S d and > 0,

P Q (n) w,± (u |X) ∈ Q ± (u |X) + S d -→ 0 as n and N → ∞,
and, for every τ ∈ (0, 1),

P C (n) ± (τ X) ⊂ C ± (τ X) + S d → 0 and P C (n) ± (τ X) ⊂ C (n) ± (τ X) + S d → 0 as n and N → ∞.
Neater convergence results-avoiding the notion of set-valued maps-are obtained if it can be assumed that Assumption (R) holds, which implies that, for any X = x and any

u ∈ S d \ {0}, the set Q ± (u |X = x) is a singleton. Then, the map u → Q ± (u |X = x)
can be seen as continuous on S d \ {0}, see Theorem 25.5 in [START_REF] Rockafellar | Convex Analysis[END_REF], hence singlevalued on S d \{0} since the gradient of a convex function is single-valued at a point if and only if it is continuous at this point. We then can state the following theorem (see Appendix A.1 for the proof), the second part of which describes the convergence of the contours in terms of the Pompeiu-Hausdorff distance d ∞ . Recall that he Pompeiu-Hausdorff distance between two sets A and B in R d is defined as [START_REF] Rockafellar | Variational Analysis[END_REF]).

d ∞ (A, B) := inf{ν ≥ 0 : A ⊂ B + νS d and B ⊂ A + νS d } (see
Theorem 3.3. Let (X, Y), (X 1 , Y 1 ), . . . , (X n , Y n ) be pairs of i.i.d. random vectors with values in R m ×R d and let w (n) be a consistent sequence of weight functions. Suppose moreover that Assumption (R) holds. Then, for every compact K ⊂ S d \ {0}, as n and N → ∞,

sup u∈K |Q (n) w,± (u |X) -Q ± (u |X)| P -→ 0
and, for every τ ∈ (0, 1) and > 0,

P d ∞ C (n) ± (τ X), C ± (τ X) > → 0.
Under the assumptions of Theorem 3.3, consistency in Pompeiu-Hausdorff distance of the quantile contours holds in case the median is a single point6 -the continuity of quantile maps then extends to the whole open unit ball. This, however, is not necessarily the case for d > 3 (see [START_REF] Figalli | On the continuity of center-outward distribution and quantile functions[END_REF]), and Pompeiu-Hausdorff consistency may fail due to the fact that our empirical version is continuous over S d while Q ± (0|x) could be a set rather than a single point: convergence then holds along subsequences of Q

(n) w,± ( • |x) towards an element of Q ± (0|x)
. This, obviously, has an impact on convergence in terms of the Pompeiu-Hausdorff distance-although it does not affect the control over the asymptotic probability contents of quantile regions. More precisely, the following corollary holds (see Appendix A.1 for the proof).

Corollary 3.4. Under the conditions of Theorem 3.3, as n and N → ∞,

P Y ∈ C (n) ± (τ X) X P -→ τ for all τ ∈ (0, 1).
(3.8)

Remark 3.5. Under the weaker conditions of Theorem 3.2, in view of (A.12), still some asymptotic control of the probability content of the empirical regions can be derived. More precisely, letting N = N (n) be such that N (n) → ∞ as n → ∞, for all τ ∈ (0, 1) and every subsequence n k , there exists a further subsequence n kj such that

lim sup m P C (n k j ) ± (τ x) X = x ≤ τ, x-a.e. in R m .
The above results, as well as the proposed regularization, are valid for any consistent sequence of weight functions. This includes-along with adequate additional assumptionsmost of the classic choices of weight functions. Here are three examples.

(i) The kernel weight function, usually defined (see chapter 5 in [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF]) as

w (n) i (x; X (n) ) := K X i -x h n / n j=1 K X j -x h n , i = 1, . . . , n
where h n is the bandwidth, and K : R m → R the kernel. Sufficient conditions for w

w (n) i (x; X (n) ) := 1 k 1 Xi∈N (n) k (x) , i = 1, . . . , n.
Sufficient conditions for this k-nearest neighbors w (n) i to form a consistent sequence of weight functions are (see [START_REF] Stone | Consistent Nonparametric Regression[END_REF]) k → +∞ and k/n → 0.

(3.9)

If a k-nearest neighbors weight function w (n) is satisfying (3.9) and Assumption (R) holds, we thus have the convergence-in probability-of the conditional quantile map described in Theorem 3.3; without this assumption we still obtain the slightly weaker result of Theorem 3.2.

The k-nearest neighbors in (ii) were understood in the classical sense of the Euclidean distance in R m , which does not take into account the distribution P X of X. An alternative k-nearest neighbors weight function can be derived from a notion of nearness based on the ordering induced by empirical center-outward distribution functions. This new weight function is obtained as follows.

(iii) An alternative k-nearest neighbors weight function. Fixing x ∈ R m , first compute, as in [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF], the empirical center-outward distribution function associated with 1 n + 1

n j=1 δ Xj + 1 n + 1 δ x ∈ P(R m ).
That distribution function is the solution T * x of the minimization problem min

T ∈Γn+1 n k=0 |X k -T (X k )| 2
where X 0 = x, Γ n+1 is the set of all bijections T between {x, X 1 , . . . , X n } and a regular grid G (n+1) of S m , of the form described in Section 3.1, consisting of (n + 1) gridpoints denoted as G 0 , G 1 , . . . , G n , obtained via a factorization of n into a product of non-negative integers of the form n + 1 = n R n S + n 0 with n 0 < min(n R , n S ). That regular grid is created as the intersection between -the rays generated by an n S -tuple u 1 , . . . , u n S ∈ S m-1 of unit vectors such that n S -1 n S j=1 δ uj converges weakly, as n S → ∞, to the uniform over S m-1 and -the n R hyperspheres with center 0 and radii j/(n R + 1), j = 1, . . . , n R , along with n 0 copies of the origin whenever n 0 > 0. Based on this grid, we define the sequence of discrete uniform measures

U (n+1) d := 1 n + 1 n+1 j=1 δ Gj ∈ P(R m ) N ∈ N
This map is defined only at the (n + 1) points x, X 1 , . . . , X n , but, as in the previous section, it can be countinuously extended (see also [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF]) to the whole space R m -call F (n)

x;± : R m → S m this extension-with the properties that F (n)

x;± coincides, on {x, X 1 , . . . , X n }, with T *

x , is the gradient of a differentiable convex function with domain R m , and satisfies F (n)

x;± (R m ) ⊂ S m . We then define the set of k-nearest center-outward neighbors of x as

K (n) k (x) := {X j : F (n) ± (X j ) ∈ N k (F (n) x;± (x))}, (3.10)
where, for each a ∈ B m and k ∈ N, N k (a) denotes the set of k-nearest neighbors (in the sense of Euclidean distance) of a. Based on this, define the center-outward nearest neighbor weight function

w (n) j (x; X (n) ) := 1 k 1 Xj ∈K (n) k (x) , j = 1, . . . , k (3.11)
and proceed as in Section 3.1 with the estimation (3.5) of the conditional quantile functions.

The next result shows that, for a suitable choice of k = k(n), center-outward nearest neighbors weight functions form a consistent sequence of weights (see the appendix for a proof).

Lemma 3.6. If k = k(n) is such that k(n) → ∞ and k(n)/n → 0 as n → ∞, the sequence of weight functions defined in (3.11) is consistent in the sense of (3.4).
This means, in particular, that Theorem 3.2 applies when the weight function (3.11) is used under the assimptions of Lemma 3.6, and that the resulting estimators are consistent.

Theorems 3.2 and 3.3 provide weak (in probability) consistency results under minimal assumptions. For sequences of weights satisfying, as n → ∞,

n j=1 w (n) j (X; X (n) )Y j -→ E[Y |X] a.s.
(3.12) (strongly consistent sequences), the conclusion in Theorem 3.2 can be upgraded to strong (almost sure) consistency. For the particular case of k-nearest neighbors, (3.12) and strong consistency hold if (3.9) is replaced with

k/log(n) → ∞ and k/n → 0, (3.13)
see [START_REF] Devroye | On the strong universal consistency of nearest neighbor regression function estimates[END_REF] and [START_REF] Devroye | Necessary and sufficient conditions for the pointwise convergence of nearest neighbor regression function estimates[END_REF].

Numerical results

This section is devoted to a numerical assessment of the performance of the estimation procedures described in Section 3. We first analyze (Section 4.1) some artificial datasetsincluding the motivating example of Hallin et al. (2015)-then turn (Section 4.2) to real-data cases. These examples showcase three important features of our estimators: their ability to capture heteroskedasticity, to deal (non-parametrically) with highly nonlinear regression, and to adapt to non-convex noise.

Simulated examples

Parabolic trend and periodic heteroskedasticity; spherical conditional densities.

We start with analyzing the motivating example given in Hallin et al. (2015). The model

(with m = 1, d = 2) is Y = Y 1 Y 2 = X X 2 + 1 + 3 2 sin π 2 X 2 e, X ∼ U [-2,2
] and e ∼ N (0, Id), (4.1)

with X and e mutually independent. In this case, the population conditional (on X = x) quantile contours are circles with radii depending on x and can be computed exactly; trend is parabolic and heteroskedasticity periodic. Figure 3 illustrates the convergence of our estimated contours to the population counterparts. Compared to Figure 1 in Hallin et al. (2015), our method produces less smooth contours, at least for smaller sample sizes. On the other hand, our method is able to capture non-convex contour shapes-something the method in Hallin et al. (2015) cannot, see Section 4.1.2. We also underline that, from a computational point of view, our method is able to handle rather large datasets (in contrast, the R packaged modQR cannot handle sample sizes over 10, 000, as explained in the documentation).

Model (4.1), as pointed out in Hallin et al. (2015), allows for testing the capacity of a method to estimate the trend while catching potential heteroskedasticity. A comparison with Figure 1 in Hallin et al. (2015) shows that both methods estimate the parabolic trend quite well, but that our method performs much better at capturing heteroskedasticity. Estimations are based on the k-nearest neighbors weights (ii) of Section 3, with N = k = 401, 1, 000, and 6, 401 and n = 3, 601, 10, 000, and 128, 020, respectively. Note that, in this univariate covariate case, the weight choices (ii) and (iii) coincide.

The performance of the Gaussian kernel (i) (with various bandwidth choices) and classical k-nearest neighbors (ii) (with various choices of k) weight functions are investigated in Figures 4 and5 x ∈ {-2, -1.6, -1.1, -0.7, -0.2, 0.2, 0.7, 1.1, 1.6, 2}.

For this sample size, the Gaussian kernel weights-due to the fact that they better exploit the information available on the x's-yield better results than the k-nearest neighbors ones. But the Gaussian kernel has a drawback for large datasets; the optimization problem (3.5) requires the whole dataset and cannot be efficiently computed. For instance, the Gaussian kernel counterpart of Figure 3 (where n = 128, 020) cannot be computed on a standard desktop computer: large-sample datasets should be handled either with nearest neighbors or compactly supported kernel weights. On the other hand, the bandwith h n and the neighborhood size k apparently have little impact on he result.

Parabolic trend and periodic heteroskedasticity; banana-shaped conditional densities.

We now consider a model in which the trend and heteroskedasticity are the same as in Model (4.1), but the quantile contours are non-convex (conditional densities are bananashaped):

Y = Y 1 Y 2 = X X 2 +   1 + 3 2 sin π 2 X 2 1.15 e 1 1 + 3 2 sin π 2 X 2 e 2 1.15 + 0.5(e 2 1 + 1.21)   , (4.3) with X ∼ U [-2,2
] and e = (e 1 , e 2 ) ∼ N (0, Id), X and e mutually independent. The dataset shown in Figure 1 was generated from that model, with X = 0 and justifies the termonology "banana-shaped." Population conditional quantile contours here cannot be computed analytically. Figure 6 shows their estimations (same method as in Section 4.1.1, with a k-nearest neighbors weight function, k = 14, 401) for orders τ = 0.2 (black), 0.4 (green), and 0.8 (yellow), along with the conditional medians (red). The sample size n = 576, 040 is very large, so that, in view of consistency results, these estimations can be considered as close numerical approximations of their theoretical counterparts. This example illustrates further the ability of our method to handle vary large sample sizes. The Hallin et al. (2015) approach produces contours that are convex by construction-hence cannot capture the "banana shape" of the conditional densities. The same comparative analysis of weight functions is performed as for Model (4.1). Figures 7 and8 show the results for Gaussian kernel weights (various bandwiths) and knearest neighbors weights (various values of k), respectively. The sample size is n = 3601. In all the examples, for the ease of computation, N = k. Note that, for the k-nearest neighbor weights, this choice creates a one-to-one (between the sample and the grid) transport map.

Still for Model (4.3), Figure 9 shows, for sample size n = 3601 and various choices of the bandwidth h and the neighborhood size k, the behavior of the empirical conditional center- outward contours in X = 0-and compares them with those of Figure 6 (considered as the population contours). The empirical conditional quantiles are computed for τ = 0.2, 0.4, and 0.8 with Gaussian kernel weights (bandwidths h = 0.1, 0.2, and 0.3) and the k-nearest neighbors weights (k = 226, 485, and 901). The influence of the choice of h and k is clearly seen here: the bigger h (the bigger k), the smoother the estimation of the shape of the contours but also, unfortunately, the worse the estimation of their location.

Some real-data examples

The CalCOFI oceanographic dataset: depth, temperature, and salinity in the oceans

The dataset "CalCOFI Over 60 Years of Oceanographic Data," available at https://www. kaggle.com/sohier/calcofi, contains the longest (1949-present) and most complete (more than 50,000 sampling stations; sample size n = 814, 247) time series of oceanographic and larval fish data worldwide. Data collected at depths down to 500 m include temperature, salinity, oxygen, phosphate, silicate, nitrate and nitrite, chlorophyll, transmissometer, PAR, C14 primary productivity, phytoplankton biodiversity, zooplankton biomass, and zooplankton biodiversity. We are focusing here on the influence of X = 'depth' (in meters) on the pair Y = ('temperature,' 'salinity') (in degrees and grams of salt per kilogram of water, respectively).

Figure 10 shows the corresponding 3D observations and the estimated conditional centeroutward quantile contours obtained from the same method as in Section 4.1 (nearest neighbors weight function with k = 6, 401); Figure 11 shows the projections of the same contours on the ('depth,' 'salinity') and ('depth,' 'temperature') axes, respectively. Inspection of these figures reveals a nonlinear center-outward conditional median; heteroskedasticity also appears as the area of the conditional quantile regions clearly decreases as a function of depth, while a positive dependence between temperatures and salinity, which is present at the surface, gradually disappears as depth increases. The projection plots of Figure 11 also provide clearer views on marginal dependencies. For example, the decrease of temperature as a function of depth is monotone and almost linear, while the dependence on depth of salinity is more complex, high at shallow depths, lower at medium depths, and higher again at greater depths. However, these marginal analyses, to some degree, are hiding the heteroskedasticity effects (in particular, the dependence on depth of the relation between salinity and temperature) which are clearly visible in Figure 10. Since the dataset is quite large, we used a nearest neighbors weight function, see the comments about the empirical performance of different weights in Section 4.1.1. 4.2.2. The Female ANSUR 2 dataset: stature, foot length and tibial height of female US Army personnel .

Our second real-data example involves a smaller sample size n. The Female Anthropometric Survey of US Army Personnel (Female ANSUR 2 or Female ANSUR II) featured in this section consists in ninety-three direct measures and 41 derived ones, as well as three-dimensional head, foot, and whole-body scans of n = 1, 986 women of the United States' army. These measurements were collected between October 4, 2010 and April 5, 2012, in May 2014 and in May 2015; they are available online at https://www.openlab.psu.edu/ansur2/.

We want to analyze the influence of the covariate X = 'stature' (in centimeters) on the variable of interest Y = ('foot length,' 'tibial height') (both in centimeters). Figure 12 provides a 3D view of the center-outward quantile contours/tubes (levels τ = 0.2, 0.4, 0.7), along with the center-outward regression median (red); Figure 13 shows the projections on the axes of the same contours. Inspecting these two figures reveals the absence of heteroskedasticity, the spherical shape of conditional distributions, and a roughly linear regression. Since the size of the model is not too large, a Gaussian kernel is convenient. The bandwidth was chosen as h = 15, which, up to scale changes, corresponds to the choice h = 0.2 in Figure 7.

Some concluding remarks

Relation to the recent literature on numerical optimal transportation

The estimation of transport maps beyond the sample points currently is a hot topic, and a fastly developing strand of literature is proposing such estimators. However, the objec- tive of most authors is to reach near-optimal convergence rates, for which they typically impose fairly strong assumptions. Some estimators (Hütter and Rigollet (2021), Manole et al. ( 2021)) are computationally quite heavy and sometimes numerically almost infeasible; others [START_REF] Pooladian | Entropic estimation of optimal transport maps[END_REF]) use a regularized version [START_REF] Cuturi | Sinkhorn distances: lightspeed computation of optimal transport[END_REF]) of the optimal transport problem to provide consistent near-optimal estimators, which needs stringent assumptions on the shape of the underlying distributions-such as being compactly supported, with densities bounded away from 0 and ∞ over their convex supports. This is redhibitory in our case, since the density of U d (a choice which plays an essential role in the interpretation of transports as quantile functions) is unbounded at 0. Other solutions (Makkuva et al. (2020) and González-Sanz et al. ( 2022)) are based on deep learning and neural network methods; they achieve excellent empirical performance, but the lack of theoretical results for the first one, the Lipschitz constraint on transport maps for the second, preclude their use in this quantile regression context.

We also could estimate the conditional quantiles through the optimal map from the uniform reference measure to P w(x) n

. From a numerical point of view, however, this would lead to the computation of a semi-discrete optimal transportation plan, which has complexity O(n d/2 ), hence is unfeasible even for moderate d. While the computational complexity of our procedure does not depend on the dimension, its statistical performance does (see [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF]) and, in that sense, we do not escape the curse of dimensionality-up to the case where P is finitely supported, see del Barrio, González-Sanz and Loubes (2021). Despite the fact that the literature on the computation of such maps is growing quite fastly (see Lévy, Mohayaee and von Hausegger (2020); Gallouët and Mérigot 2012)), the exisiting methods are restricted to dimension two, sometimes three. A further issue is that the solution of the semi-discrete problem does not produce quantile contours but creates a Voronoi tessellation of S d , each piece of which is mapped to a single sample point.

Conclusions and perspectives for further developments

Building on the concepts of center-outward quantiles recently developed in [START_REF] Shi | Center-outward sign-and rank-based quadrant, Spearman, and Kendall tests for multivariate independence[END_REF], we are proposing here a fully nonparametric solution to the problem of nonparametric multiple-output quantile regression. Contrary to earlier attempts, our solution is enjoying the quintessential property that the (conditional) probability content of its quantile regions is under control irrespective of the underlying distribution. This is only a first step into the multifarious applications of multiple-output quantile regression, though. Due to the minimality of the assumptions it requires, a completely agnostic nonparametric approach indeed is attractive, but also comes at a cost: linear or polynomial quantile regression remain justified whenever a priori knowledge of the analytical form of the regression is available and should be taken advantage of. A center-outward version of the results of Carlier, Chernozhukov and Galichon (2016), thus, is highly desirable. Single-output quantile regression has been considered in a variety of contexts: survival analysis, longitudinal data, instrumental variable regression, directional, functional, and high-dimensional data, . . . Quantile regression versions of time-series models such as the quantile autoregressive model also have been investigated (Koenker and Xiao, 2006). All these applications call for multiple-output extensions with important real-life consequences; they should and can be based on the concept of center-outward quantile, regions, and contours and are the subject of our ongoing research.

to (x, y)-exists and coincides with T and -the inner limit lim inf n T n -which is the set of (x, y) ∈ R d × R d for which there exists a sequence {(x n , y n )}, with (x n , y n ) ∈ T n , which converges to (x, y)-exists and coincides with T .

We start with some necessary properties of the population and empirical conditional quantiles and some auxiliary results. The following lemma states that the estimated conditional probability converges weakly in probability to its population counterpart. Recall from Theorem 1.12.4 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF] that weak convergence can be measured in terms of the bounded Lipschitz norm, which is defined for µ, ν ∈ P(R d ) as

d BL (µ, ν) := sup f ∈F BL (R d ) |E Z∼µ (f (Z)) -E W∼ν (f (W))| ,
where the class F BL is the class of functions f :

R d → R such that |f (z 1 ) -f (z 2 )| ≤ |z 1 -z 2 | and |f (z 1 )| ≤ 1, for all z 1 , z 2 ∈ R d .
Lemma A.1. For any > 0, under the assumptions of Theorem 3.2, we have

P d BL (P w(X) n , P Y |X ) > -→ 0.
See Section A.2 for the proof.

Let π (n) (x) be a solution of (3.5) and

π * (x) := (Id × Q ± ( • |x)) #U d .
Note that, for each value of x, (i) there exists a sequence of differentiable convex functions

ψ (n) • x : R d → R∪{+∞}, n ∈ N, such that Q (n)
w,± u j x = ∇ψ (n) u j x for j = 1, . . . , k, and n ∈ N;

(ii) there exists a convex function ψ • X = x) : R d → R ∪ {+∞} such that

Q ± u X = x = ∇ψ u X = x for U d -a.e. u ∈ S d .
Using obvious notation, it holds that π * (x) and π (n) (x), for n ∈ N, have cyclically monotone supports. Moreover, as a consequence of Corollary 3.1, we have

supp(π (n) (x)) ⊂ ∂ ψ(n) • x and supp(π * (x)) ⊂ ∂ψ • X = x , (A.1)
possibly for some other sequence ψ(n)

• x : R d → R ∪ {+∞} of convex functions such that Q (n) w,± u j x ∈ ∂ ψ(n) u j x for j = 1, . . . , k, and n ∈ N. (A.2)
The following result then follows from Lemma 9 and Corollary 14 in McCann (1995).

Lemma A.2. Let µ, ν ∈ P(R d ) be such that µ d is supported on a convex set. Let {µ n } n∈N and {ν n } n∈N ⊂ P(R d ) converge weakly as n → ∞ to µ and ν, respectively. Suppose, moreover, that there exists a sequence of probability measures {π n } n∈N ⊂ P(R d × R d ) with marginals µ n and ν n such that, for some sequence of convex functions {φ n } n∈N , it holds that, for all n ∈ N, supp(π n ) ⊂ ∂φ n . Then, (i) {π n } n converges weakly as n → ∞ to π * = (Id × ∇φ) #µ, where ∇φ is the gradient of a convex function φ pushing µ forward to ν, and (ii) there exists a sequence {a n } n∈N ⊂ R such that, µ-a.e., φ n + a n → φ and ∂φ n → ∂φ graphically as n → ∞.

See Section A.2 for the proof.

Proof of Theorem 3.2. Suppose that there exist u 0 ∈ S d and 0 > 0 such that lim sup

n→∞ P Q (n) w,± (u 0 |x) ∈ Q ± (u 0 |X) + 0 S d = δ > 0.
We can find a subsequence n k such that

lim k→∞ P Q (n k ) w,± (u 0 |w(X; X [n k ] )) ∈ Q ± (u 0 |X) + 0 S d = δ > 0. (A.3)
Note that the space of probability distributions on R d , endowed with the bounded Lipschitz metric, is separable and complete, see Theorem 1.12.4 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF]. Therefore, the convergence described in Lemma A.1 implies that, for the subsequence n k , there exists a further sub sequence n ki such that the event

Ω 0 =   sup f ∈F BL (R d ) n k i j=1 w j (X; X (n k i ) )f (Y j ) -E(f (Y)|X) -→ 0   (A.4)
has probability one. Let x = X(ω) with ω ∈ Ω 0 . Then, by Lemma A.2, we have, as i → ∞,

∂ ψ(n k i ) • x -→ Q ± ( • |X = x) graphically.
This implies, by Theorem 8.3 in [START_REF] Bagh | Convergence of set-valued mappings: equi-outer semicontinuity[END_REF], that there exists i 0 ∈ N such that, for all i > i 0

, ∂ ψ(n k i ) u 0 x ⊂ Q ± (u 0 |X = x) + 0 S d . (A.5)
From (A.2) and the fact that (A.5) holds with probability one, we deduce that

P   i0∈N i≥i0 Q (n k i ) ± (u 0 |X) ∈ Q w,± (u 0 |X) + 0 S d   = 0, (A.6) which implies that lim sup i→∞ P Q (n k i ) ± (u 0 |X) ∈ Q w,± (u 0 |X) + 0 S d ≤ 0.
This contradicts (A.3). The rest of the proof follows from compactness arguments and a refined use of Theorem 8.3. in [START_REF] Bagh | Convergence of set-valued mappings: equi-outer semicontinuity[END_REF]. Suppose that, for some q 0 ∈ (0, 1) and 0 > 0, and along a subsequence n k , we have

lim k→∞ P C (n k ) ± (q 0 X) ⊂ C ± (q 0 X) + 0 S d = δ > 0.
(A.7)

Theorem 8.3. in [START_REF] Bagh | Convergence of set-valued mappings: equi-outer semicontinuity[END_REF] and Lemma A.2 jointly imply the existence of a further subsequence n ki , such that, for every ū ∈ S d , there exists I ū ∈ N and λ ū > 0 satisfying

Q (n k i ) ± (u |x) ∈ Q ± (ū |X = x) + 0 S d for all u ∈ ū + λ ūS d and all i > I ū.
Note that q S d ⊂ ū: |ū |≤q ū + λ ūS d . Since q S d is compact, there exists a finite covering q S d ⊂ K k=1 ūk + δ ūk S d . Set I 0 = max(I ū1 , . . . , I ūK ): then for all i > I ū, we have

C (n k i ) ± (q 0 x) ⊂ C ± (q 0 x) + 0 S d ,
which holds with probability one and, using the same argument as for (A.5), contradicts (A.7).

The same reasoning holds also for the contours. w,± (u |x) are a.s. singletons, we have, as n and N tend to infinity,

P Q (n) w,± (u |x) -Q ± (u |X) > → 0.
Let K be a compact subset of S d \ {0}. In order to establish uniform convergence in K, suppose that the contrary holds. Since the space of continuous functions from K to R d (endowed with the topology of uniform convergence) is complete and separable, there exists a subsequence n k such that, for some δ > 0, the probability of

Ω = sup u∈K Q (n k ) w,± (u |X) -Q ± (u |X) → δ , (A.8)
is one. Set ω ∈ Ω and consider x = X(ω). There exists a sequence {u n k } ⊂ K such that

Q (n k ) w,± (u n k |x) -Q ± (u n k |X = x) → δ, for a possibly different δ > 0. Since the sequence {u n k } k is in K, it admits at least one point of accumulation ū ∈ K. Hence, lim inf k→∞ ( Q (n k ) w,± (u n k |x) -Q ± (ū |X = x) + Q ± (ū |X = x) -Q ± (u n k |X = x) ) ≥ δ,
where the second term tends to 0 by the continuity of u

→ Q ± (u |X = x). This implies that lim inf k→∞ Q (n k ) w,± (u n k |x) -Q ± (ū |X = x) ≥ δ. (A.9) But Lemma A.2 entails the Graphical convergence to Q ± ( • |X = x) of Q (n k )
w,± ( • |x), which contradicts (A.9). The desired uniformity over K follows.

Finally, the convergence of the contours is a consequence of the previous result on the regions.

Proof of Corollary 3.4. To prove (3.8), let us show that, for any subsequence n k , there exists a further subsequence converging a.s. To avoid repetitions, assume that N

= N (n) is such that N → ∞ as n → ∞.
Leaving aside the singular points, let us consider the empirical and population quantile rings

C (n) ± ( , τ x) := C (n) ± (τ x) \ C (n) ± ( x) and C ± ( , τ x) := C ± (τ x) \ C ± ( x), respectively. Theorem 3.3 yields, for all 0 < < τ , as n → ∞, P Y ∈ C (n) ± ( , τ X) and Y ∈ C ± ( , τ X) X P -→ 0.
(A.10) Let { j } j∈N be a monotone sequence tending to 0. For j = 1 there exists a further subsequence n 1 k , say, and a subset Ω 1 of Ω such that P(Ω 1 ) = 1 such that, for every x = X(ω)

with ω ∈ Ω 1 , P Y ∈ C (n 1 k ) ± ( 1 , τ x) and Y ∈ C ± ( 1 , τ x) X = x -→0. (A.11) By definition, P C ± ( 1 , τ x) X = x = τ -1 ; therefore, in view of (A.11), P C (n 1 k ) ± ( 1 , τ x) X = x -→ τ -1 .
Repeating the argumentat for j = 2, there exists a set Ω 2 ⊂ Ω, with P(Ω 2 ) = 1, and a subsequence (call it n 2 k ) of n 1 k such that, for every ω ∈ Ω 2 and x = X(ω),

P C (n 2 k ) ± ( 2 , τ x) X = x -→ τ -2 .
This argument can be repeated for each j ∈ N and the set Ω 0 = j∈N Ω j has probability one. Set ω ∈ Ω 0 ⊂ Ω and x = X(ω): then there exists k 1 such that

P C (n k 1 ) ± ( 1 , τ x) X = x -τ < 2 1 .
Analogously, for each m, there exist k j such that

P C (n k j ) ± ( j , τ x) X = x -τ < 2 j .
Therefore, we obtain the limit

P C (n k j ) ± ( j , τ x) X = x -→ τ
and, noticing that

P C (n k j ) ± ( j , τ x) X = x ≤ P C (n k j ) ± (τ x) X = x , also the asymptotic upper bound lim inf P C (n k j ) ± (τ x) X = x ≥ τ.
Now, Theorem 3.2 implies, for every subsequence (for which we keep the notation n) of n and j ∈ N, there existence of some n j ∈ N such that

P C (n) ± (τ X) ⊂ C ± (τ X) + 1 2 j S d ≤ 1 2 j for all n ≥ n j . The sum ∞ j=1 P C (nj ) ± (τ X) ⊂ C ± (τ X) + 1 2 j S d
thus is finite, and the Borel-Cantelli lemma yields

P   J∈N j≥J C (nj ) ± (τ X) ⊂ C ± (τ X) + 1 2 j S d   = 0,
or, equivalently, P(Ω * ) = 1, where

Ω * :=   J∈N j≥J C (nj ) ± (τ X) ⊂ C ± (τ X) + 1 2 j S d   .
Setting ω ∈ Ω * ⊂ Ω and x = X(ω), there exists J ∈ N such that

P C (nj ) ± (τ x) X = x ≤ P C ± (τ x) + 1 2 j S d X = x ,
for all j ≥ J. Since Proof of (a). Let X be a random variable independent of X (n) , with the same distribution as X i . Set u 0 = F (n) ± (X) and u j = F

(n)

± (X j ), j = 1, . . . , n. Defining

K k,(i) n (X) := {X j : i = j, F n ± (X j ) ∈ N k (F (n) ± (X i ))},
note that 1 Xi∈K k n (X) f (X i ) and 1 X∈K k,(i) n (Xi) f (X) have the same distribution. It follows that [START_REF] Stone | Consistent Nonparametric Regression[END_REF]. As for (e), it follows from the fact that k → ∞.

E   1 k n j=1 1 Xj ∈K k n (X) f (X j )   = 1 k n j=1 E 1 Xj ∈K k n (X) f (X j ) = 1 k n j=1 E 1 X∈K k,(j) n (Xj ) f (X) , which implies that E   1 k n j=1 1 Xj ∈K k n (X) f (X j )   ≤ E   f (X) 1 k n j=1 1 X∈K k,(
Proof of (c). The following lemma (see Section A.2 for a proof) is a corollary of Lemma 6.1 in [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF]. |G i -G 0 | -→ 0, a.s.

Since Q ± is a singleton with probability one, let us assume, without loss of generality, that Q ± (G 0 ) and Q This, in view of Proposition 5.33 in [START_REF] Rockafellar | Variational Analysis[END_REF], implies that sup Xi∈K k n (X)

|X i -X | -→ 0, a.s. (A.14) We are ready now to prove (c). Set a > 0. Since there are k elements in K k n (X), we have

E   1 k n j=1 1 Xj ∈K k n (X) 1 |Xj -X|>a   = E   1 k Xj ∈K k n (X) 1 |Xj -X |>a   = P sup Xi∈K k n (X)
|X j -X | > a which, owing to (A.14), tends to 0. The desired result follows as a direct consequence of Markov's inequality.

A.2. Proofs of Lemmas A.1, A.2, and A.3

Proof of Lemma A.1. Let δ > 0 and > 0 be arbitrary. Denote by K ⊂ R d a compact set such that P (Y ∈ K) ≥ 1 -δ /18. Suppose that 0 ∈ K and define F K as the class of 1-Lipschitz functions f supported on K such that f (0) = 0. Such a class, by the Arzelà-Ascoli theorem, is relatively compact for C(K) and the uniform convergence. Then, there exists a sequence f 1 , . . . , f N , such that sup

f ∈F K inf k=1•••N ||f -f k || ∞ ≤ /8.
Therefore, for every f ∈ F K , we have where the last inequality follows from the fact that the weight function is consistent, Note that every f ∈ F BL can be approximated by f 1 K . This yields

A n ≤ n j=1 w j (X; X (n) )1 R d \K (Y j ) + E(1 R d \K (Y)|X) ≤ n j=1 w j (X; X (n) )1 R d \K (Y j ) -E(1 R d \K (Y)|X) + 2E(1 R d \K (Y)|X).
Note that the bound does not depend on the function f . Consequently, Taking expectations on both sides we obtain

E   sup f ∈F BL n j=1 w j (X; X (n) )f (Y j )1 R d \K (Y j ) -E(f (Y)1 R d \K (Y)|X)   ≤ E n j=1 w j (X; X (n) )1 R d \K (Y j ) -E(1 R d \K (Y)|X) + 2P (Y / ∈ K) .
Since the weights are consistent, there exists n 1 such that

E n j=1 w j (X; X (n) )1 R d \K (Y j ) -E(1 R d \K (Y)|X) ≤ δ

Fig 1 .

 1 Fig 1. Quantile contours C ± (τ ) (left panel) and regions C ± (τ ) (right panel) for the banana-shaped Gaussian mixture of Section 4.1.2 and the quantile orders τ = 0.1, 0.2, . . . , 0.9.

  ; Hallin, La Vecchia and Liu (2020); Shi et al. (2021b); Deb and Sen (2021); Shi et al. (2021a); Hallin, Hlubinka and Hudecová (2022); Hallin, La Vecchia and Liu (2022), among others).

Fig 2 .

 2 Fig 2. Quantile regression (two-dimensional variable of interest Y; univariate regressor X), showing the conditional center-outward medians (red) and the conditional quantile contours (x, C ± (τ | x))of order τ = 0.2 (black), τ = 0.4 (green), and τ = 0.8 (yellow).

  consistent sequence of weight functions are (a) h n → 0, (b) c 1 min 1 |x |≤r, H(|x |) ≤ K(x) ≤ c 2 H(|x |), where c 1 , c 2 , r are positive constants, and H : [0, ∞) → R is bounded, decreasing, and such that H(t)t m → 0 as t → ∞, and (c) lim n→∞ n α h m n /log(n) = ∞ for any α ∈ (0, 1); when the kernel K is compactly supported, the assumptions are are much simpler, and we only need (a) and lim n→∞ h m n n = ∞, see Theorem 5.1 in[START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF]. The particular case K(x) = e -|x| 2 is known as the Gaussian kernel. (ii) The (classical) k-nearest neighbors weight function: the k-nearest neighborhood of x ∈ R m is obtained (Chapter 6 in Györfi et al. (2002)) by ordering {X 1 , . . . , X n } according to increasing values of |X j -x |. Denoting by {X (0,x) , . . . , X (n,x) } the reordered sequence, the set of k-nearest neighbors of x is N (n) k (x) := {X (j,x) : j ≤ k} and the k-nearest neighbors weight function is defined as

  , still for Model 4.1, with sample size n = 3, 601. More precisely, estimation in Figure 4 is based on the weight function w h n = 0.05, 0.1, 0.2, and 0.3; the discretization of the spherical uniform is based on a grid of size N = n. Estimation in Figure 5 is based on k-nearest neighbors weights with k = 101, 256, and 625. Here we chose N = k, which yields a one-to-one solution in (3.5). In both figures the conditional contours (of order τ = 0.2 (black) and τ = 0.4 (green)) and the estimated conditional center-outward medians (red) are shown for

Fig 3 .

 3 Fig 3. Estimated (sample sizes 3, 601 in the upper left panel, 10, 000 in the upper right panel, 128, 020 in the lower left panel) and population (lower right panel) quantile contours of order τ = 0.2 (black) and 0.4 (green) for Model (4.1); the (estimated) conditional center-outward medians are shown in red. Estimations are based on the k-nearest neighbors weights (ii) with N = k = 401, 1, 000, and 6, 401 and n = 3, 601, 10, 000, and 128, 020, respectively.

Fig 4 .

 4 Fig 4. Estimated quantile contours (Model (4.1)) based on Gaussian kernel weight functions for different choices of the bandwith h. The sample size is n = 3, 601; bandwidth values are h = 0.05 (upper left panel), h = 0.1 (upper right panel), h = 0.2 (lower left panel) and h = 0.3 (lower right panel). The estimated contour orders are τ = 0.2 (black) and 0.4 (green); the estimated conditional center-outward medians are shown in red.

Fig 5 .

 5 Fig 5. Estimated quantile contours (Model (4.1)) based on k-nearest neighbors weight functions for different choices of k. The sample size is n = 3601; k = 101 (upper left panel), k = 256 (upper right panel), k = 401 (lower left panel) and k = 625 (lower left panel). The estimated contour orders are τ = 0.2 (black) and 0.4 (green); the estimated conditional center-outward medians are shown in red.

Fig 6 .

 6 Fig 6. A numerical approximation (estimation based on a simulated sample of size n = 576, 040) of the quantile contours of order τ = 0.2 (black), 0.4 (green), and 0.8 (yellow) for Model (4.3); the conditional center-outward medians are shown in red.

Fig 7 .

 7 Fig 7. Performance, in Model (4.3), of Gaussian kernel weight functions-based estimation for various choices of the bandwidth. The sample size is n = 3, 601, and the bandwidths are h = 0.05 (upper left panel), h = 0.1 (upper right panel), h = 0.2 (lower left panel), and h = 0.3 (lower left panel). The empirical contour levels are τ = 0.2 (black) , 0.4 (green) and 0.8 (yellow); the estimated conditional center-outward medians are shown in red.

Fig 8 .

 8 Fig 8. Performance, in Model (4.3), of the k-nearest neighbors weight functions-based estimation for various choices of k. The sample size is n = 3601; k = 101 (upper left panel), k = 256 (upper right panel), k = 401 (lower left panel), and k = 625 (lower left panel). The empirical contour levels are τ = 0.2 (black), 0.4 (green), and 0.8 (yellow); the estimated conditional center-outward medians are shown in red.

Fig 9 .

 9 Fig 9. Comparison between the empirical conditional contours C (n) ± (τ 0) (levels τ = 0.2, 0.4, 0.8; sample size n = 3, 601) in Model (4.3) based on Gaussian kernel weight functions (bandwidths h = 0.1 (green), h = 0.2 (yellow), and h = 0.3 (blue)) (left panel) and those in Figure 6, based on k-nearest neighbors weight functions (k = 226 (green), k = 485 (yellow), and k = 901 (blue); sample size n = 576, 040) (right panel). The center-outward quantiles of Figure 6 (to be considered as population values) are shown in red.

Fig 10 .

 10 Fig 10. CalCOFI dataset. Left panel: the original dataset ('depth,' 'temperature,' 'salinity') for 'depth' ≤ 200 (sample size for 'depth' ≤ 200, dropping empty values, is n = 505, 829). Right panel: the empirical conditional center-outward quantile contours of orders τ = 0.2 (black), τ = 0.4 (green), and τ = 0.8 (yellow) and the empirical conditional center-outward median (red) for the multiple-output regression of (Y 1 , Y 2 ) = ('temperature,' 'salinity') with respect to X = 'depth.' Estimation based on a k-nearest neighbors weight function with k = 6, 401.

Fig 11 .

 11 Fig 11. CalCOFI dataset. Left panel: projection of the empirical center-outward quantile regions shown in Figure 10 (orders τ = 0.2, 0.4, 0.8) and median on the axes ('depth,' 'salinity'), for the multiple-output regression of (Y 1 , Y 2 ) = ('temperature,' 'salinity') with respect to X = 'depth.' Right panel: projection of the same on the axes ('depth,' 'temperature'); see Figure 10 for the color code.

Fig 12 .

 12 Fig 12. ANSUR 2 dataset (sample size n = 1, 986). Left panel: the original dataset of X = 'stature' and (Y 1 , Y 2 ) = ('foot length,' 'tibial height'). Right panel: the empirical conditional center-outward quantile contours orders τ = 0.2 (black), τ = 0.4 (green), τ = 0.7 (yellow) and the empirical conditional centeroutward median (red) for the multiple-output regression of (Y 1 , Y 2 ) with respect to X; estimation based on a Gaussian kernel weight function with bandwidth h = 15.

Fig 13 .

 13 Fig 13. ANSUR 2 dataset (sample size n = 1, 986). Left panel: projection of the empirical center-outward quantile regions (orders τ = 0.2, 0.4, 0.8) and median on the axes ('stature,' 'foot length'), for the multipleoutput regression of (Y 1 , Y 2 ) = ('foot length,' 'tibial height') with respect to X = 'stature.' Right panel: projection of the same on the axes ('stature,' 'tibial height'); see Figure 12 for the color code.

(

  2018); Meyron (2019); de Goes et al. (

Proof of Theorem 3. 3 .

 3 Under Assumption (R), it follows from del Barrio, González-Sanz and Hallin (2020) that, for ω ∈ Ω 0 ⊆ Ω where P(Ω 0 ) = 1, the center-outward quantile function Q ± (u |X = x := X(ω)) is a singleton for U d -almost all values of u ∈ S d . Therefore, we adopt here the slight abuse of notation commented before Theorem 3.3 and write Q ± (u |X = x) = {Q ± (u |X for x)}. Set u ∈ S d \ {0} and note that, since Q ± (u |X) and Q (n)

1

  To prove this lemma, we show that the following five conditions of Stone's theorem (Theorem 1 in[START_REF] Stone | Consistent Nonparametric Regression[END_REF]) are satisfied (convergence for n → ∞):(a) there exists C ≥ 1 such that, for any non-negative measurable function f , Xj ∈K k n (X) 1 |Xj -X|>a -→ 0, in probability, for all a > 0;

  k (Gj ) , we can apply Corollary 6.1. in[START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF] and conclude that there exists λ d ∈ R such that n j=1 1 X∈K k,(j) n (Xj ) ≤ kλ d . This and (A.13) complete the proof. Proof of (b), (d), and (e). Conditions (b) and (d) are direct consequences of the properties of the weight function, see Corollary 1 in

Lemma A. 3 .

 3 Let k/n → 0 as n → ∞. Then, as n → ∞ sup Gi∈N k (G0)

  Q ± graphically as n → ∞, (ii) for every X i ∈ K k n (X), there exists some G i ∈ N k (G 0 ) such that Q (n) ± (G i ) = {X i }, and (iii) Lemma A.3 yields G i → G 0 as n → ∞.

  X; X (n) )f (Y j ) -E(f (Y)|X) ≤ sup k=1•••N n j=1 w j (X; X (n) )f k (Y j ) -E(f k (Y)|X) + 4 .In consequence, there exists n 0 such that, for n ≥ n 0 , we have X;X (n) )f (Y j ) -E(f (Y)|X) X; X (n) )f (Y j ) -E(f (Y)|X) X; X (n) )f (Y j ) -E(f (Y)|X)

  X; X (n) )f (Y j ) -E(f (Y)|X) X; X (n) )f (Y j ) -E(f (Y)|X) + n j=1 w j (X; X (n) )f (Y j )1 R d \K (Y j ) -E(f (Y)1 R d \K (Y)|X) .(A.16) Inequality (A.15) provides an upper bound for the first term. The second one, denoted as A n , is bounded byn j=1 w j (X; X (n) )f (Y j )1 R d \K (Y j ) + E(f (Y)1 R d \K (Y)|X) .Since the weights are positive and sup x∈R d |f (x)| ≤ 1,

  X; X (n) )f (Y j )1 R d \K (Y j ) -E(f (Y)1 R d \K (Y)|X) ≤ n j=1 w j (X; X (n) )1 R d \K (Y j ) -E(1 R d \K (Y)|X) + 2E(1 R d \K (Y)|X).

For simplicity, in this introduction, we tacitly assume all distributions to be Lebesgue-absolutely continous.

The existence of the regular conditional probability is a direct consequence of the disintegration theorem (see, e.g., Theorem 2.5.1 inLehmann and Romano (2005)).

We are dropping the superscripts (N ) and (n) when no confusion is possible.

Note that, although N does not appear in the notation, Q (n) w,± depends on both N and n.

This is always the case for d = 2 and d = 3: see[START_REF] Figalli | On the continuity of center-outward distribution and quantile functions[END_REF].

The research of Eustasio del Barrio is partially supported by FEDER, Spanish Ministerio de Economía y Competitividad, grant MTM2017-86061-C2-1-P and Junta de Castilla y León, grants VA005P17 and VA002G18. The research of Alberto González-Sanz is partially supported by the AI Interdisciplinary Institute ANITI, which is funded by the French "Investing for the Future -PIA3" program under the Grant agreement ANR-19-PI3A-0004.

of Statistics and its Application 9 401-424. Hallin, M., Hlubinka, D. and Hudecová, Š. (2022). Efficient Fully Distribution-Free Center-Outward Rank Tests for Multiple-Output Regression and MANOVA. Journal of the American Statistical Association 0 1-43. Hallin, M., La Vecchia, D. and Liu, H. (2020). Center-outward R-estimation for semiparametric VARMA models. Journal of the American Statistical Association 0 1-14. Hallin, M., La Vecchia, D. and Liu, H. (2022). Rank-based testing for semiparametric VAR models: a measure transportation approach. Bernoulli 0 1-14. Hallin, M., Paindaveine, D. and Šiman, M. (2010). Multivariate quantiles and multipleoutput regression quantiles: from L 1 optimization to halfspace depth. The Annals of Statistics 38 635-703. Hallin, M. and Šiman, M. (2018). Multiple-output quantile regression. In Handbook of Quantile Regression (R. Koenker, V. Chernozhukov, X. He and L. Peng, eds.) 185-207. CRC Press. Hallin, M., Lu, Z., Paindaveine, D. and Šiman, M. (2015). Local bilinear multipleoutput quantile/depth regression. Bernoulli 21 1435-1466. Hallin, M., del Barrio, E., Cuesta-Albertos, J. and Matrán, C. (2021). Distribution and quantile functions, ranks and signs in dimension d: a measure transportation approach. The Annals of Statistics 49 1139-1165. Hütter, J.-C. and Rigollet, P. (2021). Minimax estimation of smooth optimal transport maps. The Annals of Statistics 49 1166-1194. Koenker, R. (2005). Quantile Regression. Econometric Society Monographs. Cambridge University Press. Koenker, R. and Bassett, G. (1978). Regression quantiles. Econometrica 46 33-50. Koenker, R. and Xiao, Z. (2006). Quantile Autoregression [with Comments, Rejoinder]. Journal of the American Statistical Association 101 980-1006. Koenker, R., Chernozhukov, V., He, X. and Peng, L., eds. (2018). Handbook of Quantile Regression. CRC Press. Konen, D. and Paindaveine, D. (2022). Multivariate ρ-quantiles: a spatial approach. Bernoulli. Kong, L. and Mizera, I. (2012). Quantile tomography: using quantiles with multivariate data. Statistica Sinica 22 1589-1610. Lehmann, E. L. and Romano, J. P. (2005). Testing Statistical Hypotheses. Springer. Lévy, B., Mohayaee, R. and von Hausegger, S. (2020). A fast semi-discrete optimal transport algorithm for a unique reconstruction of the early Universe. arXiv:2012.09074. Makkuva, A., Taghvaei, A., Oh, S. and Lee, J. (2020). Optimal transport mapping via input convex neural networks. In International Conference on Machine Learning 6672-6681. PMLR. Manole, T., Balakrishnan, S., Niles-Weed, J. and Wasserman, L. (2021). Plugin estimation of smooth optimal transport maps. arXiv:2107.12364. McCann, R. J. (1995). Existence and uniqueness of monotone measure-preserving maps. Duke Mathematical Journal 80 309-323. Merlo, L., Petrella, L., Salvati, N. and Tzavidis, N. (2022). Marginal M-quantile regression for multivariate dependent data. Computational Statistics & Data Analysis 107500. Meyron, J. (2019). Initialization procedures for discrete and semi-discrete optimal transport. Computer-Aided Design 115 13-22. Paindaveine, D. and Šiman, M. (2011). On directional multiple-output quantile regres-

Appendix

A.1. Proofs for Section 3

The convergence described of Theorem 3.2 is based on the topology of set-valued maps, in particular the Graphical convergence (or Painlevé-Kuratowski convergence of the graphs). Recall from [START_REF] Rockafellar | Variational Analysis[END_REF] that a sequence of set-valued maps {T n } n converges graphically to another set-valued map T if -the outer limit lim sup n T n -which is the set of (x, y) ∈ R d × R d for which there exists a sequence {(x n , y n )} with (x n , y n ) ∈ T n containing a subsequence which converges

Using Markov's inequality in (A.17), we obtain

for all n > n 1 . Finally using (A.15), (A.18) and (A.16), we conclude that

for all n > N = max(n 0 , n 1 ).

Proof of Lemma A.2. Due to the fact that finite second-order moments are not required in our setting, Theorem 2.8 in del Barrio and Loubes (2019) does not directly apply. Their proof, however, relies on the weak convergence of the couplings (the joint measures solving the Kantorovich problem). In our case, we can prove a similar result using Lemma 9 in McCann (1995). Indeed, since the sequences {µ n } n and {ν n } n are tight with respect to weak convergence, the same result holds also for the class Γ(µ n , ν n ) of probabilities on R d ×R d with marginals {µ n } n and {ν n } n , see Lemma 4.4 in [START_REF] Villani | Optimal Transport: Old and New[END_REF]. Note that all the measures π n , n ∈ N, belong to Γ(µ n , ν n ). Denote by π ∈ P(R d × R d ) the weak limit of {π n } n along a subsequence; for simplicity, we keep the index n for the subsequence. Lemma 9 (ii) in McCann (1995) implies that the marginals of π are µ and ν. Moreover, the support of π n is cyclically monotone: indeed, it is contained in the subdifferential of φ n . Therefore, using Lemma 9 (i) in McCann (1995), we know that also π is supported on a cyclically monotone set. Corollary 14 in McCann (1995) yields-since µ is uniformly continuous with respect to the Lebesgue measure-the existence of a unique measure with cyclically monotone support and marginals µ and U d . As a consequence, π = (Id × φ) #µ. Since this holds along all possible subsequence, we have the weak convergence of (Id × φ n ) #µ n to (Id × φ). At this point, to conclude the proof of Lemma A.2, we can repeat verbatim the rest of the proof of Theorem 2.8 in del [START_REF] Del Barrio | Central limit theorems for empirical transportation cost in general dimension[END_REF].

Proof of Lemma