1]N. N. Krasovski Institute of Mathematics and Mechanics, Yekaterinburg,
Russia

2]Moscow Institute of Physics and Technology, Dolgoprudny, Russia

3]Weierstrass Institute for Applied Analysis and Stochastics, Berlin, Germany

4]Skoltech, Moscow, Russia

5]Institute of information transmission problems, Moscow, Russia

6]Université catholique de Louvain (UCL), Louvain-la-Neuve, Belgium

# Decentralized Convex Optimization on Time-Varying Networks with Application
to Wasserstein Barycenters

Olga Yufereva<EMAIL_ADDRESS>Michael Persiianov
<EMAIL_ADDRESS>Pavel Dvurechensky  pavel.dvurechensky@wias-
berlin.de Alexander Gasnikov<EMAIL_ADDRESS>Dmitry Kovalev
<EMAIL_ADDRESS>[  [  [  [  [  [

###### Abstract

Inspired by recent advances in distributed algorithms for approximating
Wasserstein barycenters, we propose a novel distributed algorithm for this
problem. The main novelty is that we consider time-varying computational
networks, which are motivated by examples when only a subset of sensors can
make an observation at each time step, and yet, the goal is to average signals
(e.g., satellite pictures of some area) by approximating their barycenter. We
embed this problem into a class of non-smooth dual-friendly distributed
optimization problems over time-varying networks, and develop a first-order
method for this class. We prove non-asymptotic accelerated in the sense of
Nesterov convergence rates and explicitly characterize their dependence on the
parameters of the network and its dynamics. In the experiments, we demonstrate
the efficiency of the proposed algorithm when applied to the Wasserstein
barycenter problem.

###### keywords:

Distributed optimization dual oracle Wasserstein barycenter time-varying
networks consensus problem

##  1 Introduction

Optimal transport (OT) [1, 2] is getting more and more attention from the
machine learning and optimization community motivated by a long list of
applications such as unsupervised learning, semi-supervised learning,
clustering, text classification, image retrieval, and others, see [3] and
references therein. Given a basis space (e.g., pixel grid) and a
transportation cost function (e.g., squared Euclidean distance), the OT
approach defines a distance between two objects (e.g., images), modeled as two
probability measures on the basis space, as the minimal cost of transportation
of the first measure to the second. Such distances, in particular, the
Wasserstein distance, naturally capture the geometry of the data since they
are invariant to shifts and rotations. In particular, the Freché mean with
respect to the Wasserstein distance, called Wasserstein barycenter (WB) [4],
allows [5, 6] to reconstruct a template image from its random observations
obtained by random shifts and rotations.

The benefits of the use of WB in real-world applications are sometimes
outweighted by the large computational burden introduced by their definition.
The computation of the Wasserstein distance is already a large-scale
optimization problem, and to calculate the WB, one introduces a second
optimization level as the WB minimizes the sum of Wasserstein distances to a
set of probability measures. At this point, distributed optimization
algorithms turned out to be efficient to scale-up the computations of the WB
when the data is distributedly stored by a computational network [7, 8, 9, 10,
11, 12, 13, 14]. On the other hand, the nature of the data-generating process
may be distributed itself. In particular, it may be impossible to collect the
data even in one datacenter, especially if the data processing has to respect
the privacy of the individual data. Another example is a network of sensors
that measure signals following some distributions and for the analysis
purposes the whole network needs to find the WB of these distributions by
peer-to-peer communications. In this case, decentralized algorithms are
especially useful. Moreover, algorithms adapted to time-varying networks of
sensors or computing devices are required. For example, some nodes of the
network may be disconnected due to failures or, e.g., when the node is a
satellite that observes certain area, the observation is available only for a
certain period of time. At the same time, the development of decentralized
distributed algorithms for general optimization problems is important on its
own since the WB problem represents only one, yet important, example where
such algorithms are efficient.

Summarizing, the WB problem is important for applications, yet requires to
solve a large-scale optimization problem. For that problem, and other large-
scale problems, decentralized distributed optimization algorithms on time-
varying networks has to be developed. In this paper, we develop a general
decentralized accelerated gradient method on time-varying networks and apply
it to the WB problem.

###  1.1 Related work

In general, decentralized distributed optimization is an emerging and actively
developed branch of optimization, see the recent survey [15]. Our main focus
in this paper is on the setting of decentralized methods working on time-
varying networks, for which in the smooth and strongly convex case efficient
algorithms were recently proposed in [16, 17, 18] (primal oracle) and [19]
(dual oracle). Unfortunately, these methods are not directly applicable in our
case since the WB problem is not smooth, not strongly convex, and has simple
constraints. At the same time, the WB problem has tractable dual oracle, which
motivated us to extend the ADOM algorithm of [19] to the non-smooth, non-
strongly convex setting with simple constraints.

Starting with [8, 9] it was observed that decentralized methods with dual
oracle are well suited for the WB problem. In the cycle of subsequent papers
[9, 10, 11, 12, 13] different decentralized accelerated (randomized)
algorithms were proposed for dual WB problem. In [20] the authors propose to
reformulate the WB problem as a bilinear saddle-point problem (SPP).
Decentralized algorithm for this problem was proposed in [14]. Unfortunately,
all these algorithms are designed for static networks that do not change over
time.

Moreover, their extensions to time-varying networks seem to be hardly
possible. At the core of these algorithms lies the reformulation of the WB
problem as a problem with linear constraints ensuring the consensus between
the network nodes, and then solving the dual for that problem. If
communication network changes over time, then the affine-consensus constraints
also change over time, and so does the dual problem. This essentially requires
to solve a family of dual problems, which is not possible by the accelerated
gradient methods or the Mirror-Prox algorithm as in [9, 10, 11, 12, 13, 14].

The recent work [21] considers the WB problem on time-varying networks and
analyses a simple consensus method. The main difference with our paper is that
they prove asymptotic convergence rather than convergence rates, and they
consider possibly continuous measures on $\mathbb{R}$ unlike our setting of
discrete measures on $\mathbb{R}^{n}$. The paper [22] proposes Fenchel dual
gradient methods for distributed convex optimization over time-varying
networks. The main difference is that we propose an accelerated algorithm with
better complexity, yet under a stronger assumption on the network (they assume
the $B$-connectivity of the network).

###  1.2 Our contributions

Since the WB problem has an efficient dual oracle, a natural idea is to use
ADOM [19] that is an optimal decentralized algorithm for smooth strongly
convex unconstrained problems for time-varying networks with dual oracle. ADOM
can be considered as projected accelerated algorithm with inexact consensus-
based projection applied to specific dual reformulation of the distributed
optimization problem. Since the WB problem

  * •

(Smoothness) is not smooth;

  * •

(Constraints) is not unconstrained; as the space of probability measures has
simple constraints;

  * •

(Strongly convex) is not strongly convex;

the direct application of ADOM to the WB problem is not possible. Moreover,
not only WB problem, but also general non-smooth, non-strongly-convex
constrained optimization problems lack efficient algorithms on time-varying
networks.

The first main result of this paper is a generalization of ADOM for general
$\gamma$ strongly convex decentralized optimization problems with simple
constraints on time-varying networks; solving it numerically requires
$\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\frac{1}{\sqrt{\gamma\varepsilon}}\ln\frac{1}{\gamma\varepsilon}\right)$
iterations. The second contribution is the application to the WB problem on
time-varying networks; where the corresponding iteration number becomes
$\mathcal{O}\left(\frac{1}{\varepsilon}\ln\frac{1}{\varepsilon}\right)$. The
main ideas are the following. To obtain the strong convexity we use the
regularization of the primal problem by the entropy [3]. To deal with the
constraints and non-smoothness, we use special regularization of the dual
problem. This regularization goes back to [23, 24] and by infimal convolution
can be considered as the Moreau–Yosida smoothing of the primal problem [25,
26]. We emphasize that the proposed dual regularization (primal smoothing) was
earlier investigated only for non time-varying networks [27]. For time-varying
networks and problems with simple constraints the analysis is different.

##### Paper organization

Section 2 presents preliminaries and basic definitions for a general
distributed optimization problem. Section 3 states our main results, i.e. the
method, convergence rates and the parameter estimation. Section 4 introduces
the Wasserstein barycenter problem and specifies main results for the
particular case. Section 5 shows numerical experiments that illustrate and
verify the theoretical results. All the proofs are in Appendices.

##### Notation

We utilize the following dimensions:

  * •

$m$ for the number of individual devices (nodes),

  * •

$d$ for the dimension of a data in each device.

We use bold or normal font ($\mathbf{x}$ or $x$) for different spaces
$\mathbf{x}\in(\mathbb{R}^{d})^{m},\ x\in\mathbb{R}^{d}$. The $l$-th component
of a vector $x\in\mathbb{R}^{d}$ is denoted by $[x]_{l}$ and $l$-th component
of $\mathbf{x}\in(\mathbb{R}^{d})^{m}$ is denoted by $[\mathbf{x}]_{l}$ which
is the corresponding vector from $\mathbb{R}^{d}$.

Let $\mathbf{1}$ denote a column vector with all entries equal to $1$. The
$d$-dimensional simplex is denoted $S_{1}(d)$, that means
$S_{1}(d):=\\{p\in[0,1]^{d}\mid p^{\top}\mathbf{1}=1\\}$. For matrices $A$ and
$B,$ $A\circ B$ and $A/B$ stands for the element-wise product and division,
respectively. Another product we define as follows $\left\langle
M,X\right\rangle:=\sum_{i=1}^{d}\sum_{j=1}^{d}M_{ij}X_{ij}$.

Abbreviation WB means Wasserstein barycenter and ADOM refers to Accelerated
Decentralized Optimization Method proposed in [19].

##  2 Decentralized optimization

###  2.1 Decentralized computation problem

Decentralized computation simulates computation on distributed individual
devices. The devices are considered as nodes of an undirected connected graph
called a communication network. It means that each node can perform
computations based only on its local data and the data of its neighbors in
communication network. For a convex closed set $S$ and convex functions
$f_{i}$ decentralized computation of the following optimization problem

$\min\limits_{x\in S}\sum\limits_{i=1}^{m}f_{i}(x)$ (1)

requires numerical computation assuming that each function $f_{i}$ is stored
on the corresponding node $i\in[m]:=\\{1,2,\ldots,m\\}$. Such approach brings
us to an effective reformulation of the optimization problem.

####  2.1.1 Consensus condition

Since each computational node carries its own local data approximation, we can
substitute formally different variables $x_{i}$ for the mutual argument $x$ in
(1) assuming that they belong to the so-called consensus space. It means the
new variables $x_{i}$ must eventually coincide with each other and with the
wanted barycenter. We obtain an equivalent optimization problem in the
following form:

$\displaystyle\min\limits_{\mathbf{x}\in\mathcal{S}}F(\mathbf{x})=\min\limits_{\mathbf{x}\in\mathcal{S}}\sum\limits_{i=1}^{m}f_{i}([\mathbf{x}]_{i}),$
(2) $\displaystyle\text{where
}\mathcal{S}=\left\\{\mathbf{x}=([\mathbf{x}]_{1},\ldots,[\mathbf{x}]_{m})\in(S)^{m}\mid[\mathbf{x}]_{1}=\ldots=[\mathbf{x}]_{m}\right\\}.$

Here $i$-th component $[\mathbf{x}]_{i}$ is a corresponding $d$-dimensional
vector.

####  2.1.2 Time-varying communication network

We consider $m$ distributed devices that seek to reach a consensus solution of
an optimization problem. The devices are connected via an $m$-node network
that changes over time. At each time step $n$ we denote Laplacian of the
corresponding network by $\hat{W}_{n}$. In general, it suffices to take any
$\hat{W}_{n}$ satisfying the following:

  1. 1.

$\hat{W}_{n}$ is symmetric and positive semi-definite,

  2. 2.

$[\hat{W}_{n}]_{i,j}\neq 0$ if and only if $(i,j)$ are connected by the
network,

  3. 3.

$\ker\hat{W}_{n}=\\{(x_{1},\ldots,x_{m})\mid x_{1}=\ldots=x_{m}\\}$.

Further we use communication matrix that is the block matrix
$\mathbf{W}_{n}=\hat{W}_{n}\otimes I_{d}$. Hence, decentralized communication
of each vector $x_{i}$ stored on the $i$-th node at a time step $n$ can be
represented by multiplication of the $md$-dimensional vector
$(x_{1},\ldots,x_{m})$ and matrix $\mathbf{W}_{n}$: indeed, if
$\mathbf{y}=\mathbf{W}_{n}\mathbf{x}$, then it yields

$[\mathbf{y}]_{i}=\sum\limits_{j=1}^{m}[\hat{W}_{n}]_{i,j}[\mathbf{x}]_{j}=\sum\limits_{j\in\mathcal{N}_{i}}[{\hat{W}}_{n}]_{i,j}[\mathbf{x}]_{j},$

where $\mathcal{N}_{i}$ is the set of the neighboring nodes for the node $i$
according to the communication network at $n$-th iteration. Thus, for each
node $i$, vector $[\mathbf{y}]_{i}$ is a linear combination of vectors
$[\mathbf{x}]_{j}$, stored at the neighboring nodes $j\in\mathcal{N}_{i}$.

The considered algorithms require these communication matrices
$\mathbf{W}_{n}$ to have conditional numbers bounded for all
$n\in\\{0,1,2\ldots\\}$. Namely, we utilize the following assumption.

######  Assumption 1.

Let there exist constants $0<\lambda_{\min}^{+}<\lambda_{\max}$ such that

$\lambda_{\min}^{+}\leq\lambda_{\min}^{+}({\mathbf{W}}_{n})\leq\lambda_{\max}({\mathbf{W}}_{n})\leq\lambda_{\max}\qquad\forall
n,$

where $\lambda_{\min}^{+}({\mathbf{W}}_{n})$ is the smallest positive
eigenvalue of $\mathbf{W}_{n}$ and $\lambda_{\max}({\mathbf{W}}_{n})$ is the
biggest one.

A condition number of the matrix ${\mathbf{W}}_{n}$ is given as
$\frac{\lambda_{\max}({\mathbf{W}}_{n})}{\lambda_{\min}^{+}({\mathbf{W}}_{n})}$
and relates to the connectivity of a network; it appears in convergence rates
of many decentralized algorithms.

##  3 Main results

One of our main features is that general convex functions of optimization
problem can be defined on a convex set $S\subseteq\mathbb{R}^{d}$ instead of
the entire $\mathbb{R}^{d}$, that was crucial for duality of (10) and (11).
First of all, we study the case of strongly convex functions (Theorem 1). Then
we see that similarly we can approximate convex functions (Corollary 1). We
obtain an important Theorem 2 by applying Theorem 1 to the Wasserstein
barycenter problem; in this particular setup one can estimate parameters more
precisely. In the next sections we introduce necessary definitions, state
Theorem 2, and provide related numerical experiments. All the proofs are in
the Appendices.

Consider $\gamma$ strongly convex case, i.e. the following decentralized
optimization problem

$\displaystyle\mathbf{x}_{\gamma}^{*}=\operatorname*{arg\,min}\limits_{\mathbf{x}\in\mathcal{S}}F^{\gamma}(\mathbf{x})=\operatorname*{arg\,min}\limits_{\mathbf{x}\in\mathcal{S}}\sum\limits_{i=1}^{m}f_{i}^{\gamma}([\mathbf{x}]_{i})=\operatorname*{arg\,min}\limits_{x\in
S}\sum\limits_{i=1}^{m}f_{i}^{\gamma}(x),$ (3) $\displaystyle\mbox{where
}\mathcal{S}=\\{\mathbf{x}=([\mathbf{x}]_{1},\ldots,[\mathbf{x}]_{m})\in(S)^{m}\subset(\mathbb{R}^{d})^{m}\mid[\mathbf{x}]_{1}=\ldots=[\mathbf{x}]_{m}\\},$

functions $f_{i}^{\gamma}$ are $\gamma$ strongly convex, differentiable, and
defined on a convex set $S$. Recall that $i$-th component $[\mathbf{x}]_{i}$
is a corresponding $d$-dimensional vector.

######  Theorem 1.

Let $S\subset\mathbb{R}^{d}$ be a convex set, let functions
$f_{i}^{\gamma}\colon S\to\mathbb{R}$ of the problem (3) be $\gamma$ strongly
convex and differentiable, let $\mathbf{W}_{n}$ be the $n$-th communication
matrix satisfying Assumption 1 for some $\lambda_{\min}^{+},\lambda_{\max}>0$.
For any $r>0$, after $n$ iterations of Algorithm 1 we obtain
$\mathbf{x}^{n}_{r,\gamma}=\nabla H^{*}(\mathbf{z}^{n}_{g})$ that provides:

  1. 1.

consensus condition approximation: for each $i$ and $j$

$\left\|\left[\mathbf{x}^{n}_{r,\gamma}\right]_{i}-\left[\mathbf{x}^{n}_{r,\gamma}\right]_{j}\right\|^{2}_{2}\leq
C_{1}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n};$
(4)

  2. 2.

value approximation:

$F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-\min\limits_{\mathbf{x}\in\mathcal{S}}F^{\gamma}(\mathbf{x})\leq\frac{r}{2(1+r\gamma)}mK^{2}+C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2};$
(5)

where $K$ is such that $\|\nabla f^{\gamma}_{i}(x)\|_{2}<K$ for each $i$ and
for all $x$ from $\varepsilon/\gamma$-neighborhood of the solution
$\operatorname*{arg\,min}\limits_{x\in
S}\sum\limits_{i=1}^{m}f^{\gamma}_{i}(x)$. The parameters are
$C_{1}=\frac{(1+r\gamma)^{2}}{2\gamma^{2}}$,
$C_{2}=\frac{m(1+r\gamma)K}{\sqrt{2}\gamma}\sqrt{\frac{\lambda_{\max}}{\lambda_{\min}^{+}}}+\frac{m(1+r\gamma)^{2}}{4r\gamma^{2}}$.

###### Proof.

See Appendix B. ∎

######  Remark 1.

To reach $\varepsilon$ approximation of (4) and (5) it suffices to take
$r\leq\frac{\varepsilon}{2mK^{2}}$. Then the rate of the number of iterations
is

$n=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\sqrt{\frac{1+r\gamma}{r\gamma}}\ln\frac{C_{2}}{\varepsilon}\right)=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\frac{1}{\sqrt{\gamma\varepsilon}}\ln\frac{1}{\varepsilon}\right).$

###### Proof.

See Appendix B.5. ∎

Algorithm 1 Modified ADOM

1:input: $r>0$, for $i=1,\ldots,m$: $f^{\gamma}_{i}\colon S\to\mathbb{R}$,
$(f_{i}^{\gamma})^{*}(z)=\sup\\{\langle z,x\rangle-f_{i}^{\gamma}(x)\mid x\in
S\\}$

2:define $\nabla
h_{i}^{*}([\mathbf{z}]_{i})=\nabla(f_{i}^{\gamma})^{*}([\mathbf{z}]_{i})+r[\mathbf{z}]_{i}$

3:define $\nabla H^{*}(\mathbf{z})=\left(\nabla
h_{1}^{*}([\mathbf{z}]_{1}),\ldots,\nabla
h_{m}^{*}([\mathbf{z}]_{m})\right)^{\top}$

4:set $\alpha=\frac{r}{2}$,
$\eta=\frac{2\lambda_{\min}^{+}\sqrt{\gamma}}{7\lambda_{\max}\sqrt{r(1+r\gamma)}}$,
$\theta=\frac{\gamma}{\lambda_{\max}(1+r\gamma)}$,
$\sigma=\frac{1}{\lambda_{\max}}$,
$\tau=\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}$

5:set $\mathbf{z}^{0}=\mathbf{0}$, $\mathbf{z}_{f}^{0}=\mathbf{z}^{0}$,
$\mathbf{m}^{0}=\mathbf{0}$

6:for $n=0,1,2,\ldots$ do

7: $\mathbf{z}_{g}^{n}=\tau\mathbf{z}^{n}+(1-\tau)\mathbf{z}_{f}^{n}$

8: $\Delta^{n}=\sigma\mathbf{W}_{n}(\mathbf{m}^{n}-\eta\nabla
H^{*}(\mathbf{z}_{g}^{n}))$

9: $\mathbf{m}^{n+1}=\mathbf{m}^{n}-\eta\nabla
H^{*}(\mathbf{z}_{g}^{n})-\Delta^{n}$

10:
$\mathbf{z}^{n+1}=\mathbf{z}^{n}+\eta\alpha(\mathbf{z}_{g}^{n}-\mathbf{z}^{n})+\Delta^{n}$

11: $\mathbf{z}_{f}^{n+1}=\mathbf{z}_{g}^{n}-\theta\mathbf{W}_{n}\nabla
H^{*}(\mathbf{z}_{g}^{n})$

12:end for

13:output: $\mathbf{x}^{n}_{r,\gamma}=\nabla H^{*}(\mathbf{z}^{n}_{g})$

######  Corollary 1.

Let $S$ be a convex set in $\mathbb{R}^{d}$, let $f_{i}\colon S\to\mathbb{R}$
be differentiable convex functions for $i=1,\ldots,m$, and let
$\mathbf{W}_{n}$ be the $n$-th communication matrix satisfying Assumption 1
for some $\lambda_{\min}^{+},\lambda_{\max}>0$. Decentralized convex
optimization problem $\min_{x\in S}\sum_{i=1}^{m}f_{i}(x)$ over time-varying
communication networks can be $\varepsilon$ approximated numerically by
Algorithm 1 applied for $\gamma$ strongly convex regularizing
functions111e.g., one can take
$f_{i}^{\gamma}(x)=f_{i}(x)+\frac{\gamma}{2}\|x\|^{2}_{2}$ $f_{i}^{\gamma}(x)$
that satisfy

$0\leq\min\limits_{x\in S}\sum\limits_{i=1}^{m}f_{i}(x)-\min\limits_{x\in
S}\sum\limits_{i=1}^{m}f_{i}^{\gamma}(x)\leq\varepsilon/2$ (6)

if $r<\frac{\varepsilon}{4mK^{2}}$ where $K$ is such that $\|\nabla
f^{\gamma}_{i}(x)\|_{2}<K$ for each $i$ and for all $x$ from
$\varepsilon/\gamma$-neighborhood of the solution
$\operatorname*{arg\,min}\limits_{x\in S}\sum_{i=1}^{m}f^{\gamma}_{i}(x)$.
Moreover, if $\gamma=\sqrt{\varepsilon}$, then the rate of the number of
iterations is

$n=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\frac{1}{\varepsilon}\ln\frac{1}{\varepsilon}\right).$

###### Proof.

The condition $r<\frac{\varepsilon}{4mK^{2}}$ follows that the right-hand side
of the inequality (5) is less or equal $\varepsilon/2$, i.e.

$\sum_{i=1}^{m}f_{i}^{\gamma}([\mathbf{x}^{n}_{r,\gamma}]_{i})-\min_{x\in
S}\sum_{i=1}^{m}f_{i}^{\gamma}(x)\leq\varepsilon/2.$

Combining it with (6) we obtain $\varepsilon$ approximation of the convex
constrained optimization problem $\min\limits_{x\in
S}\sum\limits_{i=1}^{m}f_{i}(x)$. ∎

##  4 Wasserstein Barycenter Problem

Wasserstein barycenter problem is a motivating and challenging convex
constrained optimization problem that is convex, but not smooth and not
strongly-convex. The problem belongs to the optimal transport theory which
recently got various application in, e.g., texture mixing [28], statistical
estimation of template models [29], graphics and machine learning (for
regression, classification and generative modeling) [3]. Further we consider
WB problem as a decentralized convex optimization problem and propose a
computation method in Theorem 2.

###  4.1 Wasserstein distance

We provide here only necessary definitions and take into consideration only
finite-supported distributions since we deal with numerical experiments.
General theory, that begins with Wasserstein distance, can be found in [30].

Recall that we denote $d$-dimensional simplex by $S_{1}(d)$ and it represents
a set of possible probabilities distributions as
$S_{1}(d)=\\{p\in[0,1]^{d}\mid\sum_{i=1}^{d}[p]_{i}=1\\}$. Consider two
probability distributions $p,q\in S_{1}(d)$ with support on a finite set of
points $\\{\omega_{i}\in\mathbb{R}^{n}\\}_{i=1}^{d}$ such that
$p(\omega_{i})=p_{i}$ and $q(\omega_{i})=q_{i}$. Then, a cost (loss) matrix
$M$ is such that its element $[M]_{ij}\in\mathbb{R}_{+}$ represents the cost
of moving a unit mass from $\omega_{i}$ to $\omega_{j}$. So $M$ is a non-
negative symmetric matrix with zeros on the diagonal. It is often taken as the
Euclidean distances matrix, i.e. $[M]_{ij}=\|\omega_{i}-\omega_{j}\|_{2}^{2}$.
The set of transport plans is defined as

$\displaystyle U(p,q):=\left\\{X\in\mathbb{R}_{+}^{d\times
d}\;\big{|}\;X\mathbf{1}=p,X^{T}\mathbf{1}=q\right\\},$

i.e. the set of probabilities measures on $\mathbb{R}_{+}^{d\times d}$ with
margins $p$ and $q$. Wasserstein distance between two probability
distributions defines as the following minimum among component-wise
multiplication of the cost matrix and transport plans:

$\mathcal{W}(p,q):=\min_{X\in U(p,q)}\langle M,X\rangle.$

###  4.2 Wasserstein barycenter

Wasserstein barycenter of a set of probability distributions
$q_{1},\ldots,q_{m}$ is a probability distribution itself that is defined as
the solution to the following optimization problem

$\min_{p\in S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(p),$ (7)

where $\mathcal{W}_{q_{i}}(p):=\mathcal{W}(q_{i},p)$. In distributed approach,
each device (node) possesses its original distribution $q_{i}$ and the
corresponding function $f_{i}(\cdot)=\mathcal{W}_{q_{i}}(\cdot)$. The goal of
the whole system is, by communicating with each other, to approximate the
barycenter by an iterative algorithm. At each iteration each node computes a
new guess for the barycenter distribution using current guesses from its
neighbors. Typically, the first guess coincide with the original distribution
$q_{i}$ and the resulting distributions reach a consensus. It is known (see
e.g. [5]) that the WB captures the mean structure of given data. On example of
a dataset of hand-written digits ‘4’ of MNIST 784 [31] Figure 1 shows how
local nodes’ guesses change and tend to global barycenter, which resembles a
digit ‘4’ as well. That visually illustrates our theoretical result, Theorem
2. In this experiment communication networks are Erdős–Rényi random networks
and change every 5 iterations.

(a) initial data samples

(b) iteration 10

(c) iteration 50

(d) iteration 100

(e) iteration 200

Figure 1: Evolution of local data converging to Wasserstein barycenter of the
hand-written digit 4 of the MNIST784 dataset for a subset of 7 nodes out of 50
over Erdős–Rényi random networks varying each 5 iterations; regularization
parameters are $\gamma=0.03,r=0.001$

###  4.3 Resulting algorithm for WB problem

To apply the proposed numerical scheme we use entropy regularization of
Wasserstein distance $\mathcal{W}$ and make the following assumption on the
initial data.

######  Assumption 2.

Let vectors $q_{i}\in S_{1}(d)\subset\mathbb{R}^{d}$ be such that

$\min\limits_{i=1,\ldots,m\atop l=1,\ldots,d}[q_{i}]_{l}=\delta>0.$

It is not too restrictive as it only excludes zero probabilities of states,
which can be done by a little distortion. The entropy regularized Wasserstein
distance is defined as

$\displaystyle\mathcal{W}_{\gamma,q}(p):=\mathcal{W}_{\gamma}(q,p)=\min_{X\in
U(p,q)}\left\\{\langle M,X\rangle+\gamma\sum_{i=1}^{d}\sum_{j=1}^{d}X_{ij}\ln
X_{ij}\right\\},$ (8)

where $x\ln x$ is assumed to equal zero if $x=0$.

######  Theorem 2.

Let initial distributions $q_{i}$ satisfy Assumption 2 and let $p^{*}$ be
their Wasserstein barycenter, i.e. $p^{*}$ minimizes the problem (7). Let
communication matrices $\mathbf{W}_{n}$ satisfy Assumption 1 for some
$\lambda_{\min}^{+},\lambda_{\max}>0$. If Algorithm 1 is applied for entropy
regularized Wasserstein distance functions
$f_{i}^{\gamma}=\mathcal{W}_{\gamma,q_{i}}$, then functions $\nabla
h_{i}^{*}(z)$ are defined as

$\begin{array}[]{cc}\nabla
h_{i}^{*}(z)=\frac{r}{2}z+\sum\limits_{j=1}^{m}[q]_{j}\frac{\exp(\frac{1}{\gamma}([z]_{l}-M_{lj}))}{\sum\limits_{i=1}^{m}\exp(\frac{1}{\gamma}([z]_{i}-M_{ij}))},\\\
\nabla(H^{r,\gamma})^{*}(\mathbf{z})=(\nabla
h_{1}^{*}([\mathbf{z}]_{1}),\ldots,\nabla
h_{m}^{*}([\mathbf{z}]_{m}))^{\top},\end{array}$ (9)

and it suffices to take $\gamma=\frac{1}{8}\varepsilon\ln d$,
$K^{2}=\sum\limits_{j=1}^{d}\left(2\gamma\ln
d+\inf_{i}\sup_{l}|M_{jl}-M_{il}|-\gamma\ln\frac{\delta}{2}\\}\right)^{2}$,
and $r=\frac{\varepsilon}{4mK^{2}}$ to $\varepsilon$ approximate the solution
$p^{*}$ as follows

$\displaystyle\left|\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}([\mathbf{x}^{n}_{r,\gamma}]_{i})-\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(p^{*})\right|$
$\displaystyle\leq 2\gamma\ln
d+\frac{r}{4(1+r\gamma)}mK^{2}+C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}\leq\varepsilon.$

Thus, a sufficient number of iterations of Algorithm 1 is

$n=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\sqrt{\frac{1+r\gamma}{r\gamma}}\ln\frac{C}{\varepsilon}\right)=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\frac{1}{\varepsilon}\ln\frac{1}{\varepsilon}\right).$

##  5 Numerical Experiments

We provide numerical experiments to demonstrate performance of the proposed
method. We can fruitfully test our method on WB problem with an artificial set
of univariate, discrete and truncated Gaussian distributions. For such a
dataset, the resulting distribution (the zero-entropy Wasserstein barycenter)
is described by an analytic formula. Namely, the barycenter is a Gaussian
distribution which mean is the arithmetic average of the means of the given
Gaussians and the standard deviation of the barycenter is the arithmetic
average of the standard deviations of the given Gaussians.

For all figures of this section we generated a dataset of truncated Gaussians.
Each distribution size is $100$, while a size of a dataset (i.e. the number of
nodes) differs and is indicated at each figure. For entropy regularized
Wasserstein distance we use normalized Euclidean cost matrix and entropy
regularization parameter $\gamma=0.01$. The regularization parameter $r$ of
the method is $r=0.001$.

###  5.1 Comparison with other methods

To the best of our knowledge, we can compare our method (called here ADOM)
with local barycenters method (LB) proposed in [21] and Fenchel dual gradient
method (FDGM) proposed in [22]. They all are applicable for the WB problem on
time-varying networks. Regardless of the analytical form of the LB algorithm,
in order to implement it we need either to solve an optimization problem at
each iteration or to use approximations, e.g. methods of [32]. Note that, for
any particular setup, realization of FDGM is quite a problem, since the method
is sensitive to the step size $\alpha_{n}$. ADOM negotiates limitations
described above and reveals relatively stable convergence as it is presented
at Figures 2–3, where we test the methods on cycle networks that change every
iteration.

Figure 2:
$\frac{1}{m}\left(\sum\limits_{i=1}^{m}\mathcal{W}(q_{i},[\mathbf{x}^{n}_{r,\gamma}]_{i})-\sum\limits_{i=1}^{m}\mathcal{W}(q_{i},p^{*})\right)$-convergence
comparison on cycle networks changing every iteration Figure 3: Consensus
condition comparison of methods on cycle networks changing every iteration

###  5.2 Performance rate

We can see the difference in error value, i.e. in
$\frac{1}{m}\left(\sum\limits_{i=1}^{m}\mathcal{W}(q_{i},[\mathbf{x}^{n}_{r,\gamma}]_{i})-\sum\limits_{i=1}^{m}\mathcal{W}(q_{i},p^{*})\right)$,
in two opposite cases: Figure 4a presents evolution of errors computed for
constant networks of different topologies, while at Figure 4b networks change
every iterations within indicated topology; the exception is complete network
that cannot change. One natural way to sample a random network is to
independently sample each edge with a probability $p$. Such networks are
called Erdős–Rényi network or $(m,p)$-Erdős–Rényi network, where $m$ is the
number of nodes and $p$ is the probability of an edge. Let us notice also that
star, cycle and minimum spanning tree (of $(m,0.9)$-Erdős–Rényi) networks have
$m-1,m,$ and $m+1$ edges respectively in contrast to the complete network with
$m(m-1)/2$ edges and $(m,0.5)$-Erdős–Rényi network that has $m(m-1)/4$ edges
in average.

(a) Constant networks

(b) Time-varying networks

Figure 4: Different network topologies: error over time

For the two ‘most efficient’ topologies that prove themselves at Figures 4a–4b
we compute at Figure 5 the error evolution for different frequency of the
networks varying. We indicate the lengths of epoch, i.e. number of iteration
between network changing. Notice that the evolution for constant networks,
computed above, matches to infinite epoch length. The number of iterations
remain the same on all figures despite it is insufficient for convergence on
cycle networks. Nonetheless one can see the trends of convergence and notice
that there is no monotonicity with respect to frequency of networks varying.

(a) Cycle networks

(b) $(m,0.9)$-Erdős–Rényi networks

Figure 5: Different number of iteration between network changing: errors over
time

## Acknowledgments

The authors are grateful to Alexander Rogozin.

The work of A. Gasnikov in Section 4 of the paper was funded by Russian
Science Foundation (project 18-71-10108).

The work of O. Yufereva in the rest part of the paper was performed as part of
research conducted in the Ural Mathematical Center with the financial support
of the Ministry of Science and Higher Education of the Russian Federation
(Agreement number 075-02-2023-913).

## References

  *   * Monge [1781] Monge, G.: Mémoire sur la théorie des déblais et des remblais. Histoire de l’Académie Royale des Sciences de Paris (1781)
  * Kantorovich [1942] Kantorovich, L.: On the translocation of masses. (Doklady) Acad. Sci. URSS (N.S.) 37, 199–201 (1942)
  * Peyré et al. [2019] Peyré, G., Cuturi, M., et al.: Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning 11(5-6), 355–607 (2019)
  * Agueh and Carlier [2011] Agueh, M., Carlier, G.: Barycenters in the Wasserstein space. SIAM Journal on Mathematical Analysis 43(2), 904–924 (2011)
  * Cuturi and Doucet [2014] Cuturi, M., Doucet, A.: Fast computation of wasserstein barycenters. In: International Conference on Machine Learning, pp. 685–693 (2014). PMLR
  * del Barrio et al. [1999] Barrio, E., Gine, E., Matran, C.: Central limit theorems for the wasserstein distance between the empirical and the true distributions. The Annals of Probability 27(2), 1009–1071 (1999)
  * Staib et al. [2017] Staib, M., Claici, S., Solomon, J.M., Jegelka, S.: Parallel streaming wasserstein barycenters. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30, pp. 2647–2658. Curran Associates, Inc., ??? (2017). http://papers.nips.cc/paper/6858-parallel-streaming-wasserstein-barycenters.pdf
  * Uribe et al. [2018] Uribe, C.A., Dvinskikh, D., Dvurechensky, P., Gasnikov, A., Nedić, A.: Distributed computation of wasserstein barycenters over networks. In: 2018 IEEE Conference on Decision and Control (CDC), pp. 6544–6549 (2018). IEEE
  * Dvurechenskii et al. [2018] Dvurechenskii, P., Dvinskikh, D., Gasnikov, A., Uribe, C., Nedich, A.: Decentralize and randomize: Faster algorithm for wasserstein barycenters. Advances in Neural Information Processing Systems 31 (2018)
  * Kroshnin et al. [2019] Kroshnin, A., Tupitsa, N., Dvinskikh, D., Dvurechensky, P., Gasnikov, A., Uribe, C.: On the complexity of approximating wasserstein barycenters. In: International Conference on Machine Learning, pp. 3530–3540 (2019). PMLR
  * Dvinskikh et al. [2019] Dvinskikh, D., Gorbunov, E., Gasnikov, A., Dvurechensky, P., Uribe, C.A.: On primal and dual approaches for distributed stochastic convex optimization over networks. In: 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 7435–7440 (2019). IEEE
  * Krawtschenko et al. [2020] Krawtschenko, R., Uribe, C.A., Gasnikov, A., Dvurechensky, P.: Distributed optimization with quantization for computing wasserstein barycenters. arXiv preprint arXiv:2010.14325 (2020)
  * Dvinskikh [2021] Dvinskikh, D.: Decentralized algorithms for wasserstein barycenters. PhD thesis, Humboldt Universitaet zu Berlin (Germany) (2021)
  * Rogozin et al. [2021] Rogozin, A., Beznosikov, A., Dvinskikh, D., Kovalev, D., Dvurechensky, P., Gasnikov, A.: Decentralized distributed optimization for saddle point problems. arXiv preprint arXiv:2102.07758 (2021)
  * Gorbunov et al. [2022] Gorbunov, E., Rogozin, A., Beznosikov, A., Dvinskikh, D., Gasnikov, A.: In: Nikeghbali, A., Pardalos, P.M., Raigorodskii, A.M., Rassias, M.T. (eds.) Recent Theoretical Advances in Decentralized Distributed Convex Optimization, pp. 253–325. Springer, Cham (2022). https://doi.org/10.1007/978-3-031-00832-0_8
  * Rogozin et al. [2021] Rogozin, A., Bochko, M., Dvurechensky, P., Gasnikov, A., Lukoshkin, V.: An accelerated method for decentralized distributed stochastic optimization over time-varying graphs. In: 2021 60th IEEE Conference on Decision and Control (CDC), pp. 3367–3373 (2021). https://doi.org/10.1109/CDC45484.2021.9683110
  * Li and Lin [2021] Li, H., Lin, Z.: Accelerated gradient tracking over time-varying graphs for decentralized optimization. arXiv preprint arXiv:2104.02596 (2021)
  * Kovalev et al. [2021a] Kovalev, D., Gasanov, E., Gasnikov, A., Richtarik, P.: Lower bounds and optimal algorithms for smooth and strongly convex decentralized optimization over time-varying networks. Advances in Neural Information Processing Systems 34 (2021)
  * Kovalev et al. [2021b] Kovalev, D., Shulgin, E., Richtárik, P., Rogozin, A.V., Gasnikov, A.: ADOM: accelerated decentralized optimization method for time-varying networks. In: International Conference on Machine Learning, pp. 5784–5793 (2021). PMLR
  * Dvinskikh and Tiapkin [2021] Dvinskikh, D., Tiapkin, D.: Improved complexity bounds in Wasserstein barycenter problem. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, pp. 1738–1746 (2021). PMLR
  * Bishop and Doucet [2021] Bishop, A.N., Doucet, A.: Network consensus in the wasserstein metric space of probability measures. SIAM Journal on Control and Optimization 59(5), 3261–3277 (2021)
  * Wu and Lu [2019] Wu, X., Lu, J.: Fenchel dual gradient methods for distributed convex optimization over time-varying networks. IEEE Transactions on Automatic Control 64(11), 4629–4636 (2019) https://doi.org/10.1109/TAC.2019.2901829
  * Devolder et al. [2012] Devolder, O., Glineur, F., Nesterov, Y.: Double smoothing technique for large-scale linearly constrained convex optimization. SIAM Journal on Optimization 22(2), 702–727 (2012)
  * Gasnikov et al. [2016] Gasnikov, A.V., Gasnikova, E., Nesterov, Y.E., Chernov, A.: Efficient numerical methods for entropy-linear programming problems. Computational Mathematics and Mathematical Physics 56(4), 514–524 (2016)
  * Rockafellar [1997] Rockafellar, R.T.: Convex Analysis vol. 11. Princeton university press, Princeton (1997)
  * Lemaréchal and Sagastizábal [1997] Lemaréchal, C., Sagastizábal, C.: Practical aspects of the moreau–yosida regularization: Theoretical preliminaries. SIAM journal on optimization 7(2), 367–385 (1997)
  * Uribe et al. [2020] Uribe, C.A., Lee, S., Gasnikov, A., Nedić, A.: A dual approach for optimal algorithms in distributed optimization over networks. In: 2020 Information Theory and Applications Workshop (ITA), pp. 1–37 (2020). IEEE
  * Rabin et al. [2011] Rabin, J., Peyré, G., Delon, J., Bernot, M.: Wasserstein barycenter and its application to texture mixing. In: International Conference on Scale Space and Variational Methods in Computer Vision, pp. 435–446 (2011). Springer
  * Boissard et al. [2015] Boissard, E., Le Gouic, T., Loubes, J.-M.: Distribution’s template estimate with wasserstein metrics. Bernoulli 21(2), 740–759 (2015)
  * Villani [2009] Villani, C.: Optimal Transport: Old and New vol. 338. Springer, Cham (2009)
  * LeCun [1998] LeCun, Y.: The mnist database of handwritten digits. http://yann. lecun. com/exdb/mnist/ (1998)
  * Flamary et al. [2021] Flamary, R., Courty, N., Gramfort, A., Alaya, M.Z., Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L., Gayraud, N.T.H., Janati, H., Rakotomamonjy, A., Redko, I., Rolet, A., Schutz, A., Seguy, V., Sutherland, D.J., Tavenard, R., Tong, A., Vayer, T.: Pot: Python optimal transport. Journal of Machine Learning Research 22(78), 1–8 (2021)
  * Cuturi and Peyré [2015] Cuturi, M., Peyré, G.: A smoothed dual approach for variational Wasserstein problems. SIAM Journal on Imaging Sciences (2015)
  * Bigot et al. [2019] Bigot, J., Cazelles, E., Papadakis, N.: Data-driven regularization of Wasserstein barycenters with an application to multivariate density registration. Information and Inference: A Journal of the IMA 8(4), 719–755 (2019)

##  Appendix A ADOM and its assumptions

The state of the art numerical computation method for time-varying networks,
called ADOM, is developed in [19] and this subsection is to present its main
objects. It has natural restrictions on the class of suitable problems and,
e.g., Wasserstein barycenter problem lies beyond the requirements of this
algorithm. So we modify ADOM to solve more general optimization problems with
restrictions. For the sake of consistency, we slightly change original
notation and adduce below the results from [19].

In [19], optimization problem with the consensus condition is

$\displaystyle\min\limits_{\mathbf{x}\in\mathcal{R}}H(\mathbf{x})=\min\limits_{\mathbf{x}\in\mathcal{R}}\sum\limits_{i=1}^{m}h_{i}([\mathbf{x}]_{i}),$
(10) $\displaystyle\text{where
}\mathcal{R}=\left\\{\mathbf{x}=([\mathbf{x}]_{1},\ldots,[\mathbf{x}]_{m})\in(\mathbb{R}^{d})^{m}\mid[\mathbf{x}]_{1}=\ldots=[\mathbf{x}]_{m}\right\\},$

where functions $h_{i}\colon\mathbb{R}^{d}\to\mathbb{R}$ are assumed to be
smooth and strongly convex. Problem (10) is equivalent to the following:

$\displaystyle\min\limits_{\mathbf{z}\in\mathcal{R}^{\perp}}H^{*}(\mathbf{z}),$
(11) $\displaystyle\text{where
}\mathcal{R}^{\perp}=\left\\{\mathbf{z}=([\mathbf{z}]_{1},\ldots,[\mathbf{z}]_{m})\in(\mathbb{R}^{d})^{m}\;\Big{|}\;\sum\limits_{i=1}^{m}[\mathbf{z}]_{i}=0\right\\},$

where $H^{*}$ is the Fenchel transform of the function $H$ and
$\mathcal{R}^{\perp}$ is the orthogonal complement of $\mathcal{R}$, that
exists since $S=\mathbb{R}^{d}$ here.

######  Theorem 3 ([19, Theorem 1]).

Let functions $h_{i}\colon\mathbb{R}^{d}\to\mathbb{R}$ be $L$ smooth and $\mu$
strongly convex, $\mathbf{x}^{*}$ be the solution of the optimization problem
(10), $\mathbf{W}_{n}$ be a communication matrix at the $n$-th iteration
satisfying Assumption 1. Set parameters $\alpha,\eta,\theta,\sigma,\tau$ of
Algorithm 2 to $\alpha=\frac{1}{2L}$, $\eta=\frac{2\lambda_{\min}^{+}\sqrt{\mu
L}}{7\lambda_{\max}}$, $\theta=\frac{\mu}{\lambda_{\max}}$,
$\sigma=\frac{1}{\lambda_{\max}}$, and
$\tau=\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu}{L}}$. Then
there exists $C>0$, such that for Fenchel conjugate function
$H^{*}(\mathbf{z})$ from (11)

$\left\|\nabla H^{*}(\mathbf{z}_{g}^{n})-\mathbf{x}^{*}\right\|^{2}_{2}\leq
C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu}{L}}\right)^{n}.$
(12)

######  Remark 2.

Addressing details of the proof of Theorem 1 of [19] we see that there is a
particular choice of the constant $C$, namely

$C=\max\left\\{\frac{2\tau}{\mu^{2}},\frac{\tau(1-\tau)L}{\eta(1-\eta\alpha)\mu^{2}}\right\\}=\frac{1}{\mu^{2}}\max\left\\{\frac{2\lambda_{\min}^{+}\sqrt{\mu}}{7\lambda_{\max}\sqrt{L}},\frac{1}{2}\right\\}=\frac{1}{2\mu^{2}}.$
(13)

It means that the actual convergence rate is
$n=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\sqrt{\frac{L}{\mu}}\ln\frac{1}{\mu^{2}\varepsilon}\right)$.

Algorithm 2 ADOM: Accelerated Decentralized Optimization Method

1:input: $\nabla H^{*}\colon(\mathbb{R}^{d})^{m}\to\mathbb{R}$,
$\mathbf{z}^{0}\in\mathcal{R}^{\perp},\mathbf{m}^{0}\in(\mathbb{R}^{d})^{\mathcal{V}},\alpha,\eta,\theta,\sigma\\!>\\!0,\tau\\!\in\\!(0,1)$

2:set $\mathbf{z}_{f}^{0}=\mathbf{z}^{0}$

3:for $k=0,1,2,\ldots$ do

4: $\mathbf{z}_{g}^{n}=\tau\mathbf{z}^{n}+(1-\tau)\mathbf{z}_{f}^{n}$

5: $\Delta^{n}=\sigma\mathbf{W}_{n}(\mathbf{m}^{n}-\eta\nabla
H^{*}(\mathbf{z}_{g}^{n}))$

6: $\mathbf{m}^{n+1}=\mathbf{m}^{n}-\eta\nabla
H^{*}(\mathbf{z}_{g}^{n})-\Delta^{n}$

7:
$\mathbf{z}^{n+1}=\mathbf{z}^{n}+\eta\alpha(\mathbf{z}_{g}^{n}-\mathbf{z}^{n})+\Delta^{n}$

8: $\mathbf{z}_{f}^{n+1}=\mathbf{z}_{g}^{n}-\theta\mathbf{W}_{n}\nabla
H^{*}(\mathbf{z}_{g}^{n})$

9:end for

##  Appendix B Proof of Theorem 1

All the arguments below are applied under assumptions of Theorem 1, i.e. we
assume that $S\subset\mathbb{R}^{d}$ is a convex set,
$\mathbf{x}\in\mathcal{S}$ is equivalent to $[\mathbf{x}]_{i}\in S$ for all
$i=1,\ldots,m$, functions $f^{\gamma}_{i}\colon S\to\mathbb{R}$ are $\gamma$
strongly convex, and the output of Algorithm 1 is
$\mathbf{x}^{n}_{r,\gamma}=\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})$.
Denote also

$\displaystyle\mathbf{x}_{\gamma}^{*}=(x^{*}_{\gamma},\ldots,x^{*}_{\gamma})=\operatorname*{arg\,min}\limits_{\mathbf{x}\in\mathcal{S}}F^{\gamma}(\mathbf{x})=\operatorname*{arg\,min}\limits_{x\in
S}\sum\limits_{i=1}^{m}f_{i}^{\gamma}(x).$

###  B.1 Derivation of $(H^{r,\gamma})^{*}$

In brief, in this subsection we show that functions $h_{i}^{r,\gamma}$ from
(14) are $\frac{1}{r}$ smooth, $\frac{\gamma}{1+r\gamma}$ strongly convex, and
such that $\nabla(H^{r,\gamma})^{*}$ from Line 3 of Algorithm 1 is the
gradient of the conjugate function $(H^{r,\gamma})^{*}$ of
$H^{r,\gamma}=\sum\limits_{i=1}^{m}h_{i}^{r,\gamma}$ from (14). Then the
consensus condition (4) becomes a corollary of Theorem 3 with $L=\frac{1}{r}$
and $\mu=\frac{\gamma}{1+r\gamma}$.

From now on let functions $h_{i}^{r,\gamma}\colon\mathbb{R}^{d}\to\mathbb{R}$
and $H^{r,\gamma}\colon(\mathbb{R}^{d})^{m}\to\mathbb{R}$ be

$\begin{array}[]{cc}H^{r,\gamma}(\mathbf{x})=\sum\limits_{i=1}^{m}h_{i}^{r,\gamma}([\mathbf{x}]_{i}),\
\mbox{where}\\\ h_{i}^{r,\gamma}(x)=\inf\limits_{y\in
S}\left\\{f_{i}^{\gamma}(y)+\frac{1}{2r}\|y-x\|^{2}_{2}\right\\}.\end{array}$
(14)

Define their conjugate as $(h_{i}^{r,\gamma})^{*}$ and $(H^{r,\gamma})^{*}$.

######  Lemma 1.

If functions $h_{i}^{r,\gamma}$ and $H^{r,\gamma}$ are defined by (14), then
their Fenchel conjugate functions $(h_{i}^{r,\gamma})^{*}$ and
$(H^{r,\gamma})^{*}\colon(\mathbb{R}^{d})^{m}\to\mathbb{R}$ are

$\begin{array}[]{cc}(H^{r,\gamma})^{*}(\mathbf{z})=\sum\limits_{i=1}^{m}(h_{i}^{r,\gamma})^{*}([\mathbf{z}]_{i}),\
\mbox{where}\\\
(h_{i}^{r,\gamma})^{*}(z)=(f_{i}^{\gamma})^{*}(z)+\frac{r}{2}\|z\|^{2}_{2}.\end{array}$

Moreover, its conjugate $(H^{r,\gamma})^{**}$ coincides with $H^{r,\gamma}$.

###### Proof.

The definition (14) is similar to Moreau–Yosida smoothing, but the tricky
point is that the functions $f_{i}^{\gamma}$ are defined on a convex set $S$
instead of the $\mathbb{R}^{d}$. Let us introduce functions
$\tilde{f}_{i}^{\gamma}$ with domain $\mathbb{R}^{d}$ as follows:

$\tilde{f}_{i}^{\gamma}(x)=\left\\{\begin{array}[]{cc}&f_{i}^{\gamma}(x)\quad\mbox{
if $x\in S$}\\\ &+\infty,\quad\mbox{ otherwise}.\end{array}\right.$ (15)

Such $\tilde{f}_{i}^{\gamma}$ are $\gamma$ strongly convex as well. Moreover,
substitution $\tilde{f}_{i}^{\gamma}$ for $f_{i}^{\gamma}$ affect neither
primal $h_{i}^{r,\gamma}$:

$h_{i}^{r,\gamma}(x)=\inf\limits_{y\in
S}\left\\{f_{i}^{\gamma}(y)+\frac{1}{2r}\|y-x\|^{2}_{2}\right\\}=\inf\limits_{y\in\mathbb{R}^{d}}\left\\{\tilde{f}_{i}^{\gamma}(y)+\frac{1}{2r}\|y-x\|^{2}_{2}\right\\},$

nor $(f_{i}^{\gamma})^{*}(z)+\frac{r}{2}\|z\|^{2}_{2}$:

$\displaystyle(f^{\gamma}_{i})^{*}(z)+\frac{r}{2}\|z\|^{2}_{2}=\max\limits_{x\in
S}\left\\{\left\langle
z,x\right\rangle-f^{\gamma}_{i}(x)\right\\}+\frac{r}{2}\|z\|^{2}_{2}$
$\displaystyle=\max\limits_{x\in\mathbb{R}^{d}}\left\\{\left\langle
z,x\right\rangle-\tilde{f}^{\gamma}_{i}(x)\right\\}+\frac{r}{2}\|z\|^{2}_{2}=(\tilde{f}^{\gamma}_{i})^{*}(z)+\frac{r}{2}\|z\|^{2}_{2}.$

For each $i$ one can see that
$(h_{i}^{r,\gamma})^{*}=(f^{\gamma}_{i})^{*}(z)+\frac{r}{2}\|z\|^{2}_{2}$ is
the Fenchel conjugate of $h_{i}^{r,\gamma}$ and vice versa. Indeed, for
proper, convex and lower semicontinuous
$g_{1},g_{2}\colon\mathbb{R}^{d}\to\mathbb{R}$ we have
$(g_{1}+g_{2})^{*}(x)=g_{1}^{*}\square g_{2}^{*}$ and $(g_{1}\square
g_{2})^{*}=g_{1}^{*}+g_{2}^{*}$, where $(g_{1}\square g_{2})(x)$ means the
convolution $\inf\\{g_{1}(y)+g_{2}(x-y)\mid y\in\mathbb{R}^{d}\\}$.

Hence the Fenchel conjugate for the function $H^{r,\gamma}$ will be

$\displaystyle\sup\limits_{\mathbf{x}\in(\mathbb{R}^{d})^{m}}\left\\{\langle\mathbf{z},\mathbf{x}\rangle-H^{r,\gamma}(\mathbf{x})\right\\}$
$\displaystyle=\sup\limits_{\mathbf{x}\in(\mathbb{R}^{d})^{m}}\left\\{\sum\limits_{i=1}^{m}\left(\langle[\mathbf{z}]_{i},[\mathbf{x}]_{i}\rangle-
h_{i}^{r,\gamma}([\mathbf{x}]_{i})\right)\right\\}$ (16)
$\displaystyle=\sum\limits_{i=1}^{m}\sup\limits_{[\mathbf{x}]_{i}\in\mathbb{R}^{d}}\left\\{\langle[\mathbf{z}]_{i},[\mathbf{x}]_{i}\rangle-
h_{i}^{r,\gamma}([\mathbf{x}]_{i})\right\\}$
$\displaystyle=\sum\limits_{i=1}^{m}(h_{i}^{r,\gamma})^{*}([\mathbf{z}]_{i})=(H^{r,\gamma})^{*}(\mathbf{z}).$

In the same way one can see that $H^{r,\gamma}$ and $(H^{r,\gamma})^{**}$
coincide. ∎

######  Remark 3.

For each $i$ the function $\left(h^{r,\gamma}_{i}\right)^{*}$ from (14) is
$\left(\frac{1}{\gamma}+r\right)$ smooth and $r$ strongly convex by
definition, so we have $h^{r,\gamma}_{i}=(h^{r,\gamma}_{i})^{**}$ being
$\frac{1}{r}$ smooth and $\frac{\gamma}{1+r\gamma}$ strongly convex. In
addition

$\nabla(h^{r,\gamma}_{i})^{*}(z)=\nabla(f_{i}^{\gamma})^{*}(z)+z$

as stated in Line 3 of Algorithm 1. Then we can apply Algorithm 2 for
$L=r^{-1}$ smooth and $\mu=\frac{\gamma}{1+r\gamma}$ strongly convex functions
$h^{r,\gamma}_{i}$ and get the values of $\nabla(h^{r,\gamma}_{i})^{*}(z)$ as
output.

Thus we construct a relaxation
$\min_{\mathbf{x}\in\mathcal{R}}H^{r,\gamma}(\mathbf{x})$ of the constrained
convex optimization problem
$\min_{\mathbf{x}\in\mathcal{S}}F^{\gamma}(\mathbf{x})$.

######  Corollary 2.

Let a function $H^{r,\gamma}$ be defined in (14) and let
$\mathbf{x}^{*}_{r,\gamma}=\operatorname*{arg\,min}\limits_{\mathbf{x}\in\mathcal{R}}H^{r,\gamma}(\mathbf{x})$.
Then applying Algorithm 2 for

$\nabla(h^{r,\gamma}_{i})^{*}(z)=(f_{i}^{\gamma})^{*}(z)+rz$

we get by Theorem 3

$\left\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{n}_{r,\gamma}\right\|^{2}_{2}\leq
C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n},$
(17)

where $\mathbf{x}_{r,\gamma}^{n}=\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})$
and

$C=\frac{(1+r\gamma)^{2}}{2\gamma^{2}}.$

Moreover, since $\mathbf{x}^{*}_{r,\gamma}\in\mathcal{R}$, i.e.
$[\mathbf{x}^{*}_{r,\gamma}]_{i}=[\mathbf{x}^{*}_{r,\gamma}]_{j}$ for all $i$
and $j$, the consensus condition is approximated as follows

$\left\|\left[\mathbf{x}^{n}_{r,\gamma}\right]_{i}-\left[\mathbf{x}^{n}_{r,\gamma}\right]_{j}\right\|^{2}_{2}\leq
2C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}.$

###  B.2 Value bounds on $H^{r,\gamma}$

Despite we defined $h_{i}^{r,\gamma}$ for all $\mathbb{R}^{d}$, some
properties hold true on the initial set $S$ only.

######  Lemma 2.

Let functions $h_{i}^{r,\gamma}$ be defined in (14). If $x\in S$, then for any
$r>0$, for each $i=1,\ldots,m$ we have

$\begin{array}[]{rl}f^{\gamma}_{i}(x)-\frac{r}{2(1+r\gamma)}\left\|\nabla
f^{\gamma}_{i}(x)\right\|^{2}_{2}\leq h_{i}^{r,\gamma}(x)\leq
f^{\gamma}_{i}(x).\end{array}$ (18)

###### Proof.

The second inequality directly follows from the definition (14). To prove the
first one we recall that $f^{\gamma}_{i}$ is $\gamma$ strongly convex and the
following holds:

$\displaystyle h_{i}^{r,\gamma}(x)=\inf\limits_{y\in
S}\left\\{f^{\gamma}_{i}(y)+(2r)^{-1}\|x-y\|^{2}_{2}\right\\}$
$\displaystyle=\inf\limits_{y:\ (x-y)\in
S}\left\\{f^{\gamma}_{i}(x-y)+(2r)^{-1}\|y\|^{2}_{2}\right\\}$
$\displaystyle\geq\inf\limits_{y:\ (x-y)\in
S}\left\\{f^{\gamma}_{i}(x)+\langle\nabla
f^{\gamma}_{i}(x),-y\rangle+\gamma/2\|y\|^{2}_{2}+(2r)^{-1}\|y\|^{2}_{2}\right\\}$
$\displaystyle\geq\inf\limits_{y\in\mathbb{R}^{d}}\left\\{f^{\gamma}_{i}(x)+\langle\nabla
f^{\gamma}_{i}(x),-y\rangle+\gamma/2\|y\|^{2}_{2}+(2r)^{-1}\|y\|^{2}_{2}\right\\},$

which reaches its minimum at $y=\frac{r}{1+r\gamma}\nabla f^{\gamma}_{i}(x)$
and so equals to

$\displaystyle f^{\gamma}_{i}(x)+\frac{r}{1+r\gamma}\langle-\nabla
f^{\gamma}_{i}(x),\nabla
f^{\gamma}_{i}(x)\rangle+\frac{r}{2(1+r\gamma)}\|\nabla
f^{\gamma}_{i}(x)\|^{2}_{2}$
$\displaystyle=f^{\gamma}_{i}(x)-\frac{r}{2(1+r\gamma)}\|\nabla
f^{\gamma}_{i}(x)\|^{2}_{2}.$

∎

###  B.3 Convergence in argument

Lemma 3 shows convergence in argument in the following sense: if the
regularization parameter $r$ tends to zero, the argminimum
$\mathbf{x}^{*}_{r,\gamma}\in\mathcal{R}$ of $H^{r,\gamma}$ tends to the
argminimum $\mathbf{x}^{*}_{\gamma}\in\mathcal{S}$ of $F^{\gamma}$. By
Corollary 2 we have $\mathbf{x}^{*}_{r,\gamma}\in\mathcal{R}$ approximated by
$\mathbf{x}^{n}_{r,\gamma}\in(\mathbb{R}^{d})^{m}$ for a sufficient number of
iterations $n$.

######  Lemma 3.

Let
$\mathbf{x}^{*}_{r,\gamma}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{R}}H^{r,\gamma}(\mathbf{x})$
for $H^{r,\gamma}$ defined in (14). Let

$\|\nabla F^{\gamma}(\mathbf{x})\|_{2}^{2}\leq
mK_{\zeta}^{2}\quad\forall\mathbf{x}\in\\{\mathbf{y}\in\mathcal{S}\mid\|\mathbf{y}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\zeta\\}.$
(19)

If $r$ is such that
$\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\zeta$, then

$\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\sqrt{\frac{rm}{2\gamma}}K_{\zeta}.$
(20)

###### Proof.

Using (18) and strong convexity of $F^{\gamma}$ and $H^{r,\gamma}$ we have

$\displaystyle F^{\gamma}(\mathbf{x}^{*}_{\gamma})\geq
H^{r,\gamma}(\mathbf{x}^{*}_{\gamma})=\sum
h_{i}^{r,\gamma}([\mathbf{x}^{*}_{\gamma}]_{i})$
$\displaystyle\geq\sum\limits_{i=1}^{m}\left(h_{i}^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})+\frac{\gamma}{2(1+r\gamma)}\|[\mathbf{x}^{*}_{r,\gamma}]_{i}-[\mathbf{x}^{*}_{\gamma}]_{i}\|_{2}^{2}\right)$
$\displaystyle=H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})+\frac{\gamma}{2(1+r\gamma)}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}$
$\displaystyle\geq
F^{\gamma}(\mathbf{x}^{*}_{r,\gamma})-\frac{r}{2(1+r\gamma)}\|\nabla
F^{\gamma}(\mathbf{x}^{*}_{r,\gamma})\|^{2}_{2}+\frac{\gamma}{2(1+r\gamma)}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}$
$\displaystyle\geq
F^{\gamma}(\mathbf{x}^{*}_{r,\gamma})-\frac{r}{2(1+r\gamma)}mK_{\zeta}^{2}+\frac{\gamma}{2(1+r\gamma)}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}$
$\displaystyle\geq
F^{\gamma}(\mathbf{x}^{*}_{\gamma})+\gamma/2\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}-\frac{r}{2(1+r\gamma)}mK_{\zeta}^{2}+\frac{\gamma}{2(1+r\gamma)}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}$
$\displaystyle\geq
F^{\gamma}(\mathbf{x}^{*}_{\gamma})+\frac{\gamma}{1+r\gamma}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}-\frac{r}{2(1+r\gamma)}mK_{\zeta}^{2}.$

Then
$\frac{\gamma}{1+r\gamma}\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}^{2}-\frac{r}{2(1+r\gamma)}mK_{\zeta}^{2}\leq
0$ and hence
$\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|^{2}_{2}\leq\frac{rm}{2\gamma}K^{2}_{\zeta}$.
∎

Combining Lemma 3 with Corollary 2 we get the following.

######  Remark 4.

Let $\zeta>0$ and let $K_{\zeta}$ be such that (19) holds. If

$\displaystyle\sqrt{\frac{rm}{2\gamma}}K_{\zeta}+\sqrt{C_{1}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}\leq\zeta,$

where $C_{1}=\frac{(1+r\gamma)^{2}}{2\gamma^{2}}$, then both
$\|\mathbf{x}^{*}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\zeta$ and
$\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\zeta$ hold.

###  B.4 Value approximation

Let $\mathbf{x}^{*}_{r,\gamma}\in\mathcal{R}$ be the only argminimum of
$H^{r,\gamma}$ on the consensus space $\mathcal{R}$, i.e.

$\begin{array}[]{cc}\mathbf{x}^{*}_{r,\gamma}=\operatorname*{arg\,min}\limits_{\mathbf{x}\in\mathcal{S}}H^{r,\gamma}(\mathbf{x}).\end{array}$
(21)

In order to prove the value approximation (5) let us separate it into parts
and estimate each of them:

$\displaystyle
F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-F^{\gamma}(\mathbf{x}^{*}_{\gamma})$
(22a) $\displaystyle\leq
F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})$
(22b)
$\displaystyle+H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})$
(22c)
$\displaystyle+H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})-F^{\gamma}(\mathbf{x}^{*}_{\gamma}).$
(22d)

The last addend is negative and can be eliminated:

$H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})-F^{\gamma}(\mathbf{x}^{*}_{\gamma})\leq
H^{r,\gamma}(\mathbf{x}^{*}_{\gamma})-F^{\gamma}(\mathbf{x}^{*}_{\gamma})\leq
0.$

The rest are estimated in Lemmas 4 and 5 under additional assumptions.

######  Lemma 4.

Let $\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{\gamma}\|_{2}\leq\zeta$. If
(19) holds, then

$\displaystyle
F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})\leq\frac{r}{2(1+r\gamma)}mK^{2}_{\zeta}.$
(23)

###### Proof.

We cannot declare a uniform $K$ instead of $K_{\zeta}$ because $F^{\gamma}$ is
not smooth. Nonetheless, assuming $\mathbf{x}^{n}_{r,\gamma}$ belong to
$\zeta$-neighborhood of $\mathbf{x}^{*}_{\gamma}$, we immediately obtain from
(18) and (19) that

$F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})\leq\frac{r}{2(1+r\gamma)}\|\nabla
F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})\|^{2}_{2}\leq\frac{r}{2(1+r\gamma)}mK^{2}_{\zeta}.$

∎

######  Lemma 5.

Let (19) holds. Then

$\displaystyle
H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})$
$\displaystyle\leq
C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2},$
$\displaystyle\mbox{where }\quad C_{2}$
$\displaystyle=\frac{m(1+r\gamma)K_{\zeta}}{\sqrt{2}\gamma}\sqrt{\frac{\lambda_{\max}}{\lambda_{\min}^{+}}}+\frac{m(1+r\gamma)^{2}}{4r\gamma^{2}}.$

###### Proof.

By $\frac{m}{r}$ smoothness of $H^{r,\gamma}$

$\begin{array}[]{rl}&H^{r,\gamma}(\mathbf{x}^{n}_{r,\gamma})-H^{r,\gamma}(\mathbf{x}^{*}_{r,\gamma})\leq\langle\nabla
H^{r,\gamma}\left(\mathbf{x}^{*}_{r,\gamma}\right),\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{r,\gamma}\rangle+\frac{m}{2r}\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{r,\gamma}\|^{2}_{2}\\\
&\leq\langle\nabla
H^{r,\gamma}\left(\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{\infty})\right),\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})-\mathbf{x}^{*}_{r,\gamma}\rangle+\frac{m}{2r}\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{r,\gamma}\|^{2}_{2}\\\
&\leq\langle\mathbf{z}_{g}^{\infty},\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})-\mathbf{x}^{*}_{r,\gamma}\rangle+\frac{m}{2r}\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{r,\gamma}\|^{2}_{2},\end{array}$

where $\mathbf{z}^{\infty}_{g}$ is the limit of $\mathbf{z}_{g}^{n}$ and so it
is the argminimum of $(H^{r,\gamma})^{*}$ on $\mathcal{R}^{\perp}$. By (17) we
have

$\frac{m}{2r}\|\mathbf{x}^{n}_{r,\gamma}-\mathbf{x}^{*}_{r,\gamma}\|^{2}_{2}\leq\frac{m}{2r}C_{1}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}=\frac{m(1+r\gamma)^{2}}{4r\gamma^{2}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}$

Let us introduce an orthogonal projection matrix $\mathbf{P}$ onto the
subspace $\mathcal{R}^{\perp}$, i.e., it holds
$\mathbf{P}v=\operatorname*{arg\,min}_{z\in\mathcal{R}^{\perp}}\\{v-z\\}$ for
an arbitrary $v\in(\mathbb{R}^{d})^{n}$. Then matrix $\mathbf{P}$ is

$\mathbf{P}=\left(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{\top}\right)\otimes\mathbf{I}_{d},$
(24)

where $\mathbf{I}_{n}$ denotes $n\times n$ identity matrix,
$\mathbf{1}_{n}=(1,\ldots,1)\in\mathbb{R}^{n}$, and $\otimes$ is a Kronecker
product. Note that $\mathbf{P}^{\top}\mathbf{P}=\mathbf{P}$.

Since $\mathbf{z}_{g}^{\infty}\in\mathcal{R}^{\perp}$ and
$\mathbf{x}^{*}_{r,\gamma}\in\mathcal{R}$, the first part simplifies to
$\langle\mathbf{z}_{g}^{\infty},\mathbf{P}\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})\rangle$.
We may use Lemma 2 in [19] to get the following estimation

$\begin{array}[]{rl}\|\mathbf{P}\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})\|_{2}^{2}=\|\nabla(H^{r,\gamma})^{*}(\mathbf{z}_{g}^{n})\|_{\mathbf{P}}^{2}\leq\frac{2}{\theta\lambda_{\min}^{+}}\left((H^{r,\gamma})^{*}(\mathbf{z}^{n}_{g})-(H^{r,\gamma})^{*}(\mathbf{z}^{n+1}_{f})\right).\end{array}$

As $\mathbf{z}^{n+1}_{f}$ is a non-optimal point of Algorithm 1, this is not
greater than

$\displaystyle\frac{2}{\theta\lambda_{\min}^{+}}\left((H^{r,\gamma})^{*}(\mathbf{z}^{n}_{g})-(H^{r,\gamma})^{*}(\mathbf{z}^{*})\right)$
$\displaystyle\leq\frac{m(1+r\gamma)}{\gamma\theta\lambda_{\min}^{+}}\left\|\mathbf{z}^{n}_{g}-\mathbf{z}^{*}\right\|_{2}^{2}=\frac{m(1+r\gamma)^{2}}{\gamma^{2}}\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\left\|\mathbf{z}^{n}_{g}-\mathbf{z}^{*}\right\|_{2}^{2}$
$\displaystyle\leq\frac{m(1+r\gamma)^{2}}{2\gamma^{2}}\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}$

and the latter ones follow from the $\frac{m(1+r\gamma)}{\gamma}$ smoothness
of $(H^{r,\gamma})^{*}$ and from the fact that the proof of [19, Theorem 1]
actually covers the following chain of inequalities:

$\left\|\nabla
H^{*}(\mathbf{z}_{g}^{n})-\mathbf{x}^{*}\right\|^{2}_{2}\leq\frac{1}{\mu^{2}}\left\|\mathbf{z}^{n}_{g}-\mathbf{z}^{*}\right\|^{2}_{2}\leq
C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu}{L}}\right)^{n}=\frac{1}{2\mu^{2}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{\mu}{L}}\right)^{n}.$

By our assumption $\|\mathbf{z}^{\infty}_{g}\|_{2}=\|\nabla
H^{r,\gamma}(x_{r,\gamma}^{*})\|_{2}<\sqrt{m}K_{\zeta}$. Thus, we obtain

$\begin{array}[]{rl}(H^{r,\gamma})^{*}(\mathbf{x}^{n}_{r,\gamma})-(H^{r,\gamma})^{*}(\mathbf{x}^{*}_{r,\gamma})\\\
\leq\sqrt{m}K_{\zeta}\frac{\sqrt{m}(1+r\gamma)}{\sqrt{2}\gamma}\sqrt{\frac{\lambda_{\max}}{\lambda_{\min}^{+}}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}+\frac{m(1+r\gamma)^{2}}{4r\gamma^{2}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}\\\
\leq\left(\frac{m(1+r\gamma)K_{\zeta}}{\sqrt{2}\gamma}\sqrt{\frac{\lambda_{\max}}{\lambda_{\min}^{+}}}+\frac{m(1+r\gamma)^{2}}{4r\gamma^{2}}\right)\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}.\\\
=C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}.\end{array}$

∎

###  B.5 Final compilation

This section completes the proof of Theorem 1 and shows Remark 1.

Recall that where $C_{1}=\frac{(1+r\gamma)^{2}}{2\gamma}$ and

$C_{2}=\frac{m}{2r}C_{1}+\frac{m(1+r\gamma)K_{\zeta}}{\sqrt{2}\gamma}\frac{\lambda_{\max}}{\lambda_{\min}^{+}}=\frac{m(1+r\gamma)^{2}}{4r\gamma}+\frac{m(1+r\gamma)K_{\zeta}}{\sqrt{2}\gamma}\frac{\lambda_{\max}}{\lambda_{\min}^{+}}.$

By Remark 4 and Lemmas 4, 5 we see that
$F^{\gamma}(\mathbf{x}^{n}_{r,\gamma})-F^{\gamma}(\mathbf{x}^{*}_{\gamma})<\varepsilon$
if

$\displaystyle\forall\mathbf{x}\in\\{\mathbf{y}\in\mid\|\mathbf{y}-\mathbf{x}_{\gamma}^{*}\|_{2}<\zeta\\}\qquad\|\nabla
F^{\gamma}(\mathbf{x})\|_{2}^{2}$ $\displaystyle<mK^{2}_{\zeta},$ (25)
$\displaystyle\sqrt{\frac{rm}{2\gamma}}K_{\zeta}+\sqrt{C_{1}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}$
$\displaystyle\leq\zeta,$ (26)
$\displaystyle\frac{r}{2(1+r\gamma)}mK^{2}_{\zeta}$
$\displaystyle\leq\varepsilon/2,$ (27) $\displaystyle
C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}$
$\displaystyle\leq\varepsilon/2.$ (28)

Let $\zeta=\sqrt{\varepsilon/\gamma}$ and let
$r\leq\frac{\varepsilon}{2mK^{2}_{\zeta}}$. Then (27) holds. If (28) fulfills,
then (26) follows from (27) and (28) as
$\sqrt{\frac{rm}{2\gamma}}K_{\zeta}\leq\sqrt{\frac{\varepsilon}{2\gamma}}\leq\zeta/\sqrt{2}$
and
$\sqrt{C_{1}}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}\leq\zeta/2$
since $1\leq\sqrt{C_{1}}\leq C_{1}\leq C_{2}$ and
$\varepsilon\leq\sqrt{\varepsilon/\gamma}=\zeta$. Thus, it suffices to assume

$\displaystyle\forall i\quad\forall x\in\\{y\in
S\mid\|y-x_{\gamma}^{*}\|_{2}^{2}\leq\varepsilon/\gamma\\}\qquad\|\nabla
f_{i}^{\gamma}(x)\|_{2}$ $\displaystyle\leq K,$ $\displaystyle r$
$\displaystyle\leq\frac{\varepsilon}{2mK^{2}},$ $\displaystyle
C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}$
$\displaystyle\leq\varepsilon/2.$

So $\varepsilon$ approximation requires a number of iteration

$\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\sqrt{\frac{1+r\gamma}{r\gamma}}\ln\frac{C_{2}}{\varepsilon}\right)=\mathcal{O}\left(\frac{\lambda_{\max}}{\lambda_{\min}^{+}}\frac{1}{\sqrt{\gamma\varepsilon}}\ln\frac{1}{\varepsilon}\right).$

##  Appendix C Proof of Theorem 2

To prove Theorem 2 we combine proved Theorem 1 with features of the entropy
regularization of the Wasserstein barycenter problem.

###  C.1 Entropy regularized WB problem

Recall that for a fixed cost matrix $M$ we define the set of transport plans
as

$\displaystyle U(p,q):=\left\\{X\in\mathbb{R}_{+}^{d\times d}\mid
X\mathbf{1}=p,X^{T}\mathbf{1}=q\right\\}$

and Wasserstein distance between two probability distributions $p$ and $q$ as

$\mathcal{W}(p,q):=\min_{X\in U(p,q)}\langle M,X\rangle.$

The entropy regularized (or smoothed) Wasserstein distance is defined as

$\displaystyle\mathcal{W}_{\gamma}(p,q):=\min_{X\in U(p,q)}\left\\{\langle
M,X\rangle-\gamma E(X)\right\\},$ (29)

where $\gamma>0$ and

$\displaystyle E(X):=-\sum_{i=1}^{d}\sum_{j=1}^{d}e(X_{ij}),$
$\displaystyle\text{where }e(x)=\left\\{\begin{array}[]{ll}x\ln x&\mbox{if
}x>0\\\ 0&\mbox{if }x=0.\end{array}\right.$ (32)

So it seeks to minimize the transportation costs while maximizing the entropy.
Moreover $\mathcal{W}_{\gamma}(p,q)\to\mathcal{W}(p,q)$ as $\gamma\to 0$.

Then the convex optimization problem (7) can be relaxed to the following
$\gamma$ strongly convex optimization problem

$\displaystyle\min_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p),$ (33)

where $\mathcal{W}_{\gamma,q_{i}}(p)=\mathcal{W}_{\gamma}(q_{i},p)$. The
argminimum of (33) is called the uniform Wasserstein barycenter [4, 5] of the
family of $q_{1},\ldots,q_{m}$. Moreover, problem (33) admits a unique
solution and approximates unregularized WB problem as follows.

######  Remark 5.

Let $\gamma\leq\frac{\varepsilon}{4}\ln d$. If vectors $\hat{p}_{i}\in
S_{1}(d)$ are such that

$\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(\hat{p}_{i})-\min\limits_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p)\leq\frac{\varepsilon}{2},$

then

$\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(\hat{p}_{i})-\min\limits_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(p)\leq\varepsilon.$

Indeed, as entropy is bounded we have
$\mathcal{W}_{q_{i}}(p)\leq\mathcal{W}_{\gamma,q_{i}}(p)\leq\mathcal{W}_{q_{i}}(p)+2\gamma\ln
d$ for all $i$ and $p$. Then, for $p^{*}=\operatorname*{arg\,min}\limits_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(p)$ and
$p^{*}_{\gamma}=\operatorname*{arg\,min}\limits_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p)$ it holds that

$\displaystyle\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(\hat{p}_{i})-\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}(p^{*})$
$\displaystyle\leq\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(\hat{p})-\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p^{*})+2\gamma\ln
d$
$\displaystyle\leq\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(\hat{p})-\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p^{*}_{\gamma})+\frac{\varepsilon}{2}\leq\varepsilon.$

###  C.2 Legendre transforms

One particular advantage of entropy regularization of the Wasserstein distance
is that it yields closed-form representations for the dual function
$\mathcal{W}^{*}_{\gamma,q}(\cdot)$ and for its gradient. Recall that the
Fenchel-Legendre transform of (29) is defined as

$\displaystyle\mathcal{W}^{*}_{\gamma,q}(z)$ $\displaystyle:=\max_{p\in
S_{1}(d)}\left\\{\left\langle
z,p\right\rangle-{\it\mathcal{W}_{\gamma,q}(p)}\right\\}.$ (34)

######  Theorem 4 ([33][Theorem 2.4]).

For $\gamma>0$, the Fenchel-Legendre dual function
$\mathcal{W}^{*}_{\gamma,q}(z)$ is differentiable

$\begin{array}[]{rl}\mathcal{W}^{*}_{\gamma,q}(z)&=\gamma\left(E(q)+\left\langle
q,\ln\mathcal{K}\alpha\right\rangle\right)\\\ &=-\gamma\left\langle q,\ln
q\right\rangle+\gamma\sum\limits_{j=1}^{m}[q]_{j}\ln\left(\sum\limits_{i=1}^{m}\exp\left(\frac{1}{\gamma}\left([z]_{i}-M_{ji}\right)\right)\right)\end{array}$
(35)

and its gradient $\nabla\mathcal{W}^{*}_{\gamma,q}(z)$ is $1/\gamma$-Lipschitz
in the 2-norm with

$\begin{array}[]{rl}\nabla\mathcal{W}^{*}_{\gamma,q}(z)&=\alpha\circ\left(\mathcal{K}\cdot{q}/({\mathcal{K}\alpha})\right)\in
S_{1}(d),\\\
\left[\nabla\mathcal{W}^{*}_{\gamma,q}(z)\right]_{l}&=\sum\limits_{j=1}^{m}[q]_{j}\frac{\exp\left(\frac{1}{\gamma}([z]_{l}-M_{lj})\right)}{\sum\limits_{i=1}^{m}\exp\left(\frac{1}{\gamma}([z]_{i}-M_{ij})\right)}.\end{array}$
(36)

where $z\in\mathbb{R}^{n}$ and for brevity we denote
$\alpha=\exp({z}/{\gamma})$ and $\mathcal{K}=\exp\left({-M}/{\gamma}\right)$.

Notice that to get back and obtain the approximated barycenter we can employ
the following result (with $\lambda_{i}=1$).

######  Theorem 5 ([33][Theorem 3.1]).

The barycenter $p^{*}$ solving (33) satisfies

$\forall i=1,\ldots,m\qquad
p^{*}=\nabla\mathcal{W}^{*}_{\gamma,q_{i}}(z^{*}_{i}),$

where the set of $z^{*}_{i}$ constitutes any solution of any smoothed dual WB
problem:

$\min\limits_{z_{1},\ldots,z_{m}\in\mathbb{R}^{d}}\sum\limits_{i=1}^{m}\lambda_{i}\mathcal{W}^{*}_{\gamma,q_{i}}(z_{i})\quad\mbox{s.t.}\quad\sum\limits_{i=1}^{m}\lambda_{i}z_{i}=0.$

Thus we can apply Theorem 1 for the problem (33) with explicitly defined
$\nabla\mathcal{W}^{*}_{\gamma,q_{i}}$ and obtain $\mathbf{x}^{n}_{r,\gamma}$
that satisfies

$\displaystyle\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}([\mathbf{x}^{n}_{r,\gamma}]_{i})-\min\limits_{p\in
S_{1}(d)}\sum\limits_{i=1}^{m}\mathcal{W}_{\gamma,q_{i}}(p)$
$\displaystyle\leq\frac{r}{4(1+r\gamma)}mK^{2}+\frac{1}{2}C_{2}\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n/2}\leq\varepsilon/2.$

By Remark 5 it proves

$\displaystyle\left|\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}([\mathbf{x}^{n}_{r,\gamma}]_{i})-\sum\limits_{i=1}^{m}\mathcal{W}_{q_{i}}([\mathbf{p}^{*}]_{i})\right|$
$\displaystyle\leq 2\gamma\ln
d+\frac{r}{4(1+r\gamma)}mK^{2}+C\left(1-\frac{\lambda_{\min}^{+}}{7\lambda_{\max}}\sqrt{\frac{r\gamma}{1+r\gamma}}\right)^{n}/2\leq\varepsilon,$

for
$C=\frac{1}{2}C_{2}=\frac{(1+r\gamma)mK_{\zeta}}{2\sqrt{2}\gamma}\sqrt{\frac{\lambda_{\max}}{\lambda_{\min}^{+}}}+\frac{(1+r\gamma)^{2}}{8r\gamma^{2}}$.

###  C.3 Parameter estimation

It remains to assign $\zeta>0$ and $K=K_{\zeta}$ satisfying (25). Due to
Assumption 2 such $\zeta$ and $K$ exist.

######  Proposition 1.

Let a set $\\{q_{i}\\}_{i=1}^{m}$ satisfies Assumption 2, let $p^{*}_{\gamma}$
be the uniform Wasserstein barycenter of $\\{q_{i}\\}_{i=1}^{m}$, and let
$\zeta\in\left(0,\min\\{\frac{1}{e},\ \min_{i,l}[q_{i}]_{l}\\}\right)$. For
each $i=1,\ldots,m$ the norm of the gradient
$\|\nabla\mathcal{W}_{\gamma,q_{i}}(\cdot)\|_{2}^{2}$ is uniformly bounded
over $\\{p\in S_{1}(d)\mid\|p-p^{*}_{\gamma}\|_{2}^{2}\leq\zeta\\};$ and the
bound $K_{\rho}$ is given in (37) for $\rho\leq\min\\{\frac{1}{e},\
\min_{i,l}[q_{i}]_{l}\\}-\zeta.$

We obtain Proposition 1 as a combination of Lemma 6 from [34] and proved below
Lemma 7.

######  Lemma 6 ([34, Lemma 3.5]).

For any $\rho\in(0,1)$, $q\in S_{1}(d)$, and $p\in\\{x\in
S_{1}(d)\mid\min_{l}x_{l}\geq\rho\\}$ there is a bound:
$\|\nabla\mathcal{W}_{\gamma,q}(p)\|^{2}_{2}\leq K_{\rho}$, where

$K_{\rho}=\sum\limits_{j=1}^{d}\left(2\gamma\ln
d+\inf_{i}\sup_{l}|M_{jl}-M_{il}|-\gamma\ln\rho\right)^{2}.$ (37)

######  Lemma 7.

Let a set $\\{q_{i}\\}_{i=1}^{m}$ satisfies Assumption 2, let $p^{*}_{\gamma}$
be the uniform Wasserstein barycenter of $\\{q_{i}\\}_{i=1}^{m}$. All
components $k$ of $p^{*}_{\gamma}$ have a uniform positive lower bound:
$[p^{*}_{\gamma}]_{k}\geq\min\\{\frac{1}{e},\ \min_{i,l}[q_{i}]_{l}\\}$.

###### Proof.

Let $X^{*}_{i}$ denote the optimal transport plan between $p^{*}_{\gamma}$ and
$q_{i}$. Assume the contrary: there is $k$ such that
$[p^{*}_{\gamma}]_{k}<\min\\{\frac{1}{e},\ \min_{i,l}[q_{i}]_{l}\\}$. Then
there is another component $n$ such that
$[p^{*}_{\gamma}]_{n}>\min_{i}[q_{i}]_{n}>\min_{i,l}[q_{i}]_{l}$. Consider the
vector $p$ that consists of $[p]_{i}=[p^{*}_{\gamma}]_{i}$ except for the
components $[p]_{n}=[p^{*}_{\gamma}]_{n}+\delta$ and
$[p]_{l}=[p^{*}_{\gamma}]_{l}-\delta$, where $\delta>0$ is less than
$\min_{i,a\not=b}[X_{i}^{*}]_{a,b}$ of the optimal transport plans $X^{*}_{i}$
between $p^{*}_{\gamma}$ and $q_{n}$. Because of the entropy, all these
optimal transport plans contain only positive non-diagonal elements, so such a
$\delta$ exists.

Construct now non-optimal transport plans between $p$ and each of $q_{i}$ in
order to get the contradiction with the assumption. Initially we have
$\mathcal{W}_{\gamma,q_{i}}(p^{*}_{\gamma})=\langle C,X^{*}_{i}\rangle-\gamma
X^{*}_{i}\ln X^{*}_{i}$. Consider the matrix $X_{i}$ that differs from
$X^{*}_{i}$ only at four elements:

$\displaystyle[X_{i}]_{kk}=[X_{i}^{*}]_{kk}+\frac{1}{2}\delta,\quad[X_{i}]_{kn}=[X_{i}^{*}]_{kn}+\frac{1}{2}\delta,$
$\displaystyle[X_{i}]_{nn}=[X_{i}^{*}]_{nn}+\frac{1}{2}\delta,\quad[X_{i}]_{nk}=[X_{i}^{*}]_{nk}+\frac{1}{2}\delta.$

Then $X_{i}$ is a transport plan between $p$ and $q_{i}$ since its elements
are positive and also $X_{i}\mathbf{1}=p$ and $X_{i}^{\top}\mathbf{1}=q_{i}$.
Using the monotonicity of entropy on the interval $(0,\frac{1}{e})$ and the
assumption that diagonal elements of the cost matrix $C$ are zero, we get for
each $i$:

$\begin{array}[]{rl}\mathcal{W}_{\gamma,q_{i}}(p)\leq&\langle
C,X_{i}\rangle-\gamma X_{i}\ln X_{i}\\\ =&\langle C,X^{*}_{i}\rangle-\gamma
X^{*}_{i}\ln X^{*}_{i}+\frac{1}{2}\delta C_{kn}-\frac{1}{2}\delta C_{nk}\\\
+&([X_{i}]_{kk}\ln[X_{i}]_{kk}-[X^{*}_{i}]_{kk}\ln[X^{*}_{i}]_{kk})\\\
+&([X_{i}]_{kn}\ln[X_{i}]_{kn}-[X^{*}_{i}]_{kn}\ln[X^{*}_{i}]_{kn})\\\
+&([X_{i}]_{nk}\ln[X_{i}]_{nk}-[X^{*}_{i}]_{nk}\ln[X^{*}_{i}]_{nk})\\\
+&([X_{i}]_{nn}\ln[X_{i}]_{nn}-[X^{*}_{i}]_{nn}\ln[X^{*}_{i}]_{nn})\\\
<&\langle C,X^{*}_{i}\rangle-\gamma X^{*}_{i}\ln X^{*}_{i}+\frac{1}{2}\delta
C_{kn}-\frac{1}{2}\delta C_{nk}\\\ =&\langle C,X^{*}_{i}\rangle-\gamma
X^{*}_{i}\ln X^{*}_{i}=\mathcal{W}_{\gamma,q_{i}}(p^{*}_{\gamma}).\end{array}$

The obtained inequalities
$\mathcal{W}_{\gamma,q_{i}}(p)<\mathcal{W}_{\gamma,q_{i}}(p^{*}_{\gamma})$
contradict to the fact that $p^{*}_{\gamma}$ is the barycenter; this proves
the lemma. ∎