Penalized Barycenters in the Wasserstein Space

∗

Jérémie Bigot, Elsa Cazelles & Nicolas Papadakis
Institut de Mathématiques de Bordeaux et CNRS (UMR 5251)
Université de Bordeaux

arXiv:1606.01025v2 [math.ST] 4 Jul 2017

July 6, 2017

Abstract
A regularization of Wasserstein barycenters for random measures supported on Rd is
introduced via convex penalization. The existence and uniqueness of such barycenters is
proved for a large class of penalization functions. A stability result of regularized barycenters in terms of Bregman distance associated to the penalization term is also given. This
allows to compare the case of data made of n probability measures with the more realistic
setting where we have only access to a dataset of random variables sampled from unknown
distributions. We also analyze the convergence of the regularized empirical barycenter of
a set of n iid random probability measures towards its population counterpart, and we
discuss its rate of convergence. This approach is shown to be appropriate for the statistical analysis of discrete or absolutely continuous random measures. In this setting,
we propose efficient algorithms for the computation of penalized Wasserstein barycenters.
This approach is finally illustrated with simulated and real data sets.

1

Introduction

This paper is concerned by the statistical analysis of data sets whose elements may be modeled as random probability measures supported on Rd that are either discrete or absolutely
continuous. In the one dimensional case (d = 1), examples can be found in neuroscience [40],
biodemographic and genomics studies [42], economics [29], as well as in biomedical imaging [34], while examples for the two dimensional case (d = 2) arise in spatial statistics for
replicated point processes [27].
In this paper, we focus on first-order statistics methods for the purpose of estimating,
from such data, a population mean measure or density function.
The notion of averaging depends on the metric that is chosen to compare elements in
a given data set. In this work, we consider the Wasserstein distance W2 associated to the
quadratic cost for the comparison of probability measures. Let Ω be a convex subset of Rd and
P2 (Ω) be the set of probability measures supported on Ω with finite order second moment.
As introduced in [1], an empirical Wasserstein barycenter ν̄n of set of n probability measures
ν1 , . . . , νn (not necessarily random) in P2 (Ω) is defined as a minimizer of
µ 7→

n
1X
W 2 (µ, νi ), over µ ∈ P2 (Ω).
n i=1 2

(1.1)

∗
This work has been carried out with financial support from the French State, managed by the French
National Research Agency (ANR) in the frame of the GOTMI project (ANR-16-CE33-0010-01).

However, depending on the data at hand, such a barycenter may be irregular (and not
even uniquely defined) which is typically the case when the νi ’s are discrete measures. As
an illustrative example, we consider a dataset of the locations of reported incidents of crime
(with the exception of murders) in Chicago in 2014 which is publicly available1 and that has
been recently studied in [27]. A sample from this dataset is displayed in Figure 1. For each
P i
month 1 ≤ i ≤ 12 of 2014, we let νi = p1i pj=1
δX i,j be the discrete measure whose support
is the set of locations of reported crimes for the i-th month. As argued in [27], the locations
of crimes X i,j may be considered as a sample from random intensity functions whose values
change from one day to another as crime opportunities are not uniformly distributed in time.
In this setting, we show that a regularization of the Wasserstein barycenter of the νi ’s (as
defined below) is a meaningful way to obtain a mean distribution of crimes locations which
is absolutely continuous on the area of the city.
March

100

September

100

90

90

80

80

70

70

60

60

50

50

40

40

30

30

20

20

10

10

0

0
0

20

40

60

0

20

40

60

Figure 1: Spatial locations of reported incidents of crime in Chicago for 2 months of 2014. In
order to protect the privacy of crime victims, addresses are shown at the block level only and
specific locations are not identified.
Definition 1.1. Let Pνn = n1 ni=1 δνi where δνi is the dirac distribution at νi . A regularized
empirical barycenter µγPνn of the discrete measure Pνn on P2 (Ω) is a minimizer of
P

n
1X
µ 7→
W 2 (µ, νi ) + γE(µ) over µ ∈ P2 (Ω),
n i=1 2

(1.2)

where E : P2 (Ω) → R+ is a smooth convex penalty function, and γ > 0 is a regularization
parameter.

1.1

Related results in the literature

Statistical inference using optimal transport For d = 1, tools from optimal are used
in [32] for the registration of multiple point processes which model repeated observations
organized in samples from independent subjects or experimental units. In [32], a consistent
1

https://data.cityofchicago.org/Public-Safety/Crimes-2001-to-present/ijzp-q8t2/data

2

estimator of the population Wasserstein barycenter of multiple point processes is proposed.
It is based on a smoothed empirical Wasserstein barycenter obtained by a preliminary kernel smoothing step of the observed point processes that is followed by quantile averaging.
Therefore, the way of constructing a smoothed Wasserstein barycenter in [32] differs from the
approach followed in this paper where regularization of the empirical Wasserstein barycenter
via a penalty function is considered.
The penalized problem (1.2) is motivated by the nonparametric method introduced in
[15] for the classical density estimation problem from discrete samples based on a variational
regularization approach to optimal transport with the Wasserstein distance as a data fidelity
term. However, the adaptation of this work for the regularization of Wasserstein barycenter
has not been considered so far.
Generalized notions of Wasserstein barycenters A detailed characterization of empirical Wasserstein barycenters in terms of existence, uniqueness and regularity is given in [1].
There exists also a link between Wasserstein barycenters and the multi-marginal problem in
optimal transport as studied in [1] and [33]. Recently, the notion of Wasserstein barycenter
has been generalized in [31] for random probability measures supported on a locally compact
geodesic space. The main contributions in [31] are the proofs of existence, uniqueness and
consistency of such barycenters. The case of probability measures supported on a Riemannian
manifolds is also studied in [28]. Trimmed barycenters in the Wasserstein space P2 (Rd ) have
been introduced in [3] for the purpose of combining informations from different experimental
units in a parallelized or distributed estimation setting.
For applications in image processing, a fast approximation of the Wasserstein distance
W2 between probability measures supported on Rd with d ≥ 2 has also been proposed in
[10, 35] using a sliced framework based one dimensional Wasserstein distances along radial
projections of the input measure.
However, in all these papers, incorporating regularization (through penalization) into the
computation of Wasserstein barycenters has not been considered, which is of interest when
the data are irregular probability measures.
Regularization of the transport map Alternatively, regularized barycenters may be
obtained by adding a convex regularization on optimal transport plans when computing the
Wasserstein distance between probability measures. This approach leads to the notion of
regularized transportation problems [8, 21, 23, 25], and it has recently gained popularity
in the literature on image processing and machine learning. Recent contributions include
the fast approach in [20] to compute smooth Wasserstein barycenters of discrete measures
via entropic regularization of the transport plan, and the so-called Sinkhorn’s algorithm. In
these works, such a regularization is motivated by the need to accelerate the computation of
the Wasserstein distance between probability measures supported on Rd with d ≥ 2.

1.2

Contributions and structure of the paper

The presentation of the main results in this paper is organized as follows.
- In Section 2, we introduce various definitions and notation, and we prove a key result
called subgradient’s inequality on which a large part of the developments in the paper
lean.
3

- In Section 3, we analyze the existence, uniqueness and stability of regularized Wasserstein barycenters that are solutions of (1.2) for various of penalty functions E and any
regularization parameter γ > 0.
- In Section 4, for the Bregman distance associated to the penalization term, we derive
convergence properties of empirical regularized barycenters toward their population
counterpart in the asymptotic setting where n tends to infinity and γ = γn is let going
to zero. In this context, we demonstrate that the bias term (as classically referred
to in nonparametric statistics) converges to zero when γ → 0 in Rd . We also show
(with additional regularity assumptions for d ≥ 2) that the variance term converges to
0 when limn→∞ γn2 n = +∞. We mainly focus on penalization functions that enforce
the Wasserstein barycenter to be an absolutely continuous (a.c.) measure with a smooth
probability density function (pdf). In this case, it is natural to use the Bregman distance
to asses the quality of estimation of the pdf of the population barycenter.
- To illustrate the benefits of regularized barycenters for data analysis, we propose to
use in Section 5 efficient minimization algorithms for the computation of regularized
barycenters as well as a selection strategy for the parameter γ. This approach is finally
applied to the statistical analysis of simulated and real data sets in P2 (R) and P2 (R2 ).
Finally, a brief overview of the concepts of Bregman divergence and subgradient are gathered in the Appendix A, the proofs in a technical Appendix B and Appendix C contains
algorithmic details.

2
2.1

Definitions, notation and first results
Wasserstein distance and Kantorovich’s duality

For Ω a convex subset of Rd , we denote by M(Ω) the space of bounded Radon measures on Ω.
We recall that P2 (Ω) is the set of probability measures over (Ω, B(Ω)) with finite second order
moment, where B(Ω) is the σ-algebra of Borel subsets of Ω. In particular, P2 (Ω) ⊂ M(Ω).
Definition 2.1. The Wasserstein distance W2 (µ, ν) is defined for µ, ν ∈ P2 (Ω) by
W22 (µ, ν)

Z Z

= inf
π

|x − y|2 dπ(x, y)

(2.1)

Ω Ω

where the infimum is taken over all probability measures π on the product space Ω × Ω with
respective marginals µ and ν.
The well known Kantorovich’s duality theorem leads to another formulation of the Wasserstein distance.
Theorem 2.2 (Kantorovich’s duality). For any µ, ν ∈ P2 (Ω), one has that
W22 (µ, ν)

Z

=

φ(x)dµ(x) +

sup
(φ,ψ)∈CW

Z

Ω

ψ(y)dν(y)

(2.2)

Ω

where CW is the set of all measurable functions (φ, ψ) ∈ L1 (µ) × L1 (ν) satisfying
φ(x) + ψ(y) ≤ |x − y|2 ,

(2.3)

for µ-almost every x ∈ Ω and ν-almost every y ∈ Ω. A couple (φ, ψ) ∈ CW that attains the
supremum is called an optimal couple for (µ, ν).
4

For a detailed presentation of the Wasserstein distance and Kantorovich’s duality, we refer
to [38, 39]. For µ, ν ∈ P2 (Ω), we
denote by π µ,ν an optimal transport plan, that is a solution
RR
2
of (2.1) satisfying W2 (µ, ν) =
|x − y|2 dπ µ,ν (x, y). Likewise a pair (φµ,ν , ψ µ,ν ) ∈ L1 (dµ) ×
L1 (dν) achieving the supremum in (2.2) (under the constraint φµ,ν (x) + ψ µ,ν (y) ≤ |x − y|2 )
stands for the optimal couple in the Kantorovich duality formulation of the Wasserstein
distance between µ and ν.
For the sake of completeness, we also introduce the functional space Y := {g ∈ C(Ω) :
x 7→ g(x)/(1 + |x|2 ) is bounded endowed with the norm kgkY = supx∈Ω |g(x)|/(1 + |x|2 )
where C(Ω) is the space of continuousnfunctions from Ω to R. We finally denote
as Z the
o
closed subspace of Y given by Z = g ∈ C(Ω) : lim|x|→∞ g(x)/(1 + |x|2 ) = 0 . The space
M(Ω) of bounded Radon measures is identified with the dual of C0 (Ω) (space of continuous
functions that vanish at infinity). Finally, we denote by L1 (µ) the set of integrable functions
g : Ω → R with respect to the measure µ.

2.2

Regularized barycenters of a random measure

A probability measure ν in P2 (Ω) is said to be random if it is sampled from a distribution
P on (P2 (Ω), B (P2 (Ω)), where B (P2 (Ω)) is the Borel σ-algebra generated by the topology
induced by the distance W2 . Throughout the paper, we use bold symbols ν, X, f , . . . to
denote random objects. Then, we introduce a Wasserstein distance between distributions
of random measures (see [31] and [3] for similar concepts), and the notion of Wasserstein
barycenter of a random probability measure ν.
Definition 2.3. Let W2 (P2 (Ω)) be the space of distributions P on P2 (Ω) (endowed with the
Wasserstein distance W2 ) such that for some (thus for every) µ ∈ P2 (Ω)
W22 (δµ , P) := EP (W22 (µ, ν)) =

Z
P2 (Ω)

W22 (µ, ν)dP(ν) < +∞

where ν ∈ P2 (Ω) is a random measure with distribution P. The Wasserstein barycenter of a
random probability measure with distribution P ∈ W2 (P2 (Ω)) is defined as a minimizer of
µ ∈ P2 (Ω) 7→

W22 (δµ , P)

Z

=
P2 (Ω)

W22 (µ, ν)dP(ν) over µ ∈ P2 (Ω)

(2.4)

where δµ denotes the Dirac measure at the point µ.
Thanks to the results in [31], there exists a minimizer of (2.4), and thus the notion of
Wasserstein barycenter of a random probability measure is well defined. Throughout the
paper, the following assumptions are made on the regularizing function E.
Assumption 2.1. Let E : P2 (Ω) → R+ be a proper, lower semicontinuous and differentiable
function that is strictly convex on its domain
D(E) = {µ ∈ P2 (Ω) such that E(µ) < +∞} .

(2.5)

Regularized barycenters of a random measure are then defined as follows.
Definition 2.4. For a distribution P ∈ W2 (P2 (Ω)) and a regularization parameter γ ≥ 0,
the functional JPγ : P2 (Ω) → R+ is defined as
JPγ (µ)

Z

=
P2 (Ω)

W22 (µ, ν)dP(ν) + γE(µ), µ ∈ P2 (Ω).
5

(2.6)

If it exists, a minimizer µγP of JPγ is called a regularized Wasserstein barycenter of the random
measure ν with distribution P. In particular, if P is the discrete (or empirical) measure
P
defined by P = Pn = n1 ni=1 δνi where νi ∈ P2 (Ω) for i = 1, . . . , n, then JPγ becomes
JPγn (µ) =

n
1X
W 2 (µ, νi ) + γE(µ).
n i=1 2

(2.7)

Note that JPγ is strictly convex on D(E) by Assumption 2.1.
A typical example of a regularizing function satisfying Assumption 2.1 is the negative
entropy [15] (see e.g. Lemma 1.4.3 in [24]) defined as

E(µ) =

 R

 Rd f (x) log(f (x))dx,

 +∞

if µ admits a density f with respect to
the Lebesgue measure on Ω,
otherwise.

which enforces the barycenter to be a.c. with respect to the Lebesgue measure on Rd .

2.3

Subgradient’s inequality

In order to analyze the stability of the minimizers of JPγ with respect to the distribution P,
the notion of Bregman divergence together with the concept of subgradient will be needed.
An overview of these tools is presented in an C.
In our case, since E is supposed differentiable, for µ ∈ M(Ω) we have that ∂E(µ) =
∇E(µ). Then for µ, ν ∈ M(Ω), the (symmetric) Bregman distance is defined by
dE (µ, ν) = h∇E(µ) − ∇E(ν), µ − νi

(2.8)

R

where the linear form is understood as hf, µi = Ω f (x)dµ(x), for µ ∈ M(Ω) and f ∈ Cb (Ω)
the space of continuous bounded function from Ω to R.
Theorem 2.5 (Subgradient’s inequality). Let ν be a probability measure in P2 (Ω), and define
the functional
J : µ ∈ P2 (Ω) 7→ W22 (µ, ν) + γE(µ)
where E : P2 (Ω) → R is a proper, differentiable and convex function, and γ ≥ 0. If µ ∈ P2 (Ω)
minimizes J, then there exists φµ,ν ∈ L1 (µ) and ψ ∈ L1 (ν) verifying φµ,ν (x) + ψ(y) ≤ |x − y|2
for all x, y in Ω such that (φµ,ν , ψ) is an optimal couple of the Kantorovich’s dual problem
associated to µ, ν (Theorem 2.2). Moreover, for all η ∈ P2 (Ω),
γ h∇E(µ), µ − ηi ≤ −

Z

φµ,ν d(µ − η).

(2.9)

The proof of Theorem 2.5 is based on the two following lemmas whose proof can be found
in the Appendix B.1.
Lemma 2.6. The two following assertions are equivalent:
1. µ ∈ P2 (Ω) minimizes J over P2 (Ω),
2. there exists a subgradient φ ∈ ∂J(µ) such that hφ, η − µi ≥ 0 for all η ∈ P2 (Ω).
6

Lemma 2.7. Let µ ∈ P2 (Ω) and φ ∈ L1 (µ), then
φ ∈ ∂1 W22 (µ, ν) ⇔ ∃ ψ ∈ L1 (ν) / φ(x) + ψ(y) ≤ |x − y|2
and W22 (µ, ν) =
W22 (·, ν) at µ.

R

φdµ +

R

ψdν where ∂1 W22 (µ, ν) denote the subdifferential of the function

Proof of Theorem 2.5. Let µ ∈ P2 (Ω) be a minimizer of J. From Lemma 2.6, we know that
there exists φ a subgradient of J in µ such that hφ, η − µi ≥ 0 for all η ∈ P2 (Ω). Since
ζ 7→ E(ζ) is convex differentiable, ζ 7→ W22 (ζ, ν) is a continuous convex function and µ
minimizes J, we have by the subdifferential of the sum (Theorem 4.10 in [19]) that ∂J(µ) =
∂1 W22 (µ, ν) + γ ∇E(µ). This implies that all φ ∈ ∂J(µ) is written φ = φ1 + φ2 with φ1 = φµ,ν
optimal for the couple (µ, ν) (by Lemma 2.7) and φ2 = γ∇E(µ). Finally,
we have that
R
hφµ,ν + γ∇E(µ), η − µi ≥ 0 for all η ∈ P2 (Ω) that is γ h∇E(µ), µ − ηi ≤ − φµ,ν d(µ − η), ∀η ∈
P2 (Ω).

3

Existence, uniqueness and stability of regularized barycenters

In this section, we present some properties of the minimizers of the functional JPγ (see Definition 2.4) in terms of existence and stability.

3.1

Existence and uniqueness

In a first part, we state that the minimization problem (2.6) admits a unique minimum in the
particular setting where P is a discrete distribution on P2 (Ω) that is we study the problem
min

µ∈P2 (Ω)

JPγn (µ)

Z

=

W22 (µ, ν)dPn (ν) + γE(µ) =

n
1X
W 2 (µ, νi ) + γE(µ)
n i=1 2

(3.1)

where Pn = n1 ni=1 δνi ∈ W2 (P2 (Ω)) with ν1 , . . . , νn measures in P2 (Ω). In a second part, we
prove the existence and uniqueness of (2.6) in a general case.
P

Theorem 3.1. Suppose that Assumption 2.1 holds and that γ > 0. Then, the functional JPγn
defined by (3.1) admits a unique minimizer on P2 (Ω) which belongs to the domain D(E) of
the regularizing function E, as defined in (2.5).
The proof of Theorem 3.1 is given in the Appendix B.2. Thanks to this result, one may
impose the regularized Wasserstein barycenter µγPn to
be a.c. with respect to Lebesgue measure
R
on Ω by choosing the negative entropy E(µ) = Ω log (dµ(x)) dµ(x) for the regularization
function E. For this choice, (3.1) becomes a problem of minimization over a set of pdf with
entropy regularization. Examples of the use of the negative entropy as a regularization term
are given in Section 5 on numerical experiments.
From Theorem 3.1, it is possible to prove the general case (whose proof is also given in
Appendix B.2).
Theorem 3.2. Suppose that Assumption 2.1 holds and that γ > 0. Then, the functional JPγ
defined by (2.6) admits a unique minimizer.
7

3.2

Stability

We now study the stability of the solution with respect to the symmetric Bregman distance
dE (2.8) associated to the differentiable function E. Let ν1 , . . . , νn ∈ P2 (Ω) and η1 , . . . , ηn ∈
P
P
P2 (Ω) . We denote by Pνn (resp. Pηn ) the discrete measure n1 ni=1 δνi (resp. n1 ni=1 δηi ) in
W2 (P2 (Ω)).
Theorem 3.3. Suppose that Ω is bounded. Let µν , µη ∈ P2 (Ω) with µν minimizing JPγνn and
µη minimizing JPγη defined by (3.1). Then, the symmetric Bregman distance associated to E
n
can be upper bounded as follows
dE (µν , µη ) ≤

n
X
2
inf
W2 (νi , ησ(i) ),
γn σ∈Sn i=1

(3.2)

where Sn is the permutation group of the set {1, . . . , n}.
The proof of Theorem 3.3 is given in Appendix B.3. To better interpret the upper
bound (3.2), we need the notion of Kantorovich transport distance TW2 on the metric space
(P2 (Ω), W2 ), see [38]. For P, Q ∈ W2 (P2 (Ω)) endowed with the Wasserstein distance W2 , we
have that
Z
TW2 (P, Q) := inf
W2 (µ, ν)dΠ(µ, ν),
Π

P2 (Ω)×P2 (Ω)

where the minimum is taken over all probability measures Π on the product space P2 (Ω) ×
P2 (Ω) with marginals P and Q. Since Pνn and Pηn are discrete probability measures supported
on P2 (Ω), it follows that the upper bound (3.2) in Theorem 3.3 can also be written as (by
Birkhoff’s theorem for bi-stochastic matrices, see e.g. [38])
dE (µν , µη ) ≤

2
TW (Pν , Pη ).
γ 2 n n

The above upper bound means that the Bregman distance between the regularized Wasserstein barycenters µν and µη is controlled by the Kantorovich transport distance between the
distributions Pνn and Pηn .

3.3

Discussion

Theorem 3.3 is of particular interest in the setting where the νi ’s and ηi ’s are discrete probP
P
ability measures on Rd . If we assume that νi = p1 pj=1 δX i,j and ηi = p1 pj=1 δY i,j where
(X i,j )1≤i≤n;1≤j≤p and (Y i,j )1≤i≤n;1≤j≤p are (possibly random) vectors in Rd , then by (3.2),




1/2

p
n
1 X

X
2

dE (µν , µη ) ≤
inf
inf
|X i,j − Y σ(i),λ(j) |2 

γn σ∈Sn i=1 λ∈Sp  p j=1

.

Theorem 3.3 is also useful to compare the computation of regularized Wasserstein barycenter between the case of data made of n a.c. probability measures ν1 , . . . , νn , with the more
realistic setting where we have only access to random variables X = (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 in

8

Rd sampled from the measure νi for each 1 ≤ i ≤ n. If we denote by ν pi = p1i
usual empirical measure associated to νi , it follows from inequality (3.2) that




E d2E µγPνn , µγX



≤

Ppi

j=1 δX i,j

the

n


4 X
2
E
W
(ν
,
ν
)
,
i
p
2
i
γ 2 n i=1

where µγX is the random density satisfying


µγX



pi
n
1 X
1X
2
W2 µ,
δX  + γE(µ).
= argmin
pi j=1 i,j
µ∈P2 (Ω) n i=1

This result allows to discuss the rate of convergence (for the symmetric squared Bregman
distance) of µγX to µγPνn as a function of the rate of convergence (for the squared Wasserstein
distance) of the empirical measure ν pi to νi for each 1 ≤ i ≤ n (in the asymptotic setting
where p = min1≤i≤n pi is let going to infinity).
As an illustrative example, in the one-dimensional case (that is d = 1), one may use
Theorem 5.1 in [9], to obtain that


E

W22 (νi , ν pi )



2
≤
J2 (νi ), with J2 (νi ) =
pi + 1

Z
Ω

Fi (x)(1 − Fi (x))
dx,
fi (x)

where fi is the pdf of νi , and Fi denotes its cumulative distribution function. Therefore,
provided that J2 (νi ) is finite for each 1 ≤ i ≤ n, one obtains the following rate of convergence
of µγX to µγPνn (for d = 1)


E

d2E



µγPνn , µγX



n
J2 (νi )
8 X
8
≤ 2
≤ 2
γ n i=1 pi + 1
γ

n
1X
J2 (νi ) p−1 .
n i=1

!

(3.3)

When the measures ν1 , . . . , νn are supported on Rd with d ≥ 2, we refer to [26] for further
results on the rate of convergence of an empirical
measure
in Wasserstein distance that may


γ
γ
be used to derive rates of convergence for dE µPνn , µX .

4

Convergence properties of regularized empirical barycenters

In this section, when Ω is a compact of Rd , we study the convergence of the regularized
Wasserstein barycenter of a set ν 1 , . . . , ν n of independent random measures sampled from
a distribution P towards a minimizer of JP0 , that is a population Wasserstein barycenter of
the probability distribution P ∈ W2 (P2 (Ω)). To this end, we first introduce and recall some
notation.
Definition 4.1. For ν 1 , . . . , ν n iid random measures in P2 (Ω) sampled from a distribution
P
P ∈ W2 (P2 (Ω)), we let Pn = n1 ni=1 δν i . Moreover, we use the notation
µγPn

∈ argmin
µ∈P2 (Ω)

JPγn (µ)

Z

µγP ∈ argmin JPγ (µ) =

Z

µ0P ∈ argmin JP0 (µ) =

Z

µ∈P2 (Ω)

µ∈P2 (Ω)

W22 (µ, ν)dPn (ν) + γE(µ)

=

(4.1)

W22 (µ, ν)dP(ν) + γE(µ)

(4.2)

W22 (µ, ν)dP(ν),

(4.3)

9

that will be respectively referred as to the empirical Wasserstein barycenter (4.1), the regularized population Wasserstein barycenter (4.2) and the population Wasserstein barycenter
(4.3).
In what follows, we obtain a rate of convergence in expected squared Bregman distance
between µγPn and µγP which depends on n and γ. A general result is first stated in Section 4.1.
Complementary results are given in Section 4.2 for d = 1 and in Section 4.3 for d ≥ 2. Then,
in Section 4.4, we prove the convergence in Bregman divergence (as γ → 0) of the regularized
population Wasserstein barycenter µγP towards µ0P .

4.1

Rate of convergence of µγPn towards µγP in symmetrized Bregman distance

To compute a rate of convergence between µγPn and µγP , we will need results from the empirical
process theory. Thus, we first introduce some notions borrowed from [37].
Definition 4.2. Let F = {f : U 7→ R} be a class of real-valued functions defined on a given
set U , endowed with a norm k·k. An envelope function of F is any function u 7→ F (u) such that
|f (u)| ≤ F (u) for every u ∈ U and f ∈ F. The minimal envelope function is u 7→ supf |f (u)|.
The covering number N (, F, k · k) is the minimum number of balls {kg − f k < } of radius 
and center g needed to cover the set F. The metric entropy is the logarithm of the covering
number. Finally, we define
I(δ, F) = sup
Q

Z δq
0

1 + log N (kF kL2 (Q) , F, k · kL2 (Q) )d

(4.4)

where the supremum is taken over all discrete probability measures Q supported on U with
1/2
R
kF kL2 (Q) = |F (u)|2 dQ(u)
> 0.
Theorem 4.3. If Ω is a compact of Rd , then one has that
E(d2E (µγPn , µγP )) ≤

CI(1, H)kHkL2 (P)
γ2n

(4.5)

where C is a positive constant, and
H = {hµ : ν ∈ P2 (Ω) 7→ W22 (µ, ν) ∈ R; µ ∈ P2 (Ω)}

(4.6)

is a class of functions defined on P2 (Ω) with envelope H.
The proof of Theorem 4.3 is given in the Appendix B.4. To complete this result, one
needs to prove that I(1, H) < ∞, which depends on the rate of of convergence of the metric
entropy towards infinity as  tends to zero.

4.2

The one-dimensional case

For probability measures ν 1 , . . . , ν n supported in the real line, we can prove that the righthand side of (4.5) is finite by using existing results on the notion of bracketing number defined
below.

10

Definition 4.4. Given two real-valued functions l and r, the bracket [l, r] is the set of all
functions f with l ≤ f ≤ r. An -bracket is a bracket [l, r] with kl − rk < . The bracketing
number N[] (, F, k · k) is the minimum number of -brackets needed to cover F.
Theorem
4.5. If Ω is a compact of R, then there exists a finite constant C > 0 such that

γ
2
E(dE µPn , µγP ) ≤ γC2 n .
Proof. In what follows, C denotes a universal constant whose value may change from line to
line. We define the envelope function H : ν ∈ P2 (Ω) 7→ sup {W2 (µ, ν); W22 (µ, ν)}. Since
µ∈P2 (Ω)

for hµ ∈ H we have
|hµ (ν)| ≤ 2

Z

|x|2 dµ(x) + 2

Z

|y|2 dν(y) ≤ 4δ(Ω) for all ν ∈ P2 (Ω)

where δ(Ω) = sup|x|2 , then for all Q ∈ W2 (P2 (Ω)),
x∈Ω

Z

kHkL2 (Q) =

1/2

2

H(ν) dQ(ν)



≤ 16δ(Ω)

2

1/2

Z

≤ 4δ(Ω).

dQ(ν)

Now, it remains to control the term I(1, H) in the upper bound (4.5). By the triangle reverse
inequality, we have
|hµ (ν) − hµ0 (ν)| = |W2 (ν, µ) − W2 (ν, µ0 )| (W2 (ν, µ) + W2 (ν, µ0 ))
≤ W2 (µ, µ0 ) 2H(ν).
Then, from Theorem 2.7.11 in [37], and since Theorem 4 in [30] allows us to bound the metric
entropy by the bracket entropy, we get
log N (kHkL2 (Q) , H, k · kL2 (Q) ) ≤ log N[] (kHkL2 (Q) , H, k · kL2 (Q) )
≤ log N (, P2 (Ω), W2 ) ≤ log N[] (, P2 (Ω), W2 ).

(4.7)

Also, for d = 1, we have
0

Z 1

W2 (µ, µ ) =
0

|Fµ− (t)

−

Fµ−0 (t)|2 dt

1/2

= kFµ− − Fµ−0 kL2 ([0,1])

(4.8)

where Fµ− is the quantile function of the cumulative distribution function Fµ of µ. We denote
by G = {Fµ− , µ ∈ P2 (Ω)} the class of quantile functions of probability measures µ in P2 (Ω),
which are monotonic functions. Moreover, we can observe that Fµ− : [0, 1] → [Fµ− (0), Fµ− (1)] ⊆
Ω, where Ω is a compact included in R. Hence, G is uniformly bounded, say by a constant
M > 0. Finally, by Theorem 2.7.5. of [37] concerning the bracket entropy of the class of
monotonic functions, we obtain that log N[] (, G, L2 [0, 1]) ≤ CM
 , for some constant C > 0.
Finally, from relations (4.7) and (4.8), we can deduce that
I(1, H) =sup
Q

Z 1q
0

1 + log N (kHkL2 (Q) , H, L2 (Q))d ≤

Z 1

s

1+
0

CM
d < ∞.


Therefore, when ν 1 , . . . , ν n are iid random measures with support included in a compact
2
interval Ω, it 
follows from
 Theorem 4.5 that if γ = γn is such that limn→∞ γn n = +∞ then
limn→∞ E(d2E µγPνn , µ0P ) = 0.
11

4.3

The Rd case with additional regularization

In the case d ≥ 2, the class of functions H defined in (4.6) is too large to control the metric
entropy so that I(1, H) is finite. To solve this issue, we impose more smoothness on the
regularized Wasserstein barycenter as follows.
We assume that Ω is a smooth and uniformly convex set, and we choose
f (x) log(f (x))dx + kf k2H k (Ω) , if f = dµ
dx and f > α,
+∞
otherwise.

( R

E(µ) =

Rd

(4.9)

where k · kH k (Ω) denotes the Sobolev norm associated to the L2 (Ω) space and α > 0 is
arbitrarily small. Then, the following result holds.
Theorem 4.6. Suppose that Ω is a compact and uniformly convex set with a C 1 boundary.
Assume that the penalty function E is given by (4.9)
 for some
 α > 0 and k > d − 1. Then,
γ
γ
2
there exists a finite contant C > 0 such that E(dE µPn , µP ) ≤ γC2 n .
Proof. Supposing that Ω has a C 1 boundary, we have by the Sobolev embedding theorem
that H k (Ω) is included in the Hölder space C m,β (Ω̄) for any integer m and β ∈]0, 1] satisfying
m + β = k − d/2. Hence, the densities of µγPn and µγP given by (4.1) and (4.2) belong to
C m,β (Ω̄).
Arguing as in the proof of Theorem 4.5, we have |hµ (ν) − hµ0 (ν)| ≤ W2 (µ, µ0 ) 2H(ν) and
kHkL2 (Q) < ∞, where H(ν) = sup {W2 (µ, ν); W22 (µ, ν)} where D(E) is defined by (2.5).
µ∈D(E)

Thus, instead of controlling the metric entropy N (kHkL2 (Q) , H, k · kL2 (Q) ), it is enough to
bound the metric entropy N (, D(E), W2 ) thanks to Theorem 2.7.11 in [37].
To this end, since µ, µ0 ∈ D(E) are a.c. measures, one has that
W2 (µ, µ0 ) ≤

Z

1/2

|T (x) − T 0 (x)|2 dx

where T #λd = µ and T 0 #λd = µ0 ,

Ω

with λd denoting the Lebesgue measure on Ω. Thanks to Theorem 3.3 in [22] on the regularity
of optimal maps (results initally due to Caffarelli, [16] and [17]), the coordinates of T and T 0
are C m+1,β (Ω̄) functions λd −a.e.. Thus, we can bound N (, D(E), W2 ) by the bracket entropy
P
N[] (, C m+1,β (Ω̄), L2 (Ω)) since |T (x)−T 0 (x)|2 = dj=1 |Tj (xj )−Tj0 (xj )|2 where Tj , Tj0 : Ω → R.
Now, by Corollary 2.7.4 in [37],
log N[] (, C

m+1,β

(Ω̄), L2 (Ω)) ≤ K

 V
1



for any V ≥ d/(m + 1). Hence, as soon as V /2 < 1 (for which the condition k > d − 1 is
sufficient if V = d/(m + 1)), the upper bound in (4.5) is finite for H = {hµ : ν ∈ P2 (Ω) 7→
W22 (µ, ν) ∈ R; µ ∈ D(E)}, which yields the result of Theorem 4.6 by finally following the
arguments in the proof of Theorem 4.5.

4.4

Convergence of µγP towards µ0P in Bregman divergence

Theorem 4.7. If Ω is a compact of Rd and ∇E(µ0P ) is a bounded function on Ω then
lim DE (µγP , µ0P ) = 0,

γ→0

12

where DE denotes the Bregman divergence between two measures µ and ν defined as DE (µ, ν) =
E(µ) − E(ν) − h∇E(ν), µ − νi.
Proof. By definition (4.2) of µγP , we get that
Z

W22 (µγP , ν)dP(ν)

−

Z

W22 (µ0P , ν)dP(ν) + γ(E(µγP ) − E(µ0P )) ≤ 0.

(4.10)

By definition (4.3) of µ0P , one has W22 (µγP , ν)dP(ν) − W22 (µ0P , ν)dP(ν) ≥ 0. Therefore, by
definition of the Bregman divergence, inequality (4.10) gives
R

R

DE (µγP , µ0P ) ≤ h∇E(µ0P ), µ0P − µγP i ≤ C sup hφ, µ0P − µγP i ≤ CdBL∗ (µ0P , µγP ),
kφkBL ≤1

where dBL∗ is the bounded Lipschitz distance and kφkBL := kφk∞ + kφkLip . We denote by
k · kLip the norm define on the space of all Lipschitz functions on (Ω, d) with d(x, y) = 1x6=y .
By hypothesis, k∇E(µ0P )kBL is finite. For sequence of probability measures, convergence
in distance dBL∗ is equivalent to weak convergence (e.g. Section 1.2.1 of [38]). Hence, by
Theorem 2.1.(d) in [14], JPγ Γ-converges to JP0 . Indeed for every sequence (µγ )γ ⊂ P2 (Ω)
converging to µ ∈ P2 (Ω) in bounded Lipchitz distance,
JP0 (µ) ≤ lim inf JPγ (µγ )
γ→0

by lower semicontinuity of JPγ in W2 . Moreover, there exists a sequence (µγ )γ converging
to µ (for instance take (µγ )γ constant equals to µ) such that lim JPγ (µγ ) = lim JPγ (µ) =
γ→0

γ→0

JP0 (µ). One can also notice that JPγ : P2 (Ω) → R is equi-coercive: for all t ∈ R, the set
{ν ∈ P2 (Ω) such that JPγ (ν) ≤ t} is included in a compact Kt since it is closed in the
compact set P2 (Ω) (by compactness of Ω). Therefore, we can apply the fundamental theorem
of Γ-convergence (Theorem 2.10 in [14]) to the bounded Lipschitz metric to obtain that
dBL∗ (µ0P , µγP ) −→ 0.
γ→0

5

Algorithmic approach and numerical experiments

In this section, we first present a method to automatically choose the parameter γ. Then, we
present numerical experiments on simulated and real data sets in R and R2 . A discretization of
the minimization problem (1.2) is used to compute a numerical approximation of a regularized
Wasserstein barycenter µγPn . It consists of using a fixed grid {xk }N
of equally spaced
PNk=1 k
γ
k
d
points x ∈ R , and to approximate µPn by the discrete measure k=1 f δxk where the f k
are positive weights summing to one which minimize a discrete version of the optimisation
problem (1.2). Further algorithmic details are given in an C. In these numerical experiments,
we focus on the case where E(µ) = +∞ if µ is not a.c. to enforce the regularized Wasserstein
barycenter to have a smooth pdf (we write E(f ) = E(µf ) if µ has a density f ). In this setting,
if the grid of points is of sufficiently large size, then the weights f k yield a good approximation
of this pdf.

5.1

Choice of the parameter γ

By analogy with the work in [5] based on the Lepskii balancing principle, we use an automatic selection of the regularization parameter γ. The method proposed in [5] requires the
13

knowledge of an upper bound on the decay to zero of the variance term E(d2E (µγPn , µγP )) as
γ → 0 which is given by (4.5). To match their notation, we set λ = 1/γ. Hence λ → +∞
corresponds to γ → 0.
Without the knowledge of µP , the Lepskii balancing principle described in [5] to select an
appropriate λ works as follows:
- For σ > 1 and a threshold Λ > 0, the Look-Ahead function is defined as lΛ,σ (λ) =
min{min{κ|ρ(λ) > σρ(κ)}, Λ} where the choice ρ : λ 7→ λ1 comes from the upper bound
(4.5).
- For δ > 0, the balancing functional reads
bΛ,σ (λ) =

max

λ<κ≤lΛ,σ (λ)

{

1
1/λ
1/κ
dE (µPn , µPn )ρ(κ)}
4δ

- The data-driven choice is given by λΛ,σ, = min{λ ≤ Λ ; BΛ,σ (λ) ≤ } where BΛ,σ (λ) =
max {bΛ,σ (κ)} is the smooth balancing functional, and  > 0 is a parameter to control
λ<κ≤Λ

the stability of regularized barycenters for successive choices of λ = 1/γ.
We display in Figure 2 an example of the smooth balancing functional BΛ,σ (λ) associated
to the simulated Gaussian data from Section 5.2 for three different choices for the penalty
function E.
In this paper, instead of choosing a value for , we use the fact that BΛ,σ is a decreasing
function, and that it has the shape of a L-curve. Hence, a data-driven value λ̂Λ,σ is chosen
by determining the location where the curve BΛ,σ has an “elbow” (change of curvature). In
Figure 2, this strategy leads to a data-driven value for λ that is close to the ideal choice given
1/λ
1/λ
by the minimizer of λ 7→ dE (µP , µPn )/ minλ dE (µP , µPn ).
25

14

14

12

12
20

10

10

15
8

8

6

6
10

4

4
5

2

2

0

0
0.1

0.2

0.3

0.4

0.5

0.6

0 2

0.7

0.8

0.9

(a) E(f ) = ||f || Dirichlet

1

0
0

10

20

30

(b) E(f ) =

40

R

50

60

70

80

90

f log(f ) Entropy

100

10

20

30

40

50

60

70

80

90

100

(c) Dirichlet + Entropy

Figure 2: Simulated Gaussian data from Section 5.2. Smooth balancing functional BΛ,σ (λ)
(times 10−6 ) in solid lines for three different regularizations (with σ = 3). For these simulated
1/λ
1/λ
data, we have access to µP , and the dotted lines represent dE (µP , µPn )/ minλ dE (µP , µPn ) as
functions of λ.

14

Numerical experiments for d = 1

5.2

Simulated data We consider a simulated example where the measures ν i are discrete and
supported on a small number pi of data points. To this end, for each i = 1, . . . , n, we simulate a
sequence (X ij )1≤j≤pi of iid random variables sampled from a Gaussian distribution N (µi , σ 2i )
where the pi ’s are ranging from 5 to 10, and the µi ’s (resp. σ i ) are iid random variables such
that −2 ≤ µi ≤ 2 and 0 ≤ σ i ≤ 1 with E(µi ) = 0 and E(σ i ) = 1/2. The target measure that
we wish to estimate in these simulations is the population (or true) Wasserstein barycenter of
the random distribution N (µ1 , σ 21 ) which is N (0, 1/4) thanks to the assumptions E(µ1 ) = 0
P i
and E(σ 1 ) = 1/2. Then, we define the random discrete measure ν i = p1i pj=1
δX ij .
To illustrate the benefits of regularizing the Wasserstein barycenter of the ν i ’s, we compare
our estimator with the one obtained by the following procedure which we refer to as the
kernel method. In a preliminary step, each measure ν i is smoothed using a standard kernel
density estimator whose bandwidth hi is chosen by cross-validation. An alternative estimator
is then defined as the Wasserstein barycenter of these smoothed measures which can be
easily computed thanks to the quantile averaging formula for measures supported on R (see
e.g. Section 6.1 in [1]). This estimator corresponds to the notion of smoothed Wasserstein
barycenter of multiple point processes as considered in [32]. The density of this smoothed
Wasserstein barycenter is displayed in Figure 3. For this example, it appears a preliminary
smoothing of the ν i followed by quantile averaging is not sufficient to recover a satisfactory
Gaussian shape when the number pi of observations per unit is small.
0.09

0.09
true barycenter
kernel method
regularized barycenter

0.08

0.07

0.07

0.06

0.06

0.05

0.05

0.04

0.04

0.03

0.03

0.02

0.02

0.01

0.01

0
-4

-3

-2

-1

0

1

2

3

true barycenter
kernel method
regularized barycenter

0.08

4

(a) Dirichlet regularization

0
-4

-3

-2

-1

0

1

2

3

4

(b) Entropy regularization

Figure 3: Simulated data from Gaussian distributions with random means and variances. In
all the figures, the black curve is the density of the true Wasserstein barycenter. The blue
and dotted curve represents the pdf of the smoothed Wasserstein barycenter obtained by a
preliminary kernel smoothing step. Pdf of the regularized Wasserstein barycenter µRγPn (a) for
20 ≤ γ ≤ 50 with E(f ) = ||f 0 ||2 (Dirichlet), and (b) for 0.08 ≤ γ ≤ 14 with E(f ) = f log(f )
(negative entropy).
Alternatively, we have applied the algorithm described Section 3.1 of the C directly on
the (non-smoothed) discrete measures ν i to obtain a numerical approximation of regularized
barycenter µγPn with two different choices for the penalty function E. The results are displayed
in Figure 3 for different values of γ around the data-driven choice 1/λ̂Λ,σ from Section 5.1. For
both penalty functions and despite a small number of observations per experimental units,
15

the shape of these densities better reflects the fact that the population Wasserstein barycenter
is a Gaussian distribution.
Finally, we provide Monte-Carlo simulations to illustrate the influence of the number n
of observed measures on the convergence of these estimators. For a given 10 ≤ n0 ≤ n,
we randomly draw n0 measures ν i from the whole sample, and we compute the following
estimators: a smoothed barycenter via the kernel method, a regularized barycenter using a
data-driven choice for γ, and an ideal regularized barycenter obtained by the Lepskii method
if we had access to the true population barycenter by setting γ as the minimizer of γ 7→
dE (µP , µγPn ). For given value of n0 , this procedure is repeated 200 times, which allows to
obtain an approximation of the expected error E (d(µ̂, µP )) of each estimator µ̂, where d is
either the Bregman or the Wasserstein distance. The penalty used is a linear combinaison
of Dirichlet and negative entropy functions. The results are displayed in Figure 4. It can be
observed that our approach yields better results than the kernel method for both types of
error (using either the Bregman or Wasserstein distance).
0.5

1
Regularized barycenter
Regularized ideal barycenter
Kernel estimator

0.45

0.4

0.8

0.35

0.7

0.3

0.6

0.25

0.5

0.2

0.4

0.15

0.3

0.1

0.2

0.05
10

25

50

75

Regularized barycenter
Regularized ideal barycenter
Kernel estimator

0.9

0.1
10

(a) Error with Bregman distance

25

50

75

(b) Error with Wasserstein distance

Figure 4: Errors in terms of expected Bregman and Wasserstein distances between the population barycenter and the estimated barycenters (kernel method in dashed blue, regularized
barycenter in full red and ideal regularized barycenter in dashed red) for a sample of size
n0 = 10, 25, 50 and 75 (the whole sample is made of n = 100 discrete measures).

A real data set We consider now a real data set of neural spike trains which is publicly
available from the MBI website: http://mbi.osu.edu/2012/stwdescription.html. During
a squared-path task, the spiking activity of a movement-encoded neuron of a monkey has been
recorded during 5 seconds over n = 60 repeated trials. Each spike train is then smoothed
using a Gaussian kernel (further details on the data collection can be found in [40]). For each
trial 1 ≤ i ≤ n, we let νi be the measure with pdf proportional to the sum of these Gaussian
kernels centered at the times of spikes, see Figure 5(a). This is an example of a dataset made
of a.c. measures (histograms). The pdf of the Wasserstein barycenter ν̄n of these measures
is displayed in Figure 5(b). This approach leads to a an irregular mean density of spiking
activity.
To compute a regularized Wasserstein barycenter from this dataset, we apply the algorithm
described in Section 3.1 of the C. In Figure 6(a), we display the densities of the regularized
Wasserstein barycenters with a Dirichlet regularization obtained for 6 ≤ γ ≤ 10, where this
16

0.08

0.03

0.06
0.04

0.025

0.02
0
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0.02

5

0.06
0.015
0.04
0.02
0.01
0
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5
0.005

0.06
0.04
0
0.02
0

-0.005
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

2

(a)

2.5

3

3.5

4

4.5

5

(b)

Figure 5: Neural spike trains data. (a) A subset of 3 smoothed spikes out of n = 60 of the
neural activity of a monkey during 5 seconds. Each row represents one trial and the pdf
obtained by smoothing each spike train with a Gaussian kernel of width 50 milliseconds. (b)
Pdf of the empirical Wasserstein barycenter ν̄n for this data set.
range of values is chosen by the adaptative Lepskii balancing principle described in Section 5.1.
In Figure 6(b), we display the results obtained with a negative entropy regularization for the
data-driven range of values 0.08 ≤ γ ≤ 0.12. Comparing Figures 5 and 6, this approach allows
to clearly smooth the result obtained by quantile averaging of the νi ’s (which corresponds to
γ = 0).

0.02

0.02

0.015

0.015

0.01

0.01

0.005

0.005

0

0
0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0 2

(a) E(f ) = ||f || Dirichlet

0.5

1

1.5

(b) E(f ) =

2

R

2.5

3

3.5

4

4.5

5

f log(f ) Entropy

Figure 6: Neural spike trains data. (a) Pdf of the regularized empirical Wasserstein barycenter
with a Dirichlet regularization for 6 ≤ γ ≤ 10. (b) Pdf of the regularized empirical Wasserstein
barycenter with negative entropy regularization and 0.08 ≤ γ ≤ 0.12

17

5.3

Numerical experiments for d = 2

We consider the real dataset described previously on the locations of reported incidents of
crime in Chicago. The city of Chicago is represented as an image of size 92 × 59, and a
crime is considered as a Dirac located at some pixel (see Figure 1). To compute a regularized
Wasserstein barycenter from this dataset, we apply the algorithm described in Section 3.2 of
the C for d = 2. For a Dirichlet regularization, the adaptative Lepskii’s strategy from Section
5.1 leads us to choose λ = 2.10−3 as it can be observed from Figure 7(a). We compare our
approach with the one in [21] which consists in using a regularized barycenter associated to an
entropically regularized transportation cost. The computation of such regularized barycenters
is obtained via the so-called Sinkhorn’s algorithm, see [21] for further details. The parameter
which controls the amount of transportation plan regularization is fixed to  = 1 (following
the notation in [21]). The parameter γ controlling the amount of regularization of such
Wasserstein barycenter is again chosen in a data-driven way using the Lepskii’s strategy.
Using transportation plan regularization, it follows than the best choice is λ = 10−2 , as it can
be seen from Figure 7(b).
43.4

0.15

43.3

43.2
0.1
43.1

43

42.9
0.05
42.8

42.7

42.6
2e-4

1e-3

1.2e-3

1.4e-3

2e-3

2.9e-3

5e-3

1e-2

0
2e-5

1e-4

2e-4

1e-3

2e-3

1e-2

2e-2

0.1

0.2

1

10

(a) Our approach: regularized Wasser- (b) The approach in [21]: regularized
stein barycenter with unregularized Wasserstein barycenter with regularized
transport cost.
transport cost and Sinkhorn’s algorithm.

Figure 7: Real dataset of reported incidents of crime in Chicago. Smooth balancing functional
BΛ,σ (λ) associated to regularized Wasserstein barycenters for different values of λ.
These two choices for γ lead to the regularized Wasserstein barycenters whose pdf are displayed in Figure (8)(a) and Figure (8)(b). The main differences between the estimator is the
spreading of the mass in Figure (8)(b) due to the use of an entropically regularized transportation cost. Finally, we present in Figure (8)(c) a comparison with a standard kernel density
estimator using the whole dataset (over all the year 2014) obtained from an implementation
in Matlab2 of a bivariate Gaussian kernel density estimator which uses a data-driven choice
for the bandwidth suggested in [11]. The method of kernel density yields a much smoothed
estimator than those obtained with regularized Wasserstein barycenters.
2

https://www.mathworks.com/examples/matlab/community/20312-bivariate-kernel-density-estimation

18

(a) Our approach: pdf of regularized Wasserstein (b) The approach in [21]: regularized Wasserbarycenter with unregularized transport cost.
stein barycenter with regularized transport cost and
Sinkhorn’s algorithm.

(c) Standard Kernel density estimator of the whole
dataset with a data-driven bandwidth.

Figure 8: (a) Barycenter and (b) kernel density estimator associated to the location of crimes
in the city of Chicago during the year 2014 are represented in color scale. The red points
represent locations of all crimes during this month.

A

Bregman divergence and symmetric distance

Definition A.1 (Subdifferential). Let X be a vector space and let X 0 its dual. For a convex
functional G : X → R, we say that q ∈ X 0 is a subgradient of G at u if it satisfies the

19

inequality
G(v) ≥ G(u) + hq, v − ui for every v ∈ X

(A.1)

where h·, ·i is the duality pairing (or linear form) of X and X 0 . The set of subgradients at u
is the subdifferential of G, denoted ∂G(u).
Definition A.2 (Bregman divergence). The (generalized) Bregman divergence related to a
convex functional G : X → R is defined by
q
DG
(u, v) = G(u) − G(v) − hq, u − vi for u, v ∈ X,

where q ∈ ∂G(v). The symmetric Bregman distance is defined by
q
p
dp,q
G (u, v) = DG (u, v) + DG (v, u) = hp − q, u − vi

where p ∈ ∂G(u) and q ∈ ∂G(v).

B

Proofs of main Theorems

B.1

Proof of the subgradient’s inequality

Proof of Lemma 2.6. 2 ⇒1. Let φ ∈ ∂J(µ) such that hφ, η − µi ≥ 0 for all η ∈ P2 (Ω). By
definition of the subgradient, ∀ η ∈ P2 (Ω), we have J(η) ≥ J(µ) + hφ, η − µi which is greater
than J(µ) by assertion. Hence µ minimizes J.
1 ⇒2. Take µ ∈ int(dom J) (that is J(µ) < +∞) such that µ is a minimum of J over P2 (Ω).
Then the directional derivative of J at the point µ along (η − µ) exists (Proposition 2.22 in
[19]) and satisfies
J(µ + t(η − µ)) − J(µ)
J 0 (µ; η − µ) := lim
≥ 0.
(B.1)
t→0
t
t>0
Remark that P2 (Ω) is a convex set. By Proposition 4.3 of [19], since J is a proper convex
function and µ ∈ dom(J), we obtain the equivalence
φ ∈ ∂J(µ) ⇔ hφ, ∆i ≤ J 0 (µ; ∆) for all ∆ ∈ P2 (Ω).
Moreover, since J is proper convex and lower semi-continuous, so is J 0 (f ; ·). Given that P2 (Ω)
is a Hausdorff convex space, we get by Theorem 7.6 of [2], that for all (η−µ) ∈ P2 (Ω), J 0 (µ; η−
µ) = sup{hφ, η − µi where φ is such that hφ, ∆i ≤ J 0 (µ; ∆), ∀∆ in P2 (Ω)}. Hence by (B.1) we
get sup hφ, η − µi ≥ 0. We then define the ball B = {η + µ ∈ M(Ω) such that kηkT V ≤ },
φ∈∂J(µ)

where k · kT V is the norm of total variation. We still have
inf

sup hφ, η − µi ≥ 0.

η∈B ∩P2 (Ω) φ∈∂J(µ)

Note that ∂J(µ) in a convex set. Moreover B ∩ P2 (Ω) is compact, and (φ, η) 7→ hφ, η − µi is
bilinear. Thus we can switch the infimum and the supremum by the Ky Fan’s theorem (4.36
in [19]). In that way, there exists φ ∈ ∂J(f ) such that
inf
hφ, η − µi ≥ 0. By convexity
η∈B ∩P2 (Ω)

of P2 (Ω), any ζ ∈ P2 (Ω) can be written as t(η − µ) + µ for some t ≥ 0 and η ∈ B ∩ P2 (Ω).
This concludes the proof of the lemma.
20

Proof of Lemma 2.7. (⇐).
We first
assume that for φµ,ν ∈ L1 (µ), there exists ψ µ,ν ∈ L1 (ν)
R
R
such that W22 (µ, ν) = φµ,ν dµ + ψ µ,ν dν and φµ,ν (x) + ψ µ,ν (y) ≤ |x − y|2 . Then for all
η ∈ P2 (Ω), denoting (φη,ν , ψ η,ν ) an optimal couple for η and ν, we get
W22 (η, ν) =

Z

Z

φdη +

sup

Z

ψdν =

φη,ν dη +

Z

ψ η,ν dν

φ(x)+ψ(y)≤|x−y|2

≥ W22 (µ, ν) +

Z

φµ,ν d(η − µ).

Hence, from the definition of a subgradient, φµ,ν ∈ ∂1 W22 (µ, ν).
(⇒). We denote by F the function µ ∈ P2 (Ω) 7→ W22 (µ,
ν). Let φ∗ ∈ ∂F (µ), then by the
R ∗
∗
∗
Legendre-Fenchel theory, we have that F (φ ) + F (µ) = φ dµ, where F ∗ denote the Fenchel
conjugate of F . Recall we want to show that there exists ψ ∈ L1 (ν) verifying φ∗ (x) + ψ(y) ≤
|x − y|2 such that
Z
Z
φ∗ dµ − W22 (µ, ν) = −

ψdν,

which is equivalent to F ∗ (φ∗ ) = − ψdν. In that purpose, we first define ψ φ (·) := inf {| ·
R

y∈Ω

−y|2

R

ψ φ dν.

− φ(y)} and H(φ) := − R
By definition, H ∗ (µ) = sup { φdµ − H(φ)}. Observing that H is convex, l.s.c. on Y and
φ∈Y

proper as :
H(φ) = −

Z

ψ φ dν =

Z

sup{φ(y) − |x − y|2 } dν(x)
y∈Ω

≥

Z

(φ(y0 ) − 2|y0 | − 2|x|2 )dν(x) > −∞ by definition of ν,

where y0 ∈ Ω is such that φ(y0 ) is finite. We get H ∗∗ (φ) = H(φ) by Theorem 2.3.3. in [41].
Moreover, for µ ∈ P2 (Ω), we have by the duality formulation of Kantorovich (e.g Lemma 2.1.
of [1]) that
W22 (µ, ν)

Z

Z

φdµ +

= sup
ZΩ

Z

ZΩ
φ

We deduce that H ∗∗ (φ) =

Z

Ω



2

ψdν; φ, ψ ∈ Cb , ψ(y) ≤ inf {|x − y| − φ(x)}
x

Ω

φdµ +

= sup



Ω

φdµ +

= sup

ψdν; φ, ψ ∈ Cb , φ(x) + ψ(y) ≤ |x − y|

2

ψφ



= H ∗ (µ).

Ω

sup { φdµ−W22 (µ, ν)} = F ∗ (φ), which implies F ∗ (φ∗ ) = H(φ∗ ).
R

f ∈P2 (Ω)

∗

Thus we end up with the equality F (µ) = φ∗ dµ − F ∗ (φ∗ ) = φ∗ dµ + ψ φ dν. This exactly
∗
∗
means that Rfor φ∗ ∈ R∂1 W22 (µ, ν), there exists ψ φ such that φ∗ (x) + ψ φ (y) ≤ |x − y|2 and
∗
W22 (µ, ν) = φ∗ dµ + ψ φ dν, which concludes the proof.
R

B.2

R

R

Proof of existence, uniqueness and stability of regularized barycenters

Proof of Theorem 3.1. Let (µk )k ⊂ P2 (Ω) a minimizing sequence of probability measures of
JPγn . Hence, there exists a constant M ≥ 0 such that ∀k, JPγn (µk ) ≤ M . It follows that for

21

all k,

1
n

Pn

2 k
i=1 W2 (µ , νi )

≤ M . By Lemma 2.1 of [1] we thus have

n
n
X
1X
2 i k
W (ν , µ ) = 2
sup
n i=1 2
i=1 f ∈Z

Z

k

Z

f dµ +
Rd

i



Sf (x)dν (x) ≤ M,

Rd

1
where Sf (x) = inf { 2n
|x − y|2 − f (y)}. Since the function x 7→ |x|α (with 1 < α < 2) belongs
y∈Ω

to Z, we have that Rd |x|α dµk (x) is bounded by a constant L ≥ 0 for all k. We deduce that
(µk )k is tight (for instance, take the compact K c = {x ∈ Ω such that |x|α > L }). Since
(µk )k is tight, by Prokhorov’s theorem, there exists a subsequence of (µk )k (still denoted
(µk )) which weakly converges to a probability measure µ. Moreover, one can prove that
µ ∈ P2 (Ω). RIndeed for all lower
semicontinuous functions bounded from below by f , we have
R
that lim inf Ω f (x)dµk (x) ≥ Rd f (x)dµ(x) by weak convergence. Hence for f : x 7→ |x|2 , we
R

get

R k→∞
2
Ω |x|

dµ(x) ≤ lim inf
k→∞

R

2
k
Ω |x| dµ (x)

< +∞, and thus µ ∈ P2 (Ω).

Let (πik )1≤i≤n,1≤k be a sequence of optimal transport plans
where πik is an optimal transRR
k
2
k
port plan between µ and νi . Since supk W2 (µ , νi ) = supk Ω×Ω |x − y|2 dπik (x, y) < +∞,
we may apply Proposition 7.1.3 of [4]: (πik )k is weakly relatively compact on the probability
space over Ω × Ω and every weak limit πi Ris an optimal transport plan between µ and νi
with, for all 1 ≤ i ≤ n, W22 (µ, νi ) ≤ lim inf Rd ×Rd |x − y|2 dπik (x, y) < +∞. Since E is lower
k→∞

semicontinuous, we get that
lim inf JPγn (µk ) = lim inf
k→∞

k→∞

≥

n
1X
W 2 (µk , νi ) + γE(µk )
n i=1 2

n
1X
W 2 (µ, νi ) + γE(µ) = JPγn (µ).
n i=1 2

Hence JPγn admits at least µ ∈ P2 (Ω) as a minimizer. Finally, by the strict convexity of JPγn
on its domain, the minimizer is unique and it belongs to D(E) as defined in (2.5), which
completes the proof.
Proof of Theorem 3.2. First, let us prove the existence of a minimizer. For that purpose, we
decide to follow the sketch of the proof of the existence of a Wasserstein barycenter given by
Theorem 1 in [31]. We suppose that (Pn )n≥0 ⊆ W2 (P2 (Ω)) is a sequence of measures, such
that µn ∈ P2 (Ω) is a probability measure minimizing JPγn , for all n. Furthermore, we suppose
that there exists P ∈ W2 (P2 (Ω)) such that W2 (P, Pn ) −→ 0. We then have to prove that
n→+∞

(µn )n≥1 is precompact and that all limits minimize JPγ . We denote µ̃ a random measure with
distribution P and µ̃n a random measure with distribution Pn . Hence
W2 (µn , δx ) = W2 (δµn , δδx ) ≤ W2 (δµn , Pn ) + W2 (Pn , δδx )
= E(W22 (µn , µ̃n ))1/2 + E(W22 (µ̃n , δx ))1/2 .
Moreover, E(W22 (µn , µ̃n ))1/2 ≤ M for a constant M ≥ 0 since µn minimizes JPγn and µ̃n is of
law Pn . Then
W2 (µn , δx ) ≤ M + W2 (Pn , δδx ) ≤ M + W2 (Pn , P) + W2 (P, δδx ) ≤ L

22

since W2 (Pn , P) −→ 0 and P ∈ W2 (P2 (Ω)) by hypothesis. By Markov inequality, we have
n→+∞

for r > 0
µn (B(x, r)c ) = Pµn (|X − x|2 ≤ r2 ) ≤

Eµn (|X − x|2 )
W22 (µn , δx )
=
r2
r2

2

and µn (B(x, r)c ) ≤ Lr2 . Hence (µn )n is tight: it is possible to extract a subsequence (still
denoted (µn )) which converges weakly to a measure µ by Prokhorov’s theorem. Let us show
that µ minimizes JPγ . Let η ∈ P2 (Ω) and ν ∈ P2 (Ω) with distribution P.
JPγ (η) = EP (W22 (η, ν)) + γE(η) = W22 (δη , P) + γE(η)
= lim W22 (δη , Pn ) + γE(η) since by hypothesis W2 (Pn , P) → 0
n→+∞

≥ lim inf W22 (δµn , Pn ) + γE(µn ) since µn minimizes JPγn
n→+∞

(B.2)

Moreover, we have by the inverse triangle inequality that
lim inf W2 (δµn , Pn ) ≥ lim inf (W2 (δµ , Pn ) − W2 (δµ , δµn )) = W2 (δµ , P).
n→+∞

n→+∞

The last inequality comes from the two convergences W2 (Pn , P) → 0 and W2 (δµ , δµn ) → 0.
From (B.2) and by lower semicontinuity of E, we get JPγ (η) ≥ W22 (δµ , P) + γE(µ) = JPγ (µ).
Hence µ minimizes JPγ . To finish the proof of the existence of a minimizer, we need the
following result whose proof can be found in [31].
Theorem B.1. For all P ∈ W2 (P2 (Ω)), there is a sequence of finitely supported distributions
P
PK
2
Pn (that is Pn = K
k=1 λk δκk where
k=1 λk = 1) such that W2 (Pn , P) −→ 0.
n→+∞

Now, by Theorem B.1 it follows that for a given distribution P, one can find a sequence of
finitely supported distributions Pn such that for all n there exists a unique measure µn ∈ P2 (Ω)
minimizing JPγn using Theorem 3.1 and such that W22 (Pn , P) −→ 0 thanks to Theorem B.1.
n→+∞

Therefore there is a probability measure µ which minimizes JPγ . Let us make sure that µ is
indeed in the space P2 R(Ω). From Theorem 3.1,
we also have that µn ∈ P2 (Ω) for all n. Thus
R
by weak convergence, Ω |x|2 dµ(x) ≤ lim inf Ω |x|2 dµn (x) < +∞. Finally, the uniqueness of
n→+∞

the minimum is obtained by the strict convexity of the functional µ 7→ EP (W22 (µ, ν)) + γE(µ)
on the domain D(E), which completes the proof.

B.3

Proof of the stability’s Theorem 3.3

Proof. We denote by µ, ζ ∈ P2 (Ω) the probability measures such that µ minimizes JPγνn
and ζ minimizes JPγη . For each 1 ≤ i ≤ n, one has that θ 7→ n1 W22 (θ, νi ) is a convex,
n
proper and continuous function. Therefore, Theorem 4.10 in [19], we have that ∂JPνn (µ) =
1 Pn
∂1 W22 (µ, νi ) + γ∇E(µ). Hence by Lemma 2.7, any φ ∈ ∂JPνn (µ) is of the form
i=1 P
n
1
φ = n ni=1 φi + γ ∇E(µ) where for all i = 1, . . . , n, φi = φµ,νi is optimal in the sense
that (φµ,νi , ψ µ,νi ) is an optimal couple associated to (µ, νi ) in the Kantorovich formulation ofPthe Wasserstein distance (see Theorem 2.2). Therefore by Lemma 2.6, there exists
φ = n1 ni=1 φµ,νi + γ∇E(µ) such that hφ, θ − µi ≥ 0 for all θ ∈ P2 (Ω). Likewise, there exists
P
φ̆ = n1 ni=1 φζ,ηi + γ∇E(ζ) such that hφ̆, θ − ζi ≥ 0 for all θ ∈ P2 (Ω). Finally, we obtain that

23

!
n
1 X µ,νi
ζ,ηi
− φ ) d(µ − ζ).
(φ
n i=1

Z
γh∇E(µ) − ∇E(ζ), µ − ζi ≤ −
Ω

Following the proof of Kantorovich duality’s theorem in [38], we can restrict the supremum
over (φ, ψ) ∈ CW in Kantorovich’s duality Theorem 2.2 to the admissible pairs (φcc , φc ) where
φc (y) = inf x {|x − y|2 − φ(x)} and φcc (x) = inf y {|x − y|2 − φc (y)}. Then, we replace φµ,νi by
(φµ,νi )cc (resp. φζ,ηi by (φζ,ηi )cc ) and ψ µ,νi by (φµ,νi )c (resp. ψ ζ,ηi by (φζ,ηi )c ) and obtain
that
n

γh∇E(µ) − ∇E(ζ), µ − ζi ≤ −

1X
n i=1

Z

 µ,νi cc

(φ ) (x) − (φζ,ηi )cc (x) d(µ − ζ)(x)

Ω

n ZZ
 µ,νi cc

1X
=−
(φ ) (x) − (φζ,ηi )cc (x) d(π µ,νi − π ζ,ηi )(x, y),
n i=1 Ω×Ω

where π µ,νi is an optimal transport plan on Ω × Ω with marginals µ and νi for i ∈ {1, . . . , n}
(and π ζ,ηi optimal with marginals ζ and ηi ). Developing the right-hand size expression in the
above inequality, we get
γh∇E(µ) − ∇E(ζ), µ − ζi

ZZ
n  ZZ
1X
µ,νi cc
µ,νi
ζ,ηi cc
ζ,ηi
−
(φ ) (x)dπ
(x, y) −
(φ ) (x)dπ (x, y)
≤
n i=1

ZZ
n Z Z
1X
(φµ,νi )cc (x)dπ ζ,ηi (x, y) +
(φζ,ηi )cc (x)dπ µ,νi (x, y) .
+
n i=1

From the condition (2.3) in the Kantorovich’s dual problem, we have that (φµ,νi )cc (x) ≤
2 − (φζ,ηi )c (y) for all i ∈ {1, . . . , n}. More|x − y|2 − (φµ,νi )c (y) and (φζ,ηi )cc (x) ≤ |x −
 y| µ,ν

µ,ν
cc
µ,ν
i
i
over, we have that (φ )h (x)dπ (x, y) = −(φi i )c (y) + |x − y|2 dπ µ,νi (x, y) and likewise
(φζ,ηi )cc (x)dπ ζ,ηi (x, y) = −(φζ,ηi )c (y) + |x − y|2 dπ ζ,ηi (x, y). We therefore deduce that
n

γh∇E(µ) − ∇E(ζ), µ − ζi ≤ −

1X
n i=1

−

1X
n i=1

n

n

1X
+
n i=1
n

+
=

1
n

1X
n i=1
n Z
X
i=1

ZZ

ZZ

ZZ

ZZ



−(φµ,νi )c (y) + |x − y|2 dπ µ,νi (x, y)


−(φζ,ηi )c (y) + |x − y|2 dπ ζ,ηi (x, y)


−(φµ,νi )c (y) + |x − y|2 dπ ζ,ηi (x, y)


−(φζ,ηi )c (y) + |x − y|2 dπ µ,νi (x, y)

 µ,νi c

(φ ) (y) − (φζ,ηi )c (y) d(νi − ηi )(y).

Ω

For all 1 ≤ i ≤ n, and y, y 0 ∈ Ω, we have
(φµ,νi )c (y) − (φµ,νi )c (y 0 ) =sup{φµ,νi (x) − |x − y 0 |2 } + inf {|x − y|2 − φµ,νi (x)}
x

x

≤sup{|x − y| − |x − y | } ≤ 4cΩ |y − y 0 |
2

x

24

0 2

where cΩ = sup|x|. As a consequence,
x
1
8cΩ

h

1
µ,νi )c
4cΩ (φ

is 1-Lipschitz and so is

1
ζ,ηi )c ,
4cΩ (φ

which

i

(φµ,νi )c − (φζ,ηi )c is 1-Lipschitz for all 1 ≤ i ≤ n.
implies that
We then conclude
γh∇E(µ) − ∇E(ζ), µ − ζi ≤
=

Z

n
8cΩ X
φ d(νi − ηi ); φ ∈ ∩L1 (|νi − ηi |), kφkLip ≤ 1
sup
n i=1
n
n
8cΩ X
8cΩ X
W1 (νi , ηi ) ≤
W2 (νi , ηi ),
n i=1
n i=1

by the Kantorovich-Rubinstein theorem presented in [38], while the last inequality above
comes from Hölder inequality between the distance W2 and the distance W1 defined for θ1 , θ2
(probability measures on Ω with moment of order 1) as
Z Z

W1 (θ1 , θ2 ) = inf
π

|x − y|dπ(x, y)

Ω Ω

where π is a probability measures on Ω × Ω with respective marginals θ1 and θ2 . Since µ
and ζ are independent, we can assign to νi any ησ(i) for σ ∈ Sn the permutation group of
Pn
{1, . . . , n} and hence we obtain γh∇E(µ) − ∇E(ζ), µ − ζi ≤ n2 inf
i=1 W2 (νi , ησ(i) ), which
σ∈Sn

completes the proof.

B.4

Proof of convergence properties

Proof of Theorem 4.3. We denote by C a universal constant whose value may change from
line to line. From the subgradient’s inequality (2.9) and following the same process used in the
γ
proof of Theorem 3.3, we have that, for each ν i , i = 1, . . . , n, there exists φµPn ,ν i integrable
with respect to µγPn (x)dx such that for all η ∈ P2 (Ω):
*

n
γ
1X
φµPn ,ν i + γ∇E(µγPn ), η − µγPn
n i=1

+

≥ 0.

(B.3)

By applying once again the subgradient’s inequality, we get
µγP minimizes JPγ ⇔ ∃φ ∈ ∂JPγ (µγP ) s. t. hφ, η − µγP i ≥ 0 for all η ∈ P2 (Ω).
Let us explicit the form of a subgradient φ ∈ ∂JPγ (µγP ) using again the Theorem of the
subdifferential of a sum. We have that µ 7→ W22 (µ, ν) is continuous for all ν ∈ P2 (Ω).
Moreover
by symmetry, νR 7→ W22 (µ, ν) is
for all µ ∈ P2 (Ω) and W22 (µ, ν) ≤
RR
R measurable
2
2
2
|x − y| dµ(x)dν(y) ≤ 2 |x| dµ(x) + 2 |y| dν(y) ≤ C is integrable with respect to dP(ν)
(by compactness of Ω). Hence, by Theorem of continuity under integral sign, we deduce that
µ 7→ E[W22 (µ, ν)] is continuous. Thus we can manage the subdifferential of the following sum
and one has that ∂JPγ (µγP ) = ∂1 [E(W22 (µγP , ν))] + γ∇E(µγP ), where ν is still a random measure
with distribution P. Also the Theorem 23 in [36] implies ∂1 E[W22 (µγP , ν)] = E[∂1 W22 (µγP , ν)].
Hence, we can sum up
µγP

minimizes

JPγ

Z
⇔

φ

µγ
,ν
P

dP(ν) +

γ∇E(µγP ), η

25

−

µγP


≥ 0, ∀ η ∈ P2 (Ω).

(B.4)

In the sequel, to simplify the notation, we use µ := µγPn and η := µγP . Therefore thanks to
(B.3) and (B.4)
dE (µ, η) = h∇E(µ) − ∇E(η), µ − ηi
1
≤−
γ
=

1
γ

*

−

n
1X
φµ,ν i −
n i=1
n
1X
n i=1

ZZ

φ

+

η,ν

Z

+

Z

φ

η,ν

dP(ν), µ − η
n
1X
n i=1

φµ,ν i (x)dµ(x) +

dP(ν)dµ(x) −

ZZ

φ

η,ν

Z

φµ,ν i (x)dη(x)


dP(ν)dη(x) .

(B.5)

We would like to switch integrals of the two last terms. In that purpose, we use that
Z

W22 (η, ν)dP(ν) ≤ C + 2

Z

sup



|y|2 dν(y) < +∞.

Ω

ν∈P2 (Ω)

As 0 ≤ W22 (η, ν)dP(ν) = ( φη,ν (x)dη(x) + ψ η,ν (x)dν(y)) dP(ν), we also have that φη,ν (x)dη(x)dP(ν) <
+∞.
Since x 7→ φη,ν (x) and
ν 7→
φη,ν (x) are measurables, we obtain by Fubini’s theorem
R R
R
R
η,ν dP(ν)dη(x) =
η,ν dη(x)dP(ν). By the same tools, since
Ω P2 (Ω) φ
P2 (Ω) Ω φ
R

R R

Z

R

W22 (µ, ν)dP(ν)

Z Z

φ

=
≥

Z Z

φ

µ,ν

η,ν

RR

Z

ψ

(x)dµ(x) +
Z

ψ

(x)dµ(x) +

µ,ν

η,ν



(x)dν(y) dP(ν)


(x)dν(y) dP(ν),

we get that ( φη,ν (x)dµ(x)) dP(ν) < +∞, so Ω P2 (Ω) φη,ν dP(ν)dµ(x) =
Therefore, by the dual formulation of Kantorovich, we have that
R R

R R

−

Z

−

φµ,ν i dµ(x) =
Z

φη,ν dη(x) =

Z
Z

ψ µ,ν i (y)dν i (y) −
ψ η,ν (y)dν(y) −

ZZ

ZZ

R

R

P2 (Ω) Ω φ

|x − y|2 dπ µ,ν i (x, y)

|x − y|2 dπ η,ν (x, y)

η,ν dµ(x)dP(ν).

(B.6)
(B.7)

where π µ,ν i and π η,ν are optimal transport plans for the Wasserstein distance. Also, φµ,ν i
and φη,ν verify the Kantorovich condition, that is
φµ,ν i (x) ≤ −ψ µ,ν i (y) + |x − y|2

(B.8)

φη,ν (x) ≤ −ψ η,ν (y) + |x − y|2 .

(B.9)

Next,
the trick is to write φµ,ν i (x)dη(x) =
φµ,ν i (x)dπ η,ν i (x, y) and φη,ν (x)dµ(x) =
RR η,ν
µ,ν
φ (x)dπ (x, y). Thus, by using the equalities (B.6), (B.7) and the inequalities (B.8),
(B.9), the result (B.5) becomes
R

n
1X
γdE (µ, η) ≤−
n i=1

Z ZZ

+

RR

ZZ

2

|x − y| dπ

µ,ν i

R

n
1X
(x, y) +
n i=1

Z ZZ

|x − y|2 dπ µ,ν (x, y)dP(ν) −

26

ZZ

|x − y|2 dπ η,ν i (x, y)

|x − y|2 dπ η,ν (x, y)dP(ν).

We denote
Sµnγ
Pn

Z ZZ

2

|x − y| dπ

:=
n
1X
n i=1

Sµnγ :=
P

ZZ

µγPn ,ν

n
1X
(x, y)dP(ν) −
n i=1
γ

|x − y|2 dπ µP ,ν i (x, y) − E

ZZ

ZZ

γ

|x − y|2 dπ µPn ,ν i (x, y)


γ

|x − y|2 dπ µP ,ν (x, y) ,

and finally the last inequality writes
γdE (µγPn , µγP ) ≤ Sµnγ + Sµnγ .
Pn

(B.10)

P

γ

Remark. Since for i = 1, . . . , n the random variables
|x − y|2 dπ µP ,ν i (x, y) are independent
and identically distributed. From the law of large numbers, we can notice that Sµnγ −→ 0
P
almost surely when n → +∞.
RR

Taking the expectation with respect to the random measures, (B.10) implies
γ 2 E(d2E (µγPn , µγP )) ≤ 2E(|Sµnγ |2 ) + 2E(|Sµnγ |2 ).
Pn





• Study of E |Sµnγ |2 . Using again the fact that

RR

(B.11)

P

γ

|x − y|2 dπ µP ,ν i (x, y) are iid, we get

P



E

|Sµnγ |2
P



1
= Var
n

ZZ

2

|x − y| dπ

µγP ,ν



(x, y) =

C
.
n

(B.12)

Note that since Ω is compact, the above variance term is finite.




• Study of E |Sµnγ |2 . This term can be controlled thanks to the empirical process theory.
Pn

Define the norm associated to the class of functions H (4.6) by kGkH := sup |G(h)| where
h∈H

G : H → R. Recall that hµ ∈ H is the function hµ : ν ∈ P2 (Ω) 7→ W22 (µ, ν), hence
Sµnγ =
Pn

Z
P2 (Ω)

hµγ (ν)dP(ν) −
Pn

Z
P2 (Ω)

hµγ (ν)dPn (ν) := (P − Pn )(hµγ )
Pn

Pn

1
≤ sup |(P − Pn ) (h) | = √ kGn kH
n
h∈H
√
where Gn (h) = n(Pn − P)(h). We then obtain


E

|Sµnγ |2
Pn



≤


1 
1
E kGn k2H = k kGn kH k2L2 (P) .
n
n

(B.13)

We finally use the following Theorem 2.14.1. of [37] to control this last expression:
Theorem B.2. Let H be a Q-measurable class of measurable functions with measurable envelope function H. Then for p ≥ 1,
k kGn kH kLp (Q) ≤ CI(1, H)kHkL2∨p (Q)
with C a constant, I(1, H) defined in (4.4) and H an envelope function.
Gathering the results of (B.11), (B.12), (B.13) and (B.14), we get




E(d2E µγPn , µγP ) ≤


1 
C
+
CI(1,
H)kHk
L2 (P)
γ2n

which concludes the proof.
27

(B.14)

C

Algorithmic details

In this section, we describe how the minimization problem
min
µ

n
1X
W 2 (µ, νi ) + γE(µ) over µ ∈ P2 (Ω),
n i=1 2

(C.1)

can be solved numerically by using an appropriate discretization to compute a numerical
approximation of a regularized Wasserstein barycenter.
More precisely, given a fixed grid {xk }N
of equally spaced points xk ∈ Rd , we apk=1
PN
γ
proximate µPn by the discrete measure µf = k=1 f k δxk where the f k are positive weights
summing up to one which minimize a discrete version of the optimisation problem (C.1). In
what follows, we first describe an algorithm that is specific to the one-dimensional case, and
then we propose another algorithm that is valid for any d ≥ 1.

C.1

Discrete algorithm for d = 1 and data defined on the same grid

We first propose to compute a regularized empirical Wasserstein barycenter for a dataset
made of discrete measures ν1 , . . . , νn (or one-dimensional histograms) defined on the same
γ
grid of reals {xk }N
k=1 that the one chosen to approximate µPn . Since the grid is fixed, we
identify a discrete measure ν with the vector of weights ν = (ν(x1 ), . . . , ν(xN )) in RN
+ (with
entries that sum up to one) of its values on this grid.
The estimation of the regularized barycenter onto this grid can be formulated as:
min
f

n
X
1X
W22 (f, νi ) + γE(f ) s.t
f k = 1, and f k = f (xk ) ≥ 0,
n i=1
k

(C.2)

with the obvious abuse of notation W22 (f, νi ) = W22 (µf , νi ) and E(f ) = E(µf ).
Then, to compute a minimizer of the convex optimization problem (C.2), we perform a
subgradient descent. We denote by (f (`) )`≥1 the resulting sequence of discretized regularized
barycenters in RN along the descent. Hence, given an initial value f (1) ∈ RN
+ and for ` ≥ 1,
we thus have
"

f

(`+1)

= ΠS f

(`)

−τ

(`)

n
1X
)+
∇1 W22 (f (`) , νi )
n i=1

#!

γ∇E(f

(`)

(C.3)

where τ (`) is the `-th step time, and ΠS stands for the projection on the simplex S = {y ∈
PN
j
RN
+ such that
j=1 y = 1}. Thanks to Proposition 5 in [35], we are able to compute a subgradient of the squared Wasserstein distance W22 (f (`) , νi ) with respect to its first argument
P
(for discrete distributions). For that purpose, we denote by Rf (s) = xj ≤s f (xj ) the cdf of
P
−
k
µf = N
k=1 f (x )δxk and by Rf (t) = inf{s ∈ R : Rf (s) ≥ t} its pseudo-inverse.
Proposition C.1 ([35]). Let f = (f (x1 ), f (x2 ), . . . , f (xN )) and ν = (ν(x1 ), ν(x2 ), . . . , ν(xN ))
be two discrete distributions defined on the same grid of values x1 , . . . , xN in R. For p ≥ 1,
the subgradients of f 7→ Wpp (f, ν) can be written as
∇1 Wpp (f, ν) : xj 7→

X

|xm − x̃m |p − |xm+1 − x̃m |p

m≥j

28

(C.4)

where

(

x̃m = xk if Rg (xk−1 ) < Rf (xm ) < Rν (xk )
x̃m ∈ [xk−1 , xk ] if Rf (xm ) = Rν (xk )

Even if subgradient descent is only shown to converge with diminishing time steps [12], we
observed that using a small fixed step time (of order 10−5 ) is sufficient to obtain in practice a
convergence of the iterates (f (`) )`≥1 . Moreover, we have noticed that the principles of FISTA
(Fast Iterative Soft Thresholding, see e.g. [6]) accelerate the speed of convergence of the above
described algorithm.

C.2

Discrete algorithm for d ≥ 1 in the general case

We assume that data ν1 , . . . , νn are given in the form of n discrete probability measures
(histograms) supported on Rd (with d ≥ 1) that are not necessarily defined on the same grid.
More precisely, we assume that
νi =

pi
X
j

νi δyj

j=1

i

for 1 ≤ i ≤ n where the yij ’s are arbitrary locations in Ω ⊂ Rd , and the νij ’s are positive
weights (summing up to one for each i).
d
The estimation of the regularized barycenter onto a given grid {xk }N
k=1 of R can then be
formulated as the following minimization problem:
min
f

n
X
1X
W22 (f, νi ) + γE(f ) s.t
f k = 1, and f k ≥ 0,
n i=1
k

(C.5)

with the notation f = (f 1 , f 2 , . . . , f N ) and the convention that W22 (f, νi ) denotes the squared
P
k
Wasserstein distance between µf = N
k=1 f δxk and νi .
Problem (C.5) could be exactly solved by considering the discrete pi ×N transport matrices
Si between the barycenter µf to estimate and the data νi . Indeed, problem (C.5) is equivalent
to the convex problem
p

min min
f

S1 ···Sn

n X
N
i X
1X
||y j − xk ||2 Sij,k + γE(f )
n i=1 j=1 k=1 i

(C.6)

under the linear constraints
∀i = 1, . . . , n,

pi
X
j,k

Si

j=1

= f k,

N
X
j,k

Si

= νij , and Sij,k ≥ 0.

k=1

However, optimizing over the pi × N transport matrices Si for 1 ≤ i ≤ n involves memory
issues when using an accurate discretization grid {xk }N
k=1 with a large value of N . For this
reason, we consider subgradient descent algorithms that allow dealing directly with problem
(C.5).
To this end, we rely on the dual approach introduced in [18] and the numerical optimisation
scheme proposed in [21]. Following these works, one can show that the dual problem of (C.5)
with a regularization of the form E(Kf ) and K a discrete linear operator reads as
min

φ0 ,···φn

n
X

Hνi (φi ) + Eγ∗ (φ0 ) s.t K T φ0 +

i=1

n
X
i=1

29

φi = 0,

(C.7)

∗
where the φi ’s are dual variables (vectors in RN ) defined on the discrete grid {xk }N
k=1 , Eγ is
the Legendre transform of γE and Hνi (.) is the Legendre transform of W22 (., νi ) that reads:

Hνi (φi ) =

pi
X
j



νi min

j=1

k=1···N

1 j
||y − xk ||2 − φki .
2 i


Barycenter estimations fi can finally be recovered from the optimal dual variables φi solution
of (C.7) as:
fi ∈ ∂Hνi (φi ), for i = 1 · · · n.
(C.8)
Following [18], one value of the above subgradient can be obtained at point xk as:
∂Hνi (φi )k =

pi
X
j j,k

νi Si ,

(C.9)

j=1

where Sij,k is any row stochastic matrix of size pi × N checking:
Sij,k



6= 0 iff k ∈ argmin
k=1···N

1 j
||y − xk ||2 − φki .
2 i


i
νij Sij,k corresponds to the discrete
From the previous expressions, we see that fik = pj=1
pushforward of data νi with the transport matrix Si with the associated cost:

P

Hνi (φi ) =

pi X
N 
X
1
j=1 k=1

2



||yij − xk ||2 − φki Sij,k νij .

Numerical optimization Following [21], the dual problem (C.7), can be simplified by
P
removing one variable and thus discarding the linear constraint K T φ0 + ni=1 φi = 0. In
order to inject the regularity given by φ0 in all the reconstructed barycenters obtained by
φi , i = 1 · · · n, we modified the change of variables of [21] by setting ψi = φi + K T φ0 /n for
P
i = 1 · · · n and ψ0 = φ0 , leading to ni=1 ψi = 0. One variable, say ψn , can then be directly
Pn−1
ψi /n, we thus obtain:
obtained from the other ones. Observing that φn = −K T ψ0 − i=1
min

ψ0 ,···ψn−1

n−1
X

T

T

Hνi (ψi − K ψ0 /n) + Hνn (−K ψ0 −

i=1

n−1
X

ψi /n) + Eγ∗ (ψ0 ).

(C.10)

i=1

The subgradient (C.9) can then be used in a descent algorithm over the dual problem (C.10).
For differentiable penalizers E, we consider the L-BFGS algorithm [43, 7] that integrates a
line search method (see e.g. [13]) to select the best time step τ (`) at each iteration ` of the
subgradient descent:
(

(`+1)

(`)

(`)

= ψ0 − τ (`) (∇Eγ∗ (ψ0 ) + d`0 )
ψ0
(`+1)
(`)
ψi
= ψi − τ (`) d`i
i = 1 · · · n − 1,

where:
d`0
d`i





Pn−1 (`)  Pn−1
(`)
(`)
(`)
= K ∂Hνn −K T ψ0 /n − i=1 ψi
− i=1 ∂Hνi ψi − K T ψ0 /n



Pn−1 (`) 
(`)
(`)
(`)
= ∂Hνi ψi − K T ψ0 /n − ∂Hνn −K T ψ0 /n − i=1 ψi
.

30

(C.11)

The barycenter is finally given by (C.8), taking φi = ψi − K T ψ0 /n. Even if we only treated
differentiable functions E in the theoretical part of this paper, we can numerically consider
non differentiable penalizers E, such as Total Variation (K = ∇, E = |.|1 ). In this case, we
make use of the Fista algorithm. This just modifies the update of ψ0 in (C.11), by changing
the explicit scheme involving ∇Eγ∗ onto an implicit one through the proximity operator of Eγ∗ :
(`+1)

ψ0



(`)
= Proxτ (`) Eγ∗ ψ0 − τ (`) d`0 = argmin
ψ

1
(`)
||ψ − τ (`) d`0 − ψ||2 + Eγ∗ (ψ).
2τ (`) 0

Algorithmic issues and stabilization As detailed in [18], the computation of one subgradientqin (C.9) relies on the look for Euclidean nearest neighbors between vectors (yij , 0) and

(xk , c − φki ), with c = maxk φki . Selecting only one nearest neighbor leads to bad numerical
results in practice as subgradient descent may not be stable. For this reason, we considered
the K = 10 nearest neighbors for each j to build the row stochastic matrices Si at each
0
P
iteration as: Sij,k = wijk / k0 wijk , with wijk = exp(−( 12 kyij − xk k2 − φki )/) if k is within the
K nearest neighbors for j and data i and wijk = 0 otherwise.

References
[1] M. Agueh and G. Carlier. Barycenters in the Wasserstein space. SIAM Journal on
Mathematical Analysis, 43(2):904–924, 2011.
[2] C. D. Aliprantis and K. Border. Infinite dimensional analysis: a Hitchhiker’s guide.
Springer Science & Business Media, 2006.
[3] P. Álvarez-Esteban, E. del Barrio, J. Cuesta-Albertos, and C. Matrán. Wide consensus
for parallelized inference. arXiv e-prints 1511.05350, 2015.
[4] L. Ambrosio, N. Gigli, and G. Savaré. Gradient flows: in metric spaces and in the space
of probability measures. Springer Science & Business Media, 2008.
[5] F. Bauer and A. Munk. Optimal regularization for ill-posed problems in metric spaces.
Journal of Inverse and Ill-posed Problems, 15(2):137–148, 2007.
[6] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear
inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
[7] S. Becker.
Matlab wrapper and c implementation of L-BFGS-B-C, 2011.
https://github.com/stephenbeckr/L-BFGS-B-C.
[8] J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré. Iterative Bregman projections for regularized transportation problems. SIAM Journal on Scientific Computing,
37(2):A1111–A1138, 2015.
[9] S. Bobkov and M. Ledoux.
One-dimensional empirical measures, order
statistics and Kantorovich transport distances.
To appear in the Memoirs of the Amer. Math. Soc., 2014.
Available at http://perso.math.univtoulouse.fr/ledoux/files/2013/11/Order.statistics.10.pdf.
31

[10] N. Bonneel, J. Rabin, G. Peyré, and H. Pfister. Sliced and radon Wasserstein barycenters
of measures. Journal of Mathematical Imaging and Vision, 51(1):22–45, 2015.
[11] A. W. Bowman and A. Azzalini. Applied smoothing techniques for data analysis : the
kernel approach with S-Plus illustrations. Clarendon Press ; Oxford University Press,
1997.
[12] S. Boyd and A. Mutapcic. Subgradient methods, Winter 2006-07. Notes for EE364
Standford University.
[13] S. Boyd and L. Vandenberghe. Convex optimization. Cambridge university press, 2004.
[14] A. Braides. A handbook of γ-convergence. Handbook of Differential Equations: stationary
partial differential equations, 3:101–213, 2006.
[15] M. Burger, M. Franek, and C.-B. Schönlieb. Regularized regression and density estimation based on optimal transport. Applied Mathematics Research eXpress, 2012(2):209–
253, 2012.
[16] L. A. Caffarelli. The regularity of mappings with a convex potential. Journal of the
American Mathematical Society, 5(1):99–104, 1992.
[17] L. A. Caffarelli. Boundary regularity of maps with convex potentials–ii. Annals of
mathematics, 144(3):453–496, 1996.
[18] G. Carlier, A. Oberman, and E. Oudet. Numerical methods for matching for teams
and Wasserstein barycenters. ESAIM: Mathematical Modelling and Numerical Analysis,
49(6):1621–1642, 2015.
[19] F. Clarke. Functional Analysis, Calculus of Variations and Optimal Control, volume 264
of Graduate Texts in Mathematics. Springer - Verlag London, 2013.
[20] M. Cuturi and A. Doucet. Fast computation of Wasserstein barycenters. In International
Conference on Machine Learning 2014, JMLR W&CP, volume 32, pages 685–693, 2014.
[21] M. Cuturi and G. Peyré. A smoothed dual approach for variational Wasserstein problems.
SIAM Journal on Imaging Sciences, 9(1):320–343, 2016.
[22] G. De Philippis and A. Figalli. The monge–ampère equation and its link to optimal
transportation. Bulletin of the American Mathematical Society, 51(4):527–580, 2014.
[23] A. Dessein, N. Papadakis, and J.-L. Rouas. Regularized Optimal Transport and the Rot
Mover’s Distance. ArXiv e-prints 1610.06447, Oct. 2016.
[24] P. Dupuis and R. S. Ellis. A weak convergence approach to the theory of large deviations.
Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons,
Inc., New York, 1997. A Wiley-Interscience Publication.
[25] S. Ferradans, N. Papadakis, G. Peyré, and J.-F. Aujol. Regularized discrete optimal
transport. SIAM Journal on Imaging Sciences, 7(3):1853–1882, 2014.
[26] N. Fournier and A. Guillin. On the rate of convergence in Wasserstein distance of the
empirical measure. Probability Theory and Related Fields, 162(3-4):707, Aug. 2015.
32

[27] D. Gervini. Independent component models for replicated point processes. Spatial Statistics, 18, Part B:474 – 488, 2016.
[28] Y. H. Kim and B. Pass. Wasserstein barycenters over Riemannian manifolds. arXiv
e-prints 1412.7726, 2014.
[29] A. Kneip and K. J. Utikal. Inference for density families using functional principal
component analysis. J. Amer. Statist. Assoc., 96(454):519–542, 2001. With comments
and a rejoinder by the authors.
[30] A. N. Kolmogorov and V. M. Tikhomirov. ε-entropy and ε-capacity of sets in function
spaces. Uspekhi Matematicheskikh Nauk, 14(2):3–86, 1959.
[31] T. Le Gouic and J.-M. Loubes. Existence and Consistency of Wasserstein Barycenters.
Probability Theory and Related Fields, pages 1–17, 2016.
[32] V. M. Panaretos and Y. Zemel. Amplitude and phase variation of point processes. Annals
of Statistics, 44(2):771–812, 2016.
[33] B. Pass. Optimal transportation with infinitely many marginals. Journal of Functional
Analysis, 264(4):947–963, 2013.
[34] K. Petersen and H.-G. Müller. Functional data analysis for density functions by transformation to a hilbert space. Annals of Statistics, To be published, 2015.
[35] G. Peyré, J. Fadili, and J. Rabin. Wasserstein active contours. In IEEE International
Conference on Image Processing (ICIP), pages 2541–2544. IEEE, 2012.
[36] R. Rockafellar. Conjugate duality and optimization, volume 16. Siam, 1974.
[37] A. W. Van Der Vaart and J. A. Wellner. Weak convergence and empirical processes.
Springer, 1996.
[38] C. Villani. Topics in optimal transportation, volume 58 of Graduate Studies in Mathematics. American Mathematical Society, 2003.
[39] C. Villani. Optimal transport: old and new, volume 338. Springer Science & Business
Media, 2008.
[40] W. Wu and A. Srivastava. An information-geometric framework for statistical inferences
in the neural spike train space. Journal of Computational Neuroscience, 31(3):725–748,
Nov. 2011.
[41] C. Zalinescu. Convex analysis in general vector spaces. World Scientific, 2002.
[42] Z. Zhang and H.-G. Müller. Functional density synchronization. Computational Statistics
& Data Analysis, 55(7):2234–2249, 2011.
[43] C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw.,
23(4):550–560, Dec. 1997.

33

