Computational Guarantees for Doubly Entropic Wasserstein Barycenters
Authors: Tomas Vaskevicius
Lénaïc Chizat
Venue: NeurIPS 2023 poster (2023)
Keywords: Wasserstein barycenters
entropic penalization
optimal transport
Sinkhorn's algorithm
TL;DR: We provide computational guarantees for computing doubly regularized entropic Wasserstein barycenters.

Abstract: We study the computation of doubly regularized Wasserstein barycenters, a recently introduced family of entropic barycenters governed by inner and outer regularization strengths. Previous research has demonstrated that various regularization parameter choices unify several notions of entropy-penalized barycenters while also revealing new ones, including a special case of debiased barycenters.  In this paper, we propose and analyze an algorithm for computing doubly regularized Wasserstein barycenters. Our procedure builds on damped Sinkhorn iterations followed by exact maximization/minimization steps and guarantees convergence for any choice of regularization parameters. An inexact variant of our algorithm, implementable using approximate Monte Carlo sampling, offers the first non-asymptotic convergence guarantees for approximating Wasserstein barycenters between discrete point clouds in the free-support/grid-free setting.

Computational Guarantees for Doubly Entropic Wasserstein Barycenters via Damped Sinkhorn Iterations
Tomas Vaškevičius and Lénaïc Chizat
Institute of Mathematics
École Polytechnique Fédérale de Lausanne (EPFL) Station Z
CH-1015 Lausanne
Switzerland
Abstract
We study the computation of doubly regularized Wasserstein barycenters, a recently introduced family of entropic barycenters governed by inner and outer regular- ization strengths. Previous research has demonstrated that various regularization parameter choices unify several notions of entropy-penalize also revealing new ones, including special case of debiased barycenters. In this paper, we propose and analyze algorithm for computing doubly regularized Wasserstein barycenters. Our procedure builds on damped Sinkhorn iterations followed by exact maximization/mmnimizzaii for any choice of regularization parameters. An inexact variant of our algorithm, implementable using approximate Monte Carlo sampling. offers the first non- asymptotic convergence guarantees for approximating Wasserstein barycenters between discrete point clouds in the free-suppor//gridffree setting.
barycenters while
steps and guarantees convergence
Introduction
The Wasserstein distance between two probability distributions measures the least amount of effort needed to reconfigure one measure into the other. Unlike other notions of distances based solely on the numerical values taken by the distribution functions (e.g.. the Kullback-Leibler divergence), the Wasserstein distance incorporates an additional layer of complexity by considering pairwise distances between distinct points, measured by some predetermined cost function. As result, the Wasserstein distances can be seen to lift the geometry of the underlying space where the probability measures defined the space of the probability measures itself. This allows for more thorough and geometrically nuanced understanding of the relationships between different probability measures, which proved to be versatile tool of increasing importance
broad spectrum of areas.
Given a collection of probability measures and an associated set of positive weights that sum to one, the corresponding Wasserstein barycenter minimizes the weighted sum of Wasserstein distances to the given measures. In the special case of two measures and the squared Euclidean cost function, Wasserstein barycenters concide with the notion of McCann's displacement interpolation introduced in the seminal paper [42). The general case, encompassing an arbitrary number measures, was first studied by Agueh and Carlier [1]. where they also demonstrated a close link between Wasserstein barycenters and the multi-marginal optimal transport problem 229] Recent years have witnessed an various scientific disciplines. See,
increasing number of applications of Wasserstein barycenters
for instance, the following sample of works in economics [16] [1] statistics
[8]. image processing
[50] and machine learning [19], among other areas. For further background and references point
the interested reader to the introductory surveys 47
45 and the textbooks [58 59
54
46
28
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Despite their compelling theoretical characteristics, the computation of Wasserstein barycenters poses significant computational challenges, particularly in large-scale applications. While Wasserstein barycenters can be computed in polynomial time for fixed dimensions [3] the approximation of Wasserstein barycenters known to be NP-hard 44. Currently employed methods for approximating Wasserstein barycenters are predominantly based on space discretizations. Unfortunately, such strategies are only computationally practical for problems of relatively modest scale. Although there are a handful of grid-free techniques available for approximating Wasserstein barycenters (e.g.. 18 35 23 40), we are not aware of any existing methods that provide bounds on computational complexity. One contribution of the present paper is to introduce method that in some regimes can provably approximate Wasserstein barycenters without relying on space discretizations, but instead employing approximate Monte Carlo sampling.
More broadly, the difficulties associated with computation of the optimal transport cost has prompted the exploration of computationally efficient alternatives, leading to the consideration of regularized Wasserstein distances. Among these, the entropic penalty has emerged as one of the most successful in applications. The practical success of entropic penalization can be attributed to Sinkhorn's ' algorithm (56), which enables efficient and highly parallelizable computation, an algorithm that gained substantial traction in the machine learning community following the work of Cuturi [20]. It is worth noting that entropic Wasserstein distances are of intrinsic interest, beyond their approximation capabilities. Indeed, they hold rich historical connection to the Schrödinger bridge problem [55] 60 26), as highlighted in the recent surveys [39 [4]. Furthermore, they increasingly serve as an analytically convenient tool for studying the unregularized optimal transport problem (see, e.g..
38 31 27 15 and they underlie some favorable statistical properties that are currently under active investigation; see the works 43 30
52 49] and the references therein.
Let us now define the entropic optimal transport cost. Consider two probability measures, l and v, both supported on X' and let e X x X - [0, oo) be a cost function. The entropic Wasserstein distance with a regularization level À > 0 is defined as
(x, inf E(x,,~)(((xx vEll(u.))
H &)),
(1)
where [[, 2) denotes the set of probability measures on X' x X' with marginal distributions equal to fi and v/ and -) is the Bulllbak--eeible divergence. When À + 0, the regularized cost Tx (H, v) converges to the unregularized Wasserstein distance. Various properties of entropic optimal transport can be found in the recent lecture notes by Léonard (39]
To develop efficiently computable approximations for Wasserstein barycenters, natural approach is to replace the unregularized Wasserstein cost with the minimizer of the weighted sum of entropy- regularized costs. This method was first explored by Cuturi and Doucet [2D and it has gained additional traction in the recent years. There is some flexibility in the definition of (T). which arises from substituting the reference product measure μ with alternatives such as the Lebesgue measure. Consequently, various notions of entropic barycenters have emerged in the literature, which can be unified through the following optimization problem:
min
wy wTTx(H,, + TKT,jU, Trot).
(2)
Here pt...., w are the probability measures whose barycenter we wish to compute and Wi..... We are positive weights that sum to one. The inner regularization strength is denoted by >> while >> 0 is the outer regularization strength. The measure Tref is an arbitrary reference measure, the support of which dictates the support of the computed barycenter. For instance, if we take Tyef to be uniform measure on particular discretization of the underlying space, are dealing with a fixed- support setup. On the other hand, letting Ryet be the Lebesgue measure puts us in the fure-supportt setup. We shall refer to the minimizer of as the -baryceee,,, which exists and is unique due to the strict convexity of the outer regularization penalty; however, uniqueness may no longer holds when 5 =0. 0.
The objective 2 was recently studied in 17 it also appeared earlier in [5] for the special case r 2 X, where stochastic approximation algorithms were considered for the computation of fixed-support barycenters. In [17 Section 1. 1.3], it discussed how various choices of (λ, relate to Barycenters previously explored in the literature. To provide brief overview, (0, 0) are the unregularized

---

Wasserstein barycenters studied in M. renerrregulrizzed barycenters (A, 0) introduce shrinking bias; this can be seen already when & 1. in which case the solution computes maximum-likeliie likelihood
they introduce a blurring bias. Likewise, blurring bias is introduced by the outer-regularized barycenters (0, r), studied in 19 [2]. The only case not covered via the formulation ( appears to be the one of debiased Sinkhorn barycenters [5] 33), for which an algorithm exists but without computational guarantees. Of particular interest are the (x, X/2) barycenters: the choice = A/2 for smooth densities yields approximation bias of order X², while the choice т = À results in bias of order X, which is significantly larger than X in the regimes of interest. This is
deconvolution [53] The (X, A)-barycenters were considered in 21 7 22
I0 36
notion of entropic barycenters that was unveiled in the analysis of 17 We provide the first convergence guarantees for this type of barycenters.
The regularity, stability, approximation, and statistical sample complexity properties of (A, r)- barycenters were investigated in [17] However, the question of obtaining non-asymptotic convergence guarantees for the computation of (A, 7 -baryceeets with arbitrary regularization parameters was not addressed therein. In particular, the (X, x/2) case, which has stood out due to its compelling mathematical features, has not yet been addressed in the existing literature. This gap is addressed by the present paper, we summarize our contributions in the following section.
1.1 Contributions
The remainder of this paper is organized as follows: Section provvdess the necessary background on entropic optimal transport and a particular dual problem of the doubly regularized entropic objective (2) Section introduces damped Sinkhorn iteration scheme and complements it with convergence guarantees. An approximate version of the algorithm together with convergence results and implementation details is discussed in Section 4
We summarize key contributions:
1. Lemma presented in Section demonstrates that bounds on the dual suboptimality gap for the dual problem (8). defined in Section 2.2] can be translated into Kullback-Leible divergence bounds between the (X, -baayyceeee and the barycenters corresponding to dual-feasible variables. This translation enables us to formulate all our subsequent results in terms optimizing the dual objective
2. In Section 3 we introduce a damped Sinkhorn scheme (Algorithm I that can be employed to optimize (A, 7-)-bryycnnerrs for any choice of regularization parameters. The damping factor min(1, T 7/) accommodates the degrading smoothness properties of the dual objective (8) as a function of decreasing outer regularization parameter r. The introduced damping the Sinkhorn iterations is, in fact, necessary and it is one of our core contributions: undamped
exact scheme can be experimentally shown to diverge as soon as 7 < X/2.
3. The main result of this paper is Theorem 1 proved in Section 3 It provides convergence
guarantees for Algorithm [ with arbitrary choice of regularization parameters X, > This, in particular, results in the first algorithm with guarantees for computing (A, A/2) barycenters. For smooth densities, these barycenters incur a bias of order X² in contrast to the predominantly studied (X, X) barycenters that incur bias of order X.
4. In Section 4 we describe Algorithm an extension of Algorithm that allows us to perform inaccurate updates. We formulate sufficient conditions the inexact updates oracle under which the errors in the convergence analysis do not accumulate. Section 4.1 details an implementation of this inexact oracle, based on approximate Monte Carlo sampling.
5. Theorem proved in Section furnishes convergence guarantees for Algorithm 2 When combined with the implementation of the inexact oracle described in Section 4.1 this yields a provably convergent scheme for grid-free computation of entropic Wasserstein barycenters between discrete distributions, provided sufficient regularity on the domain X and the cost function c.
6. Appendix F complements our theoretical results with numerical experiments. Our simula- tions experimentally confirm the necessity of damping when T < A/2. They also provide experimental support for the suggested damping factor in Algorithms I and
3

---

2 Background and Notation
This section provides the background material on doubly regularized entropic Wasserstein barycenters and introduces the notation used throughout the paper. the remainder of the paper, let be compact and convex subset of R with non-empty interior. Let P(X') denote the set of probability measures on X' endowed with Borel sigma-algebra. Let e X X - [0, oo) be cost function such that Coo (x) (X) = S c(x,x") < co. We denote by KT(:,)) the Kullback-Leible divergence, || llv is the total-variation norm, and || sup, f() inf f(x') is the oscillation norm. Given two measures v, v', the notation v < v' denotes that v is absolutely continuous with respect to the measure v': in this case dv/dv' denotes the Radon-Nikodys derivative of with respect to v' Finally, throughout the paper w denotes a vector of l strictly positive elements that sum one.
2.1 Entropic Optimal Transport
For any n, (()) define the entropy regularized optimal transport problem by
Tx(,,)) = inf B(x,Y)))(((x,) Y)) + AKT( v); yell(a.))
(3)
where KL is the Kullback-Leible divergence and (T v) S P(X @ X)) is the set of probability measures such that for any 7 [[,,,)) and any Borel subset A of X' it holds that y( X) ((A)
and 7(X' x A) v(A).
Let " (I))) Li(v) + R be the function defined by
Ex (Ø. x²) = Ex~rlo(X)] + Ey-((v(())
+ > 1 / / exp
x(x) + ((y) v) u(dy)u((r))
x x
The entropic optimal transport problem 3
admits the following dual representation:
(x,, max
F. (d. 16).
0.6
For any tp define
de argmaxpe Ex (o, v)
The solution is unique p-almost everywhere to a constant; we fix particular choice
(x = log (Je
vb(y)
c(xr, w) u(dy)
exp
Likewise, we denote 2h argmax,ee E. (d) v)) with the analogous expression to the one given above, interchanging the roles of Φ and vb. Then, the maximum in (4) attained by any pair (o", si*) such that 6* = 6y- and v¹* = u such pair is said to solve the Schrödinger system and is unique up translations (6* + a, J³ a) by any constant e The optimal coupling that solves the primal problem 3 can be obtained from the pair (Ø" .b*) via the primal-dual relation
Φ (x)
+ (2) c(x, v) u(dx)v(dy))
7" (dx, dy) exp
conclude this section by listing two properties of functions of the form dy. These properties will be used repeatedly throughout this paper. First, for any ib we have
Je exp (o(()) + vbv) cxa, vl u(dy)u(do) =
A
which means, in particular, that for any 1 we have
Ex (ow, v) Exxxlox(X)) (x)m))((())
(5)

---

The second property of interest is that for any eb and any 2, x' E A' it holds that
fexp
(()))
v(dy)
(x(x) 20(x'))
-Alog
_ exp
((M) olar 20
v(ddy
I eep
(/M) (0+c(2' -/ - v(dy)
A
=ollo
f exp
((v)-c(x'' r(dy)
S sup c(l', ,11 c(xx (.) yE
partieular, for any vi we have
11 sup (x inf (hs)s()) s co0 (2).
2.2 Doubly Regularized Entropic Barycenters
Let = (v'...., w*) E ( k probability measures and w IR be a vector of positive numbers that sum to one. Given the inner regularization strength À > 0 and the outer regularization strength r > 0, the (λ, r) barycenter μX,r E P(x) of probability measures v with respect to the weights vector w is defined as the unique solution the following optimization problem:
Σ wy wyTx(H, TKL(M, Rref):
PA,r BargminJuP(X
where ret E P(2') is a reference probability measure.
We will now describe how to obtain concave dual maximization problem to the primal problem (7). following along the lines of Chizat [17 Section 2.3], where the interested reader will find comprehensive justification
junprreeosiveenssifffaaatinn of all the claims made in the rest of this section.
First, using the semi-dual formulation of entropic optimal transport problem j E {1....,k}
have, for each
Tx(u, = sup Ex~ulóps (X)] + Ey [ (Y)).
Denote vb = (+b2 ...., ()) E LI(D). Then, we may rewrite the primal problem 7 by
min max (EP(X)
wy (v) + Ex
wydes (X)] + TKL(U, Tret).
Interchanging min and max, which is justified using compactness of X' as detailed in [17 we obtain the dual optimization objective E w (vv)) + RR defined by
EX, (ub) = min pEP(X)
Σ wy ByyEy-u [2 + Ex 1=1
Σ
Wj widws + 7KT, Tret).
1=1
(8)
Σ wy Ey (V)] - f exp
Ee widde(
Tree
T
The infimum above is attained by the measure
E/=ı Oaps (x)
( -ー Pev Pvo((z)
He (dx) = "wy exp
Tree (der, Z = exp
To (dx).
To each dual variable s/ we associate the marginal measures u (dy) defined for j = 1,..., by
/ exp (ess vid (v) c(s,W) u)
v (dy) = v ('dy)
ua (dx).
(9)
A
5

---

Finally, we mention that the objective Ex7 is concave and for any 2b, 2 it holds that
Ex (eb + hab') - 1 (2b)
k
lim h-
Σ wj
- -E (w))))
h
In partieular, fixing any optimal dual variable sb", for any ub holds using concavity of Ex" that
S (x (1b))
(Bb
Wk Enw [( wi1 Eig ((00)--a01)) v (10)
This concludes our overview of the background material on (bayy-bbryeenttes.
3 Damped Sinkhorn Scheme
This section introduces a damped Sinkhorn-based optimization scheme (Algorithm D and provides guarantees for its convergence (Theorem TD. Before describing the algorithm, make quick detour to the following lemma, proved in Appendix A which shows that the sub-optimalitt gap bounds on the dual objective (8) can be transformed into corresponding bounds relative entropy between the (A, T) --baaycceteer and the barycenter associated to a given dual variable.
Lemma 1. Fix any X, > and v, w. Let up" be the maximizer of dual problem Ex and let pe be the corresponding minimizer of the primal objective
7 Then, for any ab E L1(v) we have
KL S T (EXx, (wb) Ex, (+))).
We now turn to describing an iterative scheme that ensures convergence of the dual suboptimality gap zero. Let th, be an iterate time Then, we have
& wy En Ep, + TKI(H,, Fret),
Ex (tbe) L( 0.) E
j=1
where
d argmax 1. E (d, vl() and 14 argmin
E, Σ wy + +KKL(1, Teet) }. (11)
In particular, when optimizing the dual objective EX, every time the variable zht is updated, automatically triggers the exact maximizationmniiniizztoon steps defined in (TT). is thus a natural strategy to fix Φ; and de and perform exact minimization on sb, which can be done in closed form:
w argmax, Ex (ol, v)) = v/ log duf du
(12)
where of denotes the marginal distribution vi update on each block of variables wJ Together, the update 12 followed by [I] results in the iterative Bregman projections algorithm introduced in [7]. In [6]. it was shown that this scheme converges for the (X, -JJbaryceettrs. The analysis of [36] is built upon different dual formulation from the one considered in our work; this alternative formulation is only available when 7 À [17 Section 2.3] and thus excludes the consideration of debiased barycenters (A,
vh defined in 9 The update (22 performs a Sinkhorn
We have observed empirically (see Appendix F that the iterates of the iterative Bregman projections (i.e., the scheme of updates defined (12) and (1T)) diverge whenever 7 < A/2. Indeed, decreasing the outer regularization parameter makes the minimization step in (11) less stable. As result, the
cumulative effect of performing the updates (12) and II] may result in decrease in the value of the
optimization objective "
One of the main contributions of our work is to show that this bad behaviour can be mitigated by damping the exact Sinkhorn updates (12). This leads to Algorithm | for which convergence guarantees are provided in Theorem stated below.

---

Algorithm 1:
Exact Damped Sinkhorn Scheme
Input:
regularization strengths A; O, reference measure Teef TTref, number of iterations 7 and marginal measures v' .....
: v with positive weights w1,.... wk such that Σ Wy 1.
1. Set 1 = min(1, m and initialize (v%)) = 0 for j E fl.....k).
2. For t = 0, T T- do
(a) (x(x) t looo fx exp((V/ (v) c(z, v))/A)/ for j E {l....,k) k
(b) Vi(x) -jj= wgy (z)
(c) 7, [eep(- (xxp(-V((x:)) (T)
(d) mx(dx) + Z exp( 1ap(-Vo
(c) (v) + sep exp ( (1) pe pe(dx) for j E {l....,k}.
de: ( (o-1(0) - (1) log di
(v) for j (1....,k).
3. Return (0 トー
Theorem 1. Fix any X, and v, w. Let 1h*t be the maximizer of dual problem Ex,"" Let (sbr)e2o be the sequence of iterates generated by Algorithm [ Then, for any 2 it holds that
Ex (xb")))))))))))) Ex," (vb. 2coo (2') (x))
min(>, T)
Our convergence analysis draws upon the existing analyses of Sinkhorn's algorithm [2 25 which
in turn are based on standard proof strategies in smooth convex optimization (e.g., 44 Theorem
2.1. Concerning the proof of Theorem I the main technical contribution our work lies in the following proposition proved in Appendix
Proposition 1. Consider the setup of Theorem Then, for any integer 2 0 it holds that
- EX (+b.) 2 min (T,
Σ
wy W,KKL("'
j=1
With Proposition at hand, we are ready to prove Theorem
Proof of Theorem 7 Denote 8, Ex,"" (xb*) EXP (xbe). We would like to relate the suboptimality
gap & to the increment & - -.+1-- To do this, we will first show that the iterates t/, have their oscillation norm bounded uniformly in f.. Indeed, for any j E {1.....k) k). any t 2 1, and any y E X we have
(1) (1 7) (v) + m (y).
By (6) dia has oscillation norm bounded by c co)) (X).. Because vid 0 and n E (0, 1]. by induction on t it follows that ||ebelosec S Coo (X) for any t 2 0. Combining the bound on the dual sub-optimality gap (0)) with Pinsker's inequality yields
8 < 2c (X)
Σ Ewsllw v/llvv V2c Σ wy KL
Using concavity of the square root function, Proposition Tyields for any t 2 0
min X T) 52
8 8 2 min(A,
T)
Σ wy KyKKLV [ Z 2c (x))²
1=1
By Proposition [ the sequence 6, is non-increasing. Hence, dividing the above equality by d,de+1 yields
mi) 20 (2')²
8
5,
Telescoping the left hand side completes the proof.
7

---

4 Approximate Damped Sinkhorn Scheme
In this section, we extend the analysis of Algorithm to an approximate version of the algorithm. Then, in Section 4.1 describe how inexact updates may be implemented via approximate random sampling, thus enabling the computation (A, -bbrrycettes in the free-support setting with convergence guarantees.
Algorithm 2 describes an inexact version of Algorithm TT] replaces the damped Sinkhorn iterations of Algorithm|| via approximate updates computed by an approximate Sinkhorn oracle a procedure that satisfies the properties listed Definition []
Definition (Approximate Sinkhorn Oracle). An E-approximate Sinkhorn oracle is procedure that
......,k), k). returns a Radon-Nikodym derivative d7d
given any 1/ and any index j E {L...., b
of a measure ded
v that satisfies the following properties:
do
1.
is strictly positive on the support of v²;
2. v ITV
S // (20o0
3. Ey wus
(Y) S 1 s² /2exx (2)
4. For any 1) E [0 and any E {l,...,k) it holds that 1/ 1/d mnllog(dle
/
S
(1 1-)1 ( loaee + 7/ (A)).
Algorithm 2: Approximate Damped Sinkhorn Scheme
Input: error tolerance parameter e > a function "ApproximateSmnhhoorrree properties listed in Definition]] [ regularization strengths X, T >> 0, reference measure Tref. of iterations T and marginal measures y1 ...., 5/ with positive weights
satisfying
Wi...., Wy such that Σ/= "y = 1.
1. Set 17 min(1, T/A) and initialize (v%) for j (1....,k)..
2. For / = 0,
T do
(a) di (v) + ApproxinateSinhoonOraede X, the, e, for j E (l....,k).
(b) (y) +- vi(()) - mpr di (v) for (1....,,), k].
3. Return (5r, v 1
The following theorem shows that Algorithm 2 enjoys the same convergence guarantees as Algo- rithm up to the error tolerance of the procedure used to implement the approximate updates. A noteworthy aspect of the below theorem is that the error does not accumulate over the iterations.
Theorem 2.
Fix any 7, > O and v, Let th be the maximizer of dual problem F. Let (4h) )(20 be the sequence of itterates generated by Algorithm 2 with the accuracy parameter 2 0. Let
T = [1: (xb") FX (2b.) < 2e). Then, for any 1. < T it holds that
Ex (1b*) - -x (xb.) 28 + 20 (X) 1
min (A, r) [
The proof of the above theorem can be found in Appendix C
4.1 Implementing the Approximate Sinkhorn Oracle
In this section, we show that the approximate Sinkhorn oracle (see Definition 1 can implemented using approximate random sampling when the marginal distributions v are discrete. To this end, fix

---

the regularization parameters X, 7 > 0, the weight vector w, and consider set of & discrete marginal distributions
my
v
(or (8)))
where 6 is the Dirac measure located at and 1 (vt)) is equal the probability of sampling the
point of from measure v², denote the total cardinality of the support all measures yd by
myj
Fix any 3 L1 (v).. Suppose we are given access to n iid. samples XI....,.8 X, from a probability measure 1oe that satisfies
I H S E
Then, for j == == 1,..., and l 1, my consider
DJ (v(() = (vt) 1 Σ
Dow (Xi) + v (v) c(x,
exp
A
1=1
and for any parameter Ç E (0, 1/2] define
. = -1- (1-5) Cv².
(13)
We claim that io implements the approximate Sinkhorn oracle with accuracy parameter of order ( 1/44 provided that is large enough. This is shown in the following lemma, the proof of which can be found in Appendix(] [
Lemma 2. Fix any 5 (0, 1) and consider the setup described above. With probability at least 1 8, for each j E [1.....k] it holds simultaneously that the measure 70 defined in (3)) satisfies all the properties listed in Definition 7 with accuracy parameter
1 { Ej cx (2) 25 + mneen +
1
log
2m
2/2
The above lemma shows that a step of Algorithm Zean can be implemented provided access i.i.d. sampling from some measure 1'g close to Ps in total variation norm, where vp is an arbitrary iterate of Algorithm 2 The remainder of this section is dedicated to showing that this can be achieved by sampling via Langevin Monte Carlo.
Henceforth, fix #ee to be the Lebesgue measure on X, which corresponds to the free-support barycenters setup. Then, for any 2b we have
(dx 1x exx(((VVV
where Vo(x))
Σ w
where x equal to one on X' and zero everywhere else. follows by (6) that || Ioore S (20) Further, let diamx sup,,''xx' ||x z'11 By the convexity of X. the uniform measure on X satisfies the logarithmic Sobolev inequality (LSI) with constant d 2 /4 (cf. (37). Hence, by the Holley-Stroock perturbation argument 32 the measure μe satisfies LSI with constant at most exp (2cge (x)/r) diam(X) 2// sc.
It is well-established that Langevin Monte Carlo algorithms offer convergence guarantees for ap- proximate sampling from a target measure subject to functional inequality constraints provided additional conditions hold such as the smoothness of the function V However, such guarantees do not directly apply the measure Mo due to its constrained support. Instead, it is possible to approximate Mv arbitrarily well in total variation norm by family of measures (oo.a)os0 (see Appendix
F for details) supported on all of Rd, Tuning the parameter allows us to trade-off between the approximation quality of Phb.o and its L.ST constant. Crucially, standard sampling guarantees for Langevin Monte Carlo (e.g., (57)) apply the regularized measures Pg..s which leads to provable

---

guarantees for an implementation of Algorithm 2 thus furnishing the first convergence guarantees for computation of Wasserstein barycenters in the free support setup; see Theorem stated below.
The above approximation argument applies to any cost function e that is Lipschitz on X' and exhibits quadratic growth at infinity. For the sake of simplicity,
consider the quadratic cost c(z, y)
|| - y vl.. The exact problem setup where we are able to obtain computational guarantees for free-support barycenter computation via Langevin Sampling is formalized below.
Problem Setting 1. Consider the setting described at the beginning of Section 4.1 In addition, suppose that
1. the reference measure (rer(dx) 1x is the Lebesgue measure supported on X' (free- support setup);
2. it holds that C Bn = (xr Rd | - RJ for some constant R < oo;
3. the cost function e: Rd x R + [0, oo) is defined by c(z, y) = || y1l::
4. for any &/ we have access to a stationary point Vy of V over X.
The last condition can be implemented in polynomial time using first order gradient method. For our purposes, this condition is needed to obtain a good initialization point for the Unadjusted Langevin Algorithm following the explanation in [57 Lemma I]; see Appendix E for further details.
We now proceed to the main result of this section, the proof of which can be found in Appendix E The following theorem provides the first provably convergent method for computing Wasserstein barycenters the fupe-support setting. We remark that stochastic approximation argument of a rather different flavor used to compute fixed-support Wasserstein barycenters (for T 2 X) has been previously analyzed in [5]
Then, for any confidence parameter E (0, 1) and any accuracy parameter ε >> O, we can simulate step of Algorithm 2 with success probubility least 5 in time polynomial in
Theorem 3. Consider the setup described in Problem Setting 7
-, d, R, exp /T), (Rd-1/44) -
m, log(m/5))
In particular, an e-approximation of the (X, -BBrryyyeet can be obtained within the same compula- tional complexity.
Comparing the above guarantee with the discussion following the statement of Lemma 2 we see an additional polynomial dependence on (Rd-1/1yd (note that for R s dz/1 this term disappears). We believe this term to be artefact of our analysis appearing due the approximation argument described above. Considering the setup with R < d'/4, the running time of our algorithm depends exponentially in R² /T.
We conclude with two observations. First, since approximating Wasserstein barycenters is generally NP-hard 4 an algorithm with polynomial dependence on all problem parameters does not exist P A NP. Second, notice that computing time polynomial
approximation of 7))ryBayeeeet can be done in e-1. This should be contrasted with numerical schemes based on discretizations of the set X, which would, in general, result in computational complexity of order (R/s)* to reach the same accuracy.
Conclusion
We introduced algorithms compute doubly regularized entropic Wasserstein barycenters and studied their computational complexity, both in the fixed-support and in the free-support settings. Although a naive adaptation of the usual alternate maximization scheme from [7] to our setting leads to diverging iterates (at least for small values 7), our analysis shows that it is sufficient to damp these iterations to get converging algorithm.
While we have focused on the problem barycenters measures, we note that the idea of entropic regularization is pervasive in other applications of optimal transport. There, the flexibility offered by the double entropic regularization may prove to be useful as well, and we believe that our damped algorithm could be adapted to these more general settings.
10

---

References
[1] M. Agueh and G. Carlier. Barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis, 43(2):04-924,, 2011.
[2] J. Altschuler, J. Niles-Weed, and P. Rigollet. Near-linear time approximation algorithms for optimal transport via sinkhorn iteration, Advances neural information processing systems, 30, 2017.
[3] J. M. Altschuler and E. Boix-Adsera. Wasserstein barycenters can be computed in polynomial time in fixed dimension. The Journal of Machine Learning Research, 22(1):2000-200 2021.
[4] M. Altschuler and E. Boix-Adsera. Wasserstein barycenters are np-hard to compute. SIAM Journal on Mathematics of Data Science, 4(1): 2022.
[5] M. Ballu, Q Berthet, and F. Bach. Stochastic optimization for regularized Wasserstein estima- tors. International Conference on Machine Learning, pages 602-612. PMLR, 2020.
[6] H. H. Bauschke and L. Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. Springer, 2017.
[7] 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.
[8] E. Bernton, P. E. Jacob, M. Gerber, and C. Robert. On parameter estimation with the Wasserstein distance. Information and Inference: A Journal of the IMA, 8(4):657-676, 2019.
[9] J. Bigot, E. Cazelles, and N. Papadakis. Data-driven regularization of Wasserstein barycenters with an application to multivariate density registration. Information and Inference: Journal of the IMA, 8(4):719-7555 2019.
[10] J. Bigot, E. Cazelles, and N. Papadakis. Penalization of barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis, 51(3):2261-228 2019.
[11] G. Carlier and Ekeland. Matching for teams. Economic theory, pages 397-418, 2010.
[12] G. Carlier, K. Eichinger, and A. Kroshnin. Antropie-aassense barycenters: Pde charac- terization, regularity, and clt. SIAM Journal on Mathematical Analysis, 53(5):5880-53) 5914, 2021.
[13] Cattiaux, A. Guillin, and L.-M. Wu. A note on talagrand's transportation inequality and logarithmic sobolev inequality. Probability theory and related fields, 148:285-344 1088500344444 2010.
[14] Y. Chen, T. T. Georgiou, and M. Pavon. Stochastic control liaisons: Richard sinkhorn meets gaspard monge on Schrodinger bridge. Siam Review, 63(2):249-313, 2021.
[15] S. Chewi and A-A. Pooladian. An entropic generalization of caffarelli's contraction theorem via covariance inequalities. arXiv preprint arXiv:2203, 204442200304454, 2022.
[16] P.-.. Chiappori, R. J. McCann, and L. P. Nesheim. Hedonic price equilibria, stable matching, and optimal transport: equivalence, topology, and uniqueness. Economic Theeor, pages 317-354, 2010.
[17] L. Chizat. Doubly regularized entropic Wasserstein barycenters. arXiv:2303. 1844, 2023.
arXiv preprint
[18] S. Cohen, M. Arbel, and M. P. Deisenroth. Estimating barycenters of measures in high dimensions. arXiv preprint arXiv:2007.07105, 2020.
[19] N. Courty, R. Flamary, A. Habrard, and A. Rakotomamonjy. Joint distribution optimal trans- portation for domain adaptation. Advances in neural information processing systems, 30, 2017.
[20] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.

---

[21] M. Cuturi and A. Doucet. Fast computation of Wasserstein barycenters. In International conference on machine learning, pages 685-693. PMLR, 2014.
[22] M. Cuturi and G. Peyré. Semidual regularized optimal transport. SIAM Review, 60(4):941-965, 2018.
[23] C. Daaloul, T. L. Gouic, J. Liandrat, and M. Tournus. Sampling from the Wasserstein barycenter. arXiv preprint arXiv:2105.01706 2021.
[24] E. del Barrio and J.-M. Loubes. The statistical effect of entropic regularization in optimal transportation. arXiv preprint arXiv:2006.05719, 2020.
[25] Dvurechensky, A. Gasnikov, and A. Kroshnin. Computational optimal transport: Complex- ity by accelerated gradient descent is better than by sinkhorn's 's algorithm. In International conference on machine learning, pages 1367-1376. PMLR, 2018
[26] S. Erlander and N. F. Stewart. The gravity model in transportation analysis: theory and extensions, volume 3. Vsp, 1990.
[27] M. Fathi, N. Gozlan, and M. Prod'hommes A proof the caffarelli contraction theorem via entropic regularization. Calculus of Variations and Partial Differential Equations, 59:1-18. 2020.
[28] A. Figalli and F. Glaudo. An Invitation to Optimal Transport, Wasserstein Distances, and Gradient Flows. European Mathematical Society, 2021.
[29] W. Gangbo and A. Swięch. Optimal maps for the multidimensional monge-kantorovich problem. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 51(1):23-45, 1998.
[30] A. Genevay, L. Chizat, F. Bach, M. Cuturi, and G. Peyré. Sample complexity of sinkhorn divergences. In The 22nd international conference artificial intelligence and statistics, pages 1574-1583. PMLR, 2019.
[31] I. Gentil, C. Léonard, and L. Ripani. About the analogy between optimal transport and minimal entropy. In Annales de Faculté des sciences de Toulouse: Mathématiques, volume 26, pages 569-600, 2017.
[32] R. Holley and D. W. Stroock. Logarithmic sobolev inequalities and stochastic ising models, 1986.
[33] H. Janati, M. Cuturi, and A. Gramfort. Debiased sinkhorn barycenters. In International Conference on Machine Learning, pages 4692-4701. PMLR, 2020.
[34] H. Janati, B. Muzellec, G. Peyré, and M. Cuturi. Entropic optimal transport between unbalanced gaussian measures has a closed form. Advances in neural information processing systems, 33: 10468-1000 10479, 2020.
[35] A. Korotin, L. Li, J. Solomon, and E. Burnaev. Continuous Wasserstein-2 barycenter estimation without minimax optimization. arXiv preprint arXiv:2102.01752, 2021.
[36] A. Kroshnin, N. Tupitsa, D. Dvinskikh, P. Dvurechensky, Gasnikov, and C. Uribe. On the complexity of approximating Wasserstein barycenters. In International conference on machine learning, pages 3530-3540 PMLR, 2019.
[37] J. Lehec. The langevin monte carlo algorithm in the non-smooth log-concave case. arXiv preprint arXiv:2101.10695, 2021.
[38] C. Léonard, From the Schrödinger problem to the nonge-kantorovich problem. Journal of 2012.
Functional Analysis,
[39] C. Léonard. survey of the Schrödinger problem and some of its connections with optimal transport. arXiv preprint 3rXiv:1308.0215.
2013. 2013.
12

---

[40] J. v. Lindheim. Simple approximative algorithms for free-support Wasserstein barycenters.
Computational Optimization and Applications, pages 1-34, 2023.
[41] A. Mallasto, A. Gerolin, and H. Q. Minh. Entropy-regularizeed 2-wasserstein distance between gaussian measures. Information Geometry, 5(1):289-323) 2022.
[42] R. J. McCann. A convexity principle for interacting gases. Advances in mathematics, 128(1): 153-179, 1997.
[43] G. Mena and J. Niles-Weed. Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem. Advances in Neural Information Processing Systems, 32, 2019.
[44] Y. Nesterov al. Lectures on convex optimization, volume 137. Springer, 2018.
[45] V. M. Panaretos and Y. Zemel. Statistical aspects of Wasserstein distances. Annual review of statistics and its application, 6:405-431, 2019.
[46] V. M. Panaretos and Y. Zemel. An invitation statistics in Wasserstein space. Springer Nature, 2020.
[47] G. Peyré, M. Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11 1165-6)3355607, 2019.
[48] A.-A Pooladian and J. Niles Weed. Entropic estimation optimal transport maps. arXiv preprint arXiv:2109. 12000 2021.
[49] A. Pooladian, V. Divol, and J. Niles-Weed. Minimax estimation of discontinuous optimal transport maps: The semi-discrete case. arXiv preprint urXiv:2301.11300 2023.
[50] Rabin, G. Peyré, J. Delon, and M. Bernot. Wasserstein barycenter and its application to texture mixing. In Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, I, Ein-Gedi, Israel, May 29-June 2, 2011, Revised Selected Papers 3, pages 435 446. Springer, 2012.
[51] A. Ramdas, N. Gareía Trillos, and M. Cuturi. On Wasserstein two-sample testing and related families of nonparametric tests. Entropy, 19(2):47, 2017.
[52] P. Rigollet and A. J. Stromme, On the sample complexity of entropic optimal transport. arXiv preprint arXiv:2206. 13472, 2022.
[53] Rigollet and J. Weed. Entropic optimal transport maximum-likeelhood deconvolution. Comptes Rendus Mathematique, 356(1 (-1::11111):222-22222 2018.
[54] F. Santambrogio. Optimal transport for applied mathematicians, Birkäuser, 2015.
[55] E. Schrödinger. Sur la théorie relativiste de l'électron et e l'interprétation de la mécanique quantique. In Annales de I l'institut Henri Poincaré, volume 2, pages 269-310, 1932.
[56] R. Sinkhorn. Diagonal equivalence to matrices with prescribed row and column sums. The American Mathematical Monthly, 74(4):402-405, 1967.
[57] S. Vempala and A. Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. Advances in neural information processing systems, 32, 2019.
[58] C. Villani. Topics in Optimal Transportation. American Mathematical Society, 2003.
[59] C. Villani et al. Optimal transport: old and new. Springer, 2009.
[60] A. G. Wilson. The of entropy maximising models, in the theory of trip distribution, mode split and route split. Journal of transport economics and policy, pages (08-126, 1969.
13

---

A Proof of Lemma I
To simplify the notation throughout this proof, for each j E (I,..., k) denote Φ Pvis. We have
k
7e0.
xx (2b") x (1b) =
BY [()))(()) vp (Y)] T Flg Z
(14)
j=1
Observe that for any E X' it holds that
dup dpp.b-
Zep Z
شラ (x(d"(x)- (((()))(()>()))) (dT)(z))
(x)
exp
Hence,
Zw
Zza
T log
log Ex-up-
Zao
djw
Σラ=ı ()((((x) (6 (+))(()))))))))))()())))))))((
T log E Ex-pw-
ex)
exp
djppp
7
ΣS Wy (b)()(()) (
= T log Ex-up
exp
7
k Σ j=1
sup p
E
Ex-a
(0)
-KL(H, pw))
(15)
where in the final expression we have applied the Donsker-Varadhan variational principle (i.e., convex- conjugate duality between KL-divergence and cumulant generating functions); therein, the supremum runs over probability measures μ absolutely continuous with respect to pap. and is attained by j defined as
u(dx) exp
w (x) (10))))))))))
uw (dr)
1
exp (- 1 k
(-12 Σ w (x) Tret (dx)
w (x) - (0))(()))
α exp
h
Wy ( Tree (dx) = Tap (dx)) Σ
x exp
That is, the supremum in 15)) attained by μ pap Hence, the identity (4) becomes
Ex (ab*) [xx=x(( (xb)
h
Σ
wy My [(W*)(Y)) - (o Ex
Σ
wx) - -
2=1
+ TKL(
Σ 1=1
Wi (Eywes (()))()- vs (x)] Ex- ((0))(x)) - (x))) +KL( (bb)
2 +KL(( (w))s
where the final inequality follows by noting that for each j the optimality of the pair ((d*),, (0,"),,,""" (.+)))
for the entropic optimal transport dual objective E . implies that
Eyaud [(0)) () + Ex (xx))(xx-- oi (x))
EN (6"').("") Ex (0) (bt (0))))) 200
The proof of Lemma is complete.
14

---

B Proof of Proposition] ]
Recall that for any non-negative integer we have
Σi wydl(() (z)
(+(d))= = Zi² exp
Traf (dx).
where Z is the normalizing constant defined by
E_ (yd((()
Z:
exp
Tred (dx).
X
With the notation introduced above, we have
E(bb,) = Σ
Ny [v((X]] log Zr.
wy
==1
Hence,
E(b(+1) E(bb)) Σ
wi EYnu (V)) T log
Zi+1
Ze
de der (x) des
Zi+1
W/y EYaus
log
log
Zr
5=1
Zi+1
min (A, T)
Σ
wyKL(pd vz) - - log
Zr
3=1
Therefore, to prove Proposition it suffices show that the inequality
Ze+1
log
(16)
Z
holds for any 2 O. We will complete the proof Proposition Jsiin the following lemma, the proof of which is deferred to the end of this section.
Lemma 3. Let (zbi)e2o be any sequence of the form
Vh = - + m)loo((A))
where for E (I....,kk k), (A)) is an arbitrary sequence of strictly positive functions and min(1,T/N) Then, for any t 2 0 it holds that
k Z min(A, 7) log Σ E, mor [
T log
Z
To complete the proof of Proposition
we will apply the above lemma with A dul
Indeed, we
du?
have
>log
Z+t1
< min(A,, log
Σ
wy
dus duy
(8)
Z:
j=1
min(A, T) log Σ
wi E Ey [1]
j=1
=0. 0.
By 16). the proof of Proposition is complete.
15

---

B.1 Proof of Lemma 3
We will break down the proof with the help of the following lemma, the proof of which can be found in Section B.2] Lemma 4. For any sequence (2b)) and any t 2 0 it holds that
7/A
Zibi- log
À log Ei=ı Wr Exmpr
exp
is 2 X,
Zup
log Σラ wy WyEEmm
exp
(x)
ist< <>
where Φ. ΦΦ, and (eldx)= Zap. exp( Σラ with (x:) ( Foe (dx).
Observe that the sequence (4h)) of the form stated in Lemma 3 satisfies, for any for any j {I....,k) k3 and any 2 0,
exp ( o/+1 + oi )
= exp (-
A
dpp did
dai ddp
A/r
=oxp
log
=
where
did (x) = / v(dy) exp duy,
(1+1(1) + (x) c(z, (1)
X
/ v (2x)
vi(v)) + o (x) - c(a, v)
x
Hence, by Lemma
we have
Zep.11 Zoree
log
$ 5
1 min(), r)
Σ
Wy
Ex-eee
j=1
vi(l) + (X)) max(1,A/7)
(L (dg))
exp
(17)
A
We split the remaining proof into two cases: 2 X and 5 < A.
The 2 X. yields
When 2 A, have max(1,A/7) and min(1, 1. Thus, 17
Z log
ZΦd
wy Exaus U/ w (dd)) exp (i(()) + ol x c(X, ((
< min(A, +))r)log log
Σ
j=1
k
/ Ma(da)) / i (dd))i(()) exp
( + o( c(x,
min(A,7) T) log
Σ
W/
)
1=1
k
((()))((d)) / exp
vi(r) + of (X) c(x.3) pe(dx)
= minn(,,7 log
A
1=1
= min(A, log
Σg
Wi
/ Ai()))
di dus (y)
1=1
min(A, log j=1
Σ
W.y Ey Aor
Ai(u)
This completes the proof of Lemma when
A.
16

---

The case 7 < X. For j E fl..... ... kb
and any x E X' define the measure Pe by
po(dy) v (dy) exp d( + of (x) c(x, 4)
A
By the definition of oi. we have
J (d(d))
S u(dy) exp vy(()) c(2,))
(x(x)
exp
A
1x v(dy) exp
(b(()) c(x,,)
cxp log / v (dy') exp x
(b(((') -0(o, (0)2)))))
A
=1 In particular, Ps is probability measure. Hence,
can be rewritten as
Zor11
A A/T
log
< log
Σ wy WyExxps 5=1
Ey Ex
Ze
Because A/T > 1, the function ト 2// is convex. Applying Jensen's 's inequality to the conditional expectation and using the fact that 9//T 1, it follows that
log
7 < log Zep.
Σ j=1
wy
Ex~m
E
Ex~ol
$100 x
ue(dx) / ((( (V(()) + (x) c(x.)) u(dy). (18)
log
Σ
Wy
j=1
By the definition of vd #; we have
dv db (v)
ve(i) + (((x) c(x, pe(dx).
exp
dy
x
X
Interchanging the order of integration in 18
and plugging in the above equation yields
Zepret
Jx
vi(y) (1()) c(z.v) me(de) A (c(o)
log
< log
Wy
exp
Z
x
A
1=1
dedt dus
log
W/
Ur
(v)
Af(v)
j=1
x
log
Wy
Ey
j=1
This completes the proof of Lemma
B.2 Proof of Lemma 4
To simplify the notation, denote Z, Z Let X' and 0. We have Pe & p with the Radon-Nikodym derivative dpt+1/ /dd4 given by
dpc+1
Zz
- -EE wx (x))
(x)
exp
dj4.
7e+1
Zr1
-l41(0) + ol(x)
exp
Z
j=l
17

---

Multiplying both sides by Zz+1/Ze and taking expectations with respect to br yields
zt1
Z++1 7
Ex
7
Zi+1 dpet1 (x)
Exaea
Ze duz
oy (x) (X) w/
Ex=pi
II j=1
exp
T
In the case 7 < X, the proof is complete by the -eemetrri-ee mean inequality (recall that Wy > 0 for j = l....,k and Σ Wy 1). On the other hand, if 7 2 À then + x is concave.
Hence, it follows that
"yT/A A/T
- -x)+((()++ + 2x
Zo+1 log
log Ex-mps
πI exp
Zt
7
w/7/> A/T
(x+1(X) + +x(())
I log Exrpi
II exp
À log Ex~us
K II exp
Ho (X+ + (x)
MSF/AA
1=1
k
( dx) 7/2
Ex~p
exp
where the final step follows via the Ariohnetic-Geometric mean inequality. This completes the proof of Lemma
C Proof of Theorem 2
For every / 2 and j E (1...., k), let 77 be the distribution returned by the approximate Sinkhorn oracle that satisfies the properties listed in Definition T We follow along the lines of proof of Theorem [
First, we will establish an upper bound on the oscillation norm of the iterates 3. Indeed, by the property four Definition I we have
Ilose I - (1--))l//1e ase + 1/coo (2)..
Since T6 0, for any t 2 0 we have || /losc I coo (2').
Let 8, = (1b Ex,"" (xb:) be the suboptimality gap time t. Using the concavity upper bound 100 and the property two in Definition I we have
8, < 20
Σ
Wrll ITV
j=1
+ 20 (X)
Σ
Willer D/llTV
1=1
455 /20o0 (X)
Wi KL(W,D Dz)
j=l
Set V2 (X)
Σ
W,KL(us 7T).
18

---

Combining the property three stated in the Definition I with Lemma
3 we obtain
(LL(v,,)) wo Ko,KL((' oy)) min(A, r) log Σ wi S
du
- min r)
Σ
(w)) (dy)
dis
j=1
2
my KL m2)) min(), 7) (1 + =2// 12000 (2)
min(A, (A, T)
min(A, r)
Σ my (L,KL(v', ,
8²
20 (2¹)²
2 mi)(),7)
20 (2)² max - -
min( A, T) E²
20 (x)
Provided that 8, 2 28 it holds that
min(A, (X, (ōe -25
(6, - 2z) - (61+1 2z) 2
20,0 (X)
Let T' be the first index such that 8r+1 22 and set 3o if such index exists. Then, the above equation is valid for any s 'T. In particular, repeating the proof of Theorem | for any S T we have
20 (X)
8, - 2s S
min(A, T) 1
which completes the proof of this theorem.
D Proof of Lemma 2
The first property the positivity of the probability mass function of is immediate from its definition.
To simplify the notation, denote what follows
A (2)
+ vid (1) c(xx 3)
(², exp
With this notation, recall that
P2( = 1 n n
Σ rd (m)))(X,, v/).
The above a sum of n non-negative random variables bounded by one with expectation
(v')(())) Ex-ai [wr (ur))
follows by Hoeffding's inequality and the union bound that with probability at least the following holds for any j E (1,..., k) and any I E (I....,my}:
2 log 22
(w))))) (v') (2')
In particular, the above implies that
v rr $ 2Ç + (1 5) v villav
2 + (1 - 1/11%% - -(v'P + (1- - (111))1(//) - lrv
2 + 10 - -(v))' |I + -
2log 2mm
< 2( + myep + -my/
19

---

Notice that the above bound can be made arbitrarily close to m by taking a large enough n and small enough ç. This proves the second property of Definition [ To prove the third property, observe that
(Y) #
1 (X)
1 (Y) - T (2)
Eyaud
Eyaud
+
17
(X)
2/ (2)
1-5 + v/ (Y) 1 (x)
I Ey
1 (P)
1 (X) 1 (x)
5
< By
1
+
1 1-5
(x)
< 1 +
+
|1 [11((()) 526 (1)
1 1-
1
log (2m)
S 1 + 2
+
mgew + mmy
n
This concludes the proof of the third property. It remains prove the fourth property of Definition 1
Observe that for any y: of we have
(3)
2 (v') 22 (v¹)
(b)()) nllog v2())
(1/) 17> log
( (1 4) Σi= KJ (X,,3)
(v) - (1/) + no
(+ 5 Ei= Ki(X,,
- + Σ Kd (x(():1))
(v) v + Aog
{ + Eラ= Ki(X;,
1-c
I-3 + exp
Con (20) (1'))="(1) (v)
1 K3 (x(x11)
$ (v) - -0)) + Aog
- EE KJ (Xi,y)
Now observe that for any a, b > 0 the function g: : [0, 00) + (0, xc) defined by g(x) (a+a) (x+b) is increasing if a < and decreasing if a 2 b. Thus, is maximized either at zero or at infinity. It thus follows that
C I-C
+ exp (mald)
I EE (X(),))
m log
1-5 -r= Ki((,,y)
- n(b -
w (v/)) if exp (x)) 21
0
otherwise.
This proves the claim and completes the proof of this lemma.
E Proof of Theorem 3
The purpose of this section is show how sampling via Langevin Monte Carlo algorithm yields the first provable convergence guarantees for computing barycenters in the free-support setup (cf. the discussion the end of Section 2.2 In particular, we provide computational guarantees for implementing
Algorithm 2
A measure is said to satisfy the logarithmic Sobolev inequality (LSI) with constant C' if for all sufficiently smooth functions / holds that
1f log f E, log E,, (CE,
20

---

To sample from measure ju(dz) = exp(-f(z) )d supported on R the unadjusted Langevin Monte Carlo algorithm is defined via the following recursive update rule:
2++1
Vf(x)) + 2n where Z ~ (O,,
(19)
The following Theorem is due to Vempala and Wibisono [57 Theorem 3].
Theorem 4. Let =xd=)==xp(<f((a))
be a measure on Rd. Suppose that μ satisfies LSI a with constant C and that has L-Lipschitz gradient with respect the Euclidean sequence of iterates (Tx)k20 defined via (99) and let let pk be the distribution of zi. Then, for any O, any n ST2C min (1, 32 and any k z
Consider the
2C log 2KL(po..)) it holds that
KL SE.
Thus, LSI on the measure μ provides convergence guarantees on KL(p): p). shown in (7
Lemma I] how to initialize the iterate zo so that KL(Po, (4) scales linearly with the ambient dimension up to some additional terms. The final condition described Problem Setting eensrres that (by
57 Lemma IJ) for any > 0, the initialization scheme zo (a Ia) for the Langevin algorithm
satisfies
Coo (X)
d
Le 2
KL Pub.o)
log
T 2
where L, is the smoothness constant of V + dist(x, (20²) (see Lemma 5) and Pub.o is the probability measure defined in (20).
To implement the approximate Sinkhorn oracle described in Definition I we can combine Lemmal with approximate sampling via Langevin Monte Carlo; note that by Pinsker's inequality, Kullback- Leibler divergence guarantees provide total variation guarantees which are sufficient for the applica- tion of Lemma Therefore, providing provable convergence guarantees for Algorithm amounts proving that we can do arbitrarily accurate approximate sampling from distributions of the form
Hap [xp(dx) (x(x)) (xx(-VVo(t))/))))) where Vi((x)
Σ wydh (x).
Here 1x the indicator function of X, v/ is an arbitrary iterate generated by Algorithm and we consider the free-support setup characterized via the choice (dxr = 1xda.
Notice that we cannot apply Theorem 4 directly because the measure Pu defined above has con- strained support while Theorem 4 only applies for measures supported on all Rd Nevertheless, we will show that the compactly supported measure ps can be approximated by measure Hlb.a where the parameter will trade-off LSI constant of pybto against the total variation norm between the two measures. To this end, define
Mlb. =o (xxp(-Vx(())) - dist(x, /2)/(222)))) where dist x) inf ||x - y112- (20) yEλ
The following lemma, proved at the end of this section, collects the main properties of the measure p.p.r-
Lemma 5. Consider the setup described in Problem Setting Let sb be any iterate generated by Algorithm 2 and let pap be the distribution defined in (00 Then, the measure paplo satisfies the following properties:
1. For any (0, 1/4] holds that
(8R2)
[[
(ARd d-1 +1]
Mp - Moollrr S 20 exp
2. Let Va(x) exp(-V - dist(x, x) (2/(22))); thus p = exp(-Va(x))ddr The function V, has Lipschit gradient where
Lg =
max m, + -
TX
21

---

3. The measure He,o satisfies LSI with a constant C poly exp(R2/T) Lo).
Above, the notation C = poly(x, y/ 2) denotes constant that depends polynomially on a, y and 2. With the above lemma at hand, we are ready to prove Theorem
Proof of Theorem 3
Let & be an arbitrary iterate generated via Algorithm 2 We can simulate a step of approximate Sinkhorn oracle with accuracy via Lemma (with =/ in time poly(n, n, m, d) provided access to n = poly(e- m, log(m/6) samples from any distribution Pes such that
1l Mab M || S
(21)
16m
To find a choice of Usy satisfying the above bound, consider the distribution
8R2
d-1 -1 (AR + =1)²
Hablo
with
exx
32m
Let Co and Lg be the LSI and smoothness constants the distribution Mvye provided in Lemma 5 By Theorem [ it suffices
the Langevin algorithm (19) for poly(e m, d, Ca. Tog) number
of iterations to obtain a sample from a distribution Pv such that
1pp. or l/ S 32m
particular, by the triangle inequality for the total variation norm, the choice 4/ = Tip satisfies 21 This finishes the proof.
E.1 Proof of Lemma S
To simplify the notation, denote M4b. Por Pub.or V(x) Vv(()/T,, and Va(x) V(x)/T+
dist(x, x)) /2a)))/(2)))))))))))))))>))
Total variation norm bound.
With the above shorthand notation, have
p(dr)= xx eex((-V(x))ee where Z / exp(-V(x))da
and
Ho (dx) (Z -xp()))) where Zo /
exp(-Vo exp(-Vo(z))dd.
We have
/ R
--1 expppp (VV(x)) (Z + - exp(-Vo + / -
|| MollTV =
22.
2Zg < 2 exp
(oc(d))
4R²
<
Za < 2exp
Zo.
2o
T
We thus need to upper bound 2.. Let Vol( (A) be the Lebesgue measure of the set A, let 0A denote the boundary of A, and A + {a + b: E A, E B) be the Minkowski sum sets and B. Using the facts that for each j {l...., [.....k) k} we have sup z (и) (x) < 4R² and that
22

---

X' _ BR = fz || I R)
we have
7
exp(-Vo(x))da
4R²
dist(x, 2)
2 exp
exp
de
RU
2072
4R²
x²
exp
Vol(O(X'+ + B2)) exp
dx
T
20
4R²
x²2
exp
(ol(BB exp exp
dr
T
2o
exp
4R²
π
4/2
000
(R + x:)-11
exp
x2²
dx.
r(d/2)
2x
Bounding (R + xd)=1 < 2d pd-1pd-
+ 2d dd and computing the integrals results in
8R²
1/2
VT
|/ - Ma |TV S 2exp
2d
Rd-1
+
24/2
(dd/)o
T
T(d/2)
2
88/2
(2R)a-1
< 2a exp
+ (40)
T(d/2)
Using the assumption above bound to
5 1/4 and using the bound ((d) 2 (d/2)d/2 can further simplify the
81222
[ d-1 +1|
|| Me l|v < 2o exp
(
which completes the proof of the total variation bound.
Lipschitz constant of the gradient. Recall that for any any j E (1,...,,d] .. we have
(wi)
II 12 (xr, wi'
e (x)) - ||=11 A -Alog
Σ exp
v w((rt)
Denote 2 (x) Φ (x))) - 111=1112- Fix any a, x' and define g(() / (o + (x' - x))(). Then, for any
t E 1) we have
1 Varzape [ord xx))e]=
1
g"(s)
- b||z- x' -1'1m2m,4
(22)
where
1bt (b/
May 112
+ (x + 1(x' - x))
(())) & v(yt) exp
and y E (ddxmy the matrix whose /-th column equal to the vector n/-
Because Tاه is concave, the bound 22 shows that di is 1 + $1m.4tt2ssmoo
Combining the above with the fact that the convex function dist X') has 1-Lipschitz gradient 6 Proposition 12.30] proves the desired smoothness bound on the function V.,
LSI Constant bound. The result follows, for example, by applying the sufficient Tog-Sobolev inequality criterion stated in [13 Corollary 2.1, Equation (2.3)), combined with the bound (22p. The exact constant appearing in the log-Sobolev inequality be traced from [13 Equation (3.10)).
F Numerical Experiments
In this section, we numerically validate our main theoretical results presented in Theorems and We empirically demonstrate the necessity of damping Sinkhorn iterations when < X/2. In addition.
23

---

we examine empirical convergence rates of Algorithms 1 and 2 (for a specific implementation of the approximate Sinkhorn oracle described below) simulation setup comprised of isotropic Gaussian measures
Simulation setup. Let px denote the optimal (a.aaaaaaaeeeeeeeeeee -baryceeterr; that is, M²,7 is the unique solution
to the optimization problem (2). Let (24))22 denote the iterates of either Algorithm or Algorithm Then, the dual objective sub-optimality gap bounds of Theorems I Jand combined with Lemma establish the convergence (L (4)) to zero as goes to infinity.
To numerically compute + me) we need know the true (A, r)-barycenter Mx,7 The only setup where 1 admits a known closed-form expression is when the marginal measures of are isotropic Gaussians with identical variance. In particular, it was shown in T7 Proposition 3.4] that for vd N(m,, 0² Ia) and for any non non-negative weights ( that sum to one, we have
k Σ Wim 111 Sle): where e² (0 (0
+ X²
1 N
4
Hence, in all the simulations performed in this section, we let vو = Nmm,, o²la) for some my E Rd and o² > O. While this simulation setup allows for exact computations of the divergence measure (T,()),r::14) /4), there are primary limitations in our experimental design. First, the theoretical results in this paper are proved under boundedness assumptions on the marginal measures, while Gaussian distributions are unbounded. Second, our simulations do not cover the free-support discrete point clouds setup investigated in Section 4.1 concerning the approximate Sinkhorn oracle implementable via Monte Carlo sampling.
Implementation of Algorithm isotropic) and when bi O, then Sinkhorn updates admit closed-form expression that result in quadratic Sinkhorn potentials v) o and Gaussian measure (4, In particular, the iterates of Algorithm can be written
1
When the marginals y are Gaussian measures (not necessarily
vid (v) = 2 A.y y 1/. + const, d(v)=
2 - -x
x TT const, 1t = N Ex) (23)
for some matrices A,, Ct Eg Raxd and vectors by, de, e, E Rd that can be computed using explicit recursive expressions (see, e.g., 34
We implement approximate Sinkhorn oracle (Definition) D as follows. Our updates maintain the property that the potentials v and ø² are quadratic functions the form described in Equation (23) At every iteration f, we replace each matrix A/ (for j = 1,..., .)) by performing the transformation A + A - -N/, where N is random positive-definite matrix with trace equal to In our simulations, each NO is obtained by drawing independent d x d matrix U with i.i.d. N (0, 1) entries and setting N/C approximation error the implemented approximate Sinkhorn oracle.
Implementation of Algorithm
2
U T /ttrccaeeU T U). The parameter controls the
Simulation results.
We now comment on the main findings in our numerical simulations.
In Figure I we consider the toy setup 1, v/ (0, 1). Figure la shows that undamped iterates explode when < >/2 (note the absence of the blue line due to the explosion of iterates when A/r 2.1). Figure lb investigates the critical case 7 As A/2, noticing sharp phase transition in the convergence behavior. (elureeicdemmostrraal that damping removes the bad behavior. Moreover, the damping factor suggested our work yields the fastest convergence among the tested five different damping factor choices.
we investigate undamped inexact algorithm described above. As suggested by Theorem the iterates converge up to some level governed by the accuracy parameter a of the inexact Sinkhorn oracle.
In Figure 2
'The code for reproducing the simulation results is available at https://githut com/
Tomas Conass/aalkkriiccuuu/toolllllaao http:://wwbbbyy.tttoooooeeoerrre baryyeeoppi--bbryccetetrrr
24

---

In Figure we perform the same simulations as in Figure 2 but this time with damping. We observe again that damping helps significantly in the critical case A/2, yielding much faster convergence rate. The second row in Figure 3 considers the overdamped case, where we observe slightly slower convergence than that with the optimal choice of the damping parameter.
In Figure 4 investigate the effect of different levels of damping in the noisy setup (Algorithm 2 First, observe that the undamped algorithm (1) explodes in the noisy case with T X/2. in contrast to the zero-noise setup reported in Figure Ta Moreover, we find that overdamping might have negative effect on attained accuracy: larger damping factors lead lower accuracy convergence in noisy setups.
Alportiaa I viliian daspin by 1)
Alpor withont danpins, (o
Plussi trundion aeurao Af+
Alaycillan : will desping (A 10 10-,A 200
Ind Klas Blon
10g NLL0
Number Rencicee
Number naabous
Number neatoos?
(a) Explosion for X/T > 2.
(b) Phase transition at A/T 2.
(c) Damping helps.
Figure 1: Simulations for Algorithm
with A: = 1, v² N(0, 1).
Algoniim 2 wothort damping (A 10-3)
-
lap Sto5, (ln
Nunder Zomatioss
Namtier
Nunder
1-0.00
Figure 2: Setup: l = w 1/3.1133111 = vl = = N(O, h10) The effect of noise for
the inexact algorithm.
Alpettian: 2 with damptig (A 1imlly
Aft Dule
Air 11S
Jug 00.16, [le))
t- 0.30
0.001
00001
Figure 3: Setup: A = 3, w = (1/3, 1/3), o2 = 1² - N(O, 1o). The effect of noise for the inexact algorithm. Second row overdamps.
25

---

Appeltin 2 with damping (A 10 Urall
Namber lbneiie
Namber of Ieeeeiiee
Rumber of menatiens
9-
0-48
9--385
Figure 4: Setup: 3 (1/3,1/3.1/3) v = = 2 N(0, The effect of damping for the inexact algorithm.
26


Reviews:
[{"review_id": "xgV7hjJYFH", "reviewer": "Reviewer_nAci", "summary": "This paper presents a study on the computation of doubly regularized Wasserstein barycenters, a recently introduced family of entropic barycenters with inner and outer regularization strengths. The authors build upon previous research, which has shown that different choices of regularization parameters unify various entropy-penalized barycenter concepts, including debiased barycenters.\n\nThe proposed algorithm for computing doubly regularized Wasserstein barycenters combines damped Sinkhorn iterations with exact maximization/minimization steps, ensuring convergence for any choice of regularization parameters. Additionally, the authors introduce an inexact variant of the algorithm that utilizes Sinkhorn Oracle definition, providing non-asymptotic convergence guarantees for approximating Wasserstein barycenters between discrete point clouds in the free-support/grid-free setting. The authors highlight that while a straightforward adaptation of the alternate maximization scheme leads to diverging iterates for small values of $\\tau$, their analysis demonstrates that damping these iterations is sufficient to achieve convergence.\n", "questions": "I have the following questions/comments for the authors:\n\n1. Can the authors provide more motivation for using the definition of Sinkhorn Oracle in Section 4, Approximate Damped Sinkhorn Scheme? It would be helpful to understand the rationale behind this choice and how it contributes to the proposed algorithm.\n2. In Section 2.2, Doubly Regularized Entropic Barycenters, could the optimization part (Line 172 - 184) be made more clear and detailed? This section seems to be the foundation for the updates in your main contribution algorithm. Providing additional explanations and elaboration would enhance the understanding of the optimization process. \n3. Algorithm 1 is presented without a thorough mathematical derivation. It would be beneficial to include a more detailed derivation or explanation of the algorithm to aid readers in understanding the underlying principles and steps involved.\n", "limitations": "While the paper provides theoretical guarantees and analyses, incorporating visualizations and benchmarking against alternative approaches would enhance the practical significance of the proposed algorithm\n", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Clear structure, good mathematical exposition, and solid theoretical results. \n2. Detailed proofs to support their claims, demonstrating a strong theoretical foundation\n", "weaknesses": "1. Lack of numerical experiments and empirical evaluations. \n2. While the paper discusses computational complexity, further analysis and discussion on the scalability and efficiency of the algorithm could provide a more comprehensive understanding of its practical implications.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments and questions. We respond to your questions below:\n\n1. Let us explain why the Approximate Sinkhorn Oracle (Definition 1) is defined the way it is. First, why do we need an approximate algorithm at all? We need it because, for continuous measures, we typically cannot implement the exact algorithm (Algorithm 1) due to the high-dimensional integrals that need to be computed in lines (a) and (e) of Algorithm 1. Of course, when the marginals are discrete, for any given x we can actually evaluate the line (a), as the integral inside is just a sum. Thus, the only issue is with the computation of line (e), repeated below for convenience:\n\n$$\n(e) \\quad \\frac{d\\nu\\_{t}^{j}}{d\\nu^{j}}(y) \\leftarrow \\int \\exp\\left(\\frac{\\phi_{t}^{j}(x) + \\psi_{t}^{j}(y) - c(x,y)}{\\lambda}\\right)\\mu\\_{t}(dx).\n$$\n\nWhile we cannot compute the integral above exactly, if we can generate samples from $\\mu_{t}$, then for any given $y$ we could estimate the desired quantity $\\frac{d\\nu\\_{t}^{j}}{d\\nu^{j}}(y)$ up to some error.\n\nThis leads to the question: suppose, instead of $\\frac{d\\nu_{t}^{j}}{d\\nu^{j}}(y)$, we have some other function $g^{j}\\_{t}(y)$. How close does it need to be to $\\frac{d\\nu_{t}^{j}}{d\\nu^{j}}(y)$ for us to run Algorithm 1 with $g^{j}\\_{t}$ instead of $\\frac{d\\nu_{t}^{j}}{d\\nu^{j}}(y)$? Well, the answer to this is laid out in Definition 1. To get these properties, you just need to follow the proof of Theorem 1 and put $g^{j}_{t}$ in place of  $\\frac{d\\nu\\_{t}^{j}}{d\\nu^{j}}(y)$, observing what properties need to hold for you to preserve the convergence analysis up to some tolerance. You may then see how the proof of Theorem 2 follows the lines of Theorem 1 when the properties of Definition 1 are plugged into the proof.\n\nWe hope this helps!\n\n2. and 3. Thank you for the suggestions. We will make sure to clarify these parts when we update the paper."}, {"review_id": "mCi1SFFUs0", "reviewer": "Reviewer_uehn", "summary": "This paper proposes an algorithm for solving the doubly regularized Wasserstein barycenter problem for probability measures that corresponds to adding an inner regularization based on the entropy penalty appearing in the Wasserstein distance term, and an outer regularization appearing at the level of the Wasserstein barycenter problem, based on the KL divergence. The algorithm comes with convergence guarantees, depending on inner and outer regularization parameters.The algorithm is then modified to approximately solve the doubly regularized Wasserstein barycenter problem, but with non-asymptotic convergence guarantees.", "questions": "Probability measures are supported on a compact convex subset of $\\mathbb{R}^d$, could this be relaxed?\n\nIn practice, how does the number of support points of the discrete barycenter is chosen? More generally, how do you choose $\\pi_{ref}$?", "limitations": "An experimental (and if possible theoretical) comparison of the barycenters obtained by the exact and approximate proposed algorithms, as well as a comparison with Bregman iterative projections algorithm for particular choices of regularization parameters would, in my opinion, improve this article.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "This article is very well written, proposing an algorithm for a recently introduced notion of regularization for the Wasserstein barycenter, unifying several proposals in the literature. The proposed algorithm does not rely on space discretization (prohibitive in high dimensions). The proofs are accurate, well documented and complete. The most important contribution lies in the approximate damped Sinkhorn scheme, which allows for non-asymptotic convergence guarantees on the value of the dual objective.", "weaknesses": "The idea of dumping an algorithm is not really innovative, but it stabilizes the algorithm and provides convergence guarantees for any positive values of the parameters $(\\lambda,\\tau)$.\n\nA result on the distance between barycenters obtained via the exact and approximate damped Sinkhorn scheme could have been interesting. In particular, the addition of experiments on the barycenters obtained via both algorithms would strengthen the paper.\n\nFinally, an intuition on the approximate Sinkhorn oracle (Definition 1) is important in my opinion, as the properties listed in Definition 1 are quite precise. An example of a Radon-Nikodym derivative that does not satisfy these properties could be added.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments and questions.\n\n- Regarding the compactness assumption, it is really only needed in our case because the prior work that introduced doubly entropic barycenters derived many theoretical results in the compact case. In particular, compactness was used in that context to justify a certain interchange of inf-sup [Section 2.3, https://arxiv.org/pdf/2303.11844.pdf]. Because these types of results can usually be obtained via other means (e.g., by approximation), and because it is known how to analyze Sinkhorn's algorithm for unbounded measures, we do not see any major obstacles in extending our results beyond the compact case. \n\n- Regarding the choice of $\\pi\\_{\\mathrm{ref}}$, there is, in general, no way to know where the optimal barycenter is supported or in which part of the region it is concentrated without computing it first. Thus, you would typically choose $\\pi\\_{\\mathrm{ref}}$ by gridding the space. Of course, in high dimensions, this approach is not feasible. This motivates looking into grid-free/support-free methods such as Algorithm 2 with the implementation suggested following Lemma 2."}, {"review_id": "mvhDgKSxxR", "reviewer": "Reviewer_9mYE", "summary": "The paper has proposed a computational algorithm for computing the newly developed regularized Wasserstein barycenters in [Chizat,2023] via optimizing the duality of the primal problem. The paper also characterised the convergence of algorithms on both exact and approximated algorithms which are the main contributions of the paper.", "questions": " As the doubly entropic Wasserstein barycenter may have some interesting properties for some specific type of data if the author can discuss further in terms of theoretical and experimental aspects that would be great.", "limitations": "No empirical result provided.", "rating": 4, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "The algorithms and their convergence guarantee for computing doubly entropic Wasserstein barycenters are significant contributions in terms of technical aspects. The paper is easy to follow as the related results have been well presented to make use of the derivations of the proposed algorithms. ", "weaknesses": "There are some aspects the paper can be improved:\n- A section for numerical experiments/illustrations for both synthetic and real-world data will help to demonstrate the results of the two main algorithms and their convergence algorithms.\n- Algorithms should link to the text updating equations. For instance, in Alg. 1, step 2(a) computes the left part of the equation (11) while step 2(b-d) computes the right part of that equation, authors break down into multiple lines of code without clear purposes which usually make readers harder to follow.\n- In line 199, \"We have observed empirically that the iterates of the iterative Bregman projections (i.e., the scheme of updates (12), (11)) diverge whenever \\tau < \\lamba<=2\", it is not clear where the results with the iterative Bregman projections presented.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the review.\n\nThe theoretical aspects of doubly entropic barycenters have been thoroughly investigated in https://arxiv.org/pdf/2303.11844.pdf. Our primary goal was, instead, to provide new numerical schemes for their computation and to establish their convergence, particularly covering the case $\\tau < \\lambda$ for which we obtained the first provably convergent method.\n\nIn addition, some numerical examples concerning a comparison of $(\\lambda, \\tau)$-barycenters with different parameter values are available in Section 6 of the referenced paper.\n"}, {"review_id": "GchrlT8Uvf", "reviewer": "Reviewer_ofnZ", "summary": "-----\nEDIT : 6 --> 8\n-----\n\n\nThis paper builds on [15] and considers the recently introduced model of _doubly regularized entropic Wasserstein Barycenters_ which, given a set of measures $\\nu^1,\\dots, \\nu^k$, weights $(w_j)_j$ (non-negative and sum to $1$), a reference measure $\\pi$, and two smoothing parameters $\\lambda,\\tau > 0$, consider the problem\n\n$$ \\text{ minimize } \\mu \\mapsto \\tau \\mathrm{KL}(\\mu,\\pi) + \\sum_{j=1}^k w_j T_\\lambda(\\mu, \\nu_j), $$\n\nwhere\n\n$$ T_\\lambda(\\mu,\\nu) = \\inf_\\gamma \\braket{c,\\gamma} + \\lambda \\mathrm{KL}(\\gamma,\\mu\\otimes\\nu). $$\n\nDepending on the choice of the parameters $\\lambda,\\tau$, this problem encompasses many standard formulations for Wasserstein barycenters and their variations (Schrödinger barycenters, etc.). \n\nThis problem is deeply analyzed in [15] from a theoretical perspective (existence, uniqueness, smoothness of solutions, etc.), which concludes on the need to elaborate on the computational complexity of these $(\\lambda,\\tau)$-barycenters. This is the purpose of the current paper, which has two (three) main contributions to that respect : \n1. They provide a theoretical _damped Sinkhorn algorithm_ (which can been seen as an adaptation of _Iterative Bregman Projections_ of [6]) for which they prove a $O(t^{-1})$ convergence rate toward the global minimizer of the above functional ($t$ being the number of steps). In a nutshell, the trick consists of replacing the \"learning rate\" $\\lambda$ in the standard Sinkhorn algorithm by the quantity $\\min(\\lambda,\\tau)$. \n2. They provide a  _approximated_ algorithm for which an implementation is possible, provided one has access to an \"ApproximatedSinkhornOracle\", and for which convergence rate are accessible,\n3. (2bis) They provide an example of such a practical Approximated Sinkhorn Oracle based on random sampling techniques. \n\n\n### References\n[15] : (Doubly regularized Entropic Wasserstein Barycenters, Chizat, 2023) ", "questions": "# Main questions\n\nThese are questions/suggestion I would like to see addressed during the discussion period. Given what I have written above, they naturally rely on (synthetic) numerical experiments that may hopefully shed some light on the strengths and weaknesses of the proposed approach.\n\n1. Is there a computational price to pay for the _damped_ Sinkhorn? Namely, from what is written in the paper, I assume that the \"naive\" Sinkhorn algorithm (with $\\eta = 1$ no matter $\\tau$) does converge when $\\tau \\in [\\lambda/2 , \\lambda]$ (at least empirically). If so, (i) does it provide (empirically) the same output (asymptotically in $t$) as the output of the damped algorithm? (ii) if yes, is it faster? (I would expect that it is, since the damped algorithm somewhat slow the gradient descent). \n2. Can you elaborate on $\\pi_{\\mathrm{ref}}$ be the Lebesgue measure encode \"the free support case\"? If my understanding is correct, this forces $\\mu$ to have a density, which is not exactly what \"free-support\" refers to in the work of Cuturi and Doucet (2014) (but I agree that their terminology can be discussed as well as they fix the cardinality of the support). But typically, if the $\\nu^j$ are discrete, their $(0,0)$-barycenter is discrete as well (Carlier & Ekeland, 2015) and in that case, it may be hard for $\\mu_t$ to concentrate. From [15, Thm 3.2] we may expect a $\\lambda^2$ proximity between the exact Wasserstein barycenter and the $(\\lambda,\\lambda/2)$ one, but as far as I can tell this theorem is proved only in the case where the measures have densities. So I would be interested to see if/how much things fail in the discrete setting (now that it is possible to implement it!). \n\n\n# Minor questions\n\nThese are minor questions that I'm asking because I am interested in the work, but I would not be offended if the authors do not address them during the discussion period. \n\n- From my reading of the proofs, I missed the point where using the damped algorithm is crucial. I think that I have a correct line-by-line understanding of the proof, but not a global vision; could you tell me what would fail if one tries to run the same proof but with $\\eta = 1$? I guess it appears in the proof of Lemma 4 (case $\\tau < \\lambda$), which I did not investigated; is it possible to have a global picture of why it is required?\n- In the way Theorem 2 is currently stated, nothing is said once the condition on $T$ is reached (assuming $T < \\infty$). Is it clear that once the condition is reached, the gap remains below $2\\epsilon$? As far as I can tell by checking the proof, it is not sure that the sequence $\\tilde{\\delta}(t) - \\tilde{\\delta}(t+1)$ is non-increasing (once the criterion is reached), may it happen that the gap increases again in an uncontrolled way? (i.e. at time $T$ we are below $2\\epsilon$, but at $T+1$ we suddenly get something much worse and we have to keep running the gradient descent again) I think that (if such thing can happen) this can be avoided by simply checking the variations of the objective value. ", "limitations": "I do not identify specific limitation of potential negative societal impact specific to this work. \n\nFor the former point, I guess that additional numerical experiments may showcase some limitations of the work, that would be worth discussing then (without diminishing the contributions of the work). ", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "First, the paper is overall very well written. It manages to remain concise without sacrificing the mathematical rigor, and basically every line/equation is interesting, which is highly appreciated. \n\nI also believe that its theoretical contributions are fairly significant and quite insightful and clearly speak in favor of adding more _outer regularization_ in OT-related problems. \n\nThe proofs have been investigated (with the exception of those of Lemma 2 and 3) and no major flaw was detected; the clarity of the presentation make them easy to parse while not being trivial. ", "weaknesses": "# Major issue \n\n**The paper feels somewhat incomplete.** This is my main (and almost single) issue with the work (which undoubtedly has great potential). In a nutshell, this work is dedicated to the elaboration of provably converging **practical** algorithms to compute $(\\lambda,\\tau)$-barycenters; hence having absolutely no numerical illustration of such algorithms is somewhat surprising---even more given that the implementation seems to exist ($\\ell$199, \"we have observed empirically\"; why not showcasing it?). \n\nTypically, the divergence of the standard Sinkhorn algorithm (when $\\tau < \\lambda/2$) and the convergence of the damped one would have been nice to observe (especially given that the former is only empirical). Similarly, the paper mentions both in the abstract and in the conclusion that the approach works for _both fixed and free-support_ models (which boils down to the choice of the reference measure $\\pi$), but this is (almost) not developed in the work. In my opinion, it is fascinating to have a globally converging way to approximate Wasserstein barycenter, and showcasing empirical situations where the proposed approach provides much better results than the \"naive\" one (say, IBP). \n\nSimilarly, the work mentions at several occasions the edge case $\\tau = \\lambda/2$ (because of the powerful statistical properties showcased in [15]), but never leverage this setting as far as I can tell. I do not know if there is something specific to say from the algorithmic perspective, but here as well, I feel like this could have been a nice occasion to numerically showcase the debiasing effect (for instance). \n\nDo not get me wrong: I believe that the current contributions of the paper are of interest. But in my opinion, there are few elements that are missing to make the paper complete. \n\n# Minor issues and other remark/suggestions\n\nNote: this are not actual criticism that are expected to be addressed during the discussion period immediatly, but rather comments that may hopefully be useful for the authors. \n\n- The hyperlinks are broken (probably due to the split of the pdf between the main body and the supplementary material). It may be nice to fix that for the camera ready version (which implies, I guess, to compile the supplementary material independently. Note that you can include the main body in the supplementary, which makes it \"self contained\" and quite handy for reviewers/readers). \n- [typo] In Eq. (2), I think that the $(w_j)_j$ are missing. This occurs in other places as well (Eq 8, def of $\\mu_\\psi, Z_\\phi$). \n- [typo] The parameters $\\nu,w, \\lambda,\\tau$ for the objective function $E$ are sometime missing (e.g. in Eq 10). This can be slightly confusing given the proximity with the expectation operator $\\mathbf{E}$. \n- [typo] $\\ell$193, $\\nu_j^t$ should be $\\nu_t^j$ I guess. \n-  [typo] In definition 1, the notation $\\hat{\\nu}$ and $\\tilde{\\nu}$ are both used (in my understanding, they denote the same quantity). \n- [typo] $\\ell$232, I think it should be \"the proof _of its convergence_ is deferred...\". \n- [suggestion] That's trivial but it may be handy to write in $\\ell$215 that  $\\delta_{t}-\\delta_{t+1} = E_{\\lambda,\\tau}^{\\nu,w} (\\phi_{t+1}) - E_{\\lambda,\\tau}^{\\nu,w} (\\phi_{t}) \\geq 0$ by Prop 1, making the sequence of error non-increasing (instead of saying it in $\\ell$222). \n- [suggestion] Lemma 1 is never referenced after being stated (only once in the contributions section). While it is of course obvious, it may be worth adding a short sentence after Theorem 1 saying something like \"Theorem 1 together with Lemma 1 implies the convergence of $\\mu_t$ toward $\\mu_{\\lambda,\\tau}$.\n- [typo] $\\ell$430, I think that $\\Delta_t^j$ should be defined without the $\\log$. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for such a thorough review and for your very helpful suggestions.\n\n## Answers to Main Questions\n\n### Question 1: Computational Price of Damping\n\n*Is there a computational price to pay for the damped Sinkhorn?*\n\n**Answer:** Actually, in the numerical experiments attached to the main response, we have observed that the correct amount of damping helps to speed up the convergence speed. Indeed, Figure 1 (c) shows that $\\eta = 0.5$ gives the fastest convergence (which, in this setup, is the choice suggested in the paper). Also, you may inspect Figures 2 and 3. We have three rows: no damping, theoretically suggested damping and overdamping. In all three columns, we see that the middle row (theoretical choice) gives the fastest convergence.\n\nHowever, there is a second side to this. In particular, Figure 4 suggests that damping more is worse in noisy setups. In addition, by overdamping beyond the theoretical choice specified in Algorithm 1, we do see a slowdown as expected (see Figure 3, compare the top row with the bottom one).\n\n*Namely, from what is written in the paper, I assume that the \"naive\" Sinkhorn algorithm does converge when $\\tau \\in [\\lambda/2, \\lambda]$ (at least empirically).*\n\n**Answer:** This is true, which can now be seen in Figure 1 (c). Indeed, it also answers your question (i) that we asymptotically converge to the same value regardless of the choice of $\\eta \\in (0,1]$, and we have already discussed (ii) above. \n\n### Question 2: On the Reference Measure\n\nYou are absolutely right that whenever $\\pi\\_{\\mathrm{ref}}$ has a density, it forces the optimal $(\\lambda, \\tau)$-barycenter also to have a density. And for arbitrary discrete marginals, with $\\pi\\_{\\mathrm{ref}}$ being the Lebesgue measure, the convergence via our Algorithm 2 would indeed be slow as $\\tau$ goes to $0$ (our convergence guarantee is exponential in $R/\\tau$ where $R$ is the radius of the domain). However, this is unavoidable due to the NP-hardness results for approximating $(0,0)$ barycenters due to https://arxiv.org/abs/2101.01100.  \n\nIt is true that the $\\lambda^2$ proximity to the true $(0,0)$-barycenter only holds in the continuous and smooth case. However, we can benefit from the debiaising effect of the $(\\lambda,\\lambda/2)$ barycenters -- even when we are doing computations with discrete measures -- when there is an underlying continuous structure. Consider for instance the common setting where the discrete marginals are obtained by discretizing (say, via $n$ iid samples) some continuous distributions and we would like to compute/estimate the unregularized barycenter of the continuous distributions directly. Combining the approximation ($\\lambda^2$) and estimation ($\\lambda^{-1-d/2}n^{-1/2}$) errors leads to better estimation bounds using $(\\lambda,\\lambda/2)$-barycenters compared to other entropy regularized barycenters (see [15]); even though the computation is done with discrete marginals.\n\n\n## Answers to Minor Questions\n\n- Here, I hope the performed numerical simulations will be helpful. Consider the $\\tau = \\lambda/2$ cases shown in Figure 1 (a) and Figure 2 (left-most plot). Both plots display a zig-zag-type line for the undamped algorithm ($\\eta = 1$), which essentially shows a behavior analogous to the one encountered in quadratic optimization for gradient descent with too large step size, which is when gradient descent iterates bounce back and forth between different sides of a valley. In the gradient descent case, this zig-zag behavior arises for smooth problems with a too-large step size, which is also what is happening in our simulations here. Indeed, there is a precise sense in which smoothness of the dual maximization objective for doubly entropic barycenters can be characterized [Proposition 2.6., https://arxiv.org/pdf/2303.11844.pdf], which shows that decreasing $\\tau$ (unsurprisingly) makes the dual objective less smooth, which is detrimental to the $\\eta=1$ algorithm as its convergence is essentially smoothness based. What matters is the $\\tau$ relation to $\\lambda$, and as $\\tau$ passes below $\\lambda/2$, we move beyond a critical threshold at which the exact alternating maximization/minimization stops working (see Figure 1 (b) for a display of this phase transition).\n\n- You are absolutely right. We do not say anything about what happens after the condition is reached. Indeed, numerical simulations (right-most column in Figures 2,3,4) demonstrate a curious behavior, where Algorithm 2 first reaches a good value and later gets stuck into a slightly worse one, where it remains indefinitely. "}, {"review_id": "36TSVhQQWg", "reviewer": "Reviewer_bjAF", "summary": "The paper presents an algorithm (damped Sinkhorn) and theoretical convergence guarantees for computing doubly regularized Wasserstein barycenters. The concept of doubly entropic Wasserstein barycenters extends the single entropic regularized barycenters by introducing an additional level of regularization. This addition allows for the two regularization terms to counterbalance each other, leading to a debiasing effect under the right conditions. The authors demonstrate that the damped Sinkhorn algorithm can be implemented using an approximate Sinkhorn oracle. The latter can be evaluated using Monte Carlo sampling, providing a complete computational pipeline with theoretical guarantees.", "questions": "The LSI constant derived using Holley--Stroock perturbation seems impractical due to its large value. Are there any more practical situations where this constant could be more manageable?\n\nWhat can be said about the convergence of the algorithm for an unbounded cost, e.g., the quadratic cost over $\\mathbb R^d$?\n", "limitations": "Yes.", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper tackles a complex and highly relevant problem in optimal transport and computational geometry: computing Wasserstein barycenters.\n\nThe convergence guarantees are rigorous and well-explained, enhancing the value of the proposed algorithm.\n\nThe paper connects the algorithmic problem to practical computation, demonstrating how the proposed methods can be implemented using Monte Carlo sampling. This gives a complete end-to-end view of the problem, which is valuable.\n", "weaknesses": "The paper has some minor errors and inconsistencies, such as missing weights in Equation (2) and inconsistent notations (lack of boldface in $\\psi^*$) in Lemma 1 and Theorem 1. There is a typo in line 203 (“is can”). In Definition 1, both $\\widehat\\nu$ and $\\widetilde\\nu$ are used.\n\nIn the discussion above Lemma 2, the claim about achieving accuracy close to $\\sqrt{\\varepsilon_\\mu}$ is misleading as the bound seems to be larger than $(m_j \\varepsilon_\\mu)^{1/4}$ based on the first two terms in the bound.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and for spotting some typos.\n\n- The LSI constant is indeed large, which renders the computational speed exponential in $R/\\tau$, where $R$ is the radius of the domain. This is, however, unavoidable because Wasserstein barycenters are NP-hard to compute for discrete point clouds (https://arxiv.org/abs/2101.01100). This means we must pay an exponential cost if we let $\\lambda$ and $\\tau$ go to zero.\n\nConcerning interesting practical situations where the cost could be manageable, we may consider the following:\n\n1. Use $\\lambda,\\tau = \\Theta(1)$; this way, we lose approximation properties but gain in statistical properties and in computational speed;\n\n2. Work with specific classes of marginals, e.g., it is possible to carry out exact computations for Gaussian marginals, where we can actually implement the exact scheme (Algorithm 1);\n\n3. There is some hope that some efficiency could be gained, for example, for log-concave measures, again relying on sampling approaches. However, we are not aware of any specific results in this direction yet.\n\n\n- Regarding convergence in the unbounded case, we believe this should be possible. It is known how to analyze Sinkhorn's algorithm in the unbounded case (see https://arxiv.org/abs/2212.06000), and because our scheme essentially combines Sinkhorn updates with damping, we believe it should not be too difficult to extend our results to quadratic cost in $\\mathbb{R}^{d}$."}]

Author rebuttal:
We thank all the reviewers for their feedback. We will answer minor questions raised by the reviewers individually; in this shared response to all reviewers, we will focus on the numerical simulations aspect.

Most of the reviewers pointed out the absence of numerical simulations as a significant shortcoming of our paper. We have performed numerical simulations that cover a significant fraction of the concerns raised. 

## Note: there is currently no performance benchmark for our algorithm.

Several reviewers suggested benchmarking our algorithm against others. We would like to clarify that:
- We have provided the first provably convergent method for the case $\tau < \lambda$ (both exact and inexact versions).
- Our focus on $\tau < \lambda$ is deliberate. Specifically, the recently discovered case $(\lambda, \lambda/2)$ enjoys very strong theoretical properties, and we have provided the first algorithm covering this case with convergence rates.
- When $\tau \geq \lambda$ (no damping), our algorithm reduces to the classical alternate maximization/minimization scheme, which is well understood from both practical and theoretical perspectives. However, note that we provided a unified analysis that covers all choices of $(\lambda, \tau)$.

## Simulation setup: isotropic Gaussian marginals.

However, even without a benchmark, we can perform insightful simulations to address several issues highlighted by the reviewers. These include:
- demonstrating the need for damping when $\tau < \lambda/2$;
- examining empirical convergence rates with and without damping;
- testing the inexact damped Sinkhorn scheme (Algorithm 2).

To provide meaningful computations, we need a setup for which we know what the true $(\lambda, \tau)$ barycenter is. In this response and in the attached simulations, we denote the optimal barycenter for a given problem by $\mu^{*}\_{\lambda, \tau}$, while the iterates of our algorithms (Algorithm 1 and Algorithm 2) are denoted by $\mu\_{t}$.

The only known case with closed-form $\mu^{*}\_{\lambda, \tau}$ is that of isotropic Gaussian marginals with the same variance (Proposition 3.4 in https://arxiv.org/pdf/2303.11844.pdf). Let $w \in \mathbb{R}^{k}$ denote the non-negative weights vector that sums to one, let $\sigma^{2} > 0$ be arbitrary, and for $j=1,\dots,k$ let $m\_{j} \in \mathbb{R}^{d}$ be arbitrary. If $\nu^{j} = N(m\_{j}, \sigma^{2}I\_{d})$ for $j=1,\dots,k$, then 

$$
  \mu^{*}\_{\lambda, \tau} = N(\sum\_{j=1}^{k} w\_{j}m\_{j}, \xi^{2} I\_{d})
  \text{ with } \xi^{2} = \frac{(\sigma^{2} + \sqrt{(\sigma^{2} - \lambda)^{2} + 4\sigma^{2}\tau})^{2} - \lambda^{2}}{4\sigma^{2}}.
$$

## Implementing Algorithm 1.

When the marginals are Gaussian (not necessarily isotropic) and when we initialize $\psi^{j} = 0$, then, there exist some $d\times d$ matrices $A_{t}^{j}, C_{t}^{j}, \Sigma_{t}$ and $d$-dimensional vectors $b_{t}, d_{t}, e_{t}$ such that for all $t$ the iterates of Algorithm 1 satisfy:

$$
  \psi^{j}\_{t}(y) = \frac{1}{2}y^{\top}A\_{t}^{j}y - y^{\top}b\_{t}^{j} + \mathrm{const},\quad
 \phi^{j}\_{t}(y) = \frac{1}{2}x^{\top}C\_{t}^{j}x - x^{\top}d\_{t}^{j} + \mathrm{const},\quad
\mu_{t} = N(e_{t}, \Sigma_{t}).
$$

That is, the functions $\psi$ and $\phi$ are always quadratic, and $\mu_{t}$ is always Gaussian. In this case, all the integrals involved in Algorithm 1 have closed-form solutions, and we can implement Algorithm 1. Remark: for work in this spirit, see https://arxiv.org/abs/2006.02572.

## Implementing Algorithm 2.

We implement Algorithm 2 via the transform $A_{t}^{j} \mapsto A\_{t}^{j} - N\_{t}^{j}$, where $N\_{t}^{j}$ is a random positive semi-definite matrix with trace equal to $\varepsilon$ (see attached simulations). This serves as a toy model for an inexact Sinkhorn oracle for Gaussian simulations (but note that it is different from the Monte Carlo scheme that we proposed for discrete measures).

## Limitations of our experimental design.

Let us briefly summarize some limitations of our simulation setup:
- The results in our paper are proved under boundedness assumptions, while Gaussians are unbounded. While it requires justification, extending our results to unbounded cases should not cause any major problems.
- Our simulations do not cover the free-support discrete point clouds setup.

## Comments on the uploaded simulation results.

We now comment on our simulations, referring to the figures in the attached pdf:
- Figure 1 considers the toy setup $k=1$, $\nu^{1} = N(0,1)$. We observe:
    1. figure (a) shows that undamped iterates explode for $\tau < \lambda/2$ (note the absence of a blue line, because its iterates exploded);
    2. figure (b) zooms in on the case $\tau \approx \lambda/2$, noticing a sharp phase transition in the convergence behavior at the critical case $\tau = \lambda/2$;
    3. figure (c) demonstrates that damping removes the bad behavior. Moreover, the damping factor suggested in our work yields the fastest convergence among the tested five different damping factor $\eta$ choices.
- Figure 2 investigates an undamped inexact algorithm. As suggested by Theorem 2, the iterates converge up to some level governed by the accuracy parameter of the inexact Sinkhorn oracle.
- Figure 3 performs the same simulations, but this time with damping. We observe again that damping helps significantly in the critical case $\lambda/\tau = 2$, yielding a much faster convergence. The second row in Figure 3 considers the overdamped case, where we observe slightly slower convergence than that with the optimal choice of damping parameter.
- Figure 4 investigates the effect of damping in the noisy setup (Algorithm 2). First, observe that the undamped algorithm (eta = 1) explodes in the noisy case with $\tau = \lambda/2$, while this was not the case in the zero-noise setup (Figure 1 (a)). Moreover, we find that overdamping might have a negative effect on attained accuracy: larger damping factors lead to lower accuracy at convergence in noisy setups.

Meta-review: All reviewers agree this work should appear at NeurIPS and will have a welcome audience, so this is a clear accept.  In revising the paper, please implement changes promised in the rebuttal and incorporate text from the reviewer responses.  The AC also agrees with the reviewers who suggest that additional computational experiments in this work would bolster its impact and demonstrate value to practitioners.
Decision: Accept (poster)