Improved Complexity Analysis of the Sinkhorn and Greenkhorn Algorithms for Optimal Transport⁰⁰footnotetext: JL and DY contribute equally.

Optimal transport (OT) is the problem of finding the most economic way to transfer one measure to another. The optimal transport problem dates back to Gaspard Monge [1]. Kantorovich (see [2] for historical context) provides a relaxation of this problem by allowing the mass splitting in source based on the notation of coupling, which advances its applications in areas such as economics and physics. Recently, with the increasing demand of quantifying the dissimilarity for probability measures, the optimal transport cost, which can serve as a natural notion of metric between measures, has attracted a lot of attention in the fields of statistics and machine learning. For instance, as a special case of the optimal transport cost, Wasserstein distance [3] can be used as a loss function when the output of a supervised learning method is a probability measure [4]. The application of Wasserstein-1 metric in GAN can effectively mitigate the issues of vanishing gradients and mode collapse [5]. Other applications of optimal transport include shape matching [6], color transfer [7] and object detection [8].

For ease of notation, we consider the case where the target probability vectors are of the same length, denoted $\boldsymbol{a},\boldsymbol{b}\in\Delta_{n}$. Given a cost matrix $\mathbf{C}\in\mathbb{R}_{+}^{n\times n}$, the Kantorovich relaxation attempts to compute the distance between $\boldsymbol{a}$ and $\boldsymbol{b}$ by solving the following problem:

where

As a linear programming, the Kantorovich problem can be solved by the simplex method or the interior-point method. Yet, both of them suffer from huge computational complexity (overall $O(n^{3})$) [2, 9]. A computationally efficient algorithm is developed in [10] which solves the Shannon entropic regularized Kantorovich problem

by the Sinkhorn update, where $\gamma$ and $H(\mathbf{P})=-\sum_{i,j}\mathbf{P}_{i,j}(\log\mathbf{P}_{i,j}-1)$ are the regularization parameter and discrete entropy, respectively. The easily parallel and GPU friendly features of the Sinkhorn algorithm make it popular in applications, especially in machine learning [4]. The idea of the Sinkhorn update can date back to [11, 12], and there are many studies on its convergence. The linear convergence rate of it (not for optimal transport) under Hilbert projective metric is provided in [13]. Regarding optimal transport, the computational complexity of the Sinkhorn algorithm to achive an $\varepsilon$-accuracy is shown to be $O({n^{2}\|\mathbf{C}\|_{\infty}^{3}\log n}/{\varepsilon^{3}})$ [14]. The Greenkhorn algorithm is an elementwise variant of the Sinkhorn algorithm which has the same computational complexity [14]. In [15], a variant of the Sinkhorn algorithm which utilizes modified probability vectors is developed and the $O({n^{2}\|\mathbf{C}\|_{\infty}^{2}\log n}/{\varepsilon^{2}})$ computational complexity is established. The idea has been extended to the Greenkhorn algorithm with the same improved complexity [16]. However, as can be seen later, the numerical results show that the practical performances of the original methods and the modified ones are overall indistinguishable, which motivates us to think whether we can improve the computational complexity of the vanilla Sinkhorn and Greenkhorn algorithms. The answer is affirmative.

Main contributions. Inspired by the equicontinuity for the continuous KL-divergence regularized OT problem, in this paper we have identified a discrete analogue of this property and shown that the $O({n^{2}\|\mathbf{C}\|_{\infty}^{2}\log n}/{\varepsilon^{2}})$ complexity of the vanilla Sinkhorn algorithm can be achieved based on this property. In addition, it is shown that this property can also be used to improve the complexity of the vanilla Greenkhorn algorithm to $O({n^{2}\|\mathbf{C}\|_{\infty}^{2}\log n}/{\varepsilon^{2}})$.

Notation. Let $\mathbb{R}_{+}$ and $\mathbb{R}_{++}$ denote non-negative real numbers and positive real numbers respectively. For a vector $\boldsymbol{x}\in\mathbb{R}^{n}$, we denote by $\|\boldsymbol{x}\|_{1}:=\sum_{i=1}^{n}|\boldsymbol{x}_{i}|$ the sum of absolute values its elements, and by $\|\boldsymbol{x}\|_{\infty}$ the maximum absolute value of its elements. Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, we write $\|\mathbf{A}\|_{1}:=\sum_{ij}|\mathbf{A}_{ij}|$ and $\|\mathbf{A}\|_{\infty}:=\max_{ij}|\mathbf{A}|_{ij}$. Let $\Delta_{n}:=\{\boldsymbol{x}\in\mathbb{R}^{n}_{+}:\|\boldsymbol{x}\|_{1}=1\}$ be the probability simplex in $\mathbb{R}^{n}$. For a matrix $\mathbf{A}$ and a vector $\boldsymbol{x}$, we denote by $e^{\mathbf{A}}$, $e^{\boldsymbol{x}}$, $\log\mathbf{A}$, $\log\boldsymbol{x}$ their entrywise exponentials and natural logarithms respectively. For vectors $\boldsymbol{x},\boldsymbol{y}$ of the same dimension, their product $\boldsymbol{x}\odot\boldsymbol{y}$ and division $\frac{\boldsymbol{x}}{\boldsymbol{y}}$ should be interpreted in an entrywise way. For two probability vectors $\boldsymbol{a},\boldsymbol{b}\in\Delta_{n}$, KL-divergence is defined as

The Pinsker’s inequality provides a lower bound of $KL(\boldsymbol{a}\|\boldsymbol{b})$ in terms of the $\ell_{1}$-norm:

Without loss of generality, assume $\boldsymbol{a}>0$ and $\boldsymbol{b}>0$¹¹1 The case where there are zeros in $\boldsymbol{a}$ or $\boldsymbol{b}$ is discussed in Remark 2.13. . Then the solution to the entropic regularized Kantorovich problem (2) is the unique matrix of the form

that satisfies

Here, $\mathbf{K}=\exp(-\frac{\mathbf{C}}{\gamma})$, $(\boldsymbol{u},\boldsymbol{v})\in\mathbb{R}_{++}^{n}\times\mathbb{R}_{++}^{n}$ are two (unknown) scaling vectors, and (4) is a reformulation of the constraints $\mathbf{P}\boldsymbol{1}_{n}=\boldsymbol{a},\mathbf{P}^{\mathsf{T}}\boldsymbol{1}_{n}=\boldsymbol{b}$ in terms of $\boldsymbol{u}$ and $\boldsymbol{v}$. For the details of derivation, see [2]. Based on this optimality condition, the Sinkhorn algorithm with arbitrary positive initialization vectors $\boldsymbol{u}_{0}$ and $\boldsymbol{v}_{0}$ conducts alternating projection to satisfy (4), see Algorithm 1 for a description of the algorithm. As already mentioned, the best known complexity for the vanilla Sinkhorn algorithm is $O({n^{2}\|\mathbf{C}\|_{\infty}^{3}\log n}/{\varepsilon^{3}})$ [14], while this complexity can be improved to $O({n^{2}\|\mathbf{C}\|_{\infty}^{2}\log n}/{\varepsilon^{2}})$ for the modified variant which uses lifted values of $\boldsymbol{a}$ and $\boldsymbol{b}$ [15] (more precisely, run Sinkhorn using $(1-\frac{\delta}{8})((\boldsymbol{a},\boldsymbol{b})+\frac{\delta}{n(8-\delta)}(\boldsymbol{1},\boldsymbol{1}))$ for a small mismatch $\delta$ instead of ($\boldsymbol{a}$, $\boldsymbol{b}$)). However, the numerical results in Figure 1 demonstrate that the practical performances of the vanilla and modified Sinkhorn algorithms are overall similar, which motivates us to improve the convergence analysis of the vanilla Sinkhorn algorithm. The numerical experiments are conducted on two data sets: the MNIST data set and a widely used synthetic data set[14]. In the experiments, a pair of images are randomly sampled from each data set and $\boldsymbol{a}$ and $\boldsymbol{b}$ are obtained by vectorizing the images, followed by normalization. The transport cost is computed as the $\ell_{2}$-distance of the pixel positions. The plots in Figure 1 are average results over 10 random instances.

There are several different equivalent interpretations of the Sinkhorn algorithm [2], and the convergence analysis is based on the dual perspective. First note that the dual problem of the entropic regularized Kantorovich problem (2) is given by

For the details of derivation, see for example [10]. If we define

it can be easily verified that

Moreover, when $k\geq 2$,

Therefore, in terms of $(\mathbf{f}_{k},\mathbf{g}_{k})$, the Sinkhorn algorithm can be interpreted as a block coordinate ascent on the dual problem (5), and the analysis will be primarily based on $(\mathbf{f}_{k},\mathbf{g}_{k})$.

After computing an approximate transport plan matrix $\mathbf{P}_{k}$ which satisfies $\|\mathbf{P}_{k}\boldsymbol{1}_{n}-\boldsymbol{a}\|_{1}+\|\mathbf{P}_{k}^{\mathsf{T}}\boldsymbol{1}_{n}-\boldsymbol{b}\|_{1}\leq\delta$ by the Sinkhorn algorithm, a rounding step presented in Algorithm 2 usually follows to ensure that the output matrix $\widehat{\mathbf{P}}_{k}$ satisfies $\widehat{\mathbf{P}}_{k}\boldsymbol{1}_{n}=\boldsymbol{a}$ and $\widehat{\mathbf{P}}_{k}^{\mathsf{T}}\boldsymbol{1}_{n}=\boldsymbol{b}$ simultaneously. The goal in the analysis is to study the computational complexity of the algorithm when an $\varepsilon$-approximate solution of the original Kantorovich problem (1) is obtained, that is when $\widehat{\mathbf{P}}_{k}$ satisfies

where $\mathbf{P}^{*}$ is the solution to (1).

To carry out the analysis, we first list three lemmas from [14]. For readers’ convenience, we have included the proofs in the appendix. The first lemma concerns the error bound of the rounding scheme, which is presented in a general form where $\boldsymbol{a}$ and $\boldsymbol{b}$ are not necessarily probability vectors. This will be useful in the analysis of the Greenkhorn algorithm in the next section.

The second lemma provides a bound for the error between the true OT distance $\langle\mathbf{C},\mathbf{P}^{*}\rangle$ and the approximate OT distance $\langle\mathbf{C},\widehat{\mathbf{P}}_{k}\rangle$.

The third lemma characterizes the one-iteration improvement of the Sinkhorn algorithm in terms of the function value.

The improved convergence analysis of the Sinkhorn algorithm relies heavily on the equicontinuity of the discrete regularized optimal transport problem. To state this property, we first follow the convention in the optimal transport literature[17] and define the $(\mathbf{C},\gamma)$-transform and $\overline{(\mathbf{C},\gamma)}$-transform of vectors.

The definitions of $(\mathbf{C},\gamma)$-transform and $\overline{(\mathbf{C},\gamma)}$-transform enable us to describe the Sinkhorn update as

Furthermore, we have

where $(\mathbf{f}_{*}^{\gamma},\mathbf{g}_{*}^{\gamma})$ is any solution of (5). Next we present the key property used in the analysis, which is inspired by the equicontinuity of the $(c,\gamma)$-transform and $\overline{(c,\gamma)}$-transform for the continuous KL-divergence regularized optimal transport problem.

Recall that in the Sinkhorn algorithm $\mathbf{f}_{k}=\mathbf{g}_{k}^{\overline{(\mathbf{C},\gamma)}}$ when $k\geq 2$ is odd and $\mathbf{g}_{k}=\mathbf{f}_{k}^{(\mathbf{C},\gamma)}$ when $k\geq 2$ is even. Moreover, since $\mathbf{f}_{k}$ = $\mathbf{f}_{k-1}$ if $k\geq 2$ is even and $\mathbf{g}_{k}$ = $\mathbf{g}_{k-1}$ if $k\geq 2$ is odd, the following corollary is a straightforward consequence of Lemma 2.6.

The following lemma provides a bound for the gap between $h(\mathbf{f}_{k},\mathbf{g}_{k})$ and the optimal value $h(\mathbf{f}_{*}^{\gamma},\mathbf{g}_{*}^{\gamma})$ utilizing Corollary 2.9.