#  Convergence of Batch Greenkhorn
for Regularized Multimarginal Optimal Transport

Vladimir Kostic Istituto Italiano di Tecnologia, Via Melen, 83, 16152 Genova,
Italy<EMAIL_ADDRESS>and Department of Mathematics and Informatics,
Faculty of Science, University of Novi Sad, Trg Dositeja Obradovića 4, 21000
Novi Sad, Serbia. Saverio Salzo Istituto Italiano di Tecnologia, Via Melen,
83, 16152 Genova, Italy<EMAIL_ADDRESS>and Department of Computer
Science, University College London. Massimilano Pontil Istituto Italiano di
Tecnologia, Via Melen, 83, 16152 Genova, Italy<EMAIL_ADDRESS>and
Department of Computer Science, University College London.

###### Abstract

In this work we propose a batch version of the Greenkhorn algorithm for
multimarginal regularized optimal transport problems. Our framework is general
enough to cover, as particular cases, some existing algorithms like Sinkhorn
and Greenkhorn algorithm for the bi-marginal setting, and (greedy)
MultiSinkhorn for multimarginal optimal transport. We provide a complete
convergence analysis, which is based on the properties of the iterative
Bregman projections (IBP) method with greedy control. Global linear rate of
convergence and explicit bound on the iteration complexity are obtained. When
specialized to above mentioned algorithms, our results give new insights
and/or improve existing ones.

##  1 Introduction

Over the recent years the field of optimal transport (OT) [24] has received
significant attention in machine learning and data science, as it provides
natural and powerful tools to compare probability distributions. In this paper
we study a general class of OT problems known as multimarginal optimal
transport (MOT), whereby several probability distributions are coupled
together in order to compute a measure of their association, see e.g. [20, 5].
MOT is receiving increasing interest due to its numerous applications, ranging
from density functional theory in quantum chemistry, to fluid dynamics, to
economics, to image processing, just to name a few, see [21] and references
therein. Particularly, in machine learning, MOT is important for generative
adversarial networks (GANs) [8], domain adaptation [15], Wasserstein
barycenters [1], clustering [19], and Bayesian inference of joint
distributions [14], among others.

We focus on discrete MOT, in which, given $m$ finitely supported probability
distributions, we wish to compute an optimal joint distribution which solves a
linear program, whose objective function involves a cost tensor and the
constraint set requires the joint distribution to have the given individual
ones as its marginals. It is well known that addressing directly the MOT
problem is computationally prohibitive. Furthermore, unlike the bi-marginal
case, MOT is NP-Hard for certain costs, even approximately [2]. To overcome
this issue, regularization techniques have been widely considered. The key
insight is to add a strongly convex regularizer to the MOT objective. A
popular choice in the bi-marginal setting is entropic regularization, which
leads to the well-known Sinkhorn divergence. However, the study of
computational regularized multimarginal optimal transport (RMOT) is less
developed, and the objective of this work is to devise efficient algorithms
with convergence guarantees.

### Related work.

Two popular frameworks for solving regularized OT problems are iterative
Bregman projections (IBP) and alternating dual minimization. Both have been
extensively studied in the bi-marginal case and two popular offsprings include
Sinkhorn and Greenkhorn algorithms. In contrast, the theory of general
multimarginal OT is less understood. The multimarginal version of Sinkhron
algorithm was first proposed by [4] and its convergence was established within
the framework of _cyclic IBP_ in finite dimension. More recently, the work
[12] extended the result to infinite dimension and in [7] even global linear
rate of convergence was obtained. On the other hand, the work [17] proposed a
version of _multimarginal Sinkhorn_ where, at each iteration, the marginal to
project on is chosen in a greedy way. In the same work acceleration using
Nesterov’s estimate sequences was also proposed and near-linear time
computational complexity bounds were obtained. Finally, alternating dual
minimization framework was used in [23] with Nesterov’s momentum acceleration
to build a competitive algorithm against the greedy multimarginal Sinkhorn.
Both papers [17, 23] derived sub-linear computational complexity bounds for
the considered RMOT algorithms, but no linear convergence rates were studied.

### Contribution.

After the introduction of RMOT as a Bregman projection problem in Sec. 2, in
Sec. 3 we propose a batch version of the popular Greenkhorn algorithm for RMOT
that greedily selects at each iteration a marginal and a batch of its
components of prescribed sizes. In Sec. 4, we establish its linear
convergence, providing also explicit rates depending on problem data in two
important cases: 1) batch Greenkhorn for bi-marginal optimal transport and 2)
(greedy) MultiSinkhorn of [17]. Moreover, in these two cases we also provide
iteration complexity bounds that improve state-of-the-art results. We stress
that the explicit rate we obtain for the greedy multimarginal Sinkhorn is
strictly better than the one recently derived in [7], for cyclic multimarginal
Sinkhorn. This shows the effectiveness of the greedy option for speeding up
the convergence of this type of algorithms.

### Notation.

Given $n\in\mathbb{N}$, let $[n]:=\\{1,2,\ldots,n\\}$, and let
$\Delta_{n}=\\{x\in\mathbb{R}_{+}^{n}\,|\,\lVert x\rVert_{1}=1\\}$ be the unit
simplex. Next, let $m\geq 2$ and $n_{1},\dots,n_{m}\in\mathbb{N}$. The space
of tensors (of order $m$) is denoted by
${\mathbb{X}}:=\mathbb{R}^{n_{1}\times\cdots\times n_{m}}$. For the sake of
brevity we set $\mathcal{J}=[n_{1}]\times\cdots\times[n_{m}]$ and we denote by
$j=(j_{1},\dots,j_{m})$ a general multi-index in $\mathcal{J}$. We set
${\mathbb{X}}_{+}=\\{\pi\in{\mathbb{X}}\,|\,\pi_{j}\geq 0,\text{ for every
}j\in\mathcal{J}\\}$ and
${\mathbb{X}}_{++}=\\{\pi\in{\mathbb{X}}\,|\,\pi_{j}>0,\text{ for every
}j\in\mathcal{J}\\}$. We need also to consider multi-indexes without the
$k$-th index, so we define
$\mathcal{J}_{-k}=[n_{1}]\times\dots\times[n_{k-1}]\times[n_{k+1}]\times\dots\times[n_{m}]$
and $j_{-k}$ a general multi-index in $\mathcal{J}_{-k}$. We will identify
$\mathcal{J}$ with $\mathcal{J}_{-k}\times[n_{k}]$ via the mapping $j\
\leftrightarrow\ (j_{-k},j_{k})$ with
$j_{-k}=(j_{1},\dots,j_{k-1},j_{k+1},\dots,j_{m})$. The space ${\mathbb{X}}$
is an Euclidean space when endowed with the standard scalar product and norm,
defined as

$\langle\pi,\pi^{\prime}\rangle=\sum_{j\in\mathcal{J}}\pi_{j}\quad\lVert\pi\rVert^{2}=\sum_{j\in\mathcal{J}}\pi_{j}^{2}\qquad\forall\,\pi,\pi^{\prime}\in{\mathbb{X}}.$
(1.1)

For every $\pi,\pi^{\prime}$ we denote by
$\pi\odot\pi^{\prime}\in{\mathbb{X}}$, the Hadamard product of $\pi$ and
$\pi^{\prime}$, that is, $(\pi\odot\pi^{\prime})_{j}=\pi_{j}\pi^{\prime}_{j}$.
Moreover, if
$\mathsf{v}_{1}=(v_{1,j_{1}})_{j_{1}\in[n_{1}]}\in\mathbb{R}^{n_{1}},\dots,\mathsf{v}_{m}=(v_{m,j_{m}}))_{j_{m}\in[n_{m}]}\in\mathbb{R}^{n_{m}}$,
we set $\oplus_{k=1}^{m}\mathsf{v}_{k}\in{\mathbb{X}}$ such that
$(\oplus_{k=1}^{m}\mathsf{v}_{k})_{j}=\sum_{k=1}^{m}v_{k,j_{k}}$, and
$\otimes_{k=1}^{m}\mathsf{v}_{k}\in{\mathbb{X}}$ such that
$(\otimes_{k=1}^{m}\mathsf{v}_{k})_{j}=\prod_{k=1}^{m}v_{k,j_{k}}$ For a
function $\phi\colon{\mathbb{X}}\to\left]-\infty,+\infty\right]$, the Fenchel
conjugate is denoted by $\phi^{*}$. Finally, as usual, the Dirac measure at
$x$ is denoted by $\delta_{x}$.

##  2 Entropic RMOT as the Bregman projection problem

In this section we formally introduce the discrete multimarginal optimal
transport problem and its regularized versions, emphasizing its connection
with the Bregman projection problem.

Let $\mathsf{a}_{k}\in\Delta_{n_{k}}$, $k\in[m]$ be prescribed histograms. MOT
consists in solving a linear program

$\min_{\pi\in\Pi(\mathsf{a}_{1},\dots,\mathsf{a}_{m})}\langle\mathsf{C},\pi\rangle,$
(2.1)

where $\mathsf{C}\in{\mathbb{X}}$ is a given cost tensor and
$\Pi(\mathsf{a}_{1},\dots,\mathsf{a}_{m})$, called _transport polytope_ , is
the covex set of nonnegative tensors in ${\mathbb{X}}$ whose marginals are
$\mathsf{a}_{1},\ldots,\mathsf{a}_{m}$. More specifically

$\Pi(\mathsf{a}_{1},\dots,\mathsf{a}_{m})=\big{\\{}\pi\in{\mathbb{X}}_{+}\,\big{|}\,\mathsf{R}_{k}(\pi)=\mathsf{a}_{k},\text{
for every }k\in[m]\big{\\}},$ (2.2)

where for all $k\in[m]$, we denote by
$\mathsf{R}_{k}\colon{\mathbb{X}}\to\mathbb{R}^{n_{k}}$ the _$k$ -th push-
forward projection operator_, defined as

$(\forall\,j_{k}\in[n_{k}])\quad\mathsf{R}_{k}(\pi)_{j_{k}}=\sum_{j_{-k}\in\mathcal{J}_{-k}}\pi_{(j_{-k},j_{k})},$
(2.3)

so that $\mathsf{R}_{k}(\pi)$ is the $k$-th marginal of $\pi$.

As noted in the introduction, problem (2.1) may be hard to solve. We thus
consider a regularized version using the (negative) Boltzmann-Shannon entropy
$\mathsf{H}(\pi):=\sum_{j}\pi_{j}(\log\pi_{j}-1)$, $\pi\in\text{\rm
dom}\,\mathsf{H}={\mathbb{X}}_{+}$. For a cost tensor
$\mathsf{C}\in{\mathbb{X}}_{+}$, the regularized multimarginal optimal
transport (RMOT) problem consists in computing

$\pi^{\star}=\operatorname*{\mathrm{arg\,min}}_{\pi\in\Pi(\mathsf{a}_{1},\dots,\mathsf{a}_{m})}\langle\mathsf{C},\pi\rangle+\eta\mathsf{H}(\pi),$
(2.4)

where $\eta>0$ is a regularization parameter. Now, recalling that the Kulback-
Leibler (KL) divergence
$\mathsf{KL}\colon{\mathbb{X}}\times{\mathbb{X}}\to[0,+\infty]$ is the
_Bregman distance_ associated to the (negative) entropy $\mathsf{H}$, i.e.,

$\displaystyle\mathsf{KL}(\pi,\pi^{\prime})$
$\displaystyle=\begin{cases}\mathsf{H}(\pi)-\mathsf{H}(\pi^{\prime})-\langle\pi-\pi^{\prime},\nabla\mathsf{H}(\pi^{\prime})\rangle&\text{if
}\pi^{\prime}\in\text{\rm int}(\text{\rm dom}\,\mathsf{H})\\\
+\infty&\text{otherwise}\end{cases}$
$\displaystyle=\begin{cases}\sum_{j\in\mathcal{J}}\pi_{j}\log(\pi_{j}/\pi^{\prime}_{j})-\pi_{j}+\pi^{\prime}_{j}&\text{if
}\pi^{\prime}\in{\mathbb{X}}_{++},\pi\in{\mathbb{X}}_{+}\\\
+\infty&\text{otherwise},\end{cases}$ (2.5)

and considering that $\text{\rm dom}\,\mathsf{H}={\mathbb{X}}_{+}$, it easy to
check that (2.4) can be re-written as

$\pi^{\star}=\operatorname*{\mathrm{arg\,min}}_{\mathsf{R}_{k}(\pi)=\mathsf{a}_{k},k\in[m]}\mathsf{KL}(\pi,\xi),$
(2.6)

where
$\xi=\nabla\mathsf{H}^{*}(-\mathsf{C}/\eta)=\exp(-\mathsf{C}/\eta)\in{\mathbb{X}}_{+}$
is the _Gibbs kernel_ tensor. Hence, the solution of (2.6) is nothing but the
_KL projection of $\xi$_ onto the affine set
$\\{\pi\in{\mathbb{X}}\,|\,(\forall\,k\in[m])\
\mathsf{R}_{k}(\pi)=\mathsf{a}_{k}\\}$. Indeed, in general, for an arbitrary
closed convex set ${\cal C}\subset{\mathbb{X}}$ such that
$\mathsf{C}\cap{\mathbb{X}}_{++}\neq\varnothing$ and a point
$\pi\in{\mathbb{X}}_{++}$, the KL projection of $\pi$ onto ${\cal C}$ is
defined as

${\cal P}_{\cal C}(\pi):=\operatorname*{\mathrm{arg\,min}}_{\gamma\in{\cal
C}}\mathsf{KL}(\gamma,\pi),$ (2.7)

while its KL distance to ${\cal C}$ is defined as $\mathsf{KL}_{{\cal
C}}(\pi):=\mathsf{KL}({\cal P}_{{\cal C}}(\pi),\pi)$.

##  3 Greedy KL projections for entropic RMOT

Now we focus on the entropic-regularized MOT problem (2.4), in the equivalent
form (2.6). Recalling definition (2.3), we introduce the linear operator

$\mathsf{R}\colon{\mathbb{X}}\to\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}},\qquad\mathsf{R}(\pi)=(\mathsf{R}_{1}(\pi),\dots,\mathsf{R}_{m}(\pi)).$
(3.1)

Then, the constraints in (2.6) define the affine set

$\Pi=\\{\pi\in{\mathbb{X}}\;|\;\mathsf{R}(\pi)=(\mathsf{a}_{1},\dots,\mathsf{a}_{m})\\}.$
(3.2)

We recall that the positiveness constrain embodied in problem (2.4) can be
absorbed in the entropy function $\mathsf{H}$, since $\text{\rm
dom}\,\mathsf{H}={\mathbb{X}}_{+}$, so that problem (2.4) is equivalent to the
computation of the KL projection of the Gibbs kernel $\xi$ onto the affine set
$\Pi$. However, it is possible to prove (see equation (A.11) in Appendix A)
that ${\cal P}_{\Pi}(\xi)={\cal
P}_{\Pi}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})$. Therefore, we will
equivalently target the computation of the KL projection of the normalized
Gibbs kernel $\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k}$. In order to compute
such projection, in the following, we will represent the set $\Pi$ as an
intersection of simpler affine sets and then we will use the alternating KL
projection algorithm. More precisely, we will rewrite the set (3.2) as an
intersection of affine sets, obtained via specific _sketching_ , on which KL
projections have closed forms.

We consider sketches that cover several popular algorithms (multimarginal
Sinkhorn of [4], greedy MultiSinkhorn of [17], bi-marginal Greenkhorn of [3]).
To that purpose, for each $k\in[m]$ (which refers to the $k$-th marginal) and
each batch $L\subset[n_{k}]$ we consider the _canonical injection_

$\mathsf{S}_{(k,L)}\colon\mathbb{R}^{L}\to\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}$

of $\mathbb{R}^{L}$ into
$\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}$, meaning that for
each $\mathsf{u}=(u_{j_{k}})_{j_{k}\in L}\in\mathbb{R}^{L}$,
$\mathsf{S}_{(k,L)}\mathsf{u}$ is the vector of
$\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}$ obtained from the
completion of $\mathsf{u}$ with zero entries. Then, since the adjoint operator
$\mathsf{S}^{*}_{(k,L)}\colon\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}\to\mathbb{R}^{L}$
is the standard projection, we can define

$\mathsf{R}_{(k,L)}:=\mathsf{S}^{*}_{(k,L)}\mathsf{R}\colon{\mathbb{X}}\to\mathbb{R}^{L},$
(3.3)

and the set

$\displaystyle\Pi_{(k,L)}:$
$\displaystyle=\\{\pi\in{\mathbb{X}}\,|\,\mathsf{S}_{(k,L)}^{*}\mathsf{R}(\pi)=\mathsf{S}_{(k,L)}^{*}(\mathsf{a}_{1},\dots,\mathsf{a}_{m})\\}$
(3.4)
$\displaystyle=\\{\pi\in{\mathbb{X}}\,|\,\mathsf{R}_{(k,L)}(\pi)={\mathsf{a}_{k}}_{\lvert
L}\\}$ $\displaystyle=\\{\pi\in{\mathbb{X}}\,|\,(\mathsf{R}_{k}(\pi))_{\lvert
L}={\mathsf{a}_{k}}_{\lvert L}\\}.$

Note that in the definition of $\Pi_{(k,L)}$ we ask for the $k$-th marginal of
$\pi$ to be equal to $\mathsf{a}_{k}$ only on the components in $L$. Now,
given $\tau=(\tau_{k})_{1\leq k\leq m}$ a vector of batch sizes, we set

${\cal I}(\tau)=\\{(k,L)\,|\,k\in[m],L\subset[n_{k}]\;|\;|L|\leq\tau_{k}\\}$

and obtain $\Pi=\bigcap_{(k,L)\in{\cal I}(\tau)}\Pi_{(k,L)}$.

Hence, the _iterative Bregman projections_ (IBP) algorithm, for problem (2.6),
leads to the following procedure. Let $\pi^{0}=$
$\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k}=e^{-\mathsf{C}/\eta}\odot\otimes_{k=1}^{m}\mathsf{a}_{k}\in\text{\rm
int}(\text{\rm dom}\,\mathsf{H})={\mathbb{X}}_{++}$ and define the sequence
$\pi^{t}$ recursively as follows

$\begin{array}[]{l}\text{for}\;t=0,1,\ldots\\\
\left\lfloor\begin{array}[]{l}\text{choose }(k_{t},L_{t})\in{\cal I}(\tau),\\\
\pi^{t+1}={\cal
P}_{\Pi_{(k_{t},L_{t})}}(\pi^{t}).\end{array}\right.\end{array}$ (3.5)

Since the generalized Pythagora’s theorem for Bregman projections (see
equation (A.10) in Appendix A) yields that
$\mathsf{KL}_{\Pi}(\pi^{t})=\mathsf{KL}_{\Pi_{(k,L)}}(\pi^{t})+\mathsf{KL}_{\Pi}({\cal
P}_{\Pi_{(k,L)}}(\pi^{t}))$, in (3.5) one may choose the sets in a _greedy_
manner as

$(k_{t},L_{t})=\operatorname*{\mathrm{arg\,max}}_{(k,L)\in{\cal
I}(\tau)}\mathsf{KL}_{\Pi_{(k,L)}}(\pi^{t}),$ (3.6)

so that

$(k_{t},L_{t})=\operatorname*{\mathrm{arg\,min}}_{(k,L)\in{\cal
I}(\tau)}\mathsf{KL}_{\Pi}({\cal P}_{\Pi_{(k,L)}}(\pi^{t}))\,\text{ and
}\,\mathsf{KL}_{\Pi}(\pi^{t+1})=\min_{(k,L)\in{\cal
I}(\tau)}\mathsf{KL}_{\Pi}({\cal P}_{\Pi_{(k,L)}}(\pi^{t})).$ (3.7)

This means that the next iterate is chosen, among the possible projections, as
the one which is the closest to the target set $\Pi$. Notable examples of
existing algorithms that fit in this framework are greedy multimarginal
Sinkhorn of [17] ($\tau_{k}=n_{k}$) and bi-marginal Greenkhorn of [3] ($m=2$,
$\tau_{1}=\tau_{2}=1$). We emphasize that this greedy strategy typically leads
to the best performance, provided that it can be implemented efficiently. In
the following proposition and subsequent remak we show that the projection
onto the sets $\Pi_{(k,L)}$ can be computed in a closed form and that the
greedy choice of the sets $\Pi_{(k,L)}$’s can indeed be implemented
efficiently. The proof is postponed in Appendix B.

######  Proposition 3.1.

For every $\pi\in{\mathbb{X}}_{+}$, $k\in[m]$ and $L\subset[n_{k}]$,

${\cal
P}_{\Pi_{(k,L)}}(\pi)=\nabla\mathsf{H}^{*}\big{(}\nabla\mathsf{H}(\pi)+\mathsf{R}_{(k,L)}^{*}(\bar{\mathsf{u}})\big{)}=\pi\odot\exp\big{(}\mathsf{R}^{*}_{(k,L)}(\bar{\mathsf{u}})\big{)},$
(3.8)

where

$\bar{\mathsf{u}}=\operatorname*{\mathrm{arg\,min}}_{\mathsf{u}\in\mathbb{R}^{L}}\mathsf{H}^{*}\big{(}\nabla\mathsf{H}(\pi)+\mathsf{R}^{*}_{(k,L)}(\mathsf{u})\big{)}-\mathsf{H}^{*}(\nabla\mathsf{H}(\pi))-\langle{\mathsf{a}_{k}}_{\lvert
L},\mathsf{u}\rangle=\log\frac{{\mathsf{a}_{k}}_{\lvert
L}}{{\mathsf{R}_{k}(\pi)}_{\lvert L}},$ (3.9)

and, consequently, for every $j\in\mathcal{J}$,

$({\cal
P}_{\Pi_{(k,L)}}(\pi))_{j}=\pi_{j}\times\begin{cases}\dfrac{a_{k,j_{k}}}{{\mathsf{R}_{k}(\pi)}_{j_{k}}}&\text{if
}j_{k}\in L,\\\\[8.61108pt] 1&\text{otherwise}.\end{cases}$ (3.10)

Moreover,

$\mathsf{KL}_{\Pi_{(k,L)}}(\pi)=\mathsf{KL}\big{(}{\cal
P}_{\Pi_{(k,L)}}(\pi),\pi)=\mathsf{KL}({\mathsf{a}_{k}}_{\lvert
L},\mathsf{R}_{k}(\pi)_{\lvert L}).$ (3.11)

######  Remark 3.2.

It follows from (3.11) that the greedy choice described above can be
implemented by computing $m$ vectors of sizes $n_{k}$ and then choosing
$k_{t}$ among $m$ as the index of the vector that has the maximal sum of the
largest $\tau_{k}$ components. More formally one let
$\mathsf{d}_{k}=(\mathsf{KL}(a_{k,1},(\mathsf{R}_{k}(\pi^{t}))_{1}),\dots,\mathsf{KL}(a_{k,n_{k}},(\mathsf{R}_{k}(\pi^{t}))_{n_{k}}))\in\mathbb{R}^{n_{k}}$
and consider the vector $\mathsf{d}_{k\downarrow}\in\mathbb{R}^{n_{k}}$ which
has the components of $\mathsf{d}_{k}$ arranged in a decreasing order. Then
$k_{t}=\operatorname*{\mathrm{arg\,max}}_{k\in[m]}\big{(}\sum_{j_{k}=1}^{\tau_{k}}(\mathsf{d}_{k\downarrow})_{j_{k}}\big{)}$
and $L_{t}$ corresponds to the indexes of the largest $\tau_{k_{t}}$
components of $\mathsf{d}_{k_{t}}$.

To ensure better numerical stability, especially for small regularization
parameters $\eta$, it is more convenient to work with the dual variables
$\nabla\mathsf{H}(\pi)=\log(\pi)$. More precisely, according to (3.8) we have
that each iteration of the IBP algorithm (3.5) can be parameterized as

$\pi^{t}=\exp\Big{(}-\mathsf{C}/\eta+\bigoplus_{k=1}^{m}\mathsf{v}_{k}^{t}\Big{)}\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k},\;t\in\mathbb{N},$
(3.12)

and one can implement the algorithm by updating only the dual variables
$\mathsf{v}_{k}^{t}=(v^{t}_{k,j})_{1\leq j\leq n_{k}}\in\mathbb{R}^{n_{k}}$,
$k\in[m]$, also known as _potentials_ (see Proposition B.1). Thus, in the end,
the IBP algorithm (3.5)-(3.6) can be written as the following Algorithm 1.

Input: $(\mathsf{a}_{1},\dots,\mathsf{a}_{m})$,
$\mathsf{C}\in{\mathbb{X}}_{+}$, $\eta>0$, $\tau=(\tau_{1},\dots,\tau_{m})$,
$1\leq\tau_{k}\leq n_{k}$

1 Initialization: $t=0$, $\mathsf{v}^{t}_{k}=0$, $k\in[m]$,
$\mathsf{r}_{k}^{t}=\mathsf{R}_{k}(\exp(-\mathsf{C}/\eta+\oplus_{k=1}^{m}\mathsf{v}_{k}^{t})\odot\otimes_{k=1}^{m}\mathsf{a}_{k})$

2 for  _ $t=0,1,\ldots$ _ do

3 Compute $(k_{t},L_{t})=\operatorname*{\mathrm{arg\,max}}_{(k,L)\in{\cal
I}(\tau)}\mathsf{KL}({\mathsf{a}_{k}}_{\lvert L},{\mathsf{r}^{t}_{k}}_{\lvert
L})$

4 Set $\mathsf{v}_{k}^{t+1}=\mathsf{v}_{k}^{t}$, $k\in[m]$ and update
${\mathsf{v}_{k_{t}}^{t+1}}_{\lvert
L_{t}}\leftarrow{\mathsf{v}^{t+1}_{k_{t}}}_{\\!\\!\\!\lvert
L_{t}}+\log({\mathsf{a}_{k_{t}}}_{\lvert
L_{t}})-\log({\mathsf{r}_{k_{t}}^{t}}_{\lvert L_{t}})$

5 For $k\in[m]$ compute
$\mathsf{r}_{k}^{t+1}=\mathsf{R}_{k}(\exp(-\mathsf{C}/\eta+\oplus_{k=1}^{m}\mathsf{v}_{k}^{t+1})\odot\oplus_{k=1}^{m}\mathsf{a}_{k})$

6

7 end for

Output:
$\pi^{t}=\exp(-\mathsf{C}/\eta+\oplus_{k=1}^{m}\mathsf{v}_{k}^{t})\odot\otimes_{k=1}^{m}\mathsf{a}_{k}$

Algorithm 1
$\mathtt{BatchGreenkhorn}(\mathsf{a}_{1},\dots,\mathsf{a}_{m},\mathsf{C},\eta,\tau)$

##  4 Convergence theory for Batch Greenkhorn algorithm

Results on the convergence of general IBP are typically without any rates [6,
9, 10], with the notable exception of [16] where the explicit local linear
rate for the greedy IBP was derived. In the following we prove global linear
convergence and derive the explicit dependence of the rate on the given data
in two important cases. Moreover, we also provide an analysis of the iteration
complexity of Algorithm 1.

Based on the properties of operators $\mathsf{R}$ and $\mathsf{R}_{(k,L)}$ we
can derive the main results using the properties of KL as Bregman divergence.
The proofs are given in Appendix C.

######  Theorem 4.1.

Algorithm 1 converges linearly. More precisely, if
$(\mathsf{v}_{k}^{t})_{k\in[m]}$ are generated by Algorithm 1, then the primal
iterates given by (3.12) converge linearly in KL divergence to $\pi^{\star}$
given by (2.6), i.e.,

$(\forall\,t\in\mathbb{N})\quad\mathsf{KL}(\pi^{\star},\pi^{t})\leq\bigg{(}1-\frac{e^{-(2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1})}}{b_{\tau}-1}\bigg{)}^{t}\mathsf{KL}(\pi^{\star},\pi^{0}),$
(4.1)

where $b_{\tau}=\sum_{k\in[m]}\lceil n_{k}/\tau_{k}\rceil$, and
$0<M_{1}<+\infty$ is a constant independent of the batch sizes that satisfies
$\lVert\bigoplus_{k=1}^{m}\mathsf{v}^{\star}_{k}\rVert_{\infty},\lVert\bigoplus_{k=1}^{m}\mathsf{v}^{t}_{k}\rVert_{\infty}\leq
M_{1}$ for $t\in\mathbb{N}$.

######  Theorem 4.2.

Let $\varepsilon>0$ and suppose that $\eta>\varepsilon$. For Algorithm 1, the
number of iterations required to reach the stopping criterion
$d_{t}:=\max_{k\in[m]}\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}\leq\varepsilon$
satisfies

$t\leq
2+\max_{k\in[m]}\bigg{\lceil}\frac{n_{k}}{\tau_{k}}\bigg{\rceil}\frac{5M_{2}}{\varepsilon}(2+M_{2}\eta),$
(4.2)

where $0<M_{2}<+\infty$ is a constant independent of the batch sizes such that
$\sum_{k\in[m]}\lVert\mathsf{v}_{k}^{\star}-\mathsf{v}_{k}^{t}\rVert\leq
M_{2}$, for all $t\in\mathbb{N}$.

######  Remark 4.3.

We stress that the constants $M_{1}$ and $M_{2}$ considered in the above
theorems always exist (see the proofs in the Appendix C), but we could not
have obtained general explicit expression for them depending on the problem
data, which are valid for any $m$ and $(\tau_{k})_{k\in[m]}$. On the other
hand, in the following, we show the important cases $m=2$ and
$(\tau_{k})_{k\in[m]}$ arbitrary and $m>2$ and
$(\tau_{k})_{k\in[m]}=(n_{k})_{k\in[m]}$ for which we do provide explicit
dependence on the problem data.

A natural issue for the BatchGreenkhorn algorithm is that of optimizing the
batch size. As extreme cases we have full batch ($\tau_{k}=n_{k}$), which
yields the (greedy) MultiSinkhorn algorithm proposed in [17], and
$\tau_{k}=1$, which, in the bi-marginal case, is known as the Greenkhorn
algorithm [3]. In this respect we observe that, since step 3 in Algorithm 1
can be efficiently implemented, as discussed in Remark 3.2, the largest
computational cost lies in the computation of the marginals
$\mathsf{r}_{k}^{t+1}$ in step 5. For simplicity, let us assume that $n_{k}=n$
and $\tau_{k}=\tau$, for every $k\in[m]$. Then, computing it naively yields
${\cal O}(mn^{m})$ operations, but indeed it can be done more efficiently in
${\cal O}(\tau n^{m-1})$ as we show in Appendix B. This way $n/\tau$
iterations with batch size $\tau$ have the same computational cost of one
iteration with a full batch $n$ and consequently $b_{\tau}=mn/\tau$ iterations
of BatchGreenkhorn with batch sizes $\tau$ corresponds to one cycle of cyclic
multimarginal Sinkhorn. Hence, the number of normalized cycles $T=t/b_{\tau}$
required to satisfy the stopping criterion given in Theorem 4.2 is

$T\leq 1+\frac{5M_{2}}{m\varepsilon}(2+M_{2}),$ (4.3)

which we stress is independent on the batch-sizes and the dimension $n$.
Whereas, the rate (4.1) w.r.t the normalized cycles becomes

$\mathsf{KL}(\pi^{\star},\pi^{b_{\tau}T})\leq\bigg{[}\bigg{(}1-\frac{e^{-(2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1})}}{b_{\tau}-1}\bigg{)}^{b_{\tau}}\bigg{]}^{T}\mathsf{KL}(\pi^{\star},\pi^{0}).$
(4.4)

We now introduce the analysis of the special cases discussed in Remark 4.3.
The proofs are based on bounding the potentials and are detailed in Appendix
C.

######  Theorem 4.4.

If $m=2$, then algorithm
$\mathtt{BatchGreenkhorn}(\mathsf{a}_{1},\dots,\mathsf{a}_{m},\mathsf{C},\eta,\tau)$
converges linearly with the global rate

$(\forall\,t\in\mathbb{N})\quad\mathsf{KL}(\pi^{\star},\pi^{t})\leq\left(1-\frac{e^{-20\lVert\mathsf{C}\rVert_{\infty}/\eta}}{b_{\tau}-1}\right)^{t}\mathsf{KL}(\pi^{\star},\pi^{0}).$
(4.5)

Moreover, when $\eta>\varepsilon$, the number of iterations required to reach
the stopping criterion $d_{t}\leq\varepsilon$ satisfies

$t\leq
2+\max_{k\in[m]}\bigg{\lceil}\frac{n_{k}}{\tau_{k}}\bigg{\rceil}\frac{15\lVert\mathsf{C}\rVert_{\infty}(2+3\lVert\mathsf{C}\rVert_{\infty})}{\eta\varepsilon}.$
(4.6)

The following theorem provides new insights into a known algorithm.

######  Theorem 4.5.

If for all $k\in[m]$ $\tau_{k}=n_{k}$, then algorithm
$\mathtt{BatchGreenkhorn}(\mathsf{a}_{1},\dots,\mathsf{a}_{m},\mathsf{C},\eta,\tau)$,
i.e. MultiSinkhorn algorithm of [17], converges linearly with the global rate

$(\forall\,t\in\mathbb{N})\quad\mathsf{KL}(\pi^{\star},\pi^{t})\leq\left(1-\frac{e^{-(12m-7)\lVert\mathsf{C}\rVert_{\infty}/\eta}}{m-1}\right)^{t}\mathsf{KL}(\pi^{\star},\pi^{0}).$
(4.7)

Moreover, the number of iterations required to reach stopping criterion
$d_{t}\leq\varepsilon$ satisfies

$t\leq 1+\frac{8(4m-3)\lVert\mathsf{C}\rVert_{\infty}}{\eta\varepsilon}.$
(4.8)

######  Remark 4.6.

Concerning the linear convergence rate, we notice that when the batch is full
($\tau_{k}=n_{k}$, $k\in[m]$), cyclic Sinkhorn of [4, 7] and (greedy)
MultiSinkhorn algorithm of [17] generally differ (unless $m=2$). Moreover, our
results shows (for the first time) that the rate of convergence of (greedy)
MultiSinkhorn algorithm is strictly better than that of the cyclic Sinkhorn
algorithm obtained in [7]. Indeed in [7], the following rate was provided

$(\forall\,T\in\mathbb{N})\quad\mathsf{KL}(\pi^{\star},\pi^{mT})\leq\left(1-\frac{e^{-8(2m-1)\lVert\mathsf{C}\rVert_{\infty}/\eta}}{m}\right)^{T}\mathsf{KL}(\pi^{\star},\pi^{0}),$
(4.9)

where $T$ counts the number of cycles, each one consisting of $m$
$\mathsf{KL}$ projections on the given marginals. Whereas, it follows from our
rate (4.7) that for the (greedy) MultiSinkhorn algorithm we have

$(\forall\,T\in\mathbb{N})\quad\mathsf{KL}(\pi^{\star},\pi^{mT})\leq\bigg{[}\bigg{(}1-\frac{e^{-(12m-7)\lVert\mathsf{C}\rVert_{\infty}/\eta}}{m-1}\bigg{)}^{m}\bigg{]}^{T}\mathsf{KL}(\pi^{\star},\pi^{0}),$
(4.10)

which mainly gains an $m$-root in the rate. Note that in the bi-marginal case,
the rate of the classical Sinkhorn algorithm according to [7] is
$1-(1/2)e^{-24\lVert\mathsf{C}\rVert_{\infty}/\eta}$, while we obtain
$(1-e^{-17\lVert\mathsf{C}\rVert_{\infty}/\eta})^{2}$.

######  Remark 4.7.

Concerning iteration complexity, we note that in literature the stopping
criteria concerns
$\sum_{k\in[m]}\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}$ and the
assumptions usually demand $n_{k}=n$ and $\tau_{k}=\tau$, for every $k\in[m]$.
In this setting, Theorems 4.4 and 4.5 provide the following bounds (in terms
of normalized cycles)

$T\leq
1+\frac{15\lVert\mathsf{C}\rVert_{\infty}(2+3\lVert\mathsf{C}\rVert_{\infty})}{\eta\varepsilon}\quad\text{and}\quad
T\leq 1+\frac{8(4m-3)\lVert\mathsf{C}\rVert_{\infty}}{\eta\varepsilon}.$
(4.11)

respectively. Those bounds improve existing results from [17], related to
(greedy) MultiSinkhorn, and [18], for bi-marginal Sinkhorn and Greenkhorn.
Indeed, defining $a_{\min}=\min_{k\in[m],j\in[n]}a_{k,j}$, in [13, 17], the
bi-marginal Sinkhorn and Greenkhorn algorithm are shown to feature the
following iteration complexity

$T\leq
1+\frac{2(\lVert\mathsf{C}\rVert_{\infty}/\eta+\log(a^{-1}_{\min}))}{\varepsilon}\quad\text{and}\quad
T\leq
1+\frac{56(\lVert\mathsf{C}\rVert_{\infty}/\eta+\log(n)+2\log(a^{-1}_{\min}))}{\varepsilon}$

whereas, in [18], for the MultiSinkhorn algorithm, the following bound is
provided

$T\leq
1+\frac{2m(\lVert\mathsf{C}\rVert_{\infty}/\eta+\log(1/a_{\min}))}{\varepsilon}.$
(4.12)

All those bounds contain the term $\log(a^{-1}_{\min})$ which at the best is
$\log(n)$. So they all depend on the dimension of the problem. Our results
remove this dependency.

######  Remark 4.8.

Here we look at the computational complexity in terms of arithmetic
operations. To this purpose we let ${\cal O}_{\xi}$ be the number of
arithmetic operations needed to compute a full marginal using the Gibbs kernel
$\xi$. When we factor this number out of the total computational complexity,
what remains is the number of iterations normalized with respect to the full
batch, meaning, $\bar{t}=\tau t/n$. Then according to equation (4.3) the total
computational complexity is given by

$\bigg{[}m+\frac{5M_{2}}{\varepsilon}(2+M_{2})\bigg{]}{\cal O}_{\xi}.$ (4.13)

We note that in the worst case ${\cal O}_{\xi}$ is of the order of $n^{m}$,
but this cost can be significantly reduced for structured costs [21].

##  5 Conclusions, limitations and future work

We present a new algorithm for solving multimarginal entropic regularized
optimal transport problems, called _batch Greenkhorn_ , which is an extension
of the popular Greenkhorn algorithm, that can handle multiple marginals and at
each iteration select in a greedy fashion a batch of components of a marginal.
We study the convergence of the algorithm in the framework of the iterative
Bregman projections method, providing novel linear rate of convergence as well
as iteration complexity bounds. We made a comprehensive comparison with
existing results showing the improvements over the state-of-the-art. A problem
which remains open is that of deriving bounds on the dual variables with an
explicit dependence on the given problem data for $m\geq 3$ when the batch is
not full. According to our general Theorems 4.1 and 4.2, this will allow to
have explicit linear rate and iteration complexity in all possible cases.
Additional research directions are the extension of our analysis to infinite
dimensions and general convex regularizers, implementing batch Greenkhorn with
structured costs, and analyze the impact of parallel computations.

## References

  * [1] M. Agueh and G. Carlier,  Barycenters in the Wasserstein space,  SIAM Journal on Mathematical Analysis, 43(2) (2011) pp. 904–924.
  * [2] J. Altschuler and E. Boix-Adsera,  Hardness results for multimarginal optimal transport problems,  in arXiv:2012.05398 (2020).
  * [3] J. Altschuler, J. Weedm and P. Rigollet,  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration,  in Proceedings of the 31st International Conference on Neural Information Processing Systems (2017) pp. 1961–1971.
  * [4] J-D. Benamou, G. Carlier, M. Cuturi, L. Nenna and G. Peyré,  Iterative Bregman projections for regularized transportation problems,  SIAM J. Sci. Comput., 37 (2015) pp. A1111–A1138.
  * [5] J-D. Benamou, G. Carlier, M. Cuturi and L. Nenna,  A Numerical Method to solve multimarginal optimal transport problems with Coulomb cost,  In: Glowinski R., Osher S., Yin W. (eds) Splitting Methods in Communication, Imaging, Science, and Engineering. Scientific Computation. Springer, (2017).
  * [6] L.M. Bregman,  The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,  USSR Comput. Math. & Math. Phys., 7 (1967) pp. 200–217.
  * [7] G. Carlier,  On the linear convergence of the multimarginal Sinkhorn algorithm,  in hal-03176512 (2021).
  * [8] J. Cao, L. Mo, Y. Zhang, K. Jia, C. Shen and M. Tan,  Multimarginal Wasserstein GAN,  in 33rd Conference on Neural Information Processing Systems, (2019) pp. 1776–1786.
  * [9] Y. Censor and A. Lent,  An iterative row-action method for interval convex programming,  J. Optim. Theory Appl., 34, (1981) pp. 321–353.
  * [10] Y. Censor and S. Reich,  Iteration of paracontractions and firmly nonexpansive operators with applications to feasibility optimization,  Optimization, 37, (1996) pp. 323–339.
  * [11] A. Dessein, N Papadakis and J.-L. Rouas,  Regularized optimal transport and the ROT mover’s distance,  Journal of Machine Learning Research 19, (2018) pp. 1–53.
  * [12] S. Di Marino and A. Gerolin.  An optimal transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm,  J. Sci. Comput., 85(2) (2020).
  * [13] P. Dvurechensky, A. Gasnikov, A. Kroshnin.  Computational Optimal Transport: Complexity by accelerated gradient descent is better than by Sinkhorn’s algorithm,  Proceedings of the 35th International Conference on Machine Learning, PMLR 80:1367-1376, 2018.
  * [14] C. Frogner, T. Poggio.  Fast and Flexible Inference of Joint Distributions from their Marginals  Proceedings of the 36th International Conference on Machine Learning, PMLR 97:2002-2011, 2019.
  * [15] Z. He, W. Zuo, M. Kan, S. Shan and X. Chen,  Attgan: Facial attribute editing by only changing what you want,  IEEE Transactions on Image Processing, 28(11) (2019) pp. 5464–5478.
  * [16] V. Kostić and S. Salzo,  The method of Bregman projections in deterministic and stochastic convex feasibility problems,  arXiv:2101.01704, (2021).
  * [17] T. Lin, N. Ho, M. Cuturi and M. I. Jordan,  On the complexity of approximating multimarginal optimal transport,  arXiv:1910.00152, (2020)
  * [18] T. Lin, N. Ho and M. I. Jordan,  On the Efficiency of Sinkhorn and Greenkhorn and Their Acceleration for Optimal Transport,  arXiv:1906.01437, (2021)
  * [19] L. Mi, J. Bento.  Multi-marginal optimal transport defines a generalized metric.  arXiv:2001.11114v4, 2021.
  * [20] B. Pass,  Multimarginal optimal transport: theory and applications,  ESAIM: Mathematical Modelling and Numerical Analysis, 49(6) (2015) pp. 1771–1790.
  * [21] G. Peyré and M. Cuturi,  Computational Optimal Transport: With Applications to Data Sciences,  Foundations and Trends in Machine Learning, 11(5-6) (2019) pp. 355–607.
  * [22] T. Sun and Q. Tran-Dinh.  Generalized Self-Concordant Functions: A Recipe for Newton-Type Methods.  Mathematical Programing, 178, 145–213, 2019.
  * [23] N. Tupitsa, P. Dvurechensky, A. Gasnikov, and C. Uribe,  Multimarginal optimal transport by accelerated gradient descent,  in Proceedings of the 59th IEEE Conference on Decision and Control, (2020) pp. 6132–6137.
  * [24] C. Villani,  Optimal Transport – Old and New Springer-Verlag Berlin Heidelberg, (2009).

Appendices

This supplementary material is organized as follows:

  * •

In Appendix A we provide some basic facts on Bregman divergences.

  * •

Appendix B contains the proofs of the results stated in Section 3, notably
Proposition 3.1, and gives more information on the implementation of Algorithm
1.

  * •

Finally, in Appendix C we provide the proof of the main results of Sec. 4,
concerning the linear convergence and iteration complexity of Algorithm 1.

For the reader’s convenience all results presented in the main body of the
paper are restated in this supplementary material.

##  Appendix A Bregman divergences and Bregman projections

In this section we recall few facts on Bregman distance and Bregman
projections onto affine sets. In the following $X$ is an Euclidean space and
$\phi\colon X\to\left]-\infty,+\infty\right]$ is an extended-real valued
function. The set of minimizers of the function $\phi$ is denoted by
$\operatorname*{\mathrm{arg\,min}}_{x\in X}\phi(x)$, the _domain_ of $\phi$ is
$\text{\rm dom}\,\phi:=\\{x\in X\,|\,\phi(x)<+\infty\\}$ and $\phi$ is
_proper_ when $\text{\rm dom}\,\phi\neq\varnothing$. The function $\phi$ is
_convex_ if $\phi(tx+(1-t)y)\leq t\phi(x)+(1-t)\phi(y)$ for all
$x,y\in\text{\rm dom}\,\phi$ and $t\in[0,1]$. If the above inequality is
strict when $0<t<1$ and $x\neq y$, the function is _strictly convex_. The
function $\phi$ is _closed_ if the sublevel sets $\\{x\in X\,|\,\phi(x)\leq
t\\}$ are closed in $X$ for any $t\in\mathbb{R}$. For a convex function
$\phi\colon X\to\left]-\infty,+\infty\right]$, we denote by $\phi^{*}$ its
_Fenchel conjugate_ , that is, $\phi^{*}\colon
X\to\left]-\infty,+\infty\right]$, $\phi^{*}(y):=\sup_{x\in X}\\{\langle
x,y\rangle-\phi(x)\\}$. The conjugate of a convex function is always closed
and convex, and if $\phi$ is proper closed and convex, then
$(\phi^{*})^{*}=\phi$.

A proper closed and convex function $\phi$ is _essentially smooth_ if it is
differentiable on $\mathrm{int}(\text{\rm dom}\,\phi)\neq\varnothing$, and
$\lVert\nabla\phi(x_{n})\rVert\to+\infty$ whenever $x_{n}\in\text{\rm
int}(\text{\rm dom}\,\phi)$ and $x_{n}\to x\in\mathrm{bdry}(\text{\rm
dom}\,\phi)$. The function $\phi$ is _essentially strictly convex_ if
$\text{\rm int}(\text{\rm dom}\,\phi^{*})\neq\varnothing$ and is strictly
convex on every convex subset of $\text{\rm dom}\,\partial\phi$. A _Legendre_
function is a proper closed and convex function which is also essentially
smooth and essentially strictly convex. A function is Legendre if and only if
its conjugate is so. Moreover, if $\phi$ is a Legendre function, then
$\nabla\phi\colon\text{\rm int}(\text{\rm dom}\,\phi)\to\text{\rm
int}(\text{\rm dom}\,\phi^{*})$ and $\nabla\phi^{*}\colon\text{\rm
int}(\text{\rm dom}\,\phi^{*})\to\text{\rm int}(\text{\rm dom}\,\phi)$ are
bijective, inverses of each other, and continuous. Given a Legendre function
$\phi$, the _Bregman distance_ associated to $\phi$ is the function
$D_{\phi}\colon X\times X\to[0,+\infty]$ such that

$D_{\phi}(x,y)=\begin{cases}\phi(x)-\phi(y)-\langle
x-y,\nabla\phi(y)\rangle&\text{if }y\in\text{\rm int}(\text{\rm dom}\,\phi)\\\
+\infty&\text{otherwise}.\end{cases}$ (A.1)

######  Fact A.1.

Let $\phi$ be a Legendre function on $X$. Then the following hold

  1. (i)

$(\forall\,\pi,\gamma\in\text{\rm dom}\,\phi)\quad
D_{\phi}(\pi,\gamma)+D_{\phi}(\gamma,\pi)=\langle\pi-\gamma,\nabla\phi(\pi)-\nabla\phi(\gamma)\rangle$.

  2. (ii)

$(\forall\,\pi,\gamma\in\text{\rm dom}\,\phi)\quad
D_{\phi}(\pi,\gamma)=D_{\phi^{*}}(\nabla\phi(\gamma),\nabla\phi(\pi))$.

  3. (iii)

If $\phi$ is twice differentiable, then
$(\forall\,\pi,\gamma\in\text{\rm
dom}\,\phi)(\exists\,\xi\in[\pi,\gamma])\qquad
D_{\phi}(\pi,\gamma)=\dfrac{1}{2}\langle\nabla^{2}\phi(\xi)(\pi-\gamma),\pi-\gamma\rangle$,
where $[\pi,\gamma]=\\{(1-\alpha)\pi+\alpha\gamma\,|\,\alpha\in[0,1]\\}$ is
the segment with end points $\pi$ and $\gamma$.

  4. (iv)

Suppose that $\phi$ is twice differentiable on $\text{\rm int}(\text{\rm
dom}\,\phi)$. Then

$\big{(}\forall\pi\\!\in\text{\rm int}(\text{\rm
dom}\,\phi),\,\nabla^{2}\phi(\pi)\text{ is
invertible}\big{)}\Leftrightarrow\big{(}\phi^{*}\\!\text{ is twice
differentiable}\big{)}.$ (A.2)

######  Fact A.2.

Let $\phi$ be a Legendre function on $X$ and $\phi^{*}$ be it’s Fenchel
conjugate. If $\text{\rm dom}\,\phi^{*}$ is open, then the following hold

  1. (i)

For every $\pi\in\text{\rm int}(\text{\rm dom}\,\phi)$, the sublevel sets of
$D_{\phi}(\pi,\cdot)$ are compact.

  2. (ii)

For every $\pi\in\text{\rm int}(\text{\rm dom}\,\phi)$, and every sequence
$(\gamma_{k})_{k\in\mathbb{N}}$ in $\text{\rm int}(\text{\rm dom}\,\phi)$

$D_{\phi}(\pi,\gamma_{k})\to 0\ \Rightarrow\ \gamma_{k}\to\pi.$ (A.3)

Let ${\cal C}\subset X$ be an affine set, represented as follows

$A\colon X\to Y,\qquad b\in\text{\rm Im}(A),\qquad{\cal C}:=\\{\pi\in
X\,|\,A\pi=b\\},$ (A.4)

for some linear operator $A$ between $X$ and another Euclidean space $Y$.
Given a Legendre function $\phi\colon X\to\left]-\infty,+\infty\right]$ and
$\pi\in\text{\rm int}(\text{\rm dom}\,\phi)$, the _Bregman projection_ of
$\pi$ onto ${\cal C}$ is defined as the unique solution, denoted by ${\cal
P}^{\phi}_{\cal C}(\pi)$, of the optimization problem

$\min_{\gamma\in{\cal C}}D_{\phi}(\gamma,\pi)=\min_{\gamma\in{\cal
C}}\phi(\gamma)-\phi(\pi)-\langle\gamma-\pi,\nabla\phi(\pi)\rangle$ (A.5)

and the optimal value defines the _Bregman distance from $\pi$ to ${\cal C}$_
and is denoted by $D_{{\cal C}}^{\phi}(\pi)$. The dual of the above problem is

$\min_{\lambda\in
Y}\phi^{*}(\nabla\phi(\pi)+A^{*}\lambda)-\phi^{*}(\nabla\phi(\pi))-\langle
b,\lambda\rangle$ (A.6)

and strong duality holds, meaning that

$\min_{\gamma\in{\cal C}}D_{\phi}(\gamma,\pi)=-\min_{\lambda\in
Y}\big{[}\phi^{*}(\nabla\phi(\pi)+A^{*}\lambda)-\phi^{*}(\nabla\phi(\pi))-\langle
b,\lambda\rangle\big{]}.$ (A.7)

Moreover, the following KKT conditions hold for a couple
$(\pi^{\star},\lambda^{\star})$ solving the primal and dual problem above

$\pi^{\star}\in\text{\rm int}(\text{\rm dom}\,\phi),\quad
A\pi^{\star}=b\quad\text{and}\quad\nabla\phi(\pi)+A^{*}\lambda^{\star}=\nabla\phi(\pi^{\star}).$
(A.8)

Note that the KKT conditions characterizes the projection, so that

$\pi^{\star}={\cal P}^{\phi}_{\cal C}(\pi)\ \Leftrightarrow\
\big{(}\pi^{\star}\in\text{\rm int}(\text{\rm dom}\,\phi),\ A\pi^{\star}=b,\
\text{and}\ \nabla\phi(\pi)-\nabla\phi(\pi^{\star})\in\text{\rm
Im}(A^{*})\big{)}.$ (A.9)

Finally we mention the _generalized Pythagora’s theorem_. If ${\cal C}_{1}$ is
an affine set such that ${\cal C}\subset{\cal C}_{1}$, then, for every
$\pi\in\text{\rm int}(\text{\rm dom}\,\phi)$ it holds

$D_{\cal C}^{\phi}(\pi)=D_{{\cal C}_{1}}^{\phi}(\pi)+D_{{\cal C}}^{\phi}({\cal
P}^{\phi}_{{\cal C}_{1}}(\pi)).$ (A.10)

Moreover, in this case ${\cal P}^{\phi}_{\cal C}(\pi)={\cal P}^{\phi}_{\cal
C}({\cal P}^{\phi}_{{\cal C}_{1}}(\pi))$, and

$(\forall\gamma\in\text{\rm int}(\text{\rm
dom}\,\phi))\quad\nabla\phi(\gamma)-\nabla\phi(\pi)\in\text{\rm
Im}(A^{*})\Longleftrightarrow{\cal P}_{\cal C}(\gamma)={\cal P}_{\cal
C}(\pi).$ (A.11)

In the following we let $\phi\colon\mathbb{R}\to\left]-\infty,+\infty\right[$
be the (negative) _Boltzmann-Shannon entropy_ , that is,

$\phi(t)=\begin{cases}t\log t-t&\text{if }t>0\\\ 0&\text{if }t=0\\\
+\infty&\text{if }t<0.\end{cases}$

It is clear that $\phi^{*}(s)=\exp(s)$. We define the Bregman distance
associated to $\phi$

$\displaystyle D_{\phi}(s,t)$
$\displaystyle=\begin{cases}\phi(s)-\phi(t)-\phi^{\prime}(t)(s-t)&\text{if
}t>0\\\ +\infty&\text{otherwise}.\end{cases}$
$\displaystyle=\begin{cases}s\log\frac{s}{t}-s+t&\text{if }t>0\\\
+\infty&\text{otherwise},\end{cases}$ (A.12)

which is nothing but the _Kullback-Leibler divergence_ on $\mathbb{R}$.

######  Proposition A.3.

Let $M>0$. The following hold.

  1. (i)

The function $\phi$ is strongly convex on the interval $\left]0,M\right]$ with
modulus of strong convexity equal to $1/M$. Moreover, for every $a>0$ and
$s,t\in\mathbb{R}_{++}$, $D_{\phi}(s,t)=aD_{\phi}(s/a,t/a)$.

  2. (ii)

The function $\phi^{*}$ is strongly convex on the interval
$\left[-M,+\infty\right[$ with modulus of strong convexity equal to
$\exp(-M)$.

###### Proof.

(i): It follows from the fact that the second derivative of $\phi$ is
$\phi^{\prime\prime}(t)=1/t$, which is bounded from below away from zero on
the interval $\left]0,M\right]$ by the constant $1/M$. the second part follows
directly from the definition (A.12).

(ii): It follows from the fact that the second derivative of $\exp$ is bounded
from below away from zero on the interval $\left[-M,+\infty\right[$ by
$\exp(-M)$. ∎

The negative entropy and the Kullback-Leibler divergence on ${\mathbb{X}}$ are

$\mathsf{H}(\gamma)=\sum_{j}\phi(\gamma_{j})\quad\text{and}\quad\mathsf{KL}(\gamma,\pi)=\sum_{j}D_{\phi}(\gamma_{j},\pi_{j}).$
(A.13)

######  Lemma A.4.

Let $\pi,\gamma,\alpha\in{\mathbb{X}}_{++}$ and suppose that

$0<M_{1}\leq\min_{j}\frac{\min\\{\pi_{j},\gamma_{j}\\}}{\alpha_{j}}\leq\max_{j}\frac{\max\\{\pi_{j},\gamma_{j}\\}}{\alpha_{j}}\leq
M_{2}.$

Then, setting $\mathsf{A}=\alpha\odot(\cdot)\colon{\mathbb{X}}\to{\mathbb{X}}$
(which is a positive diagonal operator), we have

$\mathsf{KL}(\pi,\gamma)\geq\max\left\\{\frac{M_{1}}{2}\lVert\log\pi-\log\gamma\rVert_{\mathsf{A}}^{2},\frac{1}{2M_{2}}\lVert\pi-\gamma\rVert_{\mathsf{A}^{-1}}^{2}\right\\}.$
(A.14)

###### Proof.

It follows from Proposition A.3(i) that, for $a>0$ and $s,t>0$ such that
$s/a,t/a\leq M$, we have $D_{\phi}(s,t)=aD_{\phi}(s/a,t/a)\geq
a(2M)^{-1}\lvert s/a-t/a\rvert^{2}=(2M)^{-1}a^{-1}\lvert s-t\rvert^{2}$. Thus,
since $\gamma_{j}/\alpha_{j},\pi_{j}/\alpha_{j}\leq M_{2}$, we have

$\mathsf{KL}(\pi,\gamma)=\sum_{j}D_{\phi}(\pi_{j},\gamma_{j})\geq\frac{1}{2M_{2}}\sum_{j}\frac{1}{\alpha_{j}}\lvert\pi_{j}-\gamma_{j}\rvert^{2}=\frac{1}{2M_{2}}\lVert\pi-\gamma\rVert^{2}_{\mathsf{A}^{-1}}.$

Now, it follows from Proposition A.3(ii) that for every $a,s,t>0$ such that
$s/a,t/a\geq e^{-M}$ we have $\log(s/a),\log(t/a)\geq-M$ and hence
$D_{\phi}(s,t)=aD_{\phi}(s/a,t/a)=aD_{\phi^{*}}(\log(s/a),\log(t/a))\geq
a(e^{-M}/2)\lvert\log s-\log t\rvert^{2}$. Therefore, since
$\pi_{j}/\alpha_{j},\gamma_{j}/\alpha_{j}\geq M_{1}$, we have

$\mathsf{KL}(\pi,\gamma)=\sum_{j}D_{\phi}(\pi_{j},\gamma_{j})\geq\frac{M_{1}}{2}\sum_{j}\alpha_{j}\lvert\log\pi_{j}-\log\gamma_{j}\rvert^{2}=\frac{M_{1}}{2}\lVert\log\pi-\log\gamma\rVert_{\mathsf{A}}^{2}.\qed$

##  Appendix B BatchGreenkhorn algorithm and its implementation

Here we provide proofs of the results in Sec. 3.

See 3.1

###### Proof.

It follows from (3.4) that

$\Pi_{(k,L)}=\Big{\\{}\pi\in{\mathbb{X}}\,\big{|}\,\mathsf{R}_{(k,L)}(\pi)=\mathsf{a}_{k\lvert
L}\Big{\\}}.$ (B.1)

Then the first equality in (3.8) follows directly form the KKT conditions
(A.8), which in this case yields

$\mathsf{R}_{(k,L)}(\bar{\pi})=\mathsf{a}_{k\lvert
L}\quad\text{and}\quad\bar{\pi}=\nabla\mathsf{H}^{*}\big{(}\nabla\mathsf{H}(\pi)+\mathsf{R}^{*}_{(k,L)}(\bar{\mathsf{u}})\big{)},$
(B.2)

where, according to (A.6), the dual parameter
$\bar{\mathsf{u}}\in\mathbb{R}^{L}$ solves the minimization problem in (3.9).
Now, since for every $j\in\mathcal{J}$,
$(\nabla\mathsf{H}(\pi))_{j}=\log(\pi_{j})$ and
$(\nabla\mathsf{H}^{*}(\gamma))_{j}=\exp(\gamma_{j})$, then (B.2) gives

$\bar{\pi}_{j}=\exp(\log(\pi_{j})+\mathsf{R}^{*}_{(k,L)}(\bar{\mathsf{u}}))_{j})=\pi_{j}\exp\big{(}(\mathsf{R}^{*}_{(k,L)}(\bar{\mathsf{u}}))_{j}\big{)}$

and the second equality in (3.8) follows.

Now, let $J^{(k)}_{L}\colon\mathbb{R}^{L}\to\mathbb{R}^{n_{k}}$ be the
canonical injection of $\mathbb{R}^{L}$ into $\mathbb{R}^{n_{k}}$. Then,
recalling the definition of $\mathsf{R}_{(k,L)}$ in (3.3), we have
$\mathsf{R}_{(k,L)}=J^{(k)*}_{L}\mathsf{R}_{k}$ and hence
$\mathsf{R}^{*}_{(k,L)}=\mathsf{R}_{k}^{*}J_{L}^{(k)}$, where
$\mathsf{R}_{k}^{*}\colon\mathbb{R}^{n_{k}}\to{\mathbb{X}}$ acts as
$(\mathsf{R}_{k}^{*}\mathsf{v})_{j}=\mathsf{v}_{j_{k}}$. Therefore, for every
$j\in\mathcal{J}$,

$(\mathsf{R}^{*}_{(k,L)}\bar{\mathsf{u}})_{j}=(\mathsf{R}_{k}^{*}J^{(k)}_{L}\bar{\mathsf{u}})_{j}=(J^{(k)}_{L}\bar{\mathsf{u}})_{j_{k}}=\begin{cases}\bar{u}_{j_{k}}&\text{if
}j_{k}\in L\\\ 0&\text{otherwise}.\end{cases}$

Hence,

$\bar{\pi}_{j}=\pi_{j}\times\begin{cases}e^{\bar{u}_{j_{k}}}&\text{if
}j_{k}\in L\\\ 1&\text{otherwise}.\end{cases}$ (B.3)

On the other hand, since $\mathsf{a}_{k\lvert
L}=J^{(k)*}_{L}\mathsf{R}_{k}\bar{\pi}$, by (B.3), we derive that, for every
$j_{k}\in L$,

$a_{k,j_{k}}=(\mathsf{R}_{k}\bar{\pi})_{j_{k}}=\sum_{j_{-k}\in\mathcal{J}_{-k}}\bar{\pi}_{(j_{-k},j_{k})}=e^{\bar{u}_{j}}(\mathsf{R}_{k}\pi)_{j},$
(B.4)

so that $e^{\bar{u}_{j_{k}}}=a_{k,j_{k}}/(\mathsf{R}_{k}\pi)_{j_{k}}$. Hence
now (3.10) follows from (B.3). Concerning the formula for the distance, by
(3.10), we have that,

$\displaystyle D^{\phi}_{\Pi_{(k,L)}}(\pi)$
$\displaystyle=D^{\phi}(\bar{\pi},\pi)$
$\displaystyle=\sum_{j\in\mathcal{J}}\bar{\pi}_{j}\log\Big{(}\frac{\bar{\pi}_{j}}{\pi_{j}}\Big{)}-\bar{\pi}_{j}+\pi_{j}$
$\displaystyle=\sum_{j_{k}\in
L}\sum_{j_{-k}\in\mathcal{J}_{-k}}\bar{\pi}_{(j_{-k},j_{k})}\log\Big{(}\frac{a_{k,j_{k}}}{(\mathsf{R}_{k}\pi)_{j_{k}}}\Big{)}-\bar{\pi}_{(j_{-k},j_{k})}+\pi_{(j_{-k},j_{k})}$
$\displaystyle=\sum_{j_{k}\in
L}a_{k,j_{k}}\log\Big{(}\frac{a_{k,j_{k}}}{(\mathsf{R}_{k}\pi)_{j_{k}}}\Big{)}-a_{k,j_{k}}+(\mathsf{R}_{k}\pi)_{j_{k}}$
$\displaystyle=\mathsf{KL}({\mathsf{a}_{k}}_{\lvert
L},\mathsf{R}_{k}(\pi)_{\lvert L}),$

which completes the proof. ∎

The following proposition justifies equation (3.12) and Algorithm 1.

######  Proposition B.1.

Let $(\pi^{t})_{t\in\mathbb{N}}$ be defined according to algorithm (3.5).
Then, we have

$(\forall\,t\in\mathbb{N})\qquad\pi^{t}=\exp\Big{(}-\mathsf{C}/\eta+\bigoplus_{k=1}^{m}\mathsf{v}^{t}_{k}\Big{)}\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k},$
(B.5)

where $\mathsf{v}^{t}_{k}=(v^{t}_{k,j})_{1\leq j\leq
n_{k}}\in\mathbb{R}^{n_{k}}$ and

$\mathsf{v}_{k}^{t+1}=\delta_{k,k_{t}}J_{L_{t}}^{(k_{t})}\mathsf{u}^{t}+\mathsf{v}_{k}^{t},\quad\mathsf{u}^{t}=\log\mathsf{a}_{k_{t}\lvert
L_{t}}-\log(\mathsf{R}_{k_{t}}\pi^{t})_{\lvert L_{t}},$ (B.6)

$J_{L_{t}}^{(k_{t})}\colon\mathbb{R}^{L_{t}}\to\mathbb{R}^{n_{k_{t}}}$ is the
canonical injection, $\delta_{k,k_{t}}$ is the Kronecker symbol, and
$\mathsf{v}^{0}_{k}$, $k\in[m]$ are arbitrary.

###### Proof.

By definition of $\pi^{t+1}$ we have

$\pi^{t+1}={\cal
P}_{\Pi_{(k_{t},L_{t})}}(\pi^{t}),\quad\Pi_{(k_{t},L_{t})}=\big{\\{}\pi\in{\mathbb{X}}\,|\,J_{L_{t}}^{(k_{t})*}\mathsf{R}_{k_{t}}\pi=\mathsf{a}_{k_{t}\lvert
L_{t}}\big{\\}}.$ (B.7)

Then it follows from Proposition 3.1 that

$\nabla\mathsf{H}(\pi^{t+1})=\nabla\mathsf{H}(\pi^{t})+\mathsf{R}_{k_{t}}^{*}J^{(k_{t})}_{L_{t}}\mathsf{u}^{t}$
(B.8)

with $\mathsf{u}^{t}\in\mathbb{R}^{L_{t}}$. Therefore, applying the above
equation recursively, we get

$\displaystyle\nabla\mathsf{H}(\pi^{t})$
$\displaystyle=\nabla\mathsf{H}(\pi^{0})+\sum_{s=0}^{t-1}\mathsf{R}_{k_{s}}^{*}J^{(k_{s})}_{L_{s}}\mathsf{u}^{s}$
$\displaystyle=\nabla\mathsf{H}(\pi^{0})+\sum_{s=0}^{t-1}\sum_{k=1}^{m}\mathsf{R}_{k}^{*}\delta_{k,k_{s}}J^{(k_{s})}_{L_{s}}\mathsf{u}^{s}$
$\displaystyle=\nabla\mathsf{H}(\pi^{0})+\sum_{k=1}^{m}\sum_{s=0}^{t-1}\mathsf{R}_{k}^{*}\delta_{k,k_{s}}J^{(k_{s})}_{L_{s}}\mathsf{u}^{s}$
$\displaystyle=\nabla\mathsf{H}(\pi^{0})+\sum_{k=1}^{m}\mathsf{R}_{k}^{*}\Big{(}\sum_{s=0}^{t-1}\delta_{k,k_{s}}J^{(k_{s})}_{L_{s}}\mathsf{u}^{s}\Big{)}.$

Then if we set
$\mathsf{v}_{k}^{t}=\sum_{s=0}^{t-1}\delta_{k,k_{s}}J^{(k_{s})}_{L_{s}}\mathsf{u}^{s}\in\mathbb{R}^{n_{k}}$,
we have (recalling that
$\pi_{0}=e^{-\mathsf{C}/\eta}\odot\otimes_{k=1}^{m}\mathsf{a}_{k}$)

$\pi^{t}_{j}=\exp\Big{(}\log(\pi^{0}_{j})+\sum_{k=1}^{m}(\mathsf{R}_{k}^{*}\mathsf{v}^{t}_{k})_{j}\Big{)}=\exp\Big{(}-\mathsf{C}_{j}/\eta+\sum_{k=1}^{m}v^{t}_{k,j_{k}}\Big{)}\prod_{k=1}^{m}a_{k,j_{k}}.$

Moreover, it is clear that

$\mathsf{v}_{k}^{t+1}=\sum_{s=0}^{t}\delta_{k,k_{s}}J_{L_{s}}^{(k_{s})}\mathsf{u}^{s}=\delta_{k,k_{t}}J_{L_{t}}^{(k_{t})}\mathsf{u}^{t}+\mathsf{v}_{k}^{t}$

Finally, it follows from (3.9) that, for every $j_{k_{t}}\in L_{t}$,
$e^{u_{j_{k_{t}}}^{t}}=a_{k_{t},j_{k_{t}}}/(\mathsf{R}_{k}\pi^{t})_{j_{k_{t}}}$,
so that

$(\forall\,j_{k_{t}}\in L_{t})\quad u^{t}_{j_{k_{t}}}=\log
a_{k_{t},j_{k_{t}}}-\log\big{(}(\mathsf{R}_{k}\pi^{t})_{j_{k_{t}}}\big{)}.$

The statement follows. ∎

Next we give more detailed implementation of the batch Greenkhorn given in
Algorithm 1.

######  Remark B.2 (Implementation details on Algorithm 1).

The most delicate part is to avoid recomputing the marginals
$\mathsf{R}_{k}(\pi^{t})=\mathsf{r}_{k}^{t}=(r_{k,j_{k}}^{t})_{j_{k}\in[n_{k}]}$,
$k\in[m]$, (step 5) that are necessary for making the greedy choice in step 3.
Now, due to equation (3.10) in Proposition 3.1, we have that

$\pi^{t+1}_{j}=\pi^{t}_{j}\times\begin{cases}\dfrac{a_{k_{t},j_{k_{t}}}}{{\mathsf{R}_{k_{t}}(\pi^{t})}_{j_{k_{t}}}}&\text{if
}j_{k_{t}}\in L_{t},\\\\[8.61108pt] 1&\text{otherwise}\end{cases}$ (B.9)

and hence

$r_{k_{t},j_{k_{t}}}^{t+1}=\begin{cases}a_{k_{t},j_{k_{t}}}&j_{k_{t}}\in
L_{t},\\\ r_{k_{t},j_{k_{t}}}^{t}&\text{otherwise}.\end{cases}$ (B.10)

To derive update formula for the other marginals, observe that for all $k\neq
k_{t}$ it follows from (B.5) that

$\displaystyle r_{k,j_{k}}^{t+1}$
$\displaystyle=\sum_{j_{-k}\in\mathcal{J}_{-k}}\exp\Big{(}\log\pi^{0}_{(j_{-k},j_{k})}+\sum_{h\neq
k}v^{t+1}_{h,j_{h}}+v^{t+1}_{k,j_{k}}\Big{)}$
$\displaystyle=\sum_{j_{-k}\in\mathcal{J}_{-k}}\exp\Big{(}\log\pi^{0}_{(j_{-k},j_{k})}+\sum_{h\not\in\\{k,k_{t}\\}}v^{t}_{h,j_{h}}+v^{t+1}_{k_{t},j_{k_{t}}}+v^{t}_{k,j_{k}}\Big{)}.$

So, since according to (B.6),
$v^{t+1}_{k_{t},j_{k_{t}}}=v^{t}_{k_{t},j_{k_{t}}}$ if $j_{k_{t}}\not\in
L_{t}$ and
$v^{t+1}_{k_{t},j_{k_{t}}}=v^{t}_{k_{t},j_{k_{t}}}+\log(a_{k_{t},j_{k_{t}}}/r^{t}_{k_{t},j_{k_{t}}})$
if $j_{k_{t}}\in L_{t}$, we have

$\displaystyle r_{k,j_{k}}^{t+1}=$
$\displaystyle\sum_{\begin{subarray}{c}j_{-k}\in\mathcal{J}_{-k}\\\
j_{k_{t}}\not\in
L_{t}\end{subarray}}\exp\Big{(}\log\pi^{0}_{(j_{-k},j_{k})}+\sum_{h\neq
k}v^{t}_{h,j_{h}}+v^{t}_{k,j_{k}}\Big{)}\;$
$\displaystyle\quad+\sum_{\begin{subarray}{c}j_{-k}\in\mathcal{J}_{-k}\\\
j_{k_{t}}\in
L_{t}\end{subarray}}\exp\Big{(}\log\pi^{0}_{(j_{-k},j_{k})}+\sum_{h\neq
k}v^{t}_{h,j_{h}}+v^{t}_{k,j_{k}}\Big{)}\frac{a_{k_{t},j_{k_{t}}}}{r^{t}_{k_{t},j_{k_{t}}}}$
$\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}j_{-k}\in\mathcal{J}_{-k}\\\
j_{k_{t}}\in
L_{t}\end{subarray}}\exp\Big{(}\log\pi^{0}_{(j_{-k},j_{k})}+\sum_{h\neq
k}v^{t}_{h,j_{h}}+v^{t}_{k,j_{k}}\Big{)}\Big{(}\frac{a_{k_{t},j_{k_{t}}}}{r^{t}_{k_{t},j_{k_{t}}}}-1\Big{)}+r_{k,j_{k}}^{t}.$

Therefore, at each iteration $t\geq 0$ we will construct an auxiliary tensor
$\widetilde{\pi}^{t}\in\mathbb{R}_{+}^{n_{1}\times\cdots\times
n_{k_{t}-1}\times 1\times n_{k_{t}+1}\times\cdots\times n_{m}}$ by

$\displaystyle\widetilde{\pi}^{t}_{j_{1},\ldots,j_{k_{t}-1},1,j_{k_{t}+1},\ldots,j_{m}}=$
$\displaystyle\displaystyle{\sum_{j_{k_{t}}\in
L_{t}}\exp\Big{(}\log\pi^{0}_{j_{1},\ldots,j_{k_{t}-1},j_{k_{t}},j_{k_{t}+1},\ldots,j_{m}}+\sum_{k\in[m]}v^{t}_{k,j_{k}}}$
(B.11) $\displaystyle\quad+\log\lvert
a_{k_{t},j_{k_{t}}}-r^{t}_{k_{t},j_{k_{t}}}\rvert-\log(r^{t}_{k_{t},j_{k_{t}}})\Big{)}\text{\rm
sgn}(a_{k_{t},j_{k_{t}}}-r^{t}_{k_{t},j_{k_{t}}}),$ (B.12)

in order to obtain that for every $k\neq k_{t}$,
$\mathsf{r}^{t+1}_{k}=\mathsf{r}^{t+1}_{k}+\mathsf{R}_{k}(\hat{\pi}^{t})$.
Hence, we can use $\hat{\pi}^{t}$ to efficiently update non-active marginals
without recomputing them from scratch. Moreover, note that using (B.11)-(B.12)
one avoids excessive numerical errors when $a_{k,j}\approx r_{k,j}^{t}$. These
observations lead us to the following implementation of BatchGreenkhorn.

Input: $\mathsf{C}\in{\mathbb{X}}_{+}$, $\eta>0$,
$(\mathsf{a}_{1},\dots,\mathsf{a}_{m})$, $(\tau_{1},\dots,\tau_{m})$,
$1\leq\tau_{k}\leq n_{k}$, $\varepsilon>0$

1 Initialization: $t=0$, $\mathsf{v}^{0}_{k}=0$, $k\in[m]$,
$\mathsf{r}_{k}^{0}=(r^{0}_{k,j})_{j\in[n_{k}]}=\mathsf{R}_{k}(\exp(-\mathsf{C}/\eta)\odot\otimes_{k=1}^{m}\mathsf{a}_{k})$

2

3while  _$\sum_{k\in[m]}\lVert\mathsf{a}_{k}-\mathsf{r}^{t}_{k}\rVert_{1}
>\varepsilon$_ do

4 for  _$k\in[m]$_ do

5 Compute vectors $\mathsf{p}_{k}$ as
$p_{k,j}:=\mathsf{KL}(a_{k,j},r^{t}_{k,j})$, for $j\in[n_{k}]$

6 Take $L_{k}^{\prime}$ to be $\tau_{k}$ largest elements of $\mathsf{p}_{k}$

7 end for

8 Choose the marginal with the best batch:
$k_{t}\leftarrow\operatorname*{\mathrm{arg\,max}}_{k\in[m]}\lVert{\mathsf{p}_{k}}_{\lvert
L_{k}^{\prime}}\rVert_{1}$ and $L_{t}=L_{k_{t}}^{\prime}$

9 Set $\mathsf{v}_{k}^{t+1}=\mathsf{v}_{k}^{t}$ and update
${\mathsf{v}_{k_{t}}^{t+1}}_{\lvert
L_{t}}\leftarrow{\mathsf{v}^{t+1}_{k_{t}}}_{\lvert
L_{t}}+\log({\mathsf{a}_{k_{t}}}_{\lvert
L_{t}})-\log({\mathsf{r}^{t}_{k_{t}}}_{\lvert L_{t}})$

10 Set $\mathsf{r}_{k_{t}}^{t+1}=\mathsf{r}_{k_{t}}^{t}$ and update
${\mathsf{r}_{k_{t}}^{t+1}}_{\lvert L_{t}}={\mathsf{a}_{k_{t}}}_{\lvert
L_{t}}$

11 For $k\in[m]\setminus\\{k_{t}\\}$ update
$\mathsf{r}^{t+1}_{k}\leftarrow\mathsf{r}^{t}_{k}+\mathsf{R}_{k}(\widetilde{\pi}^{t})$,
where $\widetilde{\pi}^{t}$ is given by (B.11)–(B.12)

12 Set $t\leftarrow t+1$

13

14 end while

Result: $\\{\mathsf{v}^{t}_{k}\\}_{k\in[m]}$

Algorithm 2 $\mathtt{BatchGreenkhorn}(\mathsf{C},\eta,\rho,\tau,\varepsilon)$

######  Remark B.3.

Let us assume that $\tau_{k}=\tau$ and $n_{k}=n$ for all $k\in[m]$ and that
$m<<n$. Then, we can conclude that the cost of one iteration of
BatchGreenkhorn is essentially determined by step 10 of Algorithm 2 which is
performed in ${\cal O}(\tau n^{m-1})$ operations. Hence, one iteration of the
MultiSinkhorn (i.e, BatchGreenkhorn with a full batch $\tau=n$) has the same
order of computational cost as $n/\tau$ iterations of BatchGreenkhorn with a
batch size $\tau$. So, we can introduce the normalized iteration counter as
$t=t_{\tau}n/\tau$, where $t_{\tau}$ is the iteration counter for the
BatchGreenkhorn with a batch size $\tau$.

##  Appendix C Convergence of Batch Greenkhorn algorithm

Here we provide proofs of the main results given in Sec. 4. We first set
notation for the rest of the section. Given $k\in[m]$ and $L\subset[n_{k}]$,
we denote by

$J_{k}\colon\mathbb{R}^{n_{k}}\to\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{k}}\times\cdots\times\mathbb{R}^{m},\quad\mathsf{v}_{k}\mapsto(0,\dots,0,\mathsf{v}_{k},0,\dots,0),$
(C.1)

the canonical injection of $\mathbb{R}^{n_{k}}$ into
$\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{k}}\times\cdots\times\mathbb{R}^{n_{m}}$
and by

$J_{L}^{(k)}\colon\mathbb{R}^{L}\to\mathbb{R}^{n_{k}}$ (C.2)

the canonical injection of $\mathbb{R}^{L}$ into $\mathbb{R}^{n_{k}}$.

We note that, referring to the operators $\mathsf{R}$ and $\mathsf{R}_{k}$
defined in (3.1) and (2.3), respectively, we have

$\mathsf{R}^{*}\colon\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}\to{\mathbb{X}},\qquad\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m})=\sum_{k=1}^{m}\mathsf{R}_{k}^{*}(\mathsf{v}_{k})=\bigotimes_{k=1}^{m}\mathsf{v}_{k},$
(C.3)

where
$(\mathsf{R}_{k}^{*}(\mathsf{v}_{k}))_{j_{1},\dots,j_{k},\dots,j_{m}}=v_{k,j_{k}}$.
Indeed the second equality in (C.3) follows from the fact that for every
$j\in\mathcal{J}$,
$(\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m}))_{j}=\sum_{k=1}^{m}(\mathsf{R}_{k}^{*}\mathsf{v}_{k})_{j}=\sum_{k=1}^{m}v_{k,j_{k}}=\big{(}\bigoplus_{k=1}^{m}\mathsf{v}_{k}\big{)}_{j}$.
We note also that $\mathsf{R}_{k}=J^{*}_{k}\circ\mathsf{R}$, since $J_{k}^{*}$
is the $k$-th canonical projection.

Then we provide a result concerning the properties of optimal potentials.

######  Lemma C.1.

Let $\pi^{\star}$ be the solution of RMOT given by (2.4). Then
$\pi^{\star}={\cal P}_{\Pi}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})$ and, for
every $k\in[m]$, there exist $\mathsf{v}_{k}^{\star}=(v^{\star}_{k,j})_{1\leq
j\leq n_{k}}\in\mathbb{R}^{n_{k}}$, such that

$\pi^{\star}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\bigoplus_{k=1}^{m}\mathsf{v}^{\star}_{k}\Big{)}\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k},$
(C.4)

and the $\mathsf{v}_{k}^{\star}$’s, can be chosen so that

$\sum_{k\in[m]}\lVert\mathsf{v}_{k}^{\star}\rVert_{\infty}\leq(4m-3)\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}.$
(C.5)

Moreover, if $m=2$, then $\mathsf{v}_{1}^{\star}$ and $\mathsf{v}_{2}^{\star}$
can be chosen such that

$\max_{k\in[m]}\lVert\mathsf{v}_{k}^{\star}\rVert_{\infty}\leq\frac{3}{2}\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}.$
(C.6)

###### Proof.

Since, by definition $\pi^{\star}={\cal P}_{\Pi}(\xi)$, it easy to see, from
the characterization of the projection given in (A.9), that

$\pi^{\star}={\cal P}_{\Pi}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})\
\Leftrightarrow\
\nabla\mathsf{H}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})-\nabla\mathsf{H}(\xi)\in\text{\rm
Im}(\mathsf{R}^{*}).$

Thus, since
$\nabla\mathsf{H}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})-\nabla\mathsf{H}(\xi)=\log\otimes_{k=1}^{m}\mathsf{a}_{k}=\oplus_{k=1}^{m}\log\mathsf{a}_{k}=\mathsf{R}^{*}(\log\mathsf{a}_{1},\ldots,\log\mathsf{a}_{m})\in\text{\rm
Im}(\mathsf{R}^{*})$, we have that ${\cal
P}_{\Pi}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})={\cal
P}_{\Pi}(\xi)=\pi^{\star}$. Now, it follows from the KKT conditions (A.8) for
the projection of $\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k}$ onto affine set
$\Pi$, that

$\pi^{\star}=\nabla\mathsf{H}^{*}\big{(}\nabla\mathsf{H}(\xi\odot\otimes_{k=1}^{m}\mathsf{a}_{k})+\mathsf{R}^{*}(\mathsf{v}_{1}^{\star},\dots,\mathsf{v}_{m}^{\star})\big{)}$

for some
$(\mathsf{v}_{1}^{\star},\dots,\mathsf{v}_{m}^{\star})\in\mathbb{R}^{n_{1}}\times\dots\times\mathbb{R}^{n_{m}}$.
Since $\nabla\mathsf{H}^{*}=\exp$ and $\nabla\mathsf{H}=\log$, (C.4) follows.
Next, observe that for every $k\in[m]$, since
$\mathsf{R}_{k}(\pi^{\star})=\mathsf{a}_{k}$, using (C.4), we obtain that for
every $j_{k}\in[n_{k}]$,

$\exp(v^{\star}_{k,j_{k}})\sum_{j_{-k}\in\mathcal{J}_{-k}}\exp(-\mathsf{C}_{(j_{-k},j_{k})}/\eta+\sum_{h\neq
k}^{m}v^{\star}_{h,j_{h}})\prod_{h\neq k}a_{h,j_{h}}=1.$ (C.7)

Hence, the vectors $\mathsf{v}^{\star}_{1},\ldots,\mathsf{v}^{\star}_{m}$
solve a (discrete) Schrödinger system, and we can apply the results from [7,
Lemma 3.1] and [12, Theorem 2.8] to obtain (C.5) and (C.6), respectively. ∎

######  Lemma C.2.

Let $\mathsf{A}\colon{\mathbb{X}}\to{\mathbb{X}}$ and
$\mathsf{A}_{k}\colon\mathbb{R}^{n_{k}}\to\mathbb{R}^{n_{k}}$ be the diagonal
and positive operators defined as
$\mathsf{A}(\pi)=\pi\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k}$ and
$\mathsf{A}_{k}\mathsf{v}_{k}=\mathsf{v}_{k}\odot\mathsf{a}_{k}$,
respectively. Let
$(\mathsf{v}_{1},\dots,\mathsf{v}_{m})\in\mathbb{R}^{n_{1}}\times\dots\times\mathbb{R}^{n_{m}}$.
Then

$(\forall\,k\in[m])\quad\lVert\mathsf{R}_{k}^{*}\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}}=\lVert\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}_{k}}.$
(C.8)

Moreover, if $\langle\mathsf{v}_{k},\mathsf{a}_{k}\rangle=0$ for every
$k=1,\dots,m-1$, then

$\lVert\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m})\rVert^{2}_{\mathsf{A}}=\sum_{k=1}^{m}\lVert\mathsf{R}_{k}^{*}\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}}=\sum_{k=1}^{m}\lVert\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}_{k}}.$
(C.9)

###### Proof.

Let $k\in[m]$. Then, recalling that
$(\mathsf{R}_{k}^{*}(\mathsf{v}_{k}))_{j}=v_{k,j_{k}}$, we have

$\displaystyle\lVert\mathsf{R}_{k}^{*}\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}}$
$\displaystyle=\sum_{j\in\mathcal{J}}v^{2}_{k,j_{k}}\prod_{\ell=1}^{m}a_{\ell,j_{\ell}}$
$\displaystyle=\sum_{j_{-k}\in\mathcal{J}_{-k}}\sum_{j_{k}=1}^{n_{k}}v^{2}_{k,j_{k}}a_{k,j_{k}}\prod_{\ell\neq
k}^{m}a_{\ell,j_{\ell}}$ $\displaystyle=\prod_{\ell\neq
k}\bigg{(}\sum_{\ell=1}^{n_{\ell}}a_{\ell,j_{\ell}}\bigg{)}\sum_{j_{k}=1}^{n_{k}}v^{2}_{k,j_{k}}a_{k,j_{k}}.$

Since $\sum_{\ell=1}^{n_{\ell}}a_{\ell,j_{\ell}}=1$, we get
$\lVert\mathsf{R}_{k}^{*}\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}}=\lVert\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}_{k}}$
and the first part of the statement follows. Concerning the second part,
equation (C.9) will follow if we prove that, for every $k,h\in[m]$ with $k\neq
h$, we have that $\mathsf{R}_{k}^{*}(\mathsf{v}_{k})$ and
$\mathsf{R}_{h}^{*}(\mathsf{v}_{h})$ are orthogonal in the metric
$\langle\cdot,\cdot\rangle_{\mathsf{A}}$. Thus, let $k,h\in[m]$ and suppose
(w.l.o.g.) that $k<h$. Then

$\displaystyle\langle\mathsf{R}_{k}^{*}(\mathsf{v}_{k}),\mathsf{R}_{h}^{*}(\mathsf{v}_{h})\rangle_{\mathsf{A}}$
$\displaystyle=\sum_{j\in\mathcal{J}}v_{k,j_{k}}v_{h,j_{h}}\prod_{\ell=1}^{m}a_{\ell,j_{\ell}}$
$\displaystyle=\sum_{\begin{subarray}{c}j_{1},\dots,j_{k-1},j_{k+1},\dots,\\\
j_{h-1},j_{h+1},\dots,j_{m}\end{subarray}}\sum_{j_{k}=1}^{n_{k}}\sum_{j_{h}=1}^{n_{h}}v_{k,j_{k}}v_{h,j_{h}}a_{k,j_{k}}a_{h,j_{h}}\prod_{\ell\neq
k,\ell\neq h}^{m}a_{\ell,j_{\ell}}$
$\displaystyle=\bigg{(}\sum_{j_{k}=1}^{n_{k}}v_{k,j_{k}}\bigg{)}\bigg{(}\sum_{j_{h}=1}^{n_{h}}v_{h,j_{h}}\bigg{)}\prod_{\ell\neq
k,\ell\neq h}\bigg{(}\sum_{\ell=1}^{n_{\ell}}a_{\ell,j_{\ell}}\bigg{)}.$

Since $\sum_{\ell=1}^{n_{\ell}}a_{\ell,j_{\ell}}=1$ and for $k<m$,
$\sum_{j_{k}=1}^{n_{k}}v_{k,j_{k}}a_{k,j_{k}}=\langle\mathsf{v}_{k},\mathsf{a}_{k}\rangle=0$,
we have
$\langle\mathsf{R}_{k}^{*}(\mathsf{v}_{k}),\mathsf{R}_{h}^{*}(\mathsf{v}_{h})\rangle_{\mathsf{A}}=0$
and hence the statement follows. ∎

See 4.1

###### Proof.

We start by recalling the two formulas

$\pi^{t}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\mathsf{V}^{t}\Big{)}\odot\alpha\quad\text{and}\quad\pi^{\star}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\mathsf{V}^{\star}\Big{)}\odot\alpha,$
(C.10)

where $\alpha:=\bigotimes_{k=1}^{m}\mathsf{a}_{k}$,
$\mathsf{V}^{t}:=\bigoplus_{k=1}^{m}\mathsf{v}^{t}_{k}$, and
$\mathsf{V}^{\star}:=\bigoplus_{k=1}^{m}\mathsf{v}_{k}^{\star}$. Moreover,
since for every $(\lambda_{k})_{k\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$
such that $\sum_{k=1}^{m}\lambda_{k}=0$, we have
$\bigoplus_{k=1}^{m}(\mathsf{v}^{t}_{k}+\lambda_{k})=\bigoplus_{k=1}^{m}\mathsf{v}^{t}_{k}$
and
$\bigoplus_{k=1}^{m}(\mathsf{v}_{k}^{\star}+\lambda_{k})=\bigoplus_{k=1}^{m}\mathsf{v}_{k}^{\star}$,
we can choose the dual variables $(\mathsf{v}^{t}_{k})_{k\in[n_{k}]}$ and
$(\mathsf{v}^{\star}_{k})_{k\in[n_{k}]}$ so that

$(\forall\,k=1,\dots,m-1)\quad\langle\mathsf{v}^{t}_{k},\mathsf{a}_{k}\rangle=0\quad\text{and}\quad\langle\mathsf{v}^{\star}_{k},\mathsf{a}_{k}\rangle=0.$
(C.11)

First, observe that Pythgora’s theorem yields that
$\mathsf{KL}_{\Pi}(\pi^{t+1})=\mathsf{KL}_{\Pi}(\pi^{t})-\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})\leq\mathsf{KL}_{\Pi}(\pi^{t})$,
which implies that, for every $t\geq 0$,
$D_{\mathsf{H}^{*}}(\log\pi^{t},\log\pi^{*})=\mathsf{KL}(\pi^{\star},\pi^{t})\leq\mathsf{KL}_{\Pi}(\pi^{0})<+\infty$.
Howevery, since $\mathsf{H}^{*}$ is a Legendre function, the sublevel sets of
$D_{\mathsf{H}^{*}}(\cdot,\log\pi^{*})$ are bounded, and hence the sequence
$(\log\pi^{t})_{t\in\mathbb{N}}$ is bounded in ${\mathbb{X}}$. Now, since the
first of (C.10) yields that
$\log\pi^{t}=-\mathsf{C}/\eta+\mathsf{V}^{t}+\log\alpha$, we have that also
the sequence $(\mathsf{V}^{t})_{t\in\mathbb{N}}$ is bounded in ${\mathbb{X}}$.
Thus let $M_{1}>0$ be such that

$\lVert\mathsf{V}^{\star}\rVert_{\infty},\lVert\mathsf{V}^{t}\rVert_{\infty}\leq
M_{1}\quad(\forall\,t\in\mathbb{N}).$

Then, recalling (C.10),

$\frac{\pi^{t}}{\alpha}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\mathsf{V}^{t}\Big{)}\geq\exp\Big{(}-\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}-\lVert\mathsf{V}^{t}\rVert_{\infty}\Big{)}\geq\exp\Big{(}-\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}-M_{1}\Big{)}$

and

$\frac{\pi^{\star}}{\alpha}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\mathsf{V}^{\star}\Big{)}\geq\exp\Big{(}-\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}-\lVert\mathsf{V}^{\star}\rVert_{\infty}\Big{)}\geq\exp\Big{(}-\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}-M_{1}\Big{)}$

and hence

$\exp\Big{(}-\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}-M_{1}\Big{)}\leq\min\Big{\\{}\frac{\pi^{t}}{\alpha},\frac{\pi^{\star}}{\alpha}\Big{\\}}\quad(\forall\,t\in\mathbb{N}).$
(C.12)

Let $t\in\mathbb{N}$, $k\in[m]$ and $L\subset[n_{k}]$. It follows from (3.10)
that

$(\forall\,j\in\mathcal{J})\quad({\cal
P}_{\Pi_{(k,L)}}(\pi^{t}))_{j}=\pi^{t}_{j}\times\begin{cases}\dfrac{a_{k,j_{k}}}{{\mathsf{R}_{k}(\pi)}_{j_{k}}}&\text{if
}j_{k}\in L,\\\\[8.61108pt] 1&\text{otherwise}\end{cases}$ (C.13)

and hence

$(\forall\,j\in\mathcal{J})\quad\frac{({\cal
P}_{\Pi_{(k,L)}}(\pi^{t}))_{j}}{\alpha_{j}}\leq\frac{\pi^{t}_{j}}{\alpha_{j}}\max\bigg{\\{}1,\dfrac{a_{k,j_{k}}}{{\mathsf{R}_{k}(\pi)}_{j_{k}}}\bigg{\\}}.$
(C.14)

Now, since $\mathsf{C}\geq 0$, we have

$\frac{\pi^{t}}{\alpha}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\mathsf{V}^{t}\Big{)}\leq\exp(\mathsf{V}^{t})\leq\exp(\lVert\mathsf{V}\rVert_{\infty})\leq\exp(M_{1})$
(C.15)

and

$\displaystyle\dfrac{{\mathsf{R}_{k}(\pi)}_{j_{k}}}{a_{k,j_{k}}}$
$\displaystyle=\dfrac{\sum_{j_{-k}\in\mathcal{J}_{-k}}\pi^{t}_{(j_{-k},j_{k})}}{a_{k,j_{k}}}$
$\displaystyle=\dfrac{\sum_{j_{-k}\in\mathcal{J}_{-k}}\exp\big{(}-\mathsf{C}_{(j_{-k},j_{k})}/\eta+\mathsf{V}^{t}_{(j_{-k},j_{k})}\big{)}\prod_{h=1}^{m}a_{h,j_{h}}}{a_{k,j_{k}}}$
$\displaystyle=\sum_{j_{-k}\in\mathcal{J}_{-k}}\exp\big{(}-\mathsf{C}_{(j_{-k},j_{k})}/\eta+\mathsf{V}^{t}_{(j_{-k},j_{k})}\big{)}\prod_{h\neq
k}^{m}a_{h,j_{h}}$
$\displaystyle\geq\exp\big{(}-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1}\big{)}\sum_{j_{-k}\in\mathcal{J}_{-k}}\prod_{h\neq k}^{m}a_{h,j_{h}}$
$\displaystyle=\exp\big{(}-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1}\big{)}\prod_{h\neq
k}^{m}\bigg{(}\sum_{j_{h}=1}^{n_{h}}a_{h,j_{h}}\bigg{)}$
$\displaystyle=\exp\big{(}-\lVert\mathsf{C}\rVert_{\infty}/\eta-M_{1}\big{)},$
(C.16)

since $\sum_{j_{h}=1}^{n_{h}}a_{h,j_{h}}=1$. Therefore, by (C.14), (C.15), and
(C),

$\frac{{\cal
P}_{\Pi_{(k,L)}}(\pi^{t})}{\alpha}\leq\exp(M_{1})\exp\big{(}\lVert\mathsf{C}\rVert_{\infty}/\eta+M_{1}\big{)}=\exp\big{(}\lVert\mathsf{C}\rVert_{\infty}/\eta+2M_{1}\big{)}$

and hence, recalling (C.15),

$\max\bigg{\\{}\frac{\pi^{t}}{\alpha},\frac{{\cal
P}_{\Pi_{(k,L)}}(\pi^{t})}{\alpha}\bigg{\\}}\leq\exp\big{(}\lVert\mathsf{C}\rVert_{\infty}/\eta+2M_{1}\big{)}.$
(C.17)

We now prove that

$\exp\big{(}-2\lVert\mathsf{C}\rVert_{\infty}-3M_{1}\big{)}b^{-1}_{\tau}\mathsf{KL}_{\Pi}(\pi^{t})=\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})=\max_{(k,L)\in\mathcal{I}(\tau)}\mathsf{KL}({\cal
P}_{\Pi_{(k,L)}}(\pi^{t}),\pi^{t}).$ (C.18)

From this inequality it will follow, using the Pythgora’s theorem
$\mathsf{KL}_{\Pi}(\pi^{t+1})+\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})=\mathsf{KL}_{\Pi}(\pi^{t})$,
that

$\exp(-2\lVert\mathsf{C}\rVert_{\infty}-3M_{1})b^{-1}_{\tau}\mathsf{KL}_{\Pi}(\pi^{t})=\mathsf{KL}_{\Pi}(\pi^{t})-\mathsf{KL}_{\Pi}(\pi^{t+1})$
(C.19)

and hence

$\mathsf{KL}_{\Pi}(\pi^{t+1})\leq\big{(}1-\exp(-2\lVert\mathsf{C}\rVert_{\infty}-3M_{1})b^{-1}_{\tau}\big{)}\mathsf{KL}_{\Pi}(\pi^{t}),$
(C.20)

which gives the statement. Thus, it remains to prove (C.18). Let
$t\in\mathbb{N}$ and, for the sake of brevity set

$\pi:=\pi^{t},\quad\pi_{k}^{h}:={\cal
P}_{\Pi_{(k,L_{k}^{h})}}(\pi^{t})\quad\mathsf{v}_{k}:=\mathsf{v}_{k}^{\star}-\mathsf{v}^{t}_{k},\quad\text{
and }\quad\mathsf{v}_{k}^{h}:=J^{(k)*}_{L^{h}_{k}}\mathsf{v}_{k}.$ (C.21)

Let, for every $k\in[m]$, $(L_{k}^{h})_{1\leq h\leq s_{k}}$ be a partition of
$[n_{k}]$ made of non empty sets of cardinality exactly $\tau_{k}$ possibly
except for the last one, such that $L_{k_{t-1}}^{1}=L_{t-1}$, where
$s_{k}=\lceil n_{k}/\tau_{k}\rceil$. Then, it follows from (C.10) that

$\frac{\pi^{\star}}{\pi}=\frac{\exp(-\mathsf{C}/\eta)\exp(\mathsf{V}^{\star})}{\exp(-\mathsf{C}/\eta)\exp(\mathsf{V}^{t})}=\exp(\mathsf{V}^{\star}-\mathsf{V}^{t}).$
(C.22)

Hence, recalling that
$\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m})=\sum_{k=1}^{m}\mathsf{R}_{k}^{*}\mathsf{v}_{k}=\bigoplus_{k=1}^{m}\mathsf{v}_{k}=\bigoplus_{k=1}^{m}(\mathsf{v}_{k}^{\star}-\mathsf{v}_{k}^{t})$,
we have

$\displaystyle\mathsf{KL}_{\Pi}(\pi)+\mathsf{KL}(\pi,\pi^{\star})$
$\displaystyle=\mathsf{KL}(\pi^{\star},\pi)+\mathsf{KL}(\pi,\pi^{\star})$
$\displaystyle=\langle\pi^{\star}-\pi,\log(\pi^{\star}/\pi)\rangle$
$\displaystyle=\langle\pi^{\star}-\pi,\mathsf{V}^{\star}-\mathsf{V}^{t}\rangle$
$\displaystyle=\langle\pi^{\star}-\pi,\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m})\rangle$
$\displaystyle=\sum_{k=1}^{m}\langle\pi^{\star}-\pi,\mathsf{R}^{*}_{k}(\mathsf{v}_{k})\rangle.$
(C.23)

Moreover, recalling that
$\mathsf{R}_{(k,L_{k}^{h})}=J^{(k)*}_{L_{k}^{h}}\circ\mathsf{R}_{k}^{*}$ and
$\mathsf{v}_{k}^{h}=J^{(k)*}_{L^{h}_{k}}(\mathsf{v}_{k})$, we have

$\mathsf{R}_{k}^{*}(\mathsf{v}_{k})=\mathsf{R}_{k}^{*}\bigg{(}\sum_{h=1}^{n_{h}}J^{(k)}_{L_{k}^{h}}\circ
J^{(k)*}_{L_{k}^{h}}(\mathsf{v}_{k})\bigg{)}=\sum_{h=1}^{n_{h}}\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})$
(C.24)

and hence

$\mathsf{KL}_{\Pi}(\pi)+\mathsf{KL}(\pi,\pi^{\star})=\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\langle\pi^{\star}-\pi,\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle.$
(C.25)

Now, recalling the general definition of $\Pi_{(k,L)}$ in (3.4), since
$\pi_{k}^{h}$ and $\pi^{\star}$ both belong to $\Pi_{(k,L_{k}^{h})}$, we have
that $\pi^{\star}-\pi_{k}^{h}\in\text{\rm
Ker}(\mathsf{R}_{(k,L_{k}^{h})})=\text{\rm
Im}(\mathsf{R}^{*}_{(k,L_{k}^{h})})^{\perp}$ and hence

$\langle\pi^{\star}-\pi,\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle=\langle\pi^{\star}-\pi_{k}^{h},\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle+\langle\pi_{k}^{h}-\pi,\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle=\langle\pi_{k}^{h}-\pi,\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle$

and hence

$\displaystyle\mathsf{KL}_{\Pi}(\pi)+\mathsf{KL}(\pi,\pi^{\star})$
$\displaystyle=\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\langle\pi_{k}^{h}-\pi,\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle$
$\displaystyle=\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\langle\mathsf{A}^{-1}(\pi_{k}^{h}-\pi),\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle_{\mathsf{A}},$
(C.26)

where $\mathsf{A}$ is the positive diagonal operator defined in Lemma C.2.
Now, it follows from (C.12), Lemma A.4, and Lemma C.2 that

$\displaystyle\mathsf{KL}(\pi^{t},\pi^{\star})$
$\displaystyle\geq(1/2)\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\lVert\log\pi^{\star}-\log\pi^{t}\rVert^{2}_{\mathsf{A}}$
$\displaystyle=(1/2)\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\lVert\mathsf{V}^{\star}-\mathsf{V}^{t}\rVert^{2}_{\mathsf{A}}$
$\displaystyle=(1/2)\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\lVert\mathsf{R}^{*}(\mathsf{v}_{1},\dots,\mathsf{v}_{m})\rVert^{2}_{\mathsf{A}}$
$\displaystyle=(1/2)\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\sum_{k=1}^{m}\lVert\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}_{k}}.$

Moreover, recalling the definition of $\mathsf{v}_{k}^{h}$ in (C.21), since
$\mathsf{v}_{k}=\sum_{h=1}^{s_{k}}J_{(k,L_{k}^{h})}\mathsf{v}^{h}_{k}$ and
$(J_{(k,L_{k}^{h})}\mathsf{v}^{h}_{k})_{h\in[s_{k}]}$ is a finite orthogonal
sequence in $\mathbb{R}^{n_{k}}$ w.r.t. the metric
$\langle\cdot,\cdot\rangle_{\mathsf{A}_{k}}$, we have

$\lVert\mathsf{v}_{k}\rVert^{2}_{\mathsf{A}_{k}}=\sum_{h=1}^{s_{k}}\lVert
J_{(k,L_{k}^{h})}\mathsf{v}^{h}_{k}\rVert^{2}_{\mathsf{A}_{k}}=\sum_{h=1}^{s_{k}}\lVert\mathsf{R}_{k}^{*}J_{(k,L_{k}^{h})}\mathsf{v}^{h}_{k}\rVert^{2}_{\mathsf{A}}=\sum_{h=1}^{s_{k}}\lVert\mathsf{R}_{(k,L_{k}^{h})}^{*}\mathsf{v}^{h}_{k}\rVert^{2}_{\mathsf{A}},$
(C.27)

where we used equation (C.8) from Lemma C.2 applied to
$J_{(k,L_{k}^{h})}\mathsf{v}^{h}_{k}$ and the fact that, by definition,
$\mathsf{R}_{(j,L_{k}^{h})}=J^{(k)*}_{L_{k}^{h}}\mathsf{R}_{k}$. Overall we
get that

$\mathsf{KL}(\pi^{t},\pi^{\star})\geq(1/2)\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\lVert\mathsf{R}_{(k,L_{k}^{h})}^{*}\mathsf{v}^{h}_{k}\rVert^{2}_{\mathsf{A}}$

and hence (C) yields

$\displaystyle\mathsf{KL}_{\Pi}(\pi)$
$\displaystyle\leq\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\langle\mathsf{A}^{-1}(\pi_{k}^{h}-\pi),\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle_{\mathsf{A}}-\mathsf{KL}(\pi,\pi^{\star})$
$\displaystyle\leq\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\langle\mathsf{A}^{-1}(\pi_{k}^{h}-\pi),\mathsf{R}^{*}_{(k,L_{k}^{h})}(\mathsf{v}_{k}^{h})\rangle_{\mathsf{A}}-\frac{1}{2}\exp(-\lVert\mathsf{C}\rVert_{\infty}/\eta-
M_{1})\lVert\mathsf{R}_{(k,L_{k}^{h})}^{*}\mathsf{v}^{h}_{k}\rVert^{2}_{\mathsf{A}}$
$\displaystyle\leq\frac{\exp(\lVert\mathsf{C}\rVert_{\infty}/\eta+M_{1})}{2}\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\lVert\mathsf{A}^{-1}(\pi_{k}^{h}-\pi)\rVert^{2}_{\mathsf{A}}$
$\displaystyle=\frac{\exp(\lVert\mathsf{C}\rVert_{\infty}/\eta+M_{1})}{2}\sum_{k=1}^{m}\sum_{h=1}^{s_{k}}\lVert(\pi_{k}^{h}-\pi)\rVert^{2}_{\mathsf{A}^{-1}},$

where in the last inequality we used the Young-Fenchel inequality $\langle
a,b\rangle_{\mathsf{A}}\leq\tfrac{\mu}{2}\lVert
a\rVert_{\mathsf{A}}^{2}+\tfrac{1}{2\mu}\lVert b\rVert_{\mathsf{A}}^{2}$. Now,
recalling that we set $\pi=\pi^{t}$ and $\pi_{k}^{h}={\cal
P}_{\Pi_{(k,L_{k}^{h})}}(\pi^{t})$, it follows from (C.17) and Lemma A.4 that
$\frac{1}{2}\lVert\pi_{k}^{h}-\pi\rVert_{\Theta^{-1}}^{2}\leq\exp\big{(}\lVert\mathsf{C}\rVert_{\infty}/\eta+2M_{1}\big{)}\mathsf{KL}(\pi_{k}^{h},\pi)$,
and consequently

$\displaystyle\mathsf{KL}_{\Pi}(\pi^{t})$
$\displaystyle\leq\exp\big{(}2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1}\big{)}\sum_{k\in[m]}\sum_{h\in[s_{k}]}\mathsf{KL}(\pi_{k}^{h},\pi)$
$\displaystyle\leq\exp\big{(}2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1}\big{)}\bigg{(}\sum_{k\in[m]}s_{k}-1\bigg{)}\max_{k\in[m]}\max_{h\in[s_{k}]}\mathsf{KL}(\pi_{k}^{h},\pi)$
$\displaystyle=\exp\big{(}2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1}\big{)}\bigg{(}\sum_{k\in[m]}\lceil
n_{k}/\tau_{k}\rceil-1\bigg{)}\max_{k\in[m]}\max_{h\in[s_{k}]}\mathsf{KL}({\cal
P}_{\Pi_{(k,L_{k}^{h})}}(\pi^{t}),\pi^{t})$
$\displaystyle\leq(b_{\tau}-1)\exp\big{(}2\lVert\mathsf{C}\rVert_{\infty}/\eta+3M_{1}\big{)}\max_{(k,L)\in{\cal
I}(\tau)}\mathsf{KL}({\cal P}_{\Pi_{(k,L)}}(\pi^{t}),\pi^{t}),$

where in the second inequality we used that, for $k=k_{t-1}$ and $h=1$,
$\pi_{k}^{h}={\cal P}_{\Pi_{(k_{t-1},L_{t-1})}}(\pi^{t})=\pi^{t}$ (since by
definition $\pi^{t}\in\Pi_{(k_{t-1},L_{t-1})}$), so that
$\mathsf{KL}(\pi_{k}^{h},\pi)=0$. This proves (C.18) and the proof is
complete. ∎

We now provide a result concerning the convergence of numerical sequences
which is critical to analyze the iteration complexity of the algorithm. This
result has been first showed implicitly in [13]. We provide here a more
explicit version together with a complete proof for the reader’s convenience.

######  Lemma C.3.

Let $M,C>0$ and let $(\delta_{t})_{t\in\mathbb{N}}$ and
$(d_{t})_{t\in\mathbb{N}}$ be two sequences of positive numbers such that, for
every $t\in\mathbb{N}$,

  1. (i)

$\delta_{t}-\delta_{t+1}\geq\bigg{(}\dfrac{d_{t}}{C}\bigg{)}^{2}$,

  2. (ii)

$\delta_{t}\leq Md_{t}$.

Let $\varepsilon>0$ and set
$\bar{t}=\min\\{t\in\mathbb{N}\,|\,d_{t}\leq\varepsilon\\}$. Then $\bar{t}\leq
1+2MC^{2}/\varepsilon$.

###### Proof.

Items (i) and (ii) imply that

$\delta_{t}-\delta_{t+1}\geq\bigg{(}\frac{\delta_{t}}{MC}\bigg{)}^{2}.$

Therefore, since $\delta_{t}\geq\delta_{t+1}$, we have

$\delta_{t}-\delta_{t+1}\geq\frac{\delta^{2}_{t}}{M^{2}C^{2}}\geq\frac{\delta_{t}\delta_{t+1}}{M^{2}C^{2}}$

and hence, dividing by $\delta_{t}\delta_{t+1}$,

$\frac{1}{\delta_{t+1}}-\frac{1}{\delta_{t}}\geq\frac{1}{M^{2}C^{2}}.$

Thus,

$\frac{1}{\delta_{t}}-\frac{1}{\delta_{0}}=\sum_{i=0}^{t-1}\bigg{(}\frac{1}{\delta_{i}}-\frac{1}{\delta_{i+1}}\bigg{)}\geq\frac{t}{M^{2}C^{2}}$

and hence we get the following rate of convergence for the sequence
$(\delta_{t})_{t\in\mathbb{N}}$

$\delta_{t}\leq\bigg{(}\frac{1}{\delta_{0}}+\frac{t}{M^{2}C^{2}}\bigg{)}^{-1}.$
(C.28)

Now, we let $\delta\in\left]0,\delta_{0}\right]$. We wish to determine the
number of iterations such that $\delta_{t}\leq\delta$. It follows from (C.28)
that

$\bigg{(}\frac{1}{\delta_{0}}+\frac{t}{M^{2}C^{2}}\bigg{)}^{-1}\leq\delta\
\Leftrightarrow\ \
\frac{1}{\delta_{0}}+\frac{t}{M^{2}C^{2}}\geq\frac{1}{\delta}\
\Leftrightarrow\ \ t\geq
M^{2}C^{2}\bigg{(}\frac{1}{\delta}-\frac{1}{\delta_{0}}\bigg{)}.$ (C.29)

This means that if we take $t\geq M^{2}C^{2}(1/\delta-1/\delta_{0})$, we have
$\delta_{t}\leq\delta$ as desired. So we set $t=\lfloor
M^{2}C^{2}(1/\delta-1/\delta_{0})\rfloor+1$. Then, we have
$\delta_{t}\leq\delta$. Now we have to cases. Suppose that $t<\bar{t}$ and let
$s\in\mathbb{N}$ be such that $t+s=\bar{t}-1$. Then, for every $i=0,\dots,s$,
since $t+i<\bar{t}$, we have $d_{t+i}>\varepsilon$ and hence, using (i),

$\delta\geq\delta_{t}-\delta_{\bar{t}}=\sum_{i=0}^{s}\bigg{(}\frac{1}{\delta_{t+i}}-\frac{1}{\delta_{t+i+1}}\bigg{)}\geq\sum_{i=0}^{s}\frac{d^{2}_{t+i}}{M^{2}C^{2}}\geq(s+1)\frac{\varepsilon^{2}}{M^{2}C^{2}},$
(C.30)

which implies that $s+1\leq C^{2}\delta/\varepsilon^{2}$. Overall we have

$\bar{t}=t+s+1\leq\bigg{\lfloor}M^{2}C^{2}\bigg{(}\frac{1}{\delta}-\frac{1}{\delta_{0}}\bigg{)}\bigg{\rfloor}+1+C^{2}\frac{\delta}{\varepsilon^{2}}\leq
1+\frac{M^{2}C^{2}}{\delta}-\frac{M^{2}C^{2}}{\delta_{0}}+\frac{C^{2}\delta}{\varepsilon^{2}}.$
(C.31)

Note that this inequality is true for any
$\delta\in\left]0,\delta_{0}\right]$. Now, suppose that
$M\varepsilon\leq\delta_{0}$. Then we have

$\bar{t}\leq
1+\min_{\delta\in\left]0,\delta_{0}\right]}\bigg{(}\frac{M^{2}C^{2}}{\delta}+\frac{C^{2}\delta}{\varepsilon^{2}}\bigg{)}=1+2\frac{MC^{2}}{\varepsilon},$
(C.32)

where the minimum is attained at
$\delta=M\varepsilon\in\left]0,\delta_{0}\right]$. On the other hand, if
$\delta_{0}<M\varepsilon$, then the minimum on the right hand side of (C.31)
is attained at $\delta=\delta_{0}$ and hence

$\bar{t}\leq 1+\frac{C^{2}\delta_{0}}{\varepsilon^{2}}\leq
1+\frac{C^{2}M\varepsilon}{\varepsilon^{2}}=1+\frac{MC^{2}}{\varepsilon}.$
(C.33)

In any case, the statement follows. ∎

######  Remark C.4.

The statement of Lemma C.3 is equivalent to the fact that the sequence
$(\min_{0\leq s<t}d_{s})_{t\in\mathbb{N}}$ converge to zero with rate ${\cal
O}(1/t)$, i.e., that for every integer $t>1$,

$\min_{0\leq s<t}d_{s}\leq\frac{2MC^{2}}{t-1}.$

Next, we prove the main result on the iteration complexity.

See 4.2

###### Proof.

For the sake of brevity let $\bar{b}=\max_{k\in[n]}\lceil
n_{k}/\tau_{k}\rceil$ and set, for every $t\in\mathbb{N}$,
$\delta_{t}:=\mathsf{KL}_{\Pi}(\pi^{t})$. Let $t\in\mathbb{N}$ be arbitrary.
Recalling (C.23), we have that

$\delta_{t}=\mathsf{KL}_{\Pi}(\pi^{t})\leq\sum_{k\in[m]}\langle\pi^{\star}-\pi^{t},\mathsf{R}_{k}^{*}(\mathsf{v}_{k}^{\star}-\mathsf{v}_{k}^{t})\rangle\\!=\\!\sum_{k\in[m]}\langle\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t}),\mathsf{v}_{k}^{\star}-\mathsf{v}_{k}^{t}\rangle,$

which, using Hoölder inequality, yields

$\delta_{t}\leq\sum_{k\in[m]}\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}\;\lVert\mathsf{v}_{k}^{\star}-\mathsf{v}_{k}^{t}\rVert_{\infty}\leq
M_{2}d_{t}.$ (C.34)

Now, we prove

$\delta_{t}-\delta_{t+1}\geq\min\Big{\\{}\frac{d_{t}^{2}}{5\bar{b}},\frac{\delta_{t}^{2}}{4M_{2}^{2}\bar{b}}\Big{\\}}\geq\frac{\delta_{t}^{2}}{5M_{2}^{2}\bar{b}}.$
(C.35)

Let for every $k\in[m]$, $(L_{k}^{h})_{1\leq h\leq[s_{k}]}$, $s_{k}:=\lceil
n_{k}/\tau_{k}\rceil$, be a partition of $[n_{k}]$ made of nonempty sets of
cardinality exactly $\tau_{k}$, except maybe for the last one, such that
$L_{k_{t}}^{1}=L_{t}$ (not that necessarily the cardinality of $L_{t}$ is
$\tau_{k_{t}}$). Then, according to the greedy choice of $(k_{t},L_{t})$ we
have that

$\bar{b}\,\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})\geq\max_{k\in[m]}s_{k}\max_{h\in[s_{k}]}\;\mathsf{KL}_{\Pi_{(k,L_{k}^{h})}}(\pi^{t})\geq\max_{k\in[m]}\sum_{h\in[s_{k}]}\mathsf{KL}_{\Pi_{(k,L_{k}^{h})}}(\pi^{t}).$

Thus, equation (3.11) of Proposition 3.1 yields

$\bar{b}\,\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})\geq\max_{k\in[m]}\sum_{h\in[s_{k}]}\mathsf{KL}({\mathsf{a}_{k}}_{\lvert
L_{k}^{h}},\mathsf{R}_{k}(\pi^{t})_{\lvert
L_{k}^{h}})=\max_{k\in[m]}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t})).$
(C.36)

Now Pinsker’s inequality guaranties that, for every $k\in[m]$

$\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))\geq\frac{\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}^{2}}{\tfrac{2}{3}\lVert\mathsf{a}_{k}\rVert_{1}+\tfrac{4}{3}\lVert\mathsf{R}_{k}(\pi^{t})\rVert_{1}}=\frac{\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}^{2}}{\tfrac{2}{3}+\tfrac{4}{3}\lVert\pi^{t}\rVert_{1}}\geq\frac{\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}^{2}}{2+\tfrac{4}{3}\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}},$
(C.37)

where in the second inequality we used that
$\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}\geq\lVert\mathsf{R}_{k}(\pi^{t})\rVert_{1}-\lVert\mathsf{a}_{k}\rVert_{1}=\lVert\pi^{t}\rVert_{1}-1$.
Thus, solving the quadratic inequality in
$\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}\geq 0$ we can conclude
that

$\lVert\mathsf{a}_{k}-\mathsf{R}_{k}(\pi^{t})\rVert_{1}\leq\tfrac{2}{3}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))+\sqrt{(\tfrac{2}{3}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t})))^{2}+2\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))}.$

Therefore, if
$\max_{k\in[m]}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))\leq 1$,
then $2+4d_{t}/3\leq 5$, and consequently,
$\max_{k\in[m]}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))\geq
d_{t}^{2}/5$, which, using Pythagoras theorem and (C.36), yields

$\delta_{t}-\delta_{t+1}=\mathsf{KL}_{\Pi}(\pi^{t})-\mathsf{KL}_{\Pi}(\pi^{t+1})=\mathsf{KL}_{\Pi_{(k_{t},L_{t})}}(\pi^{t})\geq\frac{d_{t}^{2}}{5\bar{b}}.$
(C.38)

On the other hand, if
$\max_{k\in[m]}\mathsf{KL}(\mathsf{a}_{k},\mathsf{R}_{k}(\pi^{t}))>1$, it
follows again from Pythagoras theorem and (C.36), that
$\delta_{t}-\delta_{t+1}\geq 1/\bar{b}$. Moreover, since
$\delta_{t}\leq\delta_{0}\leq M_{2}d_{0}\leq
M_{2}(1+\lVert\pi^{0}\rVert_{1})\leq 2M_{2}$, we have that
$1\geq\delta_{t}^{2}/(4M_{2}^{2})$, and (C.35) follows.

Now, similarly to what was done in the proof of Lemma C.3 we can derive from
(C.35) that

$\delta_{t}\leq\bigg{(}\frac{1}{\delta_{0}}+\frac{t}{5M_{2}^{2}\bar{b}}\bigg{)}^{-1}.$
(C.39)

Thus, if we take $r=\lfloor 5M_{2}^{2}\bar{b}\rfloor+1$ we have
$\delta_{r}\leq 1$. Then, by (3.11) with $L=[n_{k}]$ and Pythagoras theorem we
have that, for every $t\in\mathbb{N}$,

$\mathsf{KL}({\mathsf{a}_{k}},\mathsf{R}_{k}(\pi^{r+t}))=\mathsf{KL}({\cal
P}_{\Pi_{(k,[n_{k}])}}(\pi^{r+t}),\pi^{r+t})\leq\mathsf{KL}(\pi^{\star},\pi^{r+t})=\delta_{r+t}\leq\delta_{r}\leq
1.$

Thus,
$\max_{k\in[m]}\mathsf{KL}({\mathsf{a}_{k}},\mathsf{R}_{k}(\pi^{r+t}))\leq 1$
and for what we already saw,

$\delta_{r+t}-\delta_{r+t+1}\geq\frac{d^{2}_{r+t}}{5\bar{b}}.$ (C.40)

In the end the sequence $(\delta_{r+t})_{t\in\mathbb{N}}$ satisfies the two
assumptions of Lemma C.3 with $C=\sqrt{5\bar{b}}$ and $M=M_{2}$. Thus, we can
conclude that the smallest $t$ so that $d_{r+t}\leq\varepsilon$ satisfies
$t\leq 1+10M_{2}\bar{b}/\varepsilon$. Hence

$r+t\leq 2+\frac{10M_{2}\bar{b}}{\varepsilon}+5M_{2}^{2}\bar{b}\leq
2+\frac{10M_{2}\bar{b}}{\varepsilon}+\frac{5M_{2}^{2}\bar{b}\eta}{\varepsilon}=2+\frac{5M_{2}\bar{b}}{\varepsilon}(2+M_{2}\eta).\qed$

The next two results are based on novel bounds on potentials that imply
explicit dependence of constant $M>0$ in the global rate (4.1) on the given
data: $\mathsf{a}_{1},\ldots,\mathsf{a}_{m}$, $\mathsf{C}$ and $\eta$.

See 4.4

###### Proof.

Let $\mathsf{v}^{t}_{k}$, $k\in[m]$, $t\geq 0$ be given by Algorithm 1. Then
from Proposition B.1 we have that for every $t\geq 0$

$\pi^{t+1}=\exp\Big{(}-\frac{\mathsf{C}}{\eta}+\bigoplus_{k=1}^{m}\mathsf{v}^{t+1}_{k}\Big{)}\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k},$
(C.41)

with $\mathsf{v}_{k}^{0}=0$ and for $t\geq 0$, $k\in[m]$ and
$j_{k}\in[n_{k}]$,

$v_{k,j_{k}}^{t+1}=\begin{cases}v_{k,j_{k}}^{t}+\log(a_{k,j_{k}})-\log(\mathsf{R}_{k}(\pi^{t})_{j_{k}})&k=k_{t},\,j_{k}\in
L_{t},\\\ v_{k,j_{k}}^{t}&\text{otherwise}.\end{cases}$

So, to bound $\log\pi^{t}$, we will bound $\mathsf{v}^{t}_{k}$, $k\in[m]$.
Since $\mathsf{R}_{k_{t}}(\pi^{t+1})_{j_{k_{t}}}=a_{k_{t},j_{k_{t}}}$ for all
$j_{k_{t}}\in L_{t}$, (C.41) implies that

$1=\exp(v^{t+1}_{k_{t},j_{k_{t}}})\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{t}_{k,j_{k}}\Big{)}\prod_{k\neq k_{t}}a_{k,j_{k}},$

and, hence,

$\exp(-v^{t+1}_{k_{t},j_{k_{t}}})=\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{t}_{k,j_{k}}\Big{)}\prod_{k\neq k_{t}}a_{k,j_{k}}.$ (C.42)

So, using (C.7) we obtain that for every $j_{k_{t}}\in L_{t}$

$\exp(v^{t+1}_{k_{t},j_{k_{t}}}-v^{\star}_{k_{t},j_{k_{t}}})=\dfrac{\displaystyle\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{\star}_{k,j_{k}}\Big{)}}{\displaystyle\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{t}_{k,j_{k}}\Big{)}},$

while

$\exp(v^{\star}_{k_{t},j_{k_{t}}}-v^{t+1}_{k_{t},j_{k_{t}}})=\dfrac{\displaystyle\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\\!\\!\\!\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{t}_{k,j_{k}}\Big{)}}{\displaystyle\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\\!\\!\\!\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}v^{\star}_{k,j_{k}}\Big{)}}.$

However, since in general, $\alpha_{i},\beta_{i}>0$ implies that
$(\sum_{i}\alpha_{i})/(\sum_{i}\beta_{i})\leq\max_{i}\alpha_{i}/\beta_{i}$,
the last two equalities give

$\exp(v^{t+1}_{k_{t},j_{k_{t}}}-v^{\star}_{k_{t},j_{k_{t}}})\leq\max_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\exp\Big{(}\sum_{k\neq
k_{t}}(v^{\star}_{k,j_{k}}-v^{t}_{k,j_{k}})\Big{)}$

and

$\exp(v^{\star}_{k_{t},j_{k_{t}}}-v^{t+1}_{k_{t},j_{k_{t}}})\leq\max_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\exp\Big{(}\sum_{k\neq
k_{t}}(v^{t}_{k,j_{k}}-v^{\star}_{k,j_{k}})\Big{)}.$

Hence, taking the logarithm we obtain that for every $j_{k_{t}}\in L_{t}$,

$\lvert
v^{t+1}_{k_{t},j_{k_{t}}}-v^{\star}_{k_{t},j_{k_{t}}}\rvert\leq\max_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\Big{\lvert}\sum_{k\neq
k_{t}}(v^{\star}_{k,j_{k}}-v^{t}_{k,j_{k}})\Big{\rvert}\leq\max_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\sum_{k\neq
k_{t}}\lvert v^{\star}_{k,j_{k}}-v^{t}_{k,j_{k}}\rvert=\sum_{k\neq
k_{t}}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}.$

Therefore, since $v^{t+1}_{k,j_{k}}=v^{t}_{k,j_{k}}$ if $k\neq k_{t}$ or
$j_{k_{t}}\notin L_{t}$,

$\max\Big{\\{}\lVert\mathsf{v}^{t+1}_{k_{t}}-\mathsf{v}^{\star}_{k_{t}}\rVert_{\infty},\sum_{k\neq
k_{t}}\lVert\mathsf{v}^{t+1}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\Big{\\}}\leq\max\Big{\\{}\lVert\mathsf{v}^{t}_{k_{t}}-\mathsf{v}^{\star}_{k_{t}}\rVert_{\infty},\sum_{k\neq
k_{t}}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\Big{\\}},$

and, since $m=2$,

$\max_{k\in[m]}\lVert\mathsf{v}^{t+1}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\leq\max_{k\in[m]}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty},$

which implies, recalling that $\mathsf{v}^{0}=0$, that, for all $t\geq 0$,
$\max_{k\in[m]}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\leq\max_{k\in[m]}\lVert\mathsf{v}^{\star}_{k}\rVert_{\infty}$.
Now, in view of (C.6) in Lemma C.1, we have

$\max_{k\in[m]}\lVert\mathsf{v}^{\star}_{k}\rVert_{\infty}\leq\frac{3}{2}\frac{\lVert\mathsf{C}\rVert_{\infty}}{\eta}$
(C.43)

and hence, since
$\lVert\mathsf{v}^{t}_{k}\rVert_{\infty}\leq\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}+\lVert\mathsf{v}^{\star}_{k}\rVert_{\infty}$,
$\max_{k\in[m]}\lVert\mathsf{v}_{k}^{t}\rVert_{\infty}\leq
2\max_{k\in[m]}\lVert\mathsf{v}_{k}^{\star}\rVert_{\infty}\leq
3\lVert\mathsf{C}\rVert_{\infty}/\eta$. In the end, since for every
$(\mathsf{v}_{k})_{k\in[m]}\in\mathbb{R}^{n_{1}}\times\dots\times\mathbb{R}^{n_{m}}$

$\bigg{\lVert}\bigoplus_{k=1}^{m}\mathsf{v}_{k}\bigg{\rVert}_{\infty}=\max_{j\in\mathcal{J}}\bigg{\lvert}\sum_{k=1}^{m}v_{k,j_{k}}\bigg{\rvert}\leq\max_{j\in\mathcal{J}}\sum_{k=1}^{m}\lvert
v_{k,j_{k}}\rvert=\sum_{k=1}^{m}\lVert\mathsf{v}_{k}\rVert_{\infty}\leq
m\max_{k\in[m]}\lVert\mathsf{v}_{k}\rVert_{\infty},$ (C.44)

we can satisfy the boundedness assumptions on the dual variables of Theorem
4.1 with $M=6\lVert\mathsf{C}\rVert_{\infty}/\eta$ and (4.5) follows from
(4.1). Concerning the iteration complexity, again by (C.43), since
$\max_{k\in[m]}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\leq\max_{k\in[m]}\lVert\mathsf{v}^{\star}_{k}\rVert_{\infty}$,
we have
$\sum_{k\in[m]}\lVert\mathsf{v}^{t}_{k}-\mathsf{v}^{\star}_{k}\rVert_{\infty}\leq
3\lVert\mathsf{C}\rVert_{\infty}/\eta$. So, using $\eta\geq\varepsilon$, (4.8)
follows directly from (4.2) with
$M_{2}=3\lVert\mathsf{C}\rVert_{\infty}/\eta$. ∎

See 4.5

###### Proof.

Using the same notation as in the previous proof, we first show that

$(\forall\,t\in\mathbb{N})(\forall\,k\in[m])(\forall\,j_{k},\ell_{k}\in[n_{k}])\quad
v^{t}_{k,j_{k}}-v^{t}_{k,\ell_{k}}\leq 2\lVert\mathsf{C}\rVert_{\infty}/\eta.$
(C.45)

Indeed, since for $t=0$, $v^{t}_{k,j_{k}}-v^{t}_{k,\ell_{k}}=0$, we proceed by
induction assuming that (C.45) holds for $t$ and proving it for $t+1$. Noting
that for every $k\in[m]$, every $j_{-k}\in\mathcal{J}_{-k}$ and every
$j_{k},\ell_{k}\in[n_{k}]$

$\mathsf{C}_{(j_{-k},j_{k})}-\mathsf{C}_{(j_{-k},\ell_{k})}\leq
2\lVert\mathsf{C}\rVert_{\infty},$

and that $L_{t}=[n_{k_{t}}]$, from (C.42) we have that
$v^{t+1}_{k_{t},j_{k_{t}}}-v^{t+1}_{k_{t},\ell_{k_{t}}}\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta$ holds for every
$j_{k_{t}},\ell_{k_{t}}\in[n_{k_{t}}]$. On the other hand for every $k\neq
k_{t}$ $v^{t+1}_{k}=v^{t}_{k}$, which using the inductive hypothesis (C.45)
yields $v^{t+1}_{k,j_{k}}-v^{t+1}_{k,\ell_{k}}\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta$. In any case (C.45) holds for $t+1$.

Next, let $t\in\mathbb{N}$ and define the normalizing constants
$\lambda_{1}^{t+1},\ldots\lambda_{m}^{t+1}\in\mathbb{R}$ as
$\lambda_{k}^{t+1}:=-\langle\mathsf{a}_{k},\mathsf{v}^{t+1}_{k}\rangle$ for
$k\neq k_{t}$, and $\lambda_{k_{t}}^{t+1}:=-\sum_{k\neq
k_{t}}\lambda_{k}^{t+1}$. Then denoting
$\mathsf{u}_{k}^{t+1}:=\mathsf{v}_{k}^{t+1}+\lambda^{t+1}_{k}$, $k\in[m]$,
since $\sum_{k\in[m]}\lambda_{k}^{t+1}=0$, we have
$\bigoplus_{k=1}^{m}\mathsf{v}_{k}^{t+1}=\bigoplus_{k=1}^{m}\mathsf{u}_{k}^{t+1}$
and hence, recalling (B.5),

$\pi^{t+1}=\exp\Big{(}-\mathsf{C}/\eta+\bigoplus_{k=1}^{m}\mathsf{u}^{t+1}_{k}\Big{)}\odot\bigotimes_{k=1}^{m}\mathsf{a}_{k},$
(C.46)

Moreover, from (C.45) we have that for every $k\neq k_{t}$ and every
$j_{k},\ell_{k}\in[n_{k}]$

$u^{t+1}_{k,j_{k}}-u^{t+1}_{k,\ell_{k}}=v^{t+1}_{k,j_{k}}-v^{t+1}_{k,\ell_{k}}\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta,$

which, using $\sum_{j\in[n_{k}]}a_{k,j}=1$ and the fact that the
$\lambda_{k}^{t+1}$’s are chosen so that
$\langle\mathsf{a}_{k},\mathsf{u}_{k}^{t+1}\rangle=0$ for all $k\neq k_{t}$,
implies

$-u^{t+1}_{k,\ell_{k}}=\sum_{j_{k}\in[n_{k}]}a_{k,j_{k}}(u^{t+1}_{k,j_{k}}-u^{t+1}_{k,\ell_{k}})\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta,\quad\ell_{k}\in[n_{k}],$ (C.47)

and

$u^{t+1}_{k,j_{k}}=\sum_{\ell_{k}\in[n_{k}]}a_{k,\ell_{k}}(u^{t+1}_{k,j_{k}}-u^{t+1}_{k,\ell_{k}})\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta,\quad j_{k}\in[n_{k}].$ (C.48)

Therefore, we have obtained that
$\lVert\mathsf{u}^{t+1}_{k}\rVert_{\infty}\leq
2\lVert\mathsf{C}\rVert_{\infty}/\eta$ for $k\neq k_{t}$. On the other hand,
similar to what was done in the proof of Theorem 4.4 we can derive that, for
every $j_{k_{t}}\in L_{t}=[n_{k_{t}}]$,

$\exp(-u^{t+1}_{k_{t},j_{k_{t}}})=\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\\!\\!\\!\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}u^{t}_{k,j_{k}}\Big{)}\prod_{k\neq k_{t}}a_{k,j_{k}}.$

Since, recalling (C.47) and (C.48),

$\exp(-(2m-1)\lVert
C\rVert_{\infty}/\eta)\leq\exp\Big{(}-\mathsf{C}_{(j_{-k_{t}},j_{k_{t}})}/\eta+\sum_{k\neq
k_{t}}u^{t}_{k,j_{k}}\Big{)}\leq\exp((2m-1)\lVert C\rVert_{\infty}/\eta),$

and $\sum_{j_{-k_{t}}\in\mathcal{J}_{-k_{t}}}\prod_{k\neq
k_{t}}a_{k,j_{k}}=1$, we have

$\exp(-(2m-1)\lVert
C\rVert_{\infty}/\eta)\leq\exp(-u^{t+1}_{k_{t},j_{k_{t}}})\leq\exp((2m-1)\lVert
C\rVert_{\infty}/\eta).$

Therefore,

$\exp\big{(}\lvert
u^{t+1}_{k_{t},j_{k_{t}}}\rvert\big{)}=\max\big{\\{}\exp\big{(}u^{t+1}_{k_{t},j_{k_{t}}}\big{)},\exp\big{(}-u^{t+1}_{k_{t},j_{k_{t}}}\big{)}\big{\\}}\leq\exp((2m-1)\lVert
C\rVert_{\infty}/\eta)$ (C.49)

and hence

$\lVert\mathsf{u}^{t+1}_{k_{t}}\rVert_{\infty}\leq(2m-1)\lVert\mathsf{C}\rVert_{\infty}/\eta.$

Therefore, we have
$\sum_{k\in[m]}\lVert\mathsf{u}^{t}_{k}\rVert_{\infty}\leq(4m-3)\lVert\mathsf{C}\rVert_{\infty}/\eta=:M$
and due to (C.46) and the computation (C.44), we can use $M_{1}=M$ in Theorem
4.1 and get (4.7). Concerning iteration complexity, recalling (C.5) we have
that
$\sum_{k\in[m]}\lVert\mathsf{u}_{k}-\mathsf{v}_{k}^{\star}\rVert_{\infty}\leq
2(4m-3)\lVert\mathsf{C}\rVert_{\infty}/\eta$ and hence as done in (C.34) we
have

$\delta_{t}\leq\frac{2(4m-3)\lVert\mathsf{C}\rVert_{\infty}}{\eta}d_{t}.$

Moreover, since $\lVert\pi^{t}\rVert_{1}=1$ and $\bar{b}=1$, it follows from
(C.36), (C.37) and (C.38) that

$\delta_{t}-\delta_{t+1}\geq\frac{d_{t}^{2}}{2}.$ (C.50)

Thus, the statement follows from Lemma C.3. ∎