Greedy stochastic algorithms for entropy-regularized optimal transport problems

TL;DR

This paper introduces a family of fast stochastic algorithms for entropy-regularized optimal transport problems, extending Greenkhorn, with proven convergence rates and numerical evidence of improved speed under certain conditions.

Contribution

It develops a new family of stochastic algorithms for entropy-regularized OT, extending Greenkhorn, with a general convergence theorem and practical speed improvements.

Findings

Algorithms achieve convergence rates matching existing methods.

Numerical experiments demonstrate faster performance in specific regimes.

The family generalizes Greenkhorn and includes it as a special case.

Abstract

Optimal transport (OT) distances are finding evermore applications in machine learning and computer vision, but their wide spread use in larger-scale problems is impeded by their high computational cost. In this work we develop a family of fast and practical stochastic algorithms for solving the optimal transport problem with an entropic penalization. This work extends the recently developed Greenkhorn algorithm, in the sense that, the Greenkhorn algorithm is a limiting case of this family. We also provide a simple and general convergence theorem for all algorithms in the class, with rates that match the best known rates of Greenkorn and the Sinkhorn algorithm, and conclude with numerical experiments that show under what regime of penalization the new stochastic methods are faster than the aforementioned methods.

Figures

Equations

$$T^{\ast }$$

$$T_{\lambda }^{\ast }$$

$$D(u)AD(v)1=r\text{mbox}andD(v)A^{\top }D(u)1=c,$$

$$u^{k+1}$$

$$.v^{k+1}$$

$$U_{r,c}\stackrel{\text{def}}{=}\{T\in {\mathbb{R}}_{+}^{n\times n}|r(T)=r,c(T)=c\}.$$

$$dist(A,U_{r,c})\stackrel{\text{def}}{=}\Vert r(A)-r{\Vert }_1+\Vert c(A)-c{\Vert }_1$$

$$\rho (a,b)$$

$$d_{\rho }(u,v)$$

$${\rho }^r(M)={\left(\rho (r_i,r_i(M))\right)}_{i=\mathrm{1..}n}\in {\mathbb{R}}_{+}^n,$$

$${\rho }^c(M)={\left(\rho (c_i,c_i(M))\right)}_{i=\mathrm{1..}n}\in {\mathbb{R}}_{+}^n.$$

$$\rho (M)={\left({\rho }^r(M),{\rho }^c(M)\right)}_{i=\mathrm{1..}n}\in {\mathbb{R}}_{+}^{2n}.$$

$$\forall h\in {\mathbb{R}}_{+}^{2n}\Psi (h)={\left(\frac{g(h_k)}{{\sum }_{i=1}^{2n}g(h_i)}\right)}_{k=\mathrm{1..2}n}\in {\Delta }_{2n}.$$

$$\Psi (h)$$

$$\exists k\in \mathbb{N},k\le \frac{28n}{{ϵ}^2}\log \left(\frac{s}{\ell }\right),$$

$$f(x,y)=\sum _{i,j=1}^n{A}_{ij}{e}^{x_i+y_j}-⟨r,x⟩-⟨c,y⟩.$$

$$r(XAY)=r\text{mbox}andc(XAY)=c.$$

$$f(x^k,y^k)-f(x^{k+1},y^{k+1})=\rho (r_I,r_I(D(u^k)AD(v^k))),$$

$$f(x^k,y^k)-f(x^{k+1},y^{k+1})=\rho (c_J,c_J(D(u^k)AD(v^k))).$$

$$\mathbf{E}\left[f(x^k,y^k)\right]-\underset{x,y\in R}{\min }f(x,y)$$

$$f(0,0)-\mathbf{E}\left[f(x^k,y^k)\right]\le \log \left(\frac{s}{l}\right).$$

$$f(x^1,y^1)-\underset{x,y\in R}{\min }f(x,y)\le f(0,0)-\underset{x,y\in R}{\min }f(x,y)\le \log \left(\frac{s}{l}\right)$$

$$f(x^k,y^k)-\min f\le f(0,0)-\min f\le \log \left(\frac{s}{l}\right).$$

$$-\mathbf{E}\left[f(x^k,y^k)\right]$$

$$f(0,0)-\mathbf{E}\left[f(x^k,y^k)\right]$$

$$\frac{1}{n}\sum _{i=1}^n{a}_i{b}_i\ge \left(\frac{1}{n}\sum _{j=1}^n{a}_j\right)\left(\frac{1}{n}\sum _{j=1}^n{b}_j\right).$$

$$\sum _{i,j=1}^n(a_i-a_j)(b_i-b_j)$$

$$=\sum _{j=1}^n\left(\sum _{i=1}^n\left(a_i{b}_i-a_i{b}_j-a_j{a}_i\right)+n{a}_j{b}_j\right)$$

$$=2n\sum _{i=1}^n{a}_i{b}_i-2\sum _{i=1}^n{a}_i\sum _{j=1}^n{b}_j.$$

$$\Vert u-v{\Vert }_1\le \sqrt{7\Vert u{\Vert }_1{d}_{\rho }(u,v)},$$

$$||v-u|{|}_1\le \sqrt{7d_{\rho }(v,u)}$$

$$\rho (tx,ty)$$

Full text

1 Introduction

Probability distributions are the backbone of machine learning and statistics: we use them to represent a variety of objects in learning tasks, ranging from statistical models to data representations. Comparing different distributions is often done using information divergences such as Kullback-Leiber divergence, yet this discards much of structural and geometric information present in the distribution. Developing a practical measure that captures the geometry of the probability distribution is a problem to which the optimal transport (OT) distance offers an attractive solution.

First formulated by [9] then revisited by [8], OT distances have the inherent particularity of capturing the geometrical properties of the probability measures. However, they have a drawback: computing an OT distance has a typical cost of the order $O(n^{3}\log n)$ for histograms of $n$ points ([11]). This prevents the application of OT distances in large-scale machine learning problems.

The idea of entropy penalization, proposed by [3], represents a key milestone in this field. The benefits of such a regularization scheme are multiple: the regularized problem has a unique solution, greater computational stability, and can be solved efficiently using the Sinkhorn algorithm. This new family of distances has been used in a wide range of applications, such as image classification ([3]), unsupervised learning using Restricted Boltzmann Machines ([10]), learning with a Wasserstein Loss ( [5]), domain adaptation ([2]), computer graphics ([12]), and neuroimaging  ([7]). The growing interest in applications for the Sinkhorn distances has sparked the development of new and efficient algorithms for its calculation, such as stochastic gradient based algorithms by [6], and fast methods to compute Wasserstein barycenters ([4]). To this end, [1] have developed the Greenkhorn algorithm, a greedy variant of the Sinkhorn algorithm that selects columns and rows to be updated that most violate the constraints. The authors present both promising numerical results, besting the Sinkhorn algorithm, and an insightful theoretical complexity that is linear in $n$.

Our contribution: We expand on the idea of greedy column and row selection by proposing a family of algorithms that assign a probability of updating each row and column. Moreover, our family allows for any sampling so long as the probabilities are proportional to the violation of each column or row with respect to the transport polytope. We call our algorithm the Greedy Stochastic Sinkhorn. We explain the idea behind this family of methods, show how Greenkhorn is a limiting case, and propose several other instances of the algorithm. We develop an all encompassing convergence theorem that recovers the best known $O(n/\epsilon^{2})$ iteration complexity for the Greenkhorn algorithm. Finally, we exhibit some numerical experiments that explicit the relevance of Greedy Stochastic Sinkhorn in a particular regime of penalization, along with a discussion around the computational properties of the algorithms.

1.1 The Optimal Transport problem

The discrete OT problem can be seen as a problem of optimal resource allocation given by a linear program

[TABLE]

where $r,c\in\Delta_{n}\overset{\text{def}}{=}\{x\in\mathbb{R}^{n}\,|\,\sum_{i=1}^{n}x_{i}=1\}$ are respectively the initial and target distributions, $C\in\mathbb{R}_{+}^{n\times n}$ the transport cost matrix and ${\bf 1}$ is a vector of all ones of an appropriate dimension. Matrices $T\in\mathbb{R}^{n\times n}_{+}$ that satisfy the transport constraints in (1) represent valid transportation maps between $r$ and $c$, where $T_{ij}$ will represent the mass transported from $r_{i}$ to $c_{j}$. The matrix $T^{*}$ is a transportation map that minimizes the transportation cost, the computed minimum $\left<T^{*},C\right>$ is the optimal transport value and it defines a distance between $r$ and $c$  ([13]). The transportation map $T^{*}$ can be computed using the network simplex or interior point methods ([11]), but the computational cost is in both cases $O(n^{3}\log(n))$. It is this cubic cost in the dimension that makes this notion of distance infeasible in high-dimensional settings, such as in computer vision or high dimensional inference.

1.2 Entropic regularization and the Sinkhorn algorithm

An interesting approach to alleviate the computational burden was proposed by  [3] through the introduction of an entropic regularization as follows

[TABLE]

where the entropy is $E(T)=\sum_{i,j=1}^{n}-T_{ij}\log(T_{ij})$. Due to the strong convexity introduced by the entropic regularization, the problem (2) now has a unique solution. What is more, using duality theory (2) has a smooth and unconstrained dual formulation. Leveraging on the dual [3] showed that (2) can be equivalently re-written as the following matrix scaling problem: find $u,v\in\mathbb{R}^{n}_{+}$ such that

[TABLE]

where $A=e^{-\lambda C}$ with the exponential taken element-wise and $D(u)$ denotes a diagonal matrix with the elements of $u$ on the diagonal. With the $(u,v)$ solution to (3), the solution to (2) is simply given by $T_{\lambda}^{*}=D(u)AD(v).$ This matrix scaling problem can now be efficiently solved using the celebrated Sinkhorn algorithm, as proposed by [3].

The Sinkhorn algorithm is a fixed point iteration algorithm for solving (3) which alternately scales the row and column sums to match the desired marginals

[TABLE]

where we have used $x./y$ to denote elementwise division of vectors333In other words $x./y=D(y)^{-1}x$. On top of being a simple and fast algorithm, the Sinkhorn algorithm is also GPU-friendly since its highest cost is a matrix vector product which can be parallelized. The resulting distance $\left<T_{\lambda}^{*},C\right>$ defined by (2) has been dubbed the Sinkhorn distance.

Notation: For the sake of brevity, we use $r(T)=T{\bf 1}$ and $c(T)=T^{\top}{\bf 1}$ to denote the row sum and column sum vectors of $T$, respectively. Let $U_{r,c}$ be the transport polytope defined by

[TABLE]

Since we need to solve (3), in order to discuss convergence results, we need to define a distance that measures how far are the scaled iterates from the transport polytope $U_{r,c}$. We will use in all the following work the $\ell_{1}$ distance

[TABLE]

which, as argued by [1], is much more suitable to compare probability distributions than the $\ell_{2}$ distance. A simple example to see this: taking $p=(\frac{1}{n},...\frac{1}{n},0....0)\in\Delta_{2n}$ and $q=(0,0....0,\frac{1}{n},...\frac{1}{n})\in\Delta_{2n}$ as two probability distributions with disjoint supports, we see that $\lVert p-q\rVert_{1}=2$, while $\lVert p-q\rVert_{2}=\frac{1}{\sqrt{n}}$ and thus decreases as $n$ increases, despite being clearly distinct distributions for all $n$.

1.3 Greedy Sinkhorn: Greenkhorn algorithm

Greenkhorn is a greedy version of Sinkhorn proposed by [1] where at each iteration only one coordinate of $u$ or $v$ is updated in (4), picking each time the one with highest violation with respect to the corresponding marginal. These violations are computed with the following function

[TABLE]

For vectors in the simplex $u,v\in\Delta_{n}$, we have that $d_{\rho}(u,v)$ coincides with the Kullback-Leiber divergence. For this reason, $d_{\rho}$ is known as the generalized Kullback-Leiber divergence. It is not a distance because it is not symmetric, but it verifies $d_{\rho}\geq 0$ and $d_{\rho}(u,v)=0\Leftrightarrow u=v$. Therefore, if $u,v$ are two vectors of positive entries, $d_{\rho}(u,v)$ will return some measurement on how far they are from each other. Let $\rho^{r}(M)$ (resp. $\rho^{c}(M)$) be the vector of the row sum violations (resp. column sum violations) of a given matrix $M\in\mathbb{R}^{n\times n}_{+}$ with respect to $r$ (resp. $c$), where the violations are computed using $\rho$, that is

[TABLE]

We will refer to the concatenation of these two vectors as the marginal violations denoted by

[TABLE]

The marginal violations vector $\rho(M)$ measures how far the matrix $M$ is from the transport polytope $U_{r,c}$ in the sense that $M\in U_{r,c}$ if and only if all entries of $\rho(M)$ are equal to zero.

The Greenkhorn algorithm uses $\rho(M)$ to guide the choice of which row or column should be updated. As the name indicates, the algorithm chooses the row or column index greedily, that is the index of maximal value in $\rho(M)$, see Algorithm 1. This greedy variation of Sinkhorn is expected to perform better in practice, mainly because it does not update rows or columns that already match the correct marginal sum value.

[1] proved that, to reach an $\epsilon>0$ approximate solution, the Greenkhorn algorithm and the Sinkhorn algorithm converge in at most $28n\epsilon^{-2}\log(\frac{s}{l})$ and $28\epsilon^{-2}\log(\frac{s}{l})$ iterations, respectively, where $s=\lVert A\rVert_{1}$ is the total mass and $l$ the smallest entry of the matrix $A$. [1] also claimed that the Greenkhorn algorithm can be implemented in such a way that the iteration cost is linear in $n$, consequently the overall complexity of either the Sinkhorn algorithm or Greenkhorn is quadratic in $n$, which is stark contrast to the cubic dependency of the interior point type methods ([11]). Since the authors omitted the details on how such a linear iteration complexity can be achieved, we have given the details in Section 4.1. The $\epsilon^{-2}$ dependency of Greenkhorn and Sinkhorn is also in contrast with logarithmic dependency on $\epsilon$ in interior point based methods. Thus Greenkhorn and Sinkhorn are well suited for the large dimensional setting where we can tolerate an approximate solution. This is typically the case in the problems that we are interested in here, such as problems that arise in large dimensional machine learning.

2 Greedy Stochastic Sinkhorn

While the greedy strategy in the Greenkhorn algorithm is, in some sense, optimal for one step, it may not be the best strategy over a number of iterations. Here we introduce a more flexible, and less aggressive updating strategy.

At each iteration of the Greedy Stochastic Sinkhorn algorithm, instead of picking the column or row with the highest violation, as is done in the Greenkhorn algorithm, we will assign to each row and column a probability of being updated. Because we want the columns and rows with highest violation to be updated more frequently, we assign a higher probability to columns and rows with a higher violation. We do this using an increasing probability function.

Definition 1

We say that $\Psi$ is a increasing probability function if there exists an increasing positive function $g:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$ such that

[TABLE]

Several examples of an increasing probability function are given as follows

[TABLE]

where $T,\alpha>0$ are parameters. If $\rho$ is our current vector of violations, then $\Psi(\rho)=p\in\Delta_{2n}$ defines a probability distribution. Furthermore, since $\Psi$ is built on top of an increasing function, a larger violation $\rho_{i}$ will result in a larger probability $p_{i}$. See Algorithm 2 for the pseudocode of this family of stochastic algorithms. In the next section we prove that Algorithm 2 converges for any increasing probability function. This is particularly interesting when we consider that the Greenkhorn algorithm is a limiting case of the Greedy Stochastic Sinkhorn. Indeed, the selection criteria of the Greenhkorn algorithm corresponds to taking the limit over $\alpha\rightarrow\infty$ of the probability function (11).

3 Convergence analysis

We now present our main convergence theorem, discuss its consequences and proof.

Theorem 2

Consider the sequence of matrices $A^{k}\overset{\text{def}}{=}D(u^{k})AD(v^{k})$ produced by Algorithm 2 with an increasing probability function $\Psi$ as defined in (9). Then for a given $\epsilon>0,$ we have that

[TABLE]

such that $\mathbf{E}\left[dist(A^{k},U_{r,c})\right]\leq\epsilon.$

We make several interesting remarks on the consequence of this theorem.

Since $dist(A^{k},U_{r,c})$ is a positive random variable, by Markov’s inequality we have that the convergence in expectation given in Theorem 2 also proves that $dist(A^{k},U_{r,c})$ converges in probability to zero. The variance also converges to zero at a $O(\left.n\right/\epsilon)$ rate, as proven in Section 3.3. 2. 2.

The convergence rate given in Theorem 2 is exactly the same rate as given by [1] for the Greenkorn algorithm. 3. 3.

Remarkably the rate of convergence does not depend on the choice of probability function $\Psi$. Thus, in theory, a uniform selection of the coordinates gives the same asymptotic convergence as the Greenkhorn selection criteria.

Before moving onto the proof, we need several auxiliary lemmas.

3.1 Useful lemmas

Our analysis is based on the dual objective of (1) given by

[TABLE]

Let $X=D(e^{x})$ and $Y=D(e^{y})$. By writing out the first order optimality conditions of $f(x,y)$ we arrive at

[TABLE]

That is, the row sum and column sum of $XAY$ is $r$ and $c$, respectively. By denoting $u=e^{x}$ and $v=e^{y}$ the given scaling vectors, we see that (15) is the matrix scaling problem (3). Throughout this section we use $(u^{k},v^{k})$ to denote the $(u,v)$ vectors of Algorithm 2 after completing the $k$th iteration. We also denote $(x^{k},y^{k})=(\log(u^{k}),\log(v^{k})).$

The proof of Theorem 2 is based on the four next lemmas.

The first lemma was presented by [1], and it links the evolution of dual objective (14) and the marginal violations.

Lemma 3

For a given $k$, if $(u^{k+1},v^{k+1})$ was obtained by updating coordinate $I$ of $u^{k}$ then the following identity holds

[TABLE]

and if they were obtained by updating coordinate J of $v^{k}$ then

[TABLE]

Since $\rho\geq 0$ then the sequence of real numbers $(f(x^{k},y^{k}))_{k\in\textbf{N}}$ is decreasing.

Next we have an extension to the stochastic setting of another lemma by [1]. It links the expectation of the dual objective value to a type of condition number of the $A$ matrix.

Lemma 4

Let $((u^{k},v^{k}))_{k\in\textbf{N}}$ and the associated $((x^{k},y^{k}))_{k\in\textbf{N}}$ be a sequence of scaling vectors produced by the Greedy Stochastic Sinkhorn Algorithm 2. Then the following inequalities hold

[TABLE]

where $l=\min_{i,j}|A_{ij}|$ and $s=\lVert A\rVert_{1}$. As a direct consequence, we also have

[TABLE]

To prove Lemma 4 we will use this result from [1].

Proposition 5

Let $x^{0},y^{0}=0$ the initial points and $x^{1},y^{1}$ the points resulting from updating a single coordinate $u^{k}$ or a single coordinate of $v^{k}$. Then the following inequality holds

[TABLE]

where $l$ is the smallest entry of $A$ and $s=\lVert A\rVert_{1}$.

Now we write the proof of Lemma 4 of our paper.

**Proof: ** We have from Lemma 3 that the sequence of real numbers $(f(x^{k},y^{k}))_{k\in\textbf{N}}$ is decreasing for all sequences of updated indexes, in particular in the stochastic setting. Therefore the inequality of the previous proposition in fact holds for all iterations

[TABLE]

Taking expectation in the previous inequalities gives (16). To prove (17) we have that:

[TABLE]

which concludes the proof.

The next lemma is a useful inequality on ordered series of real number.

Lemma 6** (Chebyshev inequality)**

Let $a_{1}\leq a_{2}\leq\ldots\leq a_{n}\in\mathbb{R}$ and $b_{1}\leq b_{2}\leq\ldots\leq b_{n}$ be two ordered sequences. It follows that

[TABLE]

**Proof: ** Simply note the identity

[TABLE]

The proof now follows by noting that $(a_{i}-a_{j})(b_{i}-b_{j})$ is positive for every $i,j=1,\ldots,n.$.

Finally we have a lemma that is a generalization of the Pinsker inequality.

Lemma 7

The following generalized Pinsker inequality holds for $v,u\in R^{n}_{+}$

[TABLE]

where $d_{\rho}(u,v)$ is defined as in (7).

Lemma 7 is a generalization of the next proposition from [1].

Proposition 8

For $v\in\Delta_{n}$ and $u\in R^{n}_{+}$ we have that

[TABLE]

Using Proposition 8 we will now prove Lemma 7.

**Proof: ** First we note that $\rho$ is a $1-$homogeneous function, that is

[TABLE]

This implies immediatly that $d_{\rho}$ is $1-$homogeneous as well. Therefore for $v,u\in R^{n}_{+}$ and we have that

[TABLE]

Furthermore

[TABLE]

Consequently

[TABLE]

which concludes the proof.

3.2 Proof of Theorem 2

**Proof: ** Let $D_{k}\overset{\text{def}}{=}\mathbf{E}\left[dist(A^{k},U_{r,c})\right]$ and let $k^{*}\in\mathbb{N}$ be an integer such that $D_{k}>\epsilon$ for all $k<k^{*}$ (in other terms, an index such that the algorithm has not converged yet at the corresponding iteration).

Recall that $\rho(A^{k})$ is the vector of all $2n$ marginal violations for the matrix $A^{k}$, as defined in (8). We will write its components as $\rho_{i}(A^{k})$ for a given index $i$. Recall that $\Psi(\rho(A^{k}))$ is the vector of probabilities of picking each row and column, and similarly we will write its components $\Psi_{i}(\rho(A^{k}))$, which is then the probability of picking index $i$. We start the proof by showing that $D_{k}^{2}$ is upper bounded by the following conditional expectation

[TABLE]

where $I$ is the index randomly sampled at iteration $k$. Let $k<k^{*}$ and since we assume that (9) holds for some function $g$, we have that

[TABLE]

where we applied Lemma 6 in the first inequality which relies on the monotonicity of $g$, and the generalized Pinsker inequality (19) in the second inequality with $a=(r,c)$, $b=(r(A^{k}),c(A^{k}))$ and we used that $\lVert a\rVert_{1}=\lVert r\rVert_{1}+\lVert c\rVert_{1}=2.$ Taking expectation in (21), using the law of total expectation and the fact that $\mathbf{E}\left[X^{2}\right]\geq\mathbf{E}\left[X\right]^{2}$ for any random variable $X$ gives

[TABLE]

To conclude, we now use Lemma 3 to re-write $\mathbf{E}\left[\rho_{I}(A^{k})\mid A^{k}\right]$ as

[TABLE]

Thus taking expectation in (23) gives

[TABLE]

Summing over $k=0,\ldots,k^{*}-1$ in (24) and using telescopic cancellation we have that

[TABLE]

Combining the above with

[TABLE]

as proven in Lemma 4, we have that

[TABLE]

This proves that for a given integer $k^{*}$

[TABLE]

The contrapositive of the above statement is given by

[TABLE]

Choosing $k^{*}=\lceil\frac{28n}{\epsilon^{2}}log\left(\frac{s}{\ell}\right)\rceil$ concludes the proof.

Corollary 9

If we choose $\Psi$ as either (10), (11) or (12) then the Greedy Stochastic Sinkhorn algorithm converges at a $O(\frac{n}{\epsilon^{2}})$ rate according to Theorem 2.

**Proof: ** The proof follows by observing that (10), (11) or (12) are increasing probability functions. That is, the functions $x\mapsto e^{x/T}$, $x\mapsto x^{\alpha}$ and $x\mapsto 1$ are positive increasing real-valued functions.

3.3 Variance convergence

In this section we prove that the variance of the Greedy Stochastic Sinkhorn algorithm also converges to zero, and furthermore, we show that the rate of convergence is $O(n/\epsilon)$ to reach a variance of $\epsilon>0.$ Having already established the convergence of the Greedy Stochastic Sinkhorn algorithm in expectation and high probability, together with the following convergence of the variance, we can conclude that the convergence is qualitatively almost deterministic.

Theorem 10

Consider the sequence of matrices $A^{k}\overset{\text{def}}{=}D(u^{k})AD(v^{k})$ produced by Stochastin Sinkhorn Algorithm with an increasing probability function $\Psi$. Then for a given $\epsilon>0,$ we have that

[TABLE]

such that $\mathrm{Var}(dist(A^{k},U_{r,c}))\leq\epsilon.$

**Proof: ** First we note that

[TABLE]

Therefore we will prove that $\mathbf{E}\left[dist(A^{k},U_{r,c})^{2}\right]$ converges at a certain rate, which will give the same thing for the variance.

The proof is extremely similar to the proof of Theorem 2, therefore we will only highlight where we have slight differences. Now we call $D^{\prime}_{k}=\mathbf{E}\left[dist(A^{k},U_{r,c})^{2}\right]$ and consider an integer $k^{*}$ such that $\forall k<k^{*},D^{\prime}_{k}>\epsilon$ where $\epsilon>0$ is fixed.

The only modification of the proof is that now we will have

[TABLE]

Note that in the proof of Theorem 2 we had a $\epsilon^{2}$ in place of $\epsilon$ in the above lower bound. Taking this minor difference into account, the rest of the proof follows verbatim to the proof of Theorem 2.

4 Numerical experiments

In this section we provide some empirical insights into the behaviour of the Greedy Stochastic Sinkhorn Algorithm 2. We consider both real and synthetic datasets: MNIST digits and random histograms. The authors of both [3] and [1] provided numerical experiments where Sinkhorn and Greenkhorn perform considerably better than other Optimal Transport algorithms. We will show how the Greedy Stochastic Sinkhorn has a similar efficiency for monotonic probability functions (9). In particular, we will show that for (11) with $\alpha=1$, Greedy Stochastic Sinkhorn outperforms Greenkhorn in the short-term for regimes of small penalization. Finally we will discuss the computational properties of the three algorithms, in particular some drawbacks of Greenkhorn and Greedy Stochastic Sinkhorn in comparaison with Sinkhorn, and give insights on how to bypass them. But first, we explicity show how Greenkhorn and Greedy Stochastic Sinkorn have linear iteration complexities.

4.1 Updating the marginal violation

The Greenkhorn Algorithm 1 and the Greedy Stochastic Sinkhorn Algorithm 2 must re-compute marginal violations

[TABLE]

at each iteration. Calculating $\rho(A^{k})$ from scratch at each iteration would cost $O(n^{2})$, which would defeat the purpose of both algorithms of having a linear iteration complexity. Fortunately $\rho(A^{k})$ can be updated on the fly with only $O(n)$ operations. To see this, suppose we have stored $c(A^{k}),r(A^{k})$ and $\rho(A^{k})$ and now we wish to calculate $\rho(A^{k+1}).$ Suppose we sample an index $I$ in Algorithm 2 such that $I\in\{1,\ldots,n\}$, consequently we update

[TABLE]

while $v^{k+1}=v^{k}$ and $u_{i}^{k+1}=u^{k}_{i}$ for $i\neq I$ remain unaltered. We can thus calculate the $i$th component of $r(A^{k+1})$ via

[TABLE]

The column sum vector can be updated using

[TABLE]

Thus both $r(A^{k+1})$ and $c(A^{k+1})$ can be updated using $O(n)$ operations. Since $\rho(r_{i},r_{i}(A^{k+1}))=\rho(r_{i},r_{i}(A^{k}))$ for $i\neq I$, only $n+1$ components of $\rho(A^{k+1})$ need to re-computed, which costs $O(n)$ operations. The $O(n)$ cost of the case where $I\in\{n+1,\ldots 2n\}$ can be deduced in an analogous way.

4.2 Experiments

We perform experiments on MNIST dataset. We take pairs of elements from the $28\times 28$ pixels MNIST dataset, that we then vectorize into 1D arrays $r$ and $c$ (in the sense that, for example, the 2 by 2 matrix $((1,2),(3,4))$ becomes the vector $(1,2,3,4)$ ). The cost matrix $C$ is then constructed so that $C_{ij}$ equals the $\ell_{1}$ distance between pixels $i$ and $j$ in the $28\times 28$ grid. We then apply the Sinkhorn, Greenkhorn and Greedy Stochastic Sinkhorn algorithms to compute a diagonal scaling of $A=e^{-\lambda C}$. This process is then repeated $20$ times, for each time we randomly sample a pair of images from the MNIST dataset. Finally, we report the average performance of the algorithms over these $20$ experiments.

The choice of $\lambda$ defines the penalization, and we highlight the fact that regimes of low penalization, corresponding to higher values of $\lambda$ (2), are of particular interest since they change the least the solution of the original non-regularized problem (1). For this setting, Greedy Stochastic Sinkhorn with  (11) for $\alpha=1$ is clearly the best choice overall, see Figure 1.

We also compared in Figure 2 the Greedy Stochastic Sinkhorn for various choices of parameters. As expected, using the probability function (11) with $\alpha\rightarrow+\infty$ or (12) with $T\rightarrow 0$, the Greedy Stochastic Sinkhorn algorithm reduces to the Greenkhorn algorithm

Notice also that the standard deviation for Greedy Stochastic Sinkhorn (represented as errorbars) tends to [math], which is a very important property because of the stochastic nature of the algorithms. In fact, this means that not only the expectation of the distance tends to zero, but also the variance, a fact which we prove in the appendix.

4.3 Discussion and block algorithms

While Greedy Stochastic Sinkhorn and Greenkhorn empirically perform better than Sinkhorn, they also have two computational drawbacks: firstly they are not parallelizable. In fact, each iteration of Sinkhorn (4) is a rescaling of $u$ and $v$ involving a matrix-vector product that can be parallelized. Greenkhorn and Greedy Stochastic Sinkhorn update only one element per iteration which does not involve a similar product that we can parallelize. Secondly, the greedy algorithms compute at each iteration the marginal violations, which, despite only costing $O(n)$, it does represent an additional computational cost per iteration.

One solution to these issues is to re-compute these marginal violations only once every $d$ iterations: for example in Greenkhorn we compute marginal violations $\rho(A^{k})$ and we update not only the index of highest value, but the $d$ indexes of highest values. Something similar can be done in Greedy Stochastic Sinkhorn by sampling $d$ indexes without replacement instead of just one. By doing so, on the one hand we reduce computation time for computing the marginal violations by a factor $d$, and on the other hand the algorithms described are now parallelizable because updating $d$ components of $u$ and $v$ does involve a matrix-vector product.

This idea is motivated by the numerical results of [1] where the authors concluded that the efficiency of Greenkhorn is mainly due to the fact that is does not update rows and columns that already match desired sums, more than the fact that it updates indexes with highest marginal violations. This means that the procedure of updating $d$ indexes instead of just one is expected to have a similar efficiency, which is indeed the result we get in our own numerical experiments as shown in Figure 3.

5 Conclusion

We presented a family of stochastic algorithms for entropy-regularized OT problems. We were able to derive convergence rates for a very broad class of probability functions, along with numerical experiments where a simple and intuitive choice of probability functions performed the best. We also proposed and tested simple numerical solutions to the drawbacks of the greedy algorithms.

Acknowledgements

Both authors are indebted to Marco Cuturi for essentially teaching both of them about OT and many inspiring discussions. RMG acknowledges the support of the FSMP postdoctoral fund. BKA carried our this work while at an internship, and acknowledges CREST/ENSAE for funding, and the SIERRA team, Inria Paris, for hosting.