A novel notion of barycenter for probability distributions based on optimal weak mass transport
Year: 2022

Abstract:
We introduce weak barycenters of a family of probability distributions, based on the recently developed notion of optimal weak transport of mass [26], [9]. We provide a theoretical analysis of this object and discuss its interpretation in the light of convex ordering between probability measures. In particular, we show that, rather than averaging the input distributions in a geometric way (as the Wasserstein barycenter based on classic optimal transport does) weak barycenters extract common geometric information shared by all the input distributions, encoded as a latent random variable that underlies all of them. We also provide an iterative algorithm to compute a weak barycenter for a finite family of input distributions, and a stochastic algorithm that computes them for arbitrary populations of laws. The latter approach is particularly well suited for the streaming setting, i.e., when distributions are observed sequentially. The notion of weak barycenter and our approaches to compute it are illustrated on synthetic examples, validated on 2D real-world data and compared to standard Wasserstein barycenters.

heading: 1 Introduction
text: Optimal transport (OT) [41] has had a tremendous impact in the machine learning (ML) community recently, as it provides meaningful and implementable distances between probability distributions [34], thus advancing many aspects in the field, see e.g. [7, 43, 33]. The space of probability measures on Rd with finite second moment can be metrised with the Wasserstein-2 distance, the computation of which amounts to finding a transport plan that minimises the quadratic average cost of transporting mass from a source probability measure onto a target one. In this context, a natural method for averaging a finite family of probability measures is to compute their Fréchet mean, with respect to the Wasserstein-2 distance, which corresponds to the Wasserstein barycenter introduced in [1].
The goal of the present work is to explore theoretical features and potential applications to ML of barycenters of probability measures analogously defined in terms of optimal weak transport (OWT, see [26]) or more precisely quadratic barycentric transport costs. In a nutshell, for a source measure µ and a target measure ν, the OWT problem aims to transport mass so that the conditional spatial mean of target support points y, given their source support points x, is close to x in average. This amounts to finding an intermediate measure η, possibly more concentrated than ν in the sense of convex ordering of probability measures, which is close to µ with respect to the Wasserstein-2 distance.
The main motivation of our work is to investigate the effect and meaning of combining a family of probability measures using OWT instead of OT. To that end, we will define the weak barycenter of
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
this family through an optimisation problem, and discuss some of its properties. Importantly, we will see that, rather than averaging the input distributions in a metric sense, solving a weak barycenter problem corresponds to finding probability measures that encode geometric or shape information shared across all of them. In fact, the weak barycenter problem will be interpreted as finding a latent random variable common to all the input distributions. Implications of this latent variable interpretation, in terms of robustness to outliers, will also be drawn in our work.
A second motivation for our work is to develop and implement computational methods for weak barycenters, capitalising on the fact that the optimal weak coupling between any pair of distributions, with finite second moments, is always realised by a unique optimal map. This property is in sharp contrast to standard OT, where the absolute continuity with respect to the Lebesgue measure of the source or target measure is typically needed to grant the existence and uniqueness of a map—the so-called Monge map— realising the optimal coupling between them. This map is often required in different ways to compute Wasserstein barycenters (see [5], [44] or [35]).
Similarly to the Wasserstein barycenter problem, we will develop a fixed-point formulation of the weak barycenter problem, based on OWT plans. This allows us to, following [5, 44], construct an iterative procedure to compute a weak barycenter for a finite family of distributions and analyse its convergence properties. We will also define and study the so-called weak population barycenters, that deal with a population of probability measures distributed according to a given law Q supported on the Wasserstein-2 space, as in [30] for the OT case. Extending ideas from [10], we will then propose an iterative stochastic algorithm for online computation of the weak population barycenter, from a stream of probability measures sampled from Q. We will then provide numerical simulations using this proposed method, in order to illustrate the geometric meaning of the weak barycenter, and we will compare it with related objects obtained with standard OT or its entropy-regularised counterpart.
Organisation of the paper. Sec. 2 documents the background on OWT and the assumptions underlying our work. Sec. 3 analyses the weak barycenter problem, interprets it in the light of convex ordering and a latent variable model and addresses the case of an infinite population of distributions. Sec. 4 introduces two algorithms for computing the weak barycenter in the finite or population settings. Sec. 5 and 6 present the experimental setting and validation of our proposal respectively. Lastly, Sec. 7 discusses our findings and future research questions. The Appendix contains all the proofs, additional details or our simulations and the code of our experiments.

heading: 2 Background : optimal weak transport and Wasserstein barycenter
text: The optimal transport (OT) problem [42] aims to find the lowest cost to transfer the mass from one probability measure onto another. Therefore, OT is a natural way to compare two probability distributions in terms of their geometric information. In particular, the Wasserstein-p distance Wp, associated with the Euclidean cost in Rd, metrises the space Pp(Rd) of probability measures on Rd with finite p-moment. Precisely, for µ, ν ∈ Pp(Rd),
Wp(µ, ν) =
( min
π∈Π(µ,ν) ∫ Rd×Rd ‖x− y‖pdπ(x, y) )1/p , (1)
where π is a transport plan between µ and ν, that is, an element of the set Π(µ, ν) of probability measures on the product space Rd × Rd with marginals µ and ν. For p = 2 and µ absolutely continuous (a.c.), the unique optimal plan is concentrated on the graph of a measurable map called Monge map such that ν = T#µ, see eq. (14) in Appendix A.1.
Optimal weak transport. We consider here the optimal weak transport (OWT) problem introduced in [26] and in particular the special case of barycentric transport costs. The OWT problem is then defined for µ, ν ∈ P2(Rd) by
V (µ|ν) = inf π∈Π(µ,ν) ∫ Rd ‖x− ∫ Rd ydπx(y)‖2dµ(x), (2)
where πx is the disintegration of the transport plan π with respect to the first marginal µ, i.e. π(dxdy) = πx(dy)µ(dx). As our work strongly leans on OWT theory, we recall in Appendix A.2, Th. 6, that V is continuous with respect to the Wasserstein metric [9]. Additionally, the two following results from [8] (stated for our specific setting) lay the ground for our proposed weak barycenters.
Theorem 1 ([8], Theorem 1.2). The problem (2) admits a unique minimiser.
This first result strongly differs from the classical OT setting, for which the uniqueness of an optimal transport plan is not guaranteed for arbitrary measures. The optimisation problem in Eq. (2) can also be reformulated thanks to the Brenier-Strassen theorem [25], [8], through the notion of convex ordering. We denote by η ≤c ν the convex order of measures, meaning that ∫ φdη ≤ ∫ φdν for any convex function φ that is nonnegative or integrable with respect to η + ν. By Strassen’s theorem [40], two distributions are in convex order if and only if there exists a martingale coupling between them. The following theorem is a generalisation of the result originally proved in [25], Th. 1.2.
Theorem 2 ([8], Theorem 1.4). Let µ ∈ P2(Rd) and ν ∈ P1(Rd). There exists a unique η∗ ≤c ν such that
W 22 (µ, η ∗) = inf η≤cν W 22 (µ, η) = V (µ|ν). (3)
Moreover, there exists a convex function ψ : Rd → R of classC1 with∇ψ being 1-Lipschitz, such that ∇ψ#µ = η∗. Finally, the optimal coupling πµ,ν ∈ Π(µ, ν) verifies ∫ ydπµ,νx (y) = ∇ψ(x) µ-a.s.
The measurable map, or barycentric projection, Sνµ(x) := ∫ Rd ydπ µ,ν x (y) associated to the plan πµ,ν achieving the minimum in Eq. (2) is consequently uniquely defined and will be called optimal barycentric projection. From this notation, we can write the OWT cost in terms of an OT cost according to V (µ|ν) = W 22 (µ, Sνµ#µ). We emphasise that Sνµ is directly related to the optimisation problem (2), whereas applied works such as [38, 36] make use of a barycentric projection constructed from a transport plan solving an OT problem between µ and ν (often regularised) as a substitute for the Monge map, which may not exist (more details on Sνµ are displayed in Appendix A.3).
Last, let us note that OWT is somehow also related to the martingale OT problem developed in the stochastic finance community [13, 2, 27], which puts the focus on the optimal transfer of mass between distributions assumed to be in convex order themselves.
Wasserstein barycenter. The classical Wasserstein barycenter problem for a set of probability measures ν1, . . . , νn ∈ P2(Rd) with weights λ1, . . . , λn in the simplex (i.e. λi ≥ 0 and ∑n i=1 λi = 1) is defined [1] by
arg min µ∈P2(Rd) n∑ i=1 λiW 2 2 (µ, νi). (4)
The Wasserstein barycenter has been extensively studied both theoretically and numerically [30, 44, 5, 15]. Regarding the numerical part, [39] focuses on the computation of Wasserstein barycenters for a fixed number of measures and a stream of observations per measure; additionally, [32] proposed an entropy-regularised alternative via stochastic optimisation for computing the Wasserstein barycenter of a.c. distributions only from observations. Constrained by their assumption of a.c., [44] computes the Wasserstein barycenter by smoothing the observed empirical distributions. Furthermore, [21] compares the complexity of both the sample Wasserstein barycenter and a stochastic approximation to estimate a population barycenter (discrete measures and entropic regularisation). Finally, the authors of [3] recently proposed an algorithm to compute the barycenters in polynomial time.

heading: 3 Optimal weak transport barycenters and latent variable interpretation

heading: 3.1 Definition and basic properties
text: In a similar fashion, based on the weak transport cost in Eq. (2), we propose the following variant:
Definition 1. The set of weak barycenters of a finite family of measures {νi}i=1,...,n ∈ P2(Rd) with weights {λi}i=1,...,n in the simplex is defined as
arg min µ∈P2(Rd) n∑ i=1 λiV (µ|νi). (5)
Thus, a weak barycenter averages, with respect to the Wasserstein metric, an optimally chosen set of probability measures {η1, . . . , ηn} which are more concentrated than the corresponding νi, in the sense that ηi ≤c νi for each 1 ≤ i ≤ n. The existence of a solution is established as follows:
Proposition 1. The weak barycenter problem in Eq. (5) admits a minimiser µ ∈ P2(Rd).
See Sec. B of the Appendix for the proof of the above Proposition (which relies on Prokhorov’s theorem) and all the proofs for this Section. Uniqueness is in general not granted: we next show that the set of solutions is indeed an interval, with respect to the partial order of convex ordering of probability measures.
In the following, we denote by X and Yi random variables with respective laws µ and νi, for 1 ≤ i ≤ n, and δa the Dirac measure supported on a ∈ Rd. Lemma 1. If µ is a weak barycenter of {νi}i=1...,n and µ′ ≤c µ, then µ′ also is a weak barycenter. In particular, the Dirac measure supported on Eµ(X) is always a weak barycenter. Moreover, a Dirac distribution δω̄ is a weak barycenter if and only if ω̄ = ∑n i=1 λiEνi(Yi).
A consequence of the above lemma is that for any weak barycenter µ,
Eµ(X) = n∑ i=1 λiEνi(Yi) , (6)
and the value of the weak barycenter problem is given by
inf µ∈P2(Rd) n∑ i=1 λiV (µ|νi) = n∑ i=1 λi‖E(Yi)‖2 − ‖ n∑ i=1 λiE(Yi)‖2. (7)
We can also derive the following characterisation on the set of weak barycenters: Proposition 2. A measure µ ∈ P(Rd) is a weak barycenter of {νi}i=1...,n if and only if its mean satisfies (6) and µ̂ ≤c ν̂i holds for all 1 ≤ i ≤ n, where ν̂ denotes the centered version of a law ν.
For instance, in the case of one dimensional Gaussian distributions νi = N (m,σ2i ), the set of weak barycenters includes {µ = N (m,σ2) | 0 ≤ σ2 ≤ min1≤i≤n σ2i }. A natural question is whether a "maximal” weak barycenter exists, in the sense of convex ordering (up to translation by the mean). For d = 1, the answer is affirmative. When the means E(Yi) are equal, this follows from the complete lattice property of the set of probability measures with respect to the convex ordering (see [29]); the general case can then be reduced to the latter using Proposition 2. For d ≥ 2, this property is in general not true and the answer depends on the family {νi}i=1...,n. In the particular case of a.c. input measures, we can bound the distance between the Wasserstein and weak barycenters by the variances of the distributions (νi)1≤i≤n. The barycenters are then closer the more concentrated each νi is. Lemma 2. Let ν1, . . . , νn ∈ P2(Rd) be a.c., at least one of them with bounded density. Let µ̄ and µ̃ respectively denote the weak and the Wasserstein barycenters. Then
W 22 (µ̄, µ̃) ≤ 2 n∑ i=1 λi ( E‖Yi‖2 − ‖EYi‖2 ) .

heading: 3.2 Weak barycenters as latent variables
text: The weak barycenter encodes common geometric information present in all the input measures considered, therefore, it can be intuitively and rigorously interpreted as being the distribution of a latent variable underlying the realisations of random variables of laws νi for all 1 ≤ i ≤ n. Theorem 3. Let µ be a weak barycenter of {νi}i=1...,n. Then, for each 1 ≤ i ≤ n, a random variable Yi ∼ νi can be realised as
Yi = X + (EYi − EX) + Ȳi, where X ∼ µ and Ȳi = Yi − E(Yi|X) is centered conditionally on X . Moreover, one has Sνµ(X) = X + (EYi − EX) for all i = 1, . . . , n. Finally, we have E(Yi − EYi|X − EX) = X − EX or, equivalently, µ̂ ≤c ν̂i, with µ̂ and ν̂i the laws of X − EX and Yi − EYi respectively.
That is to say, each Yi ∼ νi can be realised by sampling a random variable X common to all i = 1, . . . , n and distributed according to the weak barycenter µ, translating that value by EYi − EX and adding a cluster-specific component Ȳi or idiosyncratic noise, centered conditionally on X .
Remark 1. The observations of each class (i.e. input measure) can be interpreted as outliers with respect to the (translated) law of the weak barycenter, which are statistically different and are thus left aside of its support. This way, the weak barycenter is robust to outliers, as it tends to discard them, by construction. Furthermore, this "robustness" property results in the stability of weak barycenter upon perturbation of a class with larger noise (or more scattered, outlying values). More precisely, if a class is corrupted in such a way that their observations result in a stochastically larger distribution than the original one, a weak barycenter computed in terms of the original (stochastically smaller) class will still be a weak barycenter in the new corrupted setting. An intuitive and simple way to illustrate this point follows by considering a weak barycenter µ of a one dimensional and centered family of input distributions {νi}i=1,...,n. By Proposition 2, µmust verify µ ≤c νi for all i = 1, . . . , n. In particular, from Theorem 3.A.1. in [37], we have that ∫∞ x P(X > u)du ≤ ∫∞ x
P(Yi > u)du for all x ∈ R, where X ∼ µ and Yi ∼ νi. Therefore, µ is likely to avoid outliers. Another supportive intuition in terms of robustness is that a maximal weak barycenter would be one that includes the most possible points of all classes (or distributions) in its support (all this, after re-centering) and leaves out only "outliers". A non-maximal weak barycenter is then more conservative, meaning that it counts on fewer points and leaves out more possible outliers.

heading: 3.3 Extension for the population barycenter
text: The population Wasserstein barycenter introduced in [30] and [4] extends the definition of Wasserstein barycenter for an infinite number of measures. This formulation is particularly relevant for the construction of an iterative algorithm to compute the barycenter for the streaming case, that is, when the measures are received online. The proofs are reported in Section C of the Appendix.
Let us consider a probability measure Q ∈ P2(P2(Rd)), meaning that Q is supported on a set of measures with finite moments of order 2, such that for some (and thus all) µ ∈ P2(Rd), we have that∫ P2(Rd)W 2 2 (µ, ν)dQ(ν) <∞.
Definition 2. We define the set of weak population barycenters of a distribution Q ∈ P2(P2(Rd)) as
arg min µ∈P2(Rd) ∫ P2(Rd) V (µ|ν)dQ(ν). (8)
The following lemma guarantees that the map (x, ν) 7→ Sνµ(x) appearing in Eq. (8) through V (µ|ν) =∫ ‖x− Sνµ(x)‖2dµ(x) is well defined.
Lemma 3. The function (µ, ν) ∈ (P2(Rd))2 7→ πµ,ν ∈ P2(R2d) mapping (µ, ν) to the unique optimal plan πµ,ν realising V (µ|ν) in Eq. (2) is continuous. As a consequence, for each µ ∈ P2(Rd) the function (x, ν) ∈ Rd × P2(Rd) 7→ Sνµ(x) is measurable.
Using similar arguments as those of Proposition 1 and the fact that any probability measure can be approximated by a sequence of probability measures with finite support, the following proposition confirms that the weak population barycenter problem is also well defined.
Proposition 3. The minimisation problem in Eq. (8) admits a solution.

heading: 4 Algorithms via fixed-point representations

heading: 4.1 Weak barycenter
text: For the Wasserstein barycenter problem in Eq. (4), the authors in [1] proved that if at least one of the measures ν1, . . . , νn is a.c., the Wasserstein barycenter is unique. Furthermore, if all the νi’s are a.c., and at least one of them has a bounded density, then the unique Wasserstein barycenter is also a.c. and verifies a fixed-point equation. This last property has been thoroughly studied by [5] and [44] and leveraged to compute an approximation of the barycenter via an iterative algorithm based on Monge maps, whose existence and uniqueness are guaranteed by the a.c. of the measures involved.
Akin to the fixed-point methodology in the classical Wasserstein scenario, we define an iterative procedure based on the barycentric projection computed in the optimal weak transport problem in Eq. (2), that is valid for arbitrary distributions. Therefore, we consider the following iterative rule for
probability measures ν1, . . . , νn ∈ P2(Rd):
µk+1 = G(µk), with G(µ) = ( n∑ i=1 λiS νi µ ) #µ, (9)
where for each i = 1, . . . , n, the optimal barycentric projection is given by Sνiµ : x 7→ ∫ ydπµ,νix (y), for πµ,νi ∈ Π(µ, νi) achieving the minimum in the OWT problem in Eq. (2). The proposed iterative procedure is presented in Algorithm 1.
A fundamental difference between the fixed-point computation of the Wasserstein barycenter [5] and a weak barycenter is that the optimal Monge map T νµ in the OT problem verifies T ν µ#µ = ν, whereas the pushforward measure Sνµ#µ in the OWT setting still depends on µ. We will then prove that the iterative algorithm in Eq. (9), based on the maps Sνiµ , admits converging subsequences. A convenient result is the continuity of the functional G in Eq. (9), which can be proven using Arzela-Ascoli theorem on a set of barycentric projections as well as the Skorohod’s representation theorem. Theorem 4. The function µ 7→ G(µ) defined in Eq. (9) is W2-continuous from P2(Rd) to P2(Rd).
Using an approach similar to [5] for the Wasserstein barycenter, we can state the following results for the proposed fixed-point procedure. Proposition 4. If µ is a weak-barycenter, that is a solution of problem (5), then G(µ) = µ i.e. x = ∑n i=1 λiS νi µ (x), µ(x)-a.s.
The inverse implication of Proposition 4 is not necessarily true, that is, some fixed points may not be weak barycenters. However, a Dirac delta δω, ω ∈ Rd, that meets the fixed-point condition δω = G(δω), is a weak barycenter (see Lemma 1). Proposition 5. Let (µk)k be the sequence defined by the iterative procedure µk+1 = G(µk) and starting from µ0 ∈ P2(Rd). Then (µk)k is tight and every converging subsequence must converge to a fixed point of G.
We observe that these results also hold for the classical Wasserstein barycenter of a.c. measures {νi}i=1...,n such that at least one of them has a bounded density. Moreover, the inverse implication, namely if µ is a fixed-point then it is a barycenter, is not straightforward even in the Wasserstein barycenter case, for which one considers the fixed-point equation given by µ = ( ∑n i=1 λiT νi µ )#µ, with T νi µ the Monge map verifying νi = T νi µ #µ. Indeed, [1] prove that if µ
checks x = ∑n i=1 λiT νi µ (x) for every x ∈ Rd, not only µ-almost everywhere, then µ is a Wasserstein barycenter. Also, [44, Theorem 2] provide additional conditions for this to be true by essentially invoking more smoothness on the distributions {νi}i=1...,n. Additionally, they only conjecture that under the same assumptions, the fixed-point is unique. Our method, however, includes arbitrary probability measures. Therefore, we do not expect to obtain similar results as in the Wasserstein barycenter case, for which smoothness is required.

heading: 4.2 Weak population barycenter
text: Based on [10], we construct a stochastic iterative algorithm for computing the weak population barycenter in Eq. (8). We clarify that [10] is constrained to probability measures Q supported on distributions that are a.c., whereas in our setting these distributions only need to belong to P2(Rd). Let us notice that our algorithms can be interpreted as geodesic gradient descent as in [10] and [17], however, OWT is not a metric and its potential geodesic structure is so far unknown. Therefore, the proposed algorithm only aims to mimic Riemmanian gradient descent. Our fixed-point result for the weak population barycenter problem is stated in the following Lemma: Lemma 4. If µ is a weak population barycenter of Q, then x = ∫ Sνµ(x)dQ(ν), µ(x)-a.s.
As in the finite case, the inverse implication is difficult to obtain. In particular, this has not been proven for the classical population Wasserstein barycenter in [10], where it boils down to prove the uniqueness of an absolutely continuous fixed point of µ 7→ ( ∫ T νµ dQ(ν))#µ, where T νµ is the Monge map between µ and ν. As explained in [10], the uniqueness of such fixed points has also been studied under some strong assumptions in [15] by considering parametric classes of random probability measures with compact support. This result is expected to be true by again invoking more smoothness on the distributions at hand. As our method focuses (in particular) on discrete probability measures,
the conditions under which the inverse implication holds are beyond the scope of our work. However, from the experimental results in Section 6, we believe our method presents practical advantages.
We next present an iterative scheme converging towards a distribution µ verifying the fixed-point equation in Lemma 4. This scheme is illustrated below in Algorithm 2. To prove its convergence, we will need a technical assumption on Q:
(A) There exists > 0 and R > 0 such that Q gives full measure to the set KQ := {µ ∈ P2(Rd) : ∫ |x|2+ dµ(x) ≤ R}.
Definition 3. Let µ0 ∈ KQ, νk i.i.d.∼ Q and γk > 0. We define the following iterative procedure: µk+1 = [ (1− γk)id + γkSν k
µk
] #µk, k ≥ 0, (10)
where Sν k
µk is the optimal barycentric projection between µk and νk and id is the identity operator.
The following standard conditions on the steps γk will also be assumed: ∞∑ k=1 γ2k <∞ and ∞∑ k=1 γk =∞, (11)
Theorem 5. Assume Conditions in Eq. (11), (A) and moreover that every measure verifying the fixed-point equation x = ∫ Sνµ(x)dQ(ν), µ(x)-a.s. is a weak barycenter. Then the sequence (µk)k in Eq. (10) is a.s. relatively compact w.r.t. W2 and every limit point is a weak barycenter.
The proof, provided in the supplementary material, is inspired by the standard Wasserstein barycenter case studied in [10], [35]. Assumption (A) grants that the sequence in Eq. (10) remains in some compact set, and can be replaced by more general conditions (see Remark 2 in Appendix).
Algorithm 1: Weak barycenter Input: distributions ν1, . . . , νn, # steps K; initialisation: µ0 = ν1; for k = 0, 1, . . . ,K do
for i = 1, 2, . . . , n do Solve the OWT problem between µk and νi to obtain πµk,νi ; Si = ∫ ydπµk,νix (y) end µk+1 = ( ∑n i=1 λiSi)#µk
end
Algorithm 2: Weak population barycenter Input: number of steps K; initialise distribution µ0 ∼ Q; for k = 0, 1, . . . ,K do
Sample νk ∼ Q; Update γk; Solve the OWT problem to obtain πµk,ν k
Sk = ∫ ydπµk,ν k
x (y); µk+1 = [(1− γk)id + γkSk] #µk;
end

heading: 5 Computational aspects
text: Setting and computation of OWTs. Both Algorithms 1 and 2 require the computation of the optimal barycentric projection associated to the OWT problem in Eq. (2). For two discrete measures µ =∑r i=1 aiδxi and ν = ∑m j=1 bjδyj , this boils down to solving the following quadratic programming problem
min π∈Rr×m { r∑ i=1 ai ∥∥∥xi − (πy a ) i ∥∥∥2 , πij ≥ 0, π1 = a, πT1 = b} , (12) which can be solved using a solver such as cvxpy. We also propose to solve the OWT problem in Eq. (12) with a proximal algorithm. The optimal barycentric projection is then constructed as πya . The details and examples are presented in Appendix E.1.
Comparison setting. In the next section, we compare our proposed computation for weak barycenters in Definition 2 (Algorithm 2) to the classic Wasserstein barycenter in particular for a stream of a.c. measures. Namely, we will run Algorithm 2 by, following [19, 36], replacing optimal barycentric
projections by the barycentric projections associated either to i) an optimal plan in the Kantorovich problem (1), or ii) the optimal Sinkhorn plan in the entropy regularised OT problem [18] given by
arg min π∈Π(µ,ν)
∫ ‖x− y‖2dπ(x, y) + εKL(π|µ⊗ ν), (13)
where KL denotes the Kullback-Leibler divergence. The associated barycenters will be referred to as OT barycenter and OT Sinkhorn barycenter respectively. The optimal plans for OT and regularised OT problem were computed using POT toolbox [23]. Notice that what we call OT barycenter (resp. OT Sinkhorn barycenter) is not solving a Wasserstein barycenter problem (resp. a regularised Wasserstein barycenter problem). Therefore, our method for barycentric computation differs from previous ones in the literature (see Section 2) in that it i) can process a stream of an unknown number of measures, ii) does not require the measures to be a.c., and iii) does not appeal to additional regularisation of the measures or the Wasserstein metric.

heading: 6 Experimental results
text: This section is devoted to the empirical validation of our proposal on both synthetic and real-world data. We first focused on Algorithm 2 since multiple algorithms to compute a Wasserstein barycenter for a fixed number of distributions are already available [19, 39]. We present two robustness to outliers experiments, then we validate our OWT barycenter on synthetic dataset and real-world ones. The overall conclusion of our experiments is that the weak barycenter is more likely to maintain the common (or shared) geometric features of the measures involved, as expected from Theorem 3. Additional experiments are presented in Appendix E.2, including the comparison of the energy for the computed weak barycenter in Algorithm 1 against the approximated optimal energy (using Eq.(7) and the plug-in estimator).

heading: 6.1 Robustness to outliers
text: OT’s sensitivity to outliers is a well-known problem that can be addressed e.g. with unbalanced OT [11]. We observed that OWT also allows to deal with outliers, which is coherent with the latent variable interpretation (see Remark 1). We illustrate this with two experiments. In Fig. 1 (left), we consider 50 sets of 20 − 30 observations from different 2D Gaussian measures, where each observation may be corrupted by random translations (Bernoulli p = 0.05) thus producing outliers. We show the resulting barycenters (dots), and barycenters without outliers (crosses) for Wasserstein barycenter (red) and weak barycenter (black), which shows robustness to outliers. In Fig. 1 (right), we consider two distributions supported on pair-of-ellipses, and 120 observations per distribution. Again, each observation may be corrupted by random translations (Bernoulli p = 0.05). The weak barycenter (black) shows a better preservation of the shapes than the Wasserstein barycenter (red), in particular, the red dots are more often located outside the ellipses.

heading: 6.2 Synthetic distributions
text: We implemented the proposed sequential computation of weak barycenters (Algorithm 2) on two examples of synthetic distributions: Gaussians and spirals. In each case, we sampled r observations from a random distribution at each step, and considered K steps (and thus K measures for each case).
2D Gaussians (r = 100 & K = 15). We considered distributions N (m, I), with m uniformly distributed on (−3, 3) × (−5,−5) and I the identity matrix. Fig. 2 (left) shows the empirical distributions together with the OWT and OT barycenters, the weak barycenter being the less spread out as expected. The three remaining plots illustrate the behaviour of the barycenters constructed as stated in Sec. 5. For a small regularisation parameter ε in Eq. (13), the OT and OT Sinkhorn barycenters are similar, however, as ε increases the OT Sinkhorn (OTS) barycenter becomes closer to the weak barycenter and thus even more concentrated, meaning that its samples tend to be closer to each other. Critically, for a very large ε, as the entropy tends to spread the mass in the regularised optimal plan, the associated barycentric projection will roughly move the mass to the spatial mean of the target distribution’s support.
20 10 0 10 20 30 30
20
10
0
10
20
10 0 10 20
15
10
5
0
5
10 OT Weak barycenter

heading: 6.3 Real-world dataset
text: MNIST dataset. We considered the well-known MNIST dataset [31] of grayscale images of handwritten digits. The images, of size 28 × 28 pixels, can be normalised and thus be interpreted as discrete probability measures supported on a two-dimensional grid of size 28× 28. We computed the barycenters with 30 steps of Algorithm 1 between two digits "8", that are noisy versions of the same digit with the aim to produce a more stable barycenter. To produce noisy data, we randomly (Bernoulli p = 0.1) move pixels of the prototype digit displayed in Fig. 4 (left). Fig. 4 (right) shows the barycenters using the OWT, OT, and entropic-OT (for ε = 1). This example illustrates how OWT reduces dispersion, so that weak barycenter provides the best uniformly spread results among the barycenters considered, with the two loops of the "8" well shaped.
Cytometry dataset. In biotechnology, flow cytometry is measured through intracellular markers of single cells in a biological sample with the objective of recognising common features across patients. However, these measurements are often disrupted by acquisition, rather than biological artefacts [28], thus hindering the identification of common features. To address this challenge, we compute the weak barycenter for the forward-scattered light (FSC) and side-scattered light (SSC) cell’s markers (using the flowStats package of Bioconductor [24]). We considered K = 15 patients and a variable number of cells per patient between 88 and 2185. Fig. 5 shows the 15 distributions (left) and the computed barycenters (right), thus confirming the ability of the weak barycenter to resolve the alignment of the dataset, while maintaining the expected diamond-shape. Moreover, the advantage of our proposed streaming procedure is fully exploited in this setting, since data from one or several patient can arrive
sequentially. Though this setting has been addressed with the Wasserstein barycenter in [14], also in Fig. 5, such method required a fixed grid to compute the barycenter unlike our method, thus revealing the computational simplicity of the weak barycenter.

heading: 7 Discussion
text: We have introduced the weak barycenter, which extracts common geometric information of probability measures on Rd based on optimal weak transport, and showed that it can be interpreted as a latent variable model. From the fixed-point formulation defined in terms of optimal weak transport maps, irrespective of the regularity assumptions on the measures involved, we developed practical computation via an iterative algorithm with guaranteed convergence. In particular, the proposed algorithms do not require a common grid on the sample space, when processing either observed data or samples from distributions. We have also proposed weak barycenters of a possibly infinite population of measures and developed a stochastic procedure for computing it in the streaming data regime where distributions are processes into the weak barycenter as they arrive. This has critical implications for continual-learning methods in the ML community.
Additional studies will focus on deepen the latent variable interpretation of weak barycenters, and its relationship to the aggregate information represented by the Wasserstein barycenter. Also, we identify two relevant theoretical aspects for further research: i) to exhibit general conditions on the family of input measures (or on the law of the population) for the existence of weak barycenters that are not Dirac masses; and ii) to provide conditions on those input measures for a “maximal” weak barycenter (in terms of convex ordering) to exist when d ≥ 2, among all the solutions of the weak barycenter problem (and, if possible, a way of constructing it by regularisation most probably). The statistical behaviour of the weak barycenter can also be investigated, in particular when constructed from large empirical random samples of given distributions. Lastly, the weak population barycenter could also be used to construct a predictive posterior in the context of Bayesian learning, as was done for Wasserstein barycenters in [35].
Acknowledgments. We thank Julio Backhoff-Veraguas for his valuable insight during the writing of this paper. This work was funded by ANID grants: AFB170001 & ACE210010 (CMM), FB0008 (AC3E), Fondecyt-Postdoctorado #3190926 and Fondecyt-Regular #1210606.

Review summaries:
Computing the barycenter of a finite set of probability measures is an important task in various statistical and machine learning problems, and is usually done by comparing the measures with a specific divergence. This divergence plays a significant role as it affects the properties and practical performance of the resulting barycenter, and should then be carefully chosen. Motivated by the attractive theoretical properties of the Wasserstein distance, which emerges from optimal transport (OT) theory, the Wasserstein barycenter has been introduced and extensively studied from a theoretical and practical point a view.  This paper introduces a novel barycenter of measures that is closely related to the Wasserstein barycenter, since the comparisons are made by solving OT problems: the authors introduce the "weak barycenter" (Definition 1), which is based on the recently introduced "optimal weak transport" problem, and study its theoretical properties. Specifically, they prove the existence of a solution to the optimization problem that defines the weak barycenter (Proposition 1), and they provide a characterization of the solutions (Lemma 1, Proposition 2). They also bound the deviation between the weak and Wasserstein barycenters (Lemma 2), and show that the weak barycenter is linked to a specific latent variable (Theorem 3).   Then, the authors introduce the "weak population barycenter" to average an infinite set of probability measures (Definition 2), which can be used for online settings. This barycenter correponds to an optimization problem, which is obtained by extending the definition of the weak barycenter and is well defined (Lemma 3, Proposition 3).   The authors prove that the weak and weak population barycenters are necessarily fixed points of specific operators (Proposition 4, Lemma 4), and they develop two iterative procedures (Algorithms 1 and 2, supplementary document) that are guaranteed to converge to such fixed points (Proposition 5, Theorem 5). The computational requirements of the two algorithms are investigated in Section 6. Finally, several experiments are conducted to compare the empirical performance of the weak barycenter against the Wasserstein barycenter (Section 7). 

The present manuscript introduces the natural extension of the concept of the Wasserstein Barycenter to the barycenter based on optimal weak mass transport.  Some analytical, statistical and computational properties of this object (the barycenter based on optimal weak mass transport) are discussed and shown. Some numerical experiments are shown too.

This paper studies barycenters of probability distributions using the weak transport distance, and establishes a number of results which mirror the theory of standard optimal transport barycenters that has been established over the past 10 years. In particular, they investigate existence and uniqueness, the first order optimality criterion, as well as iterative methods to compute weak barycenters. They propose a proximal sub-routine for this iterative method, and evaluate the performance of weak barycenters on a number of examples, comparing with standard barycenters and Sinkhorn barycenters.