Upper and lower risk bounds for estimating the Wasserstein barycenter of random measures on the real line

,


December 8, 2017
Abstract This paper is focused on the statistical analysis of probability measures ν 1 , . . . , ν n on R that can be viewed as independent realizations of an underlying stochastic process. We consider the situation of practical importance where the random measures ν i are absolutely continuous with densities f i that are not directly observable. In this case, instead of the densities, we have access to datasets of real random variables (X i,j ) 1≤i≤n; 1≤j≤pi organized in the form of n experimental units, such that X i,1 , . . . , X i,pi are iid observations sampled from a random measure ν i for each 1 ≤ i ≤ n. In this setting, we focus on first-order statistics methods for estimating, from such data, a meaningful structural mean measure. For the purpose of taking into account phase and amplitude variations in the observations, we argue that the notion of Wasserstein barycenter is a relevant tool. The main contribution of this paper is to characterize the rate of convergence of a (possibly smoothed) empirical Wasserstein barycenter towards its population counterpart in the asymptotic setting where both n and min 1≤i≤n p i may go to infinity. The optimality of this procedure is discussed from the minimax point of view with respect to the Wasserstein metric. We also highlight the connection between our approach and the curve registration problem in statistics. Some numerical experiments are used to illustrate the results of the paper on the convergence rate of empirical Wasserstein barycenters.

Introduction
In this paper, we are concerned with the statistical analysis of a set of absolutely continuous measures ν 1 , . . . , ν n on the real line R, with supports included in (a possibly unbounded) interval Ω ⊂ R, that can be viewed as independent copies of an underlying random measure ν. In this setting, it is of interest to define and estimate a mean measure ν 0 of the random probability measure ν. The notion of mean or averaging depends on the metric that is chosen to compare elements in a given data set. In this work, we consider the Wasserstein metric d W associated to the quadratic cost for the comparison of probability measures and we define ν 0 as the population Wasserstein barycenter of ν, given by where the above expectation is taken with respect to the distribution of ν, and W 2 (Ω) denotes the space of probability measures with support included in Ω and with finite second moment. A Wasserstein barycenter corresponds to the Fréchet mean [Fré48] that is an extension of the usual Euclidean mean to non-linear metric spaces. Throughout the paper, the population mean measure ν 0 is also referred to as the structural mean of ν, which is a terminology borrowed from curve registration (see [ZM11] and references therein). We choose to work with the Wasserstein metric because it has been shown to be successful in the presence of phase variation (we refer to Section 2 for more explanations). Data sets leading to the analysis of absolutely continuous measures appear in various research fields. Examples can be found in neuroscience [WS11], demographic and genomics studies [Del11,ZM11], economics [KU01], as well as in biomedical imaging [PM16]. Nevertheless, in such applications one does not directly observe raw data in the form of absolutely continuous measures. Indeed, we generally only have access to random observations sampled from different distributions, that represent independent subjects or experimental units.
Thus, we propose to study the estimation of the structural mean measure ν 0 (the population Wasserstein barycenter) from a data set consisting of independent real random variables (X i,j ) 1≤i≤n; 1≤j≤p i organized in the form of n experimental units, such that (conditionally on ν i ) the random variables X i,1 , . . . , X i,p i are iid observations sampled from the measure ν i with density f i , where p i denotes the number of observations for the i-th subject or experimental unit. The main purpose of this paper is to propose nonparametric estimators of the structural mean measure ν 0 and to characterize their rates of convergence with respect to the Wasserstein metric in the asymptotic setting, where both n and min 1≤i≤n p i may go to infinity. 2

Main contributions
Two types of nonparametric estimators are considered in this paper. The first one is given by the empirical Wasserstein barycenter of the set of measuresν 1 , . . . ,ν n , withν i = 1 p i p i j=1 δ X i,j for 1 ≤ i ≤ n. This estimator will be referred to as the non-smoothed empirical Wasserstein barycenter. Alternatively, since the unknown probability measures ν i are supposed to be absolutely continuous, a second estimator is based on a preliminary smoothing step which consists in using standard kernel smoothing to construct estimatorsf i of the unknown densities f i for each 1 ≤ i ≤ n. Then, an estimator of ν 0 is obtained by taking the empirical Wasserstein barycenter of the measuresν i , . . . ,ν n , withν i (A) := Af i (x)dx, A ⊂ R measurable. We refer to this class of estimators as smoothed empirical Wasserstein barycenters whose smoothness depend on the choice of the bandwidths in the preliminary kernel smoothing step.
The rates of convergence of both types of estimators are derived for their (squared) Wasserstein risks, defined as their expected (squared) Wassertein distances from ν 0 , and their optimality is discussed from the minimax point of view. Finally, some numerical experiments with simulated data are used to illustrate these results.

Related work in the literature
The notion of barycenter in the Wasserstein space, for a finite set of n probability measures supported on R d (for any d ≥ 1), has been recently introduced in [AC11] where a detailed characterization of such barycenters in terms of existence, uniqueness and regularity is given using arguments from duality and convex analysis. However, the convergence (as n → ∞) of such Wasserstein barycenters is not considered in that work.
In the one dimensional case (d = 1), computing the Wasserstein barycenter of a finite set of probability measures simply amounts to averaging (in the usual way) their quantile functions. In statistics, this approach has been referred to as quantile synchronization [ZM11]. In the presence of phase variability in the data, quantile synchronization is known to be an appropriate alternative to the usual Euclidean mean of densities to compute a structural mean density that is more consistent with the data. Various asymptotic properties of quantile synchronization are studied in [ZM11] in a statistical model and asymptotic setting similar to that of this paper with min 1≤i≤n p i ≥ n. However, other measures of risk than the one in this paper are considered in [ZM11], but the optimality of the resulting convergence rates of quantile synchronization is not discussed.
The results of this paper are very much connected with those in [PZ16] where a new framework is developed for the registration of multiple point processes on the real line for the purpose of separating amplitude and phase variation in such data. In [PZ16], consistent estimators of the structural mean of multiple point processes are obtained by the use of smoothed Wasserstein barycenters with an appropriate choice of kernel smoothing. Also, rates of convergence of such estimators are derived for the Wasserstein metric. The statistical analysis of multiple point processes is very much connected to the study of repeated observations organized in samples from independent subjects or experimental units. Therefore, some of our results in this paper on smoothed empirical Wasserstein barycenters are built upon the work in [PZ16]. Nevertheless, novel contributions include the derivation of an exact formula to compute the risk of non-smoothed Wasserstein barycenters in the case of samples of equal size, and new upper bounds on the rate of convergence of the Wasserstein risk of non-smoothed and smoothed empirical Wasserstein barycenters, together with a discussion of their optimality from the minimax point of view.
The construction of consistent estimators of a population Wasserstein barycenter for semiparametric models of random measures can also be found in [BK17] and [BLGL15], together with a discussion on their connection to the well known curve registration problem in statistics [RL01,WG97].

Organization of the paper
In Section 2, we first briefly explain why using statistics based on the Wasserstein metric is a relevant approach for the analysis of a set of random measures in the presence of phase and amplitude variations in their densities. Then, we introduce a deformable model for the registration of probability measures that is appropriate to study the statistical properties of empirical Wasserstein barycenters. The two types of nonparametric estimators described above are finally introduced at the end of Section 2. The convergence rates and the optimality of these estimators are studied in Section 3. Some numerical experiments with simulated data are proposed in Section 4 to highlight the finite sample performances of these estimators. Section 5 contains a discussion on the main contributions of this work and their potential extensions. The proofs of the main results are gathered in a technical Appendix. Finally, note that we use bold symbols f , ν, . . . to denote random objects (except real random variables).
2 Wasserstein barycenters for the estimation of the structural mean in a deformable model of probability measures 2.1 The need to account for phase and amplitude variations To estimate a mean measure from the data (X i,j ) 1≤i≤n; 1≤j≤p i , a natural approach is the following one. In a first step, one uses the X i,j 's to compute estimatorsf 1 , . . . ,f n (e.g. via kernel smoothing) of the unobserved density functions f 1 , . . . , f n of the measures ν 1 , . . . , ν n . Then, an estimator of a mean density might be defined as the usual Euclidean meanf n = 1 n n i=1f i , which is also classically referred to as the cross-sectional mean in curve registration. At the level of measures, it corresponds to computing the arithmetical mean measureν n = 1 n n i=1ν i . The Euclidean meanf n is the Fréchet mean of thef i 's with respect to the usual squared distance in the Hilbert space L 2 (Ω) of square integrable functions on Ω. Therefore, it only accounts for linear variations in amplitude in the data. However, as remarked in [ZM11], in many applications it is often of interest to also incorporate an analysis of phase variability (i.e. time warping) in such functional objects, since it may lead to a better understanding of the structure of the data. In such settings, the use of the standard squared distance in L 2 (Ω) to compare density functions ignores a possible significant source of phase variability in the data.
To better account for phase variability in the data, it has been proposed in [ZM11] to introduce the so-called method of quantile synchronization as an alternative to the cross sectional 4 meanf n . It amounts to computing the mean measure ν ⊕ n (and, if it exists, its density f ⊕ n ) whose quantile function is where F − i denotes the quantile function of the measure ν i with density f i . The statistical analysis of quantile synchronization, as studied in [ZM11], complements the quantile normalization method originally proposed in [BIAS03] to align density curves in microarray data analysis. This method is therefore appropriate for the registration of density functions and the estimation of phase and amplitude variations as explained in details in [PZ16].
Let us now assume that ν 1 , . . . , ν n are random elements taking values in the set of absolutely continuous measures contained in W 2 (Ω). In this setting, it can be checked (see e.g. Proposition 2.1 below) that quantile synchronization corresponds to computing the empirical Wasserstein barycenter of the random measures ν 1 , . . . , ν n , namely Therefore, the notion of averaging by quantile synchronization corresponds to using the Wasserstein distance d W to compare probability measures, which leads to a notion of measure averaging that may better reflect the structure of the data than the arithmetical mean in the presence of phase and amplitude variability.
To illustrate the differences between using Euclidean and Wasserstein distances to account for phase and amplitude variation, let us assume that the measures ν 1 , . . . , ν n have densities f 1 , . . . , f n obtained from the following location-scale model: we let f 0 be a density on R having a finite second moment and, for (a i , b i ) ∈ (0, ∞) × R, i = 1, . . . , n, a given set of iid random vectors, we define The sources of variability of the densities from model (2.2) are the variation in location along the x-axis, and the scaling variation. In Figure 1 In this numerical experiment, there is more variability in phase (i.e. location) than in amplitude (i.e. scaling), which can also be observed at the level of quantile functions as shown by Figure 1(b).
In the location-scale model (2.2), it can be checked, e.g. using the quantile averaging formula (2.1), that the empirical Wasserstein barycenter ν ⊕ n is the probability measure with density Hence, if we assume that E(a 1 ) = 1 and E(b 1 ) = 0, it follows that d 2 W (ν ⊕ n , ν 0 ) converges almost surely to 0 as n → ∞, meaning that ν ⊕ n is a consistent estimator of ν 0 as shown by Figure 1 6 ν n is clearly not a consistent estimator of ν 0 , as it can be observed in Figure 1(d).
Remark 2.1. It is clear that, in the above location-scale model, one may easily prove that f ⊕ n converges almost surely to f 0 as n → ∞ for various distances between density functions as illustrated by Figure 1(e). However, in this paper, we restrict our attention to the problem of how the structural mean measure ν 0 can be estimated from empirical Wasserstein barycenters with respect to the Wasserstein distance between probability measures. Showing that the density (if it exists) of such estimators converges to the density f 0 of ν 0 is not considered in this work.

Barycenters in the Wasserstein space
Let Ω be an interval of R, possibly unbounded. We let W 2 (Ω) be the set of probability measures over (Ω, B(Ω)), with finite second moment, where B(Ω) is the σ-algebra of Borel subsets of Ω. We also denote by W ac 2 (Ω) the set of measures ν ∈ W 2 (Ω) that are absolutely continuous with respect to the Lebesgue measure on R. The cumulative distribution function (cdf), the quantile function and the density function (if ν ∈ W ac 2 (Ω)) of ν are denoted respectively by F ν , F − ν and f ν .
It can be shown that W 2 (Ω) endowed with d W is a metric space, usually called Wasserstein space. For a detailed analysis of W 2 (Ω) and its connection with optimal transport theory, we refer to [Vil03].
Definition 2.2 (Square-integrability). The random measure ν is said to be square-integrable if for some (thus for every) µ ∈ W 2 (Ω).
Observe that the expectation in the definition above is well defined since ν → d W (µ, ν) is continuous and therefore, measurable. In the rest of the paper we assume that random measures ν are square integrable, in the sense of the previous definition, and so, this property is not explicitly stated in definitions or results involving ν . We use the notations F , F − and f to denote respectively the cumulative distribution function, the quantile function and the density function (if it exists) of ν.

Definition 2.3 (Population and empirical Wasserstein barycenters). A population Wasserstein barycenter of ν is defined as a minimizer of
An empirical Wasserstein barycenter of ν 1 , . . . , ν n ∈ W 2 (Ω) is defined as a minimizer of Proposition 2.1.

A deformable model of probability measures
Let ν 1 , . . . , ν n be independent copies of the random measure ν, with barycenter ν 0 . We consider a deformable model of random probability measures satisfying the following assumptions: Assumption 2.3. Conditionally on ν i , the observations X i,1 , . . . , X i,p i are iid random variables sampled from ν i , where p i ≥ 1 is a known integer, for 1 ≤ i ≤ n.
Remark 2.2. Similar assumptions are considered in [PZ16] to characterize a population barycenter in W 2 (Ω) for the purpose of estimating phase and amplitude variations from the observations of multiple point processes. For examples of parametric models satisfying the Assumptions 2.1-2.3, we refer to [BK17] and [BLGL15]. The main restriction of this deformable model is that ν 0 is assumed to be absolutely continuous.
Remark 2.3. Observe that Assumption 2.1 requires W ac 2 (Ω) to be a measurable subset of W 2 (Ω). The proof of this technical result will appear in a forthcoming paper. 8

Non-smoothed empirical barycenter
To estimate the structural mean measure ν 0 from the data (X i,j ) 1≤i≤n; 1≤j≤p i , a first approach consists in computing straightaway the barycenter of the empirical measureν 1 , . . . ,ν n wherẽ . The non-smoothed empirical barycenter is thus defined aŝ ν n,p = arg min where p = (p 1 , . . . , p n ). In the case p 1 = p 2 = . . . = p n = p, we have the following procedure for computing the non-smoothed empirical barycenter. For each 1 ≤ i ≤ n, we denote by Thanks to Proposition 2.1, the quantile function of the empirical Wasserstein barycenter is the average of the quantile functions ofν 1 , . . . ,ν n , and thus we obtain the formulâ (2.5) Note that we useν n,p instead ofν n,p to denote the non-smoothed empirical barycenter in the case p 1 = p 2 = . . . = p n = p.

Smoothed empirical barycenter
An alternative approach is to use a smoothing step to obtain estimated densities and then compute the barycenter.
1. In a first step we use kernel smoothing to obtain estimatorsf h 1 1 , . . . ,f hn n of f 1 , . . . , f n , where h 1 , . . . , h n are positive bandwidth parameters, that may be different for each subject or experimental unit. In this paper we investigate a non-standard choice for the kernel function, proposed in [PZ16], to analyze the convergence of smoothed empirical barycenter in W 2 (Ω). In Section 3 we give a precise definition of the resulting estimators. However, at this point it is not necessary to go into such details.
2. In a second step, an estimator of ν 0 is given byν h n,p , defined as the measure whose quantile function is given bŷ 3 Convergence rate for estimators of the population Wasserstein barycenter In this section we discuss the rates of convergence of the estimatorsν n,p andν h n,p , that are respectively characterized by equations (2.4) and (2.7). Some of the results presented below are using the work in [BL17], on a detailed study of the variety of rates of convergence of an empirical measure on the real line toward its population counterpart in the Wasserstein metric. Then, we discuss the optimality of these estimators from the minimax point of view following the guidelines in nonparametric statistics to derive optimal rates of convergence (see e.g. [Tsy09] for an introduction to this topic).

Non-smoothed empirical barycenter in the case of samples of equal size
Let us first characterize the rate of convergence ofν n,p , in the specific case where samples of observations per unit are of equal size, namely when p 1 = p 2 = . . . = p n = p. In what follows, we let Y 1 , . . . , Y p be iid random variables sampled from the population mean measure ν 0 (independently of the data (X i,j ) 1≤i≤n; 1≤j≤p ), and we denote by µ p = 1 p p k=1 δ Y k the corresponding empirical measure.
Theorem 3.1 provides exact formulas to compute the rate of convergence (for the expected squared Wasserstein distance) ofν n,p . Formula (3.1) relies on the computation of the variances of the order statistics of iid variables Y 1 , . . . , Y p sampled from the population mean measure ν 0 , and on the computation of the rate of convergence of E d 2 W (µ p , ν 0 ) . However, deriving a sharp rate of convergence forν n,p using inequality (3.1) requires computing the variances of the order statistics of iid random variables. To the best of our knowledge, obtaining a sharp estimate for Var Y * j for any 1 ≤ j ≤ p remains a difficult task except for specific distributions. For example, if ν 0 is assumed to be a log-concave measure, then it is possible to use the results in Section 6 in [BL17] which provide sharp bounds on the variances of order statistics for such probability measures. We discuss below some examples where equality (3.1) may be used to derive a sharp rate of convergence forν n,p .
Moreover, from Theorem 4.7 in [BL17], it follows that E d 2 W (µ p , ν 0 ) = 1 6p . Therefore, thanks to (3.1) we obtain Equality (3.2) thus shows that, when ν 0 is the uniform distribution on [0, 1], the rate of convergence ofν n,p is given by and this rate is sharp.
Beyond the specific case where ν 0 is a uniform distribution, it is in general difficult to compute E d 2 W (µ p , ν 0 ) . Nevertheless, thanks to Theorem 4.3 in [BL17], we have the following bounds for any distribution ν 0 ∈ W 2 (Ω).
The case where ν 0 is the one-sided exponential distribution. Combining inequality (3.4) with (3.1), it follows that Now (using e.g. Remark 6.13 in [BL17]) one has that if ν 0 is the one-sided exponential distribution (with density e −x , for Therefore, there exist a constant c > 0 such that for all sufficiently large p. Hence, when ν 0 is the exponential distribution the above inequalities show that the rate of convergence ofν n,p is O 1 The case where ν 0 is a Gaussian distribution. By Theorem 4.3 and Corollary 6.14 in [BL17] there exist constants c 1 , c 2 > 0 such that Therefore, combining the above upper bound with (3.5), one finally has that when ν 0 is the standard Gaussian. In this setting, the rate of convergence is thus O 1 n + 1 np + 1 p log log p .
Upper bounds in more general cases. If one is interested in deriving an upper bound on E d 2 W (ν n,p , ν 0 ) for a larger class of measures ν 0 ∈ W 2 (Ω) (e.g. beyond the log-concave case), another approach is as follows. Noting that the term 1−n pn p j=1 Var(Y * j ) in equality (3.1) is negative, a straightforward consequence of Theorem 3.1 is the following upper bound Then, thanks to inequality (3.9), to derive the rate of convergence ofν n,p , it remains to control the rate of convergence of the empirical measure µ p to ν 0 for the expected squared Wasserstein distance. This issue is discussed in detail in [BL17]. In particular, the work in [BL17] describes a variety of rates for the expected distance E d 2 W (µ p , ν 0 ) , from the standard one O 1 p to slower rates. For example, by Theorem 5.1 in [BL17], the following upper bound holds where, the so-called J 2 -functional is defined as The J 2 functional is shown to be measurable in Proposition A.1 of the Appendix. Therefore, provided that J 2 (ν 0 ) is finite, the empirical measure µ p converges to ν 0 at the rate O 1 p . Hence, using inequality (3.10), we have: By Corollary 3.1, if J 2 (ν 0 ) < ∞, then it follows thatν n,p converges to ν 0 at the rate O 1 n + 1 p . Hence, in the setting where p ≥ n and J 2 (ν 0 ) < ∞,ν n,p converges at the classical parametric rate O 1 n . The case p ≥ n, usually refereed to as the dense case in the literature on functional data analysis (see e.g. [LH10] and references therein), corresponds to the situation where the number of observations per unit/subject is larger than the sample size n of functional objects. In the sparse case (when p < n), the non-smoothed Wasserstein barycenter converges at the rate O 1 p , if J 2 (ν 0 ) < ∞. Remark 3.1. When ν 0 is the uniform distribution on [0, 1] one has that J 2 (ν 0 ) < ∞, but we have shown that E d 2 W (ν n,p , ν 0 ) 1 n + 1 np + 1 p 2 . Hence, in this setting,ν n,p converges at the parametric rate O 1 n provided that p ≥ √ n, which is a dense regime condition weaker than p ≥ n.
To conclude this discussion on the rate of convergence of the non-smoothed Wasserstein barycenter in the case of samples of equal size, we study in more detail the control of the rate of convergence of the term E d 2 W (µ p , ν 0 ) in inequality (3.9). As pointed out in many works (see for example [dBGU05,BL17] and the references therein) the finiteness of J 2 (ν 0 ) is the key point to control the convergence of the empirical measure µ p to the population measure ν 0 in the Wasserstein space. Some known facts concerning this issue are the following.
1. If J 2 (ν 0 ) < ∞ then ν 0 is supported on an interval of R and its density is a.e. strictly positive on this interval.
2. If ν 0 is compactly supported with a density bounded away from zero or with a log-concave density then J 2 (ν 0 ) < ∞.
3. If the density of ν 0 is of the form C α e −|x| α then J 2 (ν 0 ) is finite if and only if α > 2. In particular, J 2 (ν 0 ) = ∞ for the Gaussian distribution.
Some further comments can be made in the case where ν 0 is Gaussian. In this setting, one has that J 2 (ν 0 ) = ∞ and the rate of convergence of E d 2 W (µ p , ν 0 ) to zero is slower than O 1 p . Indeed, from Corollary 6.14 in [BL17], if ν 0 is the standard Gaussian, then the rate , which leads to the upper bound (3.8) for the non-smoothed Wasserstein barycenterν n,p in the Gaussian case. Hence, thanks to (3.8), one has that if p is sufficiently large with respect to n (namely when p ≥ n log log p), thenν n,p also converges at the classical parametric rate O 1 n when ν 0 is the standard Gaussian distribution. Remark 3.2. Following the work in [BL17], if ν 0 has a log-concave distribution, then one may obtain rates of convergence for E d 2 W (ν n,p , ν 0 ) that are slower than the standard O 1 p rate 13 (e.g. for beta or exponential distributions). Moreover, it is also possible to considerer for any q ≥ 1 and for any absolutely continuous probability ν on Ω, the functional in order to control the rate of convergence of the empirical measure to ν for the q-Wasserstein distance.

Non-smoothed empirical barycenter in the general case
Let us now consider the general situation where the p i 's are possibly different. The result below gives an upper bound on the rate of convergence ofν n,p where p = (p 1 , . . . , p n ).
For the random measure ν, we consider the extended random variable J 2 (ν) (see Proposition A.1). Since the ν i 's are independent copies of ν, by applying inequality (3.10) it follows that Hence, from Theorem 3.2, we finally obtain the following upper bound on the rate of convergence for the non-smoothed empirical barycenter Remark 3.3. Knowing if J 2 (ν) has a finite expectation is in general a difficult task. But, if we assume that the density f of ν is bounded below by a non-random positive constant then (obviously) E [J 2 (ν)] < ∞.

The case of smoothed empirical barycenters
In this section, we assume that Ω = [0, 1] and we discuss the rate of convergence of smoothed empirical barycentersν h n,p (note that the following results hold if Ω is any compact interval). To choose an appropriate kernel function to study the convergence rate of the estimatorν h n,p , we follow the proposal made in [PZ16]. We let ψ be a positive, smooth and symmetric density on the real line, such that R x 2 ψ(x)dx = 1. We also denote by Ψ the cdf of the density ψ and, for a bandwidth parameter h > 0, we let ψ h (x) = 1 h ψ x h . Then, for any y ∈ [0, 1] and h > 0, we denote by µ y h the measure supported on [0, 1] whose density f µ y h is defined as where h i > 0 is a bandwidth parameter depending on i. For a discussion on the intuition for this choice of kernel smoothing, we refer to [PZ16]. A key property to analyze the convergence rate ofν h n,p is the following lemma which relates the Wasserstein distance betweenν h i i and the empirical measureν i = 1 Remark 3.4. The results in [PZ16] strongly depend on the assumption that Ω is compact. To go beyond this assumption, one should be able to extend (in an appropriate way) the density f µ y h (x) to a non-compact setting, which we believe to be a difficult task. Nevertheless, it should be remarked that our results on non-smoothed empirical Wasserstein barycenter hold in the general case where Ω = R.
Lemma 3.1. Let 1 ≤ i ≤ n. Suppose that 0 < h i ≤ 1/4, then one has the following upper bound for h i small enough and some constant C ψ > 0 depending only on ψ.
Proof. The upper bound (3.16) follows immediately from Lemma 1 in [PZ16] and the symmetry of ψ. Then, by applying inequality (3.17) and since ψ is symmetric, it follows that for h small enough Hence, the second part of Lemma 3.1 is a consequence of the above inequality, the fact that α ≥ 5, and the upper bound (3.16), which completes the proof.
The result below gives a rate of convergence for the estimatorν h n,p .
Theorem 3.3. Suppose that Assumptions 2.1, 2.2 and 2.3 are satisfied, and that the density ψ, used to define kernel smoothing in (3.15), satisfies inequality (3.17). If J 2 (ν) has a finite expectation, and the bandwidth parameters h i are small enough, then we have Theorem 3.3 can then be used to discuss choices of bandwidth parameters that may lead to a parametric rate of convergence. For example, if 0 < h i ≤ n −1/2 for all 1 ≤ i ≤ n and min 1≤i≤n p i ≥ n (dense case), then Theorem 3.3 implies that (for sufficiently large n to ensure that max 1≤i≤n {h i } is small enough) Remark 3.5. In the dense case (namely min 1≤i≤n p i ≥ n) it can be seen, by comparing the upper bounds (3.13) and (3.19), that a preliminary smoothing step of the data (namely kernel smoothing the empirical measuresν i = 1 does not improve the parametric rate of convergence n −1/2 . Moreover, the bandwidth values have to be small to ensure the rate of convergence n −1/2 for E d W (ν h n,p , ν 0 ) . This result comes from the fact that we evaluate the risk of empirical barycenters at the level of measures in W 2 (Ω), and that we do not aim to control an estimation of the density f 0 of the population mean measure ν 0 .
Remark 3.6. Theorem 3.3 shares similarities with the results from Theorem 2 in [PZ16], which gives the rate of convergence for smoothed Wasserstein barycenters, computed from the realizations of multiple Poisson processes in a deformable model of measures similar to that of this paper. The main difference in [PZ16] is that the number p i of observations for each experimental unit are independent Poisson random variables with expectation E(p i ) = τ n for each 1 ≤ i ≤ n (they are not deterministic integers). From such observations and under similar assumptions, it is proved in [PZ16] that the following upper bound holds (in probability) (3.20) The upper bound (3.20) is very similar to (3.18), proposed in this paper. Both approaches essentially split the distance d W (ν h n,p , ν 0 ) into three terms: 1. a parametric term of the order √ n coming from the observation of a sample of size n of the random measure ν, 2. a term involving kernel smoothing and 16 3. a term involving the Wasserstein distance between ν 0 and its empirical counterpart Under the conditions that τ n ≥ O(n 2 ) and max 1≤i≤n h i ≤ O P n −1/2 , it follows from Theorem 2 in [PZ16] thatν h n,p converges at the parametric rate O n −1/2 , for the Wasserstein distance. The quantity τ n represents the averaged number of points observed for each Poisson process. As remarked in [PZ16] the condition τ n ≥ O(n 2 ) corresponds to a dense sampling regime where the number n of observed Poisson processes should not grow too fast with respect to the expected number of points observed for each process. Comparing the upper bounds (3.18) and (3.20), the main difference in the control of the risk ofν h n,p between our approach and the one in [PZ16] is that we use the condition E [J 2 (ν)] < ∞. Under such an assumption, the smoothed Wasserstein barycenter (for the model considered in this paper) may be shown to converge at the rate O n −1/2 , for the expected Wasserstein distance, under the dense case setting p := min{p i , 1 ≤ i ≤ n} ≥ n which is somehow a weaker condition than E(p i ) ≥ n 2 for all 1 ≤ i ≤ n, as in [PZ16]. Therefore, it might be argued that there is a sort of "optimality gap" in [PZ16], and that the results in this paper are a first step towards closing this gap. But, on the other hand, the result in [PZ16] is more general because nothing is assumed about the finiteness of the functional J 2 . In particular, the upper bound (3.20) given in [PZ16] also holds when J 2 is infinite, which is not the case for the upper bound (3.18) in this paper.

A lower bound on the minimax risk
In the rest of this section, we show that, in the dense case and for the expected squared Wasserstein distance, the rate of convergence O n −1 for non-smoothed empirical Wasserstein barycenters is optimal from the minimax point of view over a large class of random measures ν satisfying the deformable model defined in Section 2.3 through Assumptions 2.1, 2.2 and 2.3.
Definition 3.2. Let A > 0. We denote by F(R, A) ⊆ W ac 2 (R) a given set of measures with variance bounded by A, which contains at least all Gaussian distributions with variance bounded by A.
Then, by inequality (3.9), we obtain the following corollary giving a uniform rate of convergence for the non-smoothed empirical barycenter in the case of samples of equal size.
it follows that The condition (3.21) may be interpreted as the generalization of the dense case setting that has been discussed in the previous sections as it is valid only if p is sufficiently large with respect to n. As an example, let A ≥ 0 and suppose that the set F(R, A) can be partitioned as where F 0 (R, A) denotes a set of measures ν 0 ∈ W ac 2 (R) with variance bounded by A satisfying while G(R, A) denotes the set of Gaussian distributions with variance bounded by A. For this example, it follows from inequalities (3.4), (3.7) and (3.10) in Section 3.1 (with samples of equal size) that provided that log log p ≥ 1, where c 2 is a constant from inequality (3.7). Hence, if p is such that p ≥ n log log p then condition (3.21) is satisfied with The following theorem shows that the upper bound (3.22) in Corollary 3.3 is optimal (in term of rate of convergence) from the minimax point of view in nonparametric statistics. whereν =ν ((X i,j ) 1≤i≤n; 1≤j≤p i ) denotes any estimator taking values in (W 2 (R), B (W 2 (R))) withν denoting a measurable function of the data (X i,j ) 1≤i≤n; 1≤j≤p i sampled from the deformable model defined in Section 2.3. Now, by using inequalities (3.13) and (3.19) and Definitions 3.1 and 3.2 introduced above, we also obtain the following corollary giving uniform rates of convergence for the non-smoothed Wasserstein barycenter in the general situation where the p i s are possibly different.
Corollary 3.4. Let A > 0 and σ > 0. Suppose that the assumptions of Corollary 3.2 are satisfied, and that p i ≥ n, for all 1 ≤ i ≤ n. Then, the following upper bound holds Hence, under the assumptions made in Corollary 3.4, the estimatorν n,p converges at the optimal rate of convergence n −1/2 provided that We conclude this discussion by a few remarks on the rate of convergence that may be obtained in the sparse case.
Remark 3.7. In the case of samples of equal size, the results above show that the rate of convergence n −1 is optimal in the dense case (for the risk E d 2 W (ν n,p , ν 0 ) ), namely when the number p = p 1 = . . . = p n of observations per units is sufficiently large with respect to n. We believe that deriving a lower bound on the minimax risk depending on p in the sparse case (e.g. when p < n) is more involved. Indeed, from the discussion in Section 3.1 on the rate of convergence of the non-smoothed empirical barycenter, it appears that the exact decay of E d 2 W (ν n,p , ν 0 ) as a function of p is difficult to establish as it depends on ν 0 . Indeed, from Section 3.1, one has that -if ν 0 is the uniform distribution on [0, 1], then E d 2 W (ν n,p , ν 0 ) From Theorem 3.1, one has that the risk of the non-smoothed empirical barycenter may be bounded from below as follows The quantity p j=1 j/p dα may be interpreted as a bias term when estimating the unknown measure by the nonparametric estimator µ p = 1 p p j=1 δ Y j . Therefore, for samples of equal size and in the sparse case (when p < n), the lower bound (3.24) may be used to control (as a function of p) the best rate of convergence forν n,p that may be obtained over the class of measures ν 0 ∈ F(R, A).
Remark 3.8. Finally, we remark that better rates of convergence may be obtained if one assumes a parametric model for the random measure ν. Indeed, suppose that µ 0 ∈ W ac 2 (Ω) denotes a known probability measure with expectation m 0 and variance σ 2 0 and consider that the data (X i,j ) 1≤i≤n; 1≤j≤p i are sampled from iid random measures ν 1 , . . . , ν n satisfying the location model where a 1 , . . . , a n are iid random variables with unknown expectationā and variance γ 2 . In this model, the population Wasserstein barycenter is the measure ν 0 with quantile function F − ν 0 (·) = F − µ 0 (·) +ā. Since, the measure µ 0 is assumed to be known, a natural estimator for ν 0 is to take the measureν 0 with quantile function F − ν 0 (·) = F − µ 0 (·) +â, witĥ Then, it is clear that In the case where all the p i 's are equal to p, then the above equality simplifies to and thus the parametric estimatorν 0 converges at the rate O 1 n + 1 np . Therefore, either in the dense (p ≥ n) or sparse case (p < n), the parametric estimatorν 0 converges at the rate O 1 n in the location model (3.25) when the "reference measure" µ 0 is known. Moreover, in the sparse case (p < n) the parametric estimatorν 0 converges faster than the non-smoothed empirical Wasserstein barycenterν n,p , thanks to the results in Section 3.1.

Numerical experiments
In this simulation study we perform Monte Carlo experiments to compare the decay of the squared Wassertein risks E d 2 W (ν h n,p , ν 0 ) and E d 2 W (ν n,p , ν 0 ) , of the smoothed and non-smoothed empirical Wasserstein barycentersν h n,p andν n,p , as a function of the number n of units and the sample size p. The theoretical results in this paper indicate that, in the dense case, both estimators converge at the optimal parametric rate O 1 n . However, in the sparse case it remains unclear if a preliminary smoothing step may improve the quality of estimation of the population Wasserstein barycenter. The purpose of these numerical experiments is thus to compare the behavior of smoothed and non-smoothed empirical Wasserstein barycenters, in these two settings, and analyze the influence of the number n of measures and the sample size p.
We analyze the case of random samples (X i,j ) 1≤i≤n; 1≤j≤p , with 10 ≤ n ≤ 200 and 10 ≤ p ≤ 200. Data are generated from densities supported on a compact interval Ω that are sampled from the following model, accounting for vertical and horizontal variations  2, 2]). This setting corresponds to the the simulation study conducted in [PM16]. For each choice of f (either the Gaussian or Uniform case), the interval Ω is taken such that each random function f i has a compact support included in Ω. Therefore, the population Wasserstein barycenter in model (4.1) is the measure with density f thanks to the fact that E(b i ) = 0 and E(a i ) = 1. The Gaussian case (resp. Uniform case) corresponds to the estimation of a Wasserstein barycenter having smooth (resp. non-differentiable) density f . For given values of n and p, we evaluate the Wasserstein risk ofν h n,p by repeating M = 100 times the following experiment. First, data are simulated from model (4.1). Then, for each 1 ≤ i ≤ n, we use kernel smoothing to compute the densityf h i i and its associated measurê ν h i i . We slightly deviate from the analysis carried out in Section 3, as we use a Gaussian kernel to smooth the data (X i,j ) 1≤j≤p , with bandwidth h i chosen by cross validation, instead of the specific kernel defined in (3.14), that has been proposed for the convergence analysis ofν h n,p . We found that this modification has no substantial effect on the finite sample performance of the procedure, and a similar choice has been made in the numerical experiments in [PZ16]. In Figure 2(a) (resp. Figure 3(a)) we display an example of densities estimated from realizations of the model (4.1) with n = p = 100 and f the truncated Gaussian density (resp. f the Uniform density). After computing the quantile function F − ν h n,p of the empirical smoothed Wasserstein barycenterν h n,p , we approximate d 2 ) 2 dα by discretizing the integral over a fine grid of values for α ∈]0, 1[. This approximated value of d 2 W (ν h n,p , ν 0 ) is then averaged over the M = 100 repeated experiments to approximate E d 2 W (ν h n,p , ν 0 ) . Thanks to the explicit expression (2.5) of the non-smoothed empirical Wasserstein barycenter ν n,p , its quantile function F − νn,p is straightforward to compute on a grid of values for α, and the Wasserstein risk E d 2 W (ν n,p , ν 0 ) is then approximated in the same way by using Monte Carlo repetitions.
For values of n and p ranging from 10 to 200, we display in Figure 2(c) and 2(d) (resp. Figure  3(c) and 3(d)) these approximations of E d 2 W (ν h n,p , ν 0 ) and E d 2 W (ν n,p , ν 0 ) (in logarithmic scale) for f the truncated Gaussian density (resp. f the Uniform density). For both estimators, it appears that the Wasserstein risk is clearly a decreasing function of the number n of units. To the contrary, increasing p does not lead to a significant decay of this risk. This suggest that 1 n Var(ν) is the most significant term in the upper bound (3.12) of the Wasserstein risk ofν n,p .
In Figure 2(b) and Figure 3(b), we also display the logarithm of the ratio It can be observed that: -When the population Wasserstein barycenter has a smooth density f (Gaussian case) then, for values of p larger than 100, both estimators (smoothed and non-smoothed empirical Wasserstein barycenters) appear to have squared Wasserstein risks of approximately the same magnitude. This tends to confirm the results on convergence rates obtained in Section 3, in the dense case (when p is sufficiently large with respect to n), which show that a preliminary smoothing is not necessary in this setting. For smaller values of p (between 10 and 50), the smoothed empirical Wasserstein barycenter has a smaller Wasserstein risk. This suggests that introducing a smoothing step through kernel smoothing of the data, in each experimental unit, may improve the quality of the estimation of ν 0 , when the sample size p is small (which corresponds to the sparse case), in the setting where the population Wasserstein barycenter is a distribution with a smooth density. In this example, Gaussian kernel smoothing is particularly well suited, which may explain the better performances obtained with a smoothed empirical Wasserstein barycenter in the sparse case.
-When the population Wasserstein barycenter has a non-smooth density f (Uniform case) then, except for very small values of p ≤ 10, the non-smoothed empirical Wasserstein barycenter has always a lower squared Wasserstein risk, both in the sparse and dense cases. A preliminary smoothing with a Gaussian kernel does not improve the estimation of a Wasserstein barycenter having piecewise constant density and, in this setting, such a step is thus not necessary.

Conclusion and perspectives
In this paper we have studied the rate of convergence for the (squared) Wasserstein distance of (possibly smoothed) empirical barycenters in a deformable model of measures. The main contributions of this work can be summarized as follows. In the case of samples of equal size, we have derived a closed-form expression for the risk of non-smooth empirical barycenter, as a function of n and p, which allows to derive sharp rates of convergence whose decay in p depends on the population mean measure ν 0 . A second conclusion of the paper is that, in the dense case (when the minimal number min 1≤i≤n p i ≥ n of observations per unit is sufficiently large with respect to the number n of observed measures), the non-smoothed empirical barycenter converges at the parametric rate of convergence n −1 . Moreover, this rate is shown to be a lower bound on the decay of a novel notion of minimax risk, in the deformable model of measures introduced in this paper. In the dense case, the numerical experiments that have been carried out are in agreement with the theoretical results which show that, in this setting, one may only consider the non-smoothed empirical Wasserstein barycenter and that a preliminary smoothing step is not necessary to obtain an optimal estimator. A first perspective would be to find a lower bound on the minimax risk depending on p in the sparse case. However, to this end, we believe that one has to first obtain sharper rates of convergence as a function of p, for the non-smoothed empirical barycenter.
Finally, a natural perspective is to ask how these results can be extended to higher dimensional settings for measures supported on R d , with d > 1. However, we believe that this is far from being obvious as the results in this paper rely heavily on the closed-form formula of Wasserstein barycenters, in the one-dimensional setting though quantile averaging. Such results do not hold in higher-dimension, for data sets consisting of iid random vectors, sampled from unknown random measures supported on R 2 or R 3 , for example.     Var Y * j .
Proof. Let 1 ≤ j ≤ p and 1 ≤ i ≤ n. Thanks to the expression (A.3) for the density of X * i,j , one has that where we used the change of variable α = F i (x) to obtain the last equality. By Proposition 2.1, Therefore, using Fubini's theorem, it follows that where we used the change of variable y = F − 0 (α) and (A.2) to obtain the last equality above. Given that E X * j = 1 n n i=1 E X * i,j , the first statement of Lemma A.1 follows from equality (A.4). Now, let us prove the second statement of Lemma A.1. Thanks to formula (A.3) for the density of X * i,j , one has that, for each 1 ≤ j ≤ p and 1 ≤ i ≤ n, where, again, we use the change of variable α = F i (x), and Fubini's theorem to obtain the last equality. Similarly, from (A.2) it follows that, for each 1 ≤ j ≤ p, Hence, using equalities (A.4), (A.5) and (A.6), and the fact that for each 1 ≤ i ≤ n, we obtain where, for the last inequalities, we used that E F − = F − 0 by Proposition 2.1. Therefore, from the above equality, one finally obtains which completes the proof of Lemma A.1.
A.2 Proof of Theorem 3.1 By Definition 2.1 of the Wasserstein distance, and sinceν n,p = 1 p p j=1 δX * j , it follows by using Fubini's theorem that From Lemma A.1, one has that E X * j = E Y * j . Therefore, by combining (A.7) with (A.1), we obtain where the last equalities also follow from Lemma A.1 and (A.1), which completes the proof of Theorem 3.1.

A.3 Proof of Theorem 3.2
We recall that ν ⊕ n denotes the measure with quantile function given by equation (2.1). By the triangle inequality, we have that Thanks to Definition 2.1 of the Wasserstein distance, it follows by Fubini's theorem that By Assumption 2.3, one has that E F − i (α) = F − 0 (α) for any 1 ≤ i ≤ n, and thus, by independence of the random variables F − i (α), one obtains Hence, by (A.10) and the inequality Now, let us remark that

28
where F − ν i denotes the quantile function of the measureν i = 1 p i p i j=1 δ X i,j for each 1 ≤ i ≤ n, and · denotes the usual norm in L 2 ([0, 1], dx). Hence, the above inequality leads to the following upper bound Therefore, Theorem 3.2 follows from inequality (A.9) combined with (A.11) and (A.12), which completes its proof.
A.4 Proof of Theorem 3.3 The proof follows the same lines as the proof of Theorem 3.2. By the triangle inequality, we have that where ν ⊕ n is the measure with quantile function given by equation (2.1). The expectation of the second term in the right-hand size of inequality (A.13) is controlled by inequality (A.11). Then, to control the first term, it suffices to remark that denotes the quantile function of the measureν h i i defined in (3.15), and · denotes the usual norm in L 2 ([0, 1], dx). Therefore, one has that j=1 δ X i,j for each 1 ≤ i ≤ n. Hence, the above inequalities lead to the following upper bound Finally, by applying Lemma 3.1 and inequality (3.10), we obtain that (A.14) Therefore, Theorem 3.3 follows from inequality (A.13) combined with (A.11) and (A.14), which completes the proof.
A.5 Proof of Theorem 3.4 Let A > 0 and σ > 0. To derive Theorem 3.4, we follow the classical scheme in nonparametric statistics to obtain optimal rates of convergence (see Chapter 2 in [Tsy09]). To this end we introduce appropriate random measures in W 2 (Ω), satisfying the deformable model defined in Section 2.3, that will serve as the basic hypotheses to obtain a lower bound. Let m (1) and m (2) be two real numbers such that where C is a positive constant to be specified later on. For k = 1, 2, we let a (k) be independent Gaussian random variables with E a (k) = m (k) and Var(a (k) ) = γ 2 with γ = min(A 1/2 , σ). We also let H (k) denote the hypothesis that the data are sampled according the following deformable model: where a i,j 's are iid random variables sampled from the Gaussian distribution with zero mean and variance γ 2 , that are independent of the a where e i is the vector in R p i with all entries equal to one, the notation Var (X) denotes the covariance matrix of a random vector X, and I i is the identity p i × p i matrix. For each k = 1, 2, if we denote by ν in W 2 (R) of the random measure ν (k) is the Gaussian distribution with mean m (k) and variance γ 2 , and that Var a (k) dα = γ 2 ≤ σ 2 .
Then, for k = 1, 2, we let P (k) be the probability measure of the data in model (A.16) under the hypothesis H (k) . From our remark above, one has that P (k) is the product of n Gaussian measures P (k) i on R p i with mean and covariance given by (A.17) for 1 ≤ i ≤ n. Hence, the Kullback divergence K P (1) , P (2) between P (1) and P (2) can be decomposed as follows where the last inequality follows from (A.15) and the fact that γ 2 = min(A, σ 2 ).