Overrelaxed Sinkhorn–Knopp Algorithm for Regularized Optimal Transport

Abstract

This article describes a set of methods for quickly computing the solution to the regularized optimal transport problem. It generalizes and improves upon the widely used iterative Bregman projections algorithm (or Sinkhorn–Knopp algorithm). We first proposed to rely on regularized nonlinear acceleration schemes. In practice, such approaches lead to fast algorithms, but their global convergence is not ensured. Hence, we next proposed a new algorithm with convergence guarantees. The idea is to overrelax the Bregman projection operators, allowing for faster convergence. We proposed a simple method for establishing global convergence by ensuring the decrease of a Lyapunov function at each step. An adaptive choice of the overrelaxation parameter based on the Lyapunov function was constructed. We also suggested a heuristic to choose a suitable asymptotic overrelaxation parameter, based on a local convergence analysis. Our numerical experiments showed a gain in convergence speed by an order of magnitude in certain regimes.

1. Introduction

Optimal transport (OT) is an efficient and flexible tool to compare two probability distributions, which has been popularized in the computer vision community in the context of discrete histograms [1]. The introduction of entropic regularization of the optimal transport problem in [2] has made possible the use of the fast Sinkhorn–Knopp algorithm [3], scaling with high-dimensional data. Regularized optimal transport has thus been intensively used in machine learning with applications such as geodesic PCA [4], domain adaptation [5], data fitting [6], training of a Boltzmann machine [7], or dictionary learning [8,9].

The computation of optimal transport between two data relies on the estimation of an optimal transport matrix, the entries of which represent the quantity of mass transported between data locations. Regularization of optimal transport with strictly convex regularization [2,10], nevertheless, involves a spreading of the mass. Hence, for particular purposes such as color interpolation [11] or gradient flow [12], it is necessary to consider small parameters $\epsilon$ for the entropic regularization term. The Sinkhorn–Knopp (SK) algorithm is a state-of-the-art algorithm to solve the regularized transport problem. The SK algorithm performs alternated projections, and the sequence of generated iterates converges to a solution of the regularized transport problem. Unfortunately, the lower $\epsilon$ is, the slower the SK algorithm converges. To improve the convergence rate of the SK algorithm, several acceleration strategies have been proposed in the literature, based for example on mixing or overrelaxation.

1.1. Accelerations of the Sinkhorn–Knopp Algorithm

In the literature, several accelerations of the Sinkhorn–Knopp algorithm have been proposed, using for instance greedy coordinate descent [13] or screening strategies [14]. In another line of research, the introduction of relaxation variables through heavy ball approaches [15] has recently gained popularity to speed up the convergence of algorithms optimizing convex [16] or non-convex [17,18] problems. In this context, the use of regularized nonlinear accelerations (RNAs) [19,20,21] based on Anderson mixing has led to important numerical improvements, although the global convergence is not guaranteed with such approaches, as is shown further. In this paper, we also investigated another approach related to the successive overrelaxation (SOR) algorithm [22], which is a classical way to solve linear systems. Similar schemes have been empirically considered to accelerate the SK algorithm in [9,23]. The convergence of these algorithms has nevertheless not been studied yet in the context of regularized optimal transport.

Recent progress has been made on computational complexity guarantees for the Sinkhorn–Knopp algorithm and accelerated versions [13,24,25,26]. Since the methods we discuss in this paper are based on asymptotic acceleration techniques, it is challenging to show their efficiency via global computational complexity guarantee, and we do not cover these aspects here.

1.2. Overview and Contributions

The contribution of this paper is twofold. First, the numerical efficiency of the RNA methods applied to the SK algorithm to solve the regularized transport problem is shown. Second, a new extrapolation and relaxation technique for accelerating the Sinkhorn–Knopp (SK) algorithm while ensuring convergence is given. The numerical efficiency of this new algorithm is demonstrated, and a heuristic rule is also proposed to improve the rate of the algorithm.

Section 2 is devoted to the Sinkhorn–Knopp algorithm. In Section 3, we propose to apply regularized nonlinear acceleration (RNA) schemes to the SK algorithm. We experimentally show that such methods lead to impressive accelerations for low values of the entropic regularization parameter. In order to have a globally converging method, we then propose a new overrelaxed algorithm: Sinkhorn–Knopp with successive overrelaxation (SK-SOR). In Section 4, we show the global convergence of this algorithm and analyze its local convergence rate to justify the acceleration. We finally demonstrate numerically in Section 5 the interest of our method. Larger accelerations are indeed observed for decreasing values of the entropic regularization parameter.

Remark 1.

This paper is an updated version of an unpublished work [27] presented at the NIPS 2017 Workshop on Optimal Transport & Machine Learning. In the meantime, complementary results on the global convergence of our method presented in Section 4 were provided in [28]. The authors showed the existence of a parameter ${\theta }_0$ such that both global convergence and local acceleration were ensured for overrelaxation parameters $\omega \in (1,{\theta }_0)$. This result was nevertheless theoretical, and the numerical estimation of ${\theta }_0$ is still an open question. With respect to our unpublished work [27], the current article presents an original contribution in Section 3: the application of RNA methods to accelerate the convergence of the SK algorithm.

2. Sinkhorn Algorithm

Before going into further details, we now briefly introduce the main notations and concepts used throughout this article.

2.1. Discrete Optimal Transport

We considered two discrete probability measures ${\mu }_k\in {\mathbb{R}}_{+*}^{n_k}$. Let us define the two following linear operators:

$$\begin{array}{cccc}\hfill A_1& :\left\{\begin{array}{c}{\mathbb{R}}^{n_1{n}_2}\to {\mathbb{R}}^{n_1}\hfill \\ {\left(A_{1}x\right)}_i=\sum _j{x}_{i,j}\hfill \end{array}\right.\hfill & \hfill A_2& :\left\{\begin{array}{c}{\mathbb{R}}^{n_1{n}_2}\to {\mathbb{R}}^{n_2}\hfill \\ {\left(A_{2}x\right)}_j=\sum _i{x}_{i,j},\hfill \end{array}\right.\hfill \end{array}$$

as well as the affine constraint sets:

$$\begin{array}{cc}\hfill {𝒞}_k& =\left\{\gamma \in {\mathbb{R}}^{n_1{n}_2}\mid A_k\gamma ={\mu }_k\right\}.\hfill \end{array}$$

Given a cost matrix c with nonnegative coefficients, where $c_{i,j}$ represents the cost of moving mass ${\left({\mu }_1\right)}_i$ to ${\left({\mu }_2\right)}_j$, the optimal transport problem corresponds to the estimation of an optimal transport matrix $\gamma$ solution of:

$$\underset{\gamma \in {𝒞}_1\cap {𝒞}_2\cap {\mathbb{R}}_{+}^{n_1{n}_2}}{min}\langle c,\gamma \rangle :=\sum _{i,j}{c}_{i,j}{\gamma }_{i,j}.$$

This is a linear programming problem whose resolution becomes intractable for large problems.

2.2. Regularized Optimal Transport

In [2], it was proposed to regularize this problem by adding a strictly convex entropy regularization:

$$\underset{\gamma \in {𝒞}_1\cap {𝒞}_2\cap {\mathbb{R}}_{+}^{n_1{n}_2}}{min}{K}^{\epsilon }\left(\gamma \right):=〈c,\gamma 〉+\epsilon KL(\gamma ,\mathbf{1}),$$

(1)

with $\epsilon >0$, $\mathbf{1}$ the matrix of size $n_1\times n_2$ full of ones, and the Kullback–Leibler divergence is:

$$KL(\gamma ,\zeta )=\sum _{i,j}{\gamma }_{i,j}\left(log\left(\frac{{\gamma }_{i,j}}{{\zeta }_{i,j}}\right)-1\right)+\sum _{i,j}{\zeta }_{i,j}$$

(2)

with the convention $0log0=0$. It was shown in [29] that the regularized optimal transport matrix ${\gamma }^{*}$, which is the unique minimizer of Problem (1), is the Bregman projection of ${\gamma }^0=e^{-c/\epsilon }$ (here and in the sequel, exponentiation is meant entry-wise) onto ${𝒞}_1\cap {𝒞}_2$:

$${\gamma }^{*}=\underset{{𝒞}_1\cap {𝒞}_2}{argmin}{K}^{\epsilon }\left(\gamma \right)=P_{{𝒞}_1\cap {𝒞}_2}\left(e^{-c/\epsilon }\right),$$

(3)

where $P_{𝒞}$ is the Bregman projection onto $𝒞$ defined as:

$$P_{𝒞}\left(\zeta \right):=\underset{\gamma \in 𝒞}{argmin}KL(\gamma ,\zeta ).$$

2.3. Sinkhorn–Knopp Algorithm

Iterative Bregman projections onto ${𝒞}_1$ and ${𝒞}_2$ converge to a point in the intersection ${𝒞}_1\cap {𝒞}_2$ [30]. Hence, the so-called Sinkhorn–Knopp (SK) algorithm [3] that performs alternate Bregman projections can be considered to compute the regularized transport matrix:

$$\begin{array}{ccc}\hfill {\gamma }^0& =e^{-c/\epsilon }\hfill & \hfill {\gamma }^{\ell +1}=P_{{𝒞}_2}\left(P_{{𝒞}_1}\left({\gamma }^{\ell }\right)\right),\end{array}$$

and we have ${lim}_{\ell \to +\infty }{\gamma }^{\ell }=P_{{𝒞}_1\cap {𝒞}_2}\left({\gamma }^0\right)={\gamma }^{*}.$

In the discrete setting, these projections correspond to diagonal scalings of the input:

$$\begin{array}{cccc}\hfill P_{{𝒞}_1}\left(\gamma \right)& =diag\left(u\right)\gamma \hfill & \hfill \mathrm{with}u& ={\mu }_1\oslash A_1\gamma \hfill \\ \hfill P_{{𝒞}_2}\left(\gamma \right)& =\gamma diag\left(v\right)\hfill & \hfill \mathrm{with}v& ={\mu }_2\oslash A_2\gamma \hfill \end{array}$$

(4)

where ⊘ is the pointwise division. To compute the solution numerically, one simply has to store $(u^{\ell },v^{\ell })\in {\mathbb{R}}^{n_1}\times {\mathbb{R}}^{n_2}$ and to iterate:

$$\begin{array}{cccc}\hfill u^{\ell +1}& ={\mu }_1\oslash {\gamma }^0{v}^{\ell }\hfill & \hfill v^{\ell +1}& ={\mu }_2\oslash {}^t{\gamma }^0{u}^{\ell +1}.\hfill \end{array}$$

(5)

We then have ${\gamma }^{\ell }=diag\left(u^{\ell }\right){\gamma }^{0}diag\left(v^{\ell }\right).$

Another way to interpret the SK algorithm is as an alternate maximization algorithm on the dual of the regularized optimal transport problem; see [31], Remark 4.24. The dual problem of (1) is:

$$\underset{\begin{array}{c}\alpha \in {\mathbb{R}}^{n_1}\\ \beta \in {\mathbb{R}}^{n_2}\end{array}}{max}E(\alpha ,\beta ):=\langle \alpha ,{\mu }_1\rangle +\langle \beta ,{\mu }_2\rangle -\epsilon \sum _{i,j}{e}^{({\alpha }_i+{\beta }_j-c_{i,j})/\epsilon }.$$

(6)

As the function E is concave, continuously differentiable, and admits a maximizer, the following alternate maximization algorithm converges to a global optimum:

$$\begin{array}{cc}\hfill {\alpha }^{\ell +1}& =\underset{\alpha }{argmax}E(\alpha ,{\beta }^{\ell })\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\beta }^{\ell +1}& =\underset{\beta }{argmax}E({\alpha }^{\ell +1},\beta ).\hfill \end{array}$$

(8)

The explicit solutions of these problems are:

$$\begin{array}{ccc}\hfill {\alpha }_i^{\ell +1}& =\epsilon log\left(\sum _{j}exp\left(log{\left({\mu }_1\right)}_i-\left({\beta }_j^{\ell }-c_{i,j}\right)/\epsilon \right)\right)\hfill & \hfill i=1\cdots n_1\end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\beta }_j^{\ell +1}& =\epsilon log\left(\sum _{i}exp\left(log{\left({\mu }_2\right)}_j-\left({\alpha }_i^{\ell +1}-c_{i,j}\right)/\epsilon \right)\right)\hfill & \hfill j=1\cdots n_2\end{array}$$

(10)

and we recover the SK algorithm (5) by taking $u_i=e^{{\alpha }_i/\epsilon }$, $v_j=e^{{\beta }_j/\epsilon }$, and ${\gamma }_{i,j}^0=e^{-c_{i,j}/\epsilon }$.

Efficient parallel computations can be considered [2], and one can almost reach real-time computation for large-scale problems for certain classes of cost matrices c, allowing the use of separable convolutions [32]. For low values of the parameter $\epsilon$, numerical issues can arise, and the log stabilization version of the algorithm presented in Relations (9) and (10) is necessary [12]. Above all, the linear rate of convergence degrades as $\epsilon \to 0$ (see for instance Chapter 4 in [31]). In the following sections, we introduce different numerical schemes that accelerate the convergence in the regime $\epsilon \to 0$.

3. Regularized Nonlinear Acceleration of the Sinkhorn–Knopp Algorithm

In order to accelerate the SK algorithm for low values of the regularization parameter $\epsilon$, we propose to rely on regularized nonlinear acceleration (RNA) techniques. In Section 3.1, we first introduce RNA methods. The application to SK is then detailed in Section 3.2.

3.1. Regularized Nonlinear Acceleration

To introduce RNA, we first rewrite the SK algorithm (7) and (8) as:

$${\beta }^{\ell +1}=\mathrm{SK}\left({\beta }^{\ell }\right):=\underset{\beta }{argmax}E\left(\underset{\alpha }{argmax}E(\alpha ,{\beta }^{\ell }),\beta \right).$$

(11)

The goal of this algorithm is to build a sequence ${\left({\beta }^{\ell }\right)}_{l⩾1}$ converging to a fixed point of SK, i.e., to a point ${\beta }^{*}$ satisfying:

$${\beta }^{*}=SK\left({\beta }^{*}\right).$$

(12)

Many optimization problems can be recast as fixed-point problems. The Anderson acceleration or Anderson mixing is a classical method to build a sequence ${\left(x^n\right)}_n$ that converges numerically fast to a fixed point of any operator T from ${\mathbb{R}}^N$ to ${\mathbb{R}}^N$. This method defines at each step a linear (but not necessarily convex) combination of some previous values of ${\left(x^k\right)}_k$ and ${\left(T{x}^k\right)}_k$ to provide a value of $x^n$ such that $∥x^n-T{x}^n∥$ is as low as possible.

Numerically, fast local convergence rates can be observed when the operator T is smooth. This method can nevertheless be unstable even in the favorable setting where T is affine. Such case arises for instance when minimizing a quadratic function F with the descent operator $T=I-h\mathsf{\nabla }F$, with time step $h>0$. Unfortunately, there are no convergence guarantees, in a general setting or in the case $T=SK$, that the RNA sequence ${\left(x^n\right)}_n$ converges for any starting point $x^0$.

RNA is an algorithm that can be seen as a generalization of the Anderson acceleration. It can also be applied to any fixed-point problem. The RNA method [21] applied to algorithm (11) using at each step the N previous iterates is:

$$\begin{array}{cc}\hfill {\left(s_l^{\ell +1}\right)}_{l=0}^{N-1}& =\underset{\mathbf{s};{\sum }_{l=0}^{N-1}{s}_l=1}{argmin}{\mathcal{P}}_{\lambda }\left(\mathbf{s},{({\beta }^{\ell -l}-y^{\ell -l})}_{l=0}^{N-1}\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill y^{\ell +1}& =\sum _{l=0}^{N-1}{s}_l^{\ell +1}(y^{\ell -l}+\omega ({\beta }^{\ell -l}-y^{\ell -l}))\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\beta }^{\ell +1}& =\mathrm{SK}\left(y^{\ell +1}\right),\hfill \end{array}$$

(15)

where extrapolation weights $\mathbf{s}={\left(s_l\right)}_{l=0}^{N-1}$ are the solution to the optimization problem (16) defined below, and $\omega$ is a relaxation parameter. RNA uses the memory of the past trajectory when $N>1$. We can remark that for $N=1$, $s_0^{\ell }=1$ for all $\ell ⩾0$, and RNA is reduced to a simple relaxation parameterized by $\omega$:

$$y^{\ell +1}=y^{\ell }+\omega ({\beta }^{\ell }-y^{\ell }).$$

Let us now discuss the role of the different RNA parameters $\omega$ and $\mathbf{s}$.

Relaxation: Taking origins from Richardson’s method [33], relaxation leads to numerical convergence improvements in gradient descent schemes [34]. Anderson suggested underrelaxing with $\omega \in (0;1]$, while the authors of [21] proposed to take $\omega =1$.

Extrapolation: Let us define the residual $r\left(y\right)=\mathrm{SK}\left(y\right)-y$. As the objective is to estimate the fixed point of SK, the extrapolation step builds a vector y such that $\left|\right|r\left(y\right)\left|\right|$ is minimal. A relevant guess of such y is obtained by looking at a linear combination of previous iterates that reaches this minimum. More precisely, RNA methods estimate the weight vector ${\left(s_l^{\ell +1}\right)}_{l=0}^{N-1}$ as the unique solution of the regularized problem:

$$\begin{array}{cc}\hfill {\left(s_l^{\ell +1}\right)}_{l=0}^{N-1}& =\underset{\mathbf{s};{\sum }_{l=0}^{N-1}{s}_l=1}{argmin}{\mathcal{P}}_{\lambda }(\mathbf{s},R):={\left|\right|R\mathbf{s}\left|\right|}^2+{\lambda \left|\right|\mathbf{s}\left|\right|}^2\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & =\frac{{{(}^{t}RR+\lambda {\mathrm{Id}}_N)}^{-1}{\mathbf{1}}_N}{\langle {{(}^{t}RR+\lambda {\mathrm{Id}}_N)}^{-1}{\mathbf{1}}_N,{\mathbf{1}}_N\rangle }.\hfill \end{array}$$

(17)

where the columns of $R:=[r\left(y^{\ell }\right),\cdots ,r\left(y^{\ell +1-N}\right)]$ are the N previous residuals. The regularization parameter $\lambda >0$ generalizes the original Anderson acceleration [19] introduced for $\lambda =0$. Taking $\lambda >0$ indeed leads to a more stable numerical estimation of the extrapolation parameters.

3.2. Application to SK

We now detail the whole Sinkhorn–Knopp algorithm using regularized nonlinear acceleration, which is presented in Algorithm 1.

In all our experiments corresponding to $N>1$, we considered a regularization $\lambda =1\times {10}^{-10}$ for the weight estimation (16) within the RNA scheme. For the SK algorithm, we considered the log-stabilization implementation proposed in [9,12,35] to avoid numerical errors for low values of $\epsilon$. This algorithm acts on the dual variables $\alpha ,\beta$. We refer to the aforementioned papers for more details.

As the SK algorithm successively projects the matrix ${\gamma }^{\ell }$ onto the set of linear constraints ${𝒞}_k$, $k=1,2$, we took as the convergence criterion the error realized on the first marginal of the transport matrix $\gamma =diag(exp(\alpha /\epsilon \left)\right)exp(-c/\epsilon )diag(exp(\beta /\epsilon \left)\right)$, that is ${\sum }_i|{\sum }_j{\gamma }_{i,j}^{\ell }-{\left({\mu }_1\right)}_i|<\eta =1e-9$. Note that the variable $\alpha$ is introduced into the algorithm only for computing the convergence criterion.

Algorithm 1 RNA SK Algorithm in the Log Domain.

Require: ${\mu }_1\in {\mathbb{R}}^{n_1}$, ${\mu }_2\in {\mathbb{R}}^{n_2}$, $c\in {\mathbb{R}}_{+}^{n_1\times n_2}$  Set $\ell =0$, ${\beta }^0,y^0={\mathbf{0}}_{n_2}$, ${\gamma }^0=exp-c/\epsilon$, $\omega \in \mathbb{R}$, $N>0$ and $\eta >0$  Set $x_i=-{max}_j\left(c_{i,j}-y_j^0\right)$ and ${\alpha }_i=-{max}_j\left(c_{i,j}-{\beta }_j^0\right)$,         $i=1\cdots n_1$  while $\left|\right|exp(\alpha /\epsilon )\otimes ({\gamma }^{0}exp({\beta }^{\ell }/\epsilon ))-{\mu }_1\left|\right|>\eta$ do   $\tilde{N}=min(N,\ell +1)$   $R=[{\beta }^{\ell }-y^{\ell },\cdots ,{\beta }^{\ell +1-\tilde{N}}-y^{\ell +1-\tilde{N}}]$   $\mathbf{w}={{(}^{t}RR+\lambda {\mathrm{Id}}_{\tilde{N}})}^{-1}{\mathbf{1}}_{\tilde{N}}/\langle {{(}^{t}RR+\lambda {\mathrm{Id}}_{\tilde{N}})}^{-1}{\mathbf{1}}_{\tilde{N}},{\mathbf{1}}_{\tilde{N}}\rangle$   $y^{\ell +1}={\sum }_{l=0}^{\tilde{N}-1}{w}_l((1-\omega )y^{\ell -l}+\omega {\beta }^{\ell -l})$   ${\tilde{x}}_i=-{max}_j(c_{i,j}-y_j^{\ell +1})$,                       $i=1\cdots n_1$   $x_i={\tilde{x}}_i-\epsilon log({\sum }_{j}exp((-c_{i,j}+{\tilde{x}}_i+y_j^{\ell +1})/\epsilon -log{\left({\mu }_1\right)}_i))$,      $i=1\cdots n_1$   ${\beta }_j^{\ell +1}=y^{\ell +1}-\epsilon log({\sum }_{i}exp((-c_{i,j}+x_i+y_j^{\ell +1})/\epsilon -log{\left({\mu }_2\right)}_j))$,      $j=1\cdots n_2$   ${\alpha }_i=-{max}_j(c_{i,j}-{\beta }_j^{\ell +1})$,                       $i=1\cdots n_1$   $\ell \leftarrow \ell +1$  end while  return${\gamma }_{i,j}=diag(exp(\alpha /\epsilon )){\gamma }^{0}diag(exp(\beta /\epsilon ))$

We now present numerical results obtained with random cost matrices of size $100\times 100$ with entries uniform in $[0,1]$ and uniform marginals ${\mu }_1$ and ${\mu }_2$. All convergence plots are mean results over 20 realizations.

We first considered the relaxation parameter $\omega =1$, in order to recover the original SK algorithm for $N=1$. In Figure 1, we show the convergence results obtained with RNA orders $N\in \{1,2,4,8\}$ on four regularized transport problems corresponding to entropic parameters $\epsilon \in \{0.003,0.01,0.03,0.1\}$. Figure 1 first illustrates that the convergence is improved with higher RNA orders N.

The acceleration is also greater for low values of the regularization parameter $\epsilon$. This is an important behavior, as many iterations are required to have an accurate estimation of these challenging regularized transport problems. In the settings $\epsilon \in \{0.003,0.01\}$, a speedup of more than $\times 100$ in terms of the iteration number was observed between RNA orders $N=8$ and $N=1$ (SK) to reach the same convergence threshold. We did not observe a significant improvement by considering higher RNA orders such as $N=16$.

Next, we focus on the influence of the relaxation parameter $\omega \in \{0.5,1,1.5,1.9\}$ onto the behavior of RNA schemes of orders $N\in \{1,2,8\}$. We restricted our analysis to the more challenging settings $\epsilon =0.01$ and $\epsilon =0.003$. As illustrated in Figure 2, increasing $\omega$ systematically led to improvements in the case $N=1$. For other RNA orders satisfying $N>1$, we did not observe clear tendencies. Taking $\omega \approx 1.5$ generally allowed accelerating the convergence.

To sum up, we recommend the use of the parametrization $N=8$ and $\omega =1.5$ when applying the RNA method to the SK algorithm.

We recall that the convergence of such approaches is not ensured. This last experiment nevertheless suggested that in the case $N=1$, there is room to accelerate the original SK algorithm ($\omega =1$), while keeping its global convergence guarantees, by looking at overrelaxed schemes with parameters $\omega >1$.

3.3. Discussion

The nonlinear regularized acceleration algorithm involves relevant numerical accelerations without convergence guarantees. To build an algorithm that ensures the convergence of iterates, but also improves the numerical behavior of the SK algorithm, we now propose to follow a different approach using Lyapunov sequences, which is a classical tool to study optimization algorithms. The new scheme proposed here uses the specific form of the SK algorithm with the two variables $\alpha$ and $\beta$. It performs two successive overrelaxations (SOR) at each step, one for the update of $\alpha$ and one for $\beta$. The algorithm does not use any mixing scheme, but the simple structure allows defining a sequence, called a Lyapunov sequence, which decreases at each step. This Lyapunov approach allows ensuring the convergence of the algorithm for a suitable choice of the overrelaxation parameter.

The algorithm can be summarized as follows:

$$\begin{array}{cc}\hfill {\alpha }^{\ell +1}& =(1-\omega ){\alpha }^{\ell }+\omega \underset{\alpha }{argmax}E(\alpha ,{\beta }^{\ell })\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\beta }^{\ell +1}& =(1-\omega ){\beta }^{\ell }+\omega \underset{\beta }{argmax}E({\alpha }^{\ell +1},\beta ).\hfill \end{array}$$

(19)

Our convergence analysis relied on an online adaptation of the overrelaxation parameter $\omega \in [1,2)$. As illustrated by Figure 3, in the case for $\epsilon =0.01$, the proposed SK-SOR method was not as performant as high RNA orders $N>1$ with $\omega =1.5$. It nevertheless gave an important improvement with respect to the original SK method, while being provably convergent.

4. Sinkhorn–Knopp with Successive Overrelaxation

In this section, we propose a globally convergent overrelaxed SK algorithm. Different from the RNA point of view of the previous section, our algorithm SK-SOR relies on successive overrelaxed (SOR) projections.

As illustrated in Figure 4a,b, the original SK algorithm (5) performs alternate Bregman projections (4) onto the affine sets ${𝒞}_1$ and ${𝒞}_2$. In practice, the convergence may degrade when $\epsilon \to 0$. The idea developed in this section is to perform overrelaxed projections in order to accelerate the process, as displayed in Figure 4c.

In what follows, we first define overrelaxed Bregman projections. We then propose a Lyapunov function that is used to show the global convergence of our proposed algorithm in Section 4.3. The local convergence rate is then discussed in Section 4.4.

4.1. Overrelaxed Projections

We recall that $P_{{𝒞}_k}$ are operators realizing the Bregman projection of matrices $\gamma \in {\mathbb{R}}^{n_1{n}_2}$ onto the affine sets ${𝒞}_k$, $k=1,2$, as defined in (4). For $\omega \ge 0$, we now define the $\omega$-relaxed projection operator $P_{{𝒞}_k}^{\omega }$ as:

$$log{P}_{{𝒞}_k}^{\omega }\left(\gamma \right)=(1-\omega )log\gamma +\omega log{P}_{{𝒞}_k}\left(\gamma \right),$$

(20)

where the logarithm is taken coordinate-wise. Note that $P_{{𝒞}_k}^0$ is the identity, $P_{{𝒞}_k}^1=P_{{𝒞}_k}$ is the standard Bregman projection, and $P_{{𝒞}_k}^2$ is an involution (in particular because ${𝒞}_k$ is an affine subspace). In the following, we will consider overrelaxations corresponding to $\omega \in [1;2)$. A naive algorithm would then consist of iteratively applying $P_{{𝒞}_2}^{\omega }\circ P_{{𝒞}_1}^{\omega }$ for some choice of $\omega$. While it often behaves well in practice, this algorithm may sometimes not converge even for reasonable values of $\omega$. Our goal in this section is to make this algorithm robust and to guarantee its global convergence.

Remark 2.

Duality gives another point of view on the iterative overrelaxed Bregman projections: they indeed correspond to a successive overrelaxation (SOR) algorithm on the dual objective E given in (6). This is a procedure that, starting from $({\alpha }^0,{\beta }^0)=(\mathbf{0},\mathbf{0})$, defines for $\ell \in {\mathbb{N}}^{*}$,

$$\begin{array}{cc}\hfill {\alpha }^{\ell +1}& =(1-\omega ){\alpha }^{\ell }+\omega \underset{\alpha }{argmax}E(\alpha ,{\beta }^{\ell })\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\beta }^{\ell +1}& =(1-\omega ){\beta }^{\ell }+\omega \underset{\beta }{argmax}E({\alpha }^{\ell +1},\beta ),\hfill \end{array}$$

(22)

From the definition of the projections in (4) and using again the relationships $u_i=e^{{\alpha }_i/\epsilon }$, $v_j=e^{{\beta }_j/\epsilon }$ and ${\gamma }_{i,j}^0=e^{-c_{i,j}/\epsilon }$, Expressions (21) and (22) correspond to the overrelaxed projections (20).

4.2. Lyapunov Function

Convergence of the successive overrelaxed projections is not guaranteed in general. In order to derive a robust algorithm with provable convergence, we introduced the Lyapunov function:

$$F\left(\gamma \right)=KL({\gamma }^{*},\gamma ),$$

(23)

where ${\gamma }^{*}$ denotes the solution of the regularized OT problem. We used this function to enforce the strict descent criterion $F\left({\gamma }^{\ell +1}\right)<F\left({\gamma }^{\ell }\right)$ as long as the process has not converged.

The choice of (23) as a Lyapunov function is of course related to the fact that Bregman projections are used throughout the algorithm. Further, we show (Lemma 1) that its decrease is simple to compute, and this descent criterion still allows enough freedom in the choice of the overrelaxation parameter.

Crucial properties of this Lyapunov function are gathered in the next lemma.

Lemma 1.

For any $M\in {\mathbb{R}}_{+}^{*}$, the sublevel set $\left\{\gamma \mid F\left(\gamma \right)\le M\right\}$ is compact. Moreover, for any γ in ${\mathbb{R}}_{+*}^{n_1{n}_2}$, the decrease of the Lyapunov function after an overrelaxed projection can be computed as:

$$F\left(\gamma \right)-F(P_{{𝒞}_k}^{\omega }\left(\gamma \right))=〈{\mu }_k,{\phi }_{\omega }\left(\left(A_k\gamma \right)\oslash {\mu }_k\right)〉,$$

(24)

where

$${\phi }_{\omega }\left(x\right)=x(1-x^{-\omega })-\omega logx$$

(25)

is a real function, applied coordinate-wise.

Proof. 

The fact that the Kullback–Leibler divergence is jointly lower semicontinuous implies in particular that K is closed. Moreover, $K\subset {\mathbb{R}}_{+}^{n_1\times n_2}$ is bounded because F is the sum of nonnegative, coercive functions of each component of its argument $\gamma$.

Formula (24) comes from the expression:

$$F\left({\gamma }^1\right)-F\left({\gamma }^2\right)=\sum _{i,j}\left({\gamma }_{i,j}^{*}log({\gamma }_{i,j}^2/{\gamma }_{i,j}^1)+{\gamma }_{i,j}^1-{\gamma }_{i,j}^2\right)$$

and Relations (4) and (20).    □

It follows from Lemma 1 that the decrease of F for an overrelaxed projection is very inexpensive to estimate, since its computational cost is linear with respect to the dimension of data ${\mu }_k$. In Figure 5, we display the function ${\phi }_{\omega }\left(x\right)$. Notice that for the Sinkhorn–Knopp algorithm, which corresponds to $\omega =1$, the function ${\phi }_{\omega }$ is always nonnegative. For other values $1\le \omega <2$, it is nonnegative for x close to one.

4.3. Proposed Algorithm

We first give a general convergence result that later serves as a basis to design an explicit algorithm.

Theorem 1.

Let ${\mathsf{\Theta }}_1$ and ${\mathsf{\Theta }}_2$ be two continuous functions of γ such that:

$$\forall \gamma \in {\mathbb{R}}_{+*}^{n_1{n}_2},F\left(P_{{𝒞}_k}^{{\mathsf{\Theta }}_k\left(\gamma \right)}\left(\gamma \right)\right)\le F\left(\gamma \right),$$

(26)

where the inequality is strict whenever $\gamma \notin {𝒞}_k$. Consider the sequence defined by ${\gamma }^0=e^{-c/\epsilon }$ and:

$$\begin{array}{cc}\hfill {\tilde{\gamma }}^{\ell +1}& =P_{{𝒞}_1}^{{\mathsf{\Theta }}_1\left({\gamma }^{\ell }\right)}\left({\gamma }^{\ell }\right)\hfill \\ \hfill {\gamma }^{\ell +1}& =P_{{𝒞}_2}^{{\mathsf{\Theta }}_2\left({\tilde{\gamma }}^{\ell +1}\right)}\left({\tilde{\gamma }}^{\ell +1}\right).\hfill \end{array}$$

Then, the sequence $\left({\gamma }^{\ell }\right)$ converges to ${\gamma }^{*}$.

Lemma 2.

Let us take ${\gamma }^0$ in ${\mathbb{R}}_{+*}^{n_1{n}_2}$, and denote:

$$S=\left\{diag\left(u\right){\gamma }^{0}diag\left(v\right),(u,v)\in {\mathbb{R}}_{+*}^{n_1+n_2}\right\}$$

the set of matrices that are diagonally similar to ${\gamma }^0$. Then, the set $S\cap {𝒞}_1\cap {𝒞}_2$ contains exactly one element ${\gamma }^{*}=P_{{𝒞}_1\cap {𝒞}_2}\left({\gamma }^0\right)$.

Proof. 

We refer to [2] for a proof of this lemma.    □

Proof of the Theorem.

First of all, notice that the operators $P_{{𝒞}_k}^{\theta }$ apply a scaling to lines or columns of matrices. All $\left({\gamma }^{\ell }\right)$ are thus diagonally similar to ${\gamma }^0$:

$$\forall \ell \ge 0,{\gamma }^{\ell }\in S$$

By construction of the functions ${\mathsf{\Theta }}_k$, the sequence of values of the Lyapunov function $\left(F\left({\gamma }^{\ell }\right)\right)$ is non-increasing. Hence, $\left({\gamma }^{\ell }\right)$ is precompact. If $\zeta$ is a cluster point of $\left({\gamma }^{\ell }\right)$, let us define:

$$\begin{array}{cc}\hfill \tilde{\zeta }& =P_{{𝒞}_1}^{{\mathsf{\Theta }}_1\left(\zeta \right)}\left(\zeta \right)\hfill \\ \hfill {\zeta }^{\prime }& =P_{{𝒞}_2}^{{\mathsf{\Theta }}_2\left(\tilde{\zeta }\right)}\left(\tilde{\zeta }\right).\hfill \end{array}$$

Then, by the continuity of the applications, $F\left(\zeta \right)=F\left(\tilde{\zeta }\right)=F\left({\zeta }^{\prime }\right)$. From the hypothesis made on ${\mathsf{\Theta }}_1$ and ${\mathsf{\Theta }}_2$, it can be deduced that $\zeta$ is in ${𝒞}_1$ and is thus a fixed point of $P_{{𝒞}_1}$, while $\tilde{\zeta }$ is in ${𝒞}_2$. Therefore, ${\zeta }^{\prime }=\tilde{\zeta }=\zeta$ is in the intersection ${𝒞}_1\cap {𝒞}_2$. By Lemma 2, $\zeta ={\gamma }^{*}$, and the whole sequence $\left({\gamma }^{\ell }\right)$ converges to the solution.    □

We can construct explicitly functions ${\mathsf{\Theta }}_k$ as needed in Theorem 1 using the following lemma.

Lemma 3.

Let $1\le \omega <\theta$. Then, for any $\gamma \in {\mathbb{R}}_{+*}^{n_1{n}_2}$, one has:

$$F(P_{{𝒞}_k}^{\omega }\left(\gamma \right))\le F(P_{{𝒞}_k}^{\theta }\left(\gamma \right)).$$

(27)

Moreover, equality occurs if and only if $\gamma \in {𝒞}_k$.

Proof. 

Thanks to Lemma 1, one knows that:

$$F(P_{{𝒞}_k}^{\omega }\left(\gamma \right))-F(P_{{𝒞}_k}^{\theta }\left(\gamma \right))=〈{\mu }_k,({\phi }_{\theta }-{\phi }_{\omega })\left(\left(A_k\gamma \right)\oslash {\mu }_k\right)〉.$$

The function that maps $t\in [1,\infty )$ to ${\phi }_t\left(x\right)$ is non-increasing since ${\partial }_t{\phi }_t\left(x\right)$ = $(x^{1-t}-1)logx.$ Moreover, for $x\ne 1$, it is strictly decreasing. Thus, Inequality (27) is valid, with equality iff $A_k\gamma ={\mu }_k$.    □

We now argue that a good choice for the functions ${\mathsf{\Theta }}_k$ may be constructed as follows. Pick a target parameter ${\theta }_0\in [1;2)$, which will act as an upper bound for the overrelaxation parameter $\omega$, and a small security distance $\delta >0$. Define the functions ${\mathsf{\Theta }}^{*}$ and $\mathsf{\Theta }$ as:

$$\begin{array}{cc}\hfill {\mathsf{\Theta }}^{*}\left(w\right)& =sup\left\{\omega \in [1;2]\mid {\phi }_{\omega }\left(minw\right)\ge 0\right\},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \mathsf{\Theta }\left(w\right)& =min(max(1,{\mathsf{\Theta }}^{*}\left(w\right)-\delta ),{\theta }_0),\hfill \end{array}$$

(29)

where $minw$ denotes the smallest coordinate of the vector w.

Proposition 1.

The function:

$${\mathsf{\Theta }}_k\left(\gamma \right)=\mathsf{\Theta }\left(\left(A_k\gamma \right)\oslash {\mu }_k\right)$$

(30)

is continuous and verifies the descent condition (26).

Proof. 

Looking at Figure 5 can help understand this proof. Since the partial derivative of ${\partial }_{\omega }{\phi }_{\omega }\left(x\right)$ is nonzero for any $x<1$, the implicit function theorem proves the continuity of ${\mathsf{\Theta }}^{*}$. The function ${\mathsf{\Theta }}^{*}\left(\left(A_k\gamma \right)\oslash {\mu }_k)\right)$ is such that every term in Relation (24) is non-negative. Therefore, by Lemma 3, using this parameter minus $\delta$ ensures the strong decrease (26) of the Lyapunov function. Constraining the parameter to $[1,{\theta }_0]$ preserves this property.   □

This construction, which is often an excellent choice in practice, has several advantages:

• it allows choosing arbitrarily the parameter ${\theta }_0$ that will be used eventually when the algorithm is close to convergence (we motivate what are good choices for ${\theta }_0$ in Section 4.4);
• it is also an easy approach to having an adaptive method, as the approximation of ${\mathsf{\Theta }}^{*}$ has a negligible cost (it only requires solving a one-dimensional problem that depends on the smallest value of $\left(A_k\gamma \right)\oslash {\mu }_k$, which can be done in a few iterations of Newton’s method).

The resulting algorithm, which is proven to be convergent by Theorem 1, is written in pseudocode in Algorithm 2. The implementation in the log domain is also given in Algorithm 3. Both processes use the function $\mathsf{\Theta }$ defined implicitly in (29). The evaluation of $\mathsf{\Theta }$ is approximated in practice with a few iterations of Newton’s method on the function $\omega \mapsto {\phi }_{\omega }(minu)$, which is decreasing, as can be seen in Figure 5. With the choice ${\theta }_0=1$, one recovers exactly the original SK algorithm.

Algorithm 2 Overrelaxed SK Algorithm (SK-SOR).

Require: ${\mu }_1\in {\mathbb{R}}^{n_1}$, ${\mu }_2\in {\mathbb{R}}^{n_2}$, $c\in {\mathbb{R}}_{+}^{n_1\times n_2}$  Set $u={\mathbf{1}}_{n_1}$, $v={\mathbf{1}}_{n_2}$, ${\gamma }^0=e^{-c/\epsilon }$, ${\theta }_0\in [1;2)$ and $\eta >0$  while $\left|\right|u\otimes {\gamma }^{0}v-{\mu }_1\left|\right|>\eta$ do   $\tilde{u}={\mu }_1\oslash \left({\gamma }^{0}v\right)$,   $\omega =\mathsf{\Theta }(u\oslash \tilde{u})$   $u=u^{1-\omega }\otimes {\tilde{u}}^{\omega }$   $\tilde{v}={\mu }_2\oslash {(}^t{\gamma }^{0}u)$   $\omega =\mathsf{\Theta }(v\oslash \tilde{v})$   $v=v^{1-\omega }\otimes {\tilde{v}}^{\omega }$  end while  return$\gamma =diag\left(u\right){\gamma }^{0}diag\left(v\right)$

Algorithm 3 Overrelaxed SK Algorithm (SK-SOR) in the Log Domain.

Require: ${\mu }_1\in {\mathbb{R}}^{n_1}$, ${\mu }_2\in {\mathbb{R}}^{n_2}$, $c\in {\mathbb{R}}_{+}^{n_1\times n_2}$  Set $\alpha ={\mathbf{0}}_{n_1}$, $\beta ={\mathbf{0}}_{n_2}$, ${\gamma }^0=e^{-c/\epsilon }$, ${\theta }_0\in [1;2)$ and $\eta >0$  while $\left|\right|exp(\alpha /\epsilon )\otimes ({\gamma }^{0}exp(\beta /\epsilon ))-{\mu }_1\left|\right|>\eta$ do   $r_i={\sum }_{j}exp((-c_{i,j}+{\alpha }_i+{\beta }_j)/\epsilon -log{\left({\mu }_1\right)}_i)$,           $i=1\cdots n_1$   $\tilde{\alpha }=\alpha -\epsilon logr$   $\omega =\mathsf{\Theta }\left(r\right)$   $\alpha =(1-\omega )\alpha +\omega \tilde{\alpha }$   $s_j={\sum }_{i}exp((-c_{i,j}+{\alpha }_i+{\beta }_j)/\epsilon -log{\left({\mu }_2\right)}_j)$,         $j=1\cdots n_2$   $\tilde{\beta }=\beta -\epsilon logs$   $\omega =\mathsf{\Theta }\left(s\right)$   $\beta =(1-\omega )\beta +\omega \tilde{\beta }$  end while  return $\gamma =diag(exp(\alpha /\epsilon )){\gamma }^{0}diag(exp(\beta /\epsilon ))$

4.4. Acceleration of the Local Convergence Rate

In order to justify the acceleration of convergence that is observed in practice, we now study the local convergence rate of the overrelaxed algorithm, which follows from the classical convergence analysis of the linear SOR method. Our result involves the second largest eigenvalue of the matrix:

$$M_1=diag(1\oslash {\mu }_1){\gamma }^{*}diag(1\oslash {\mu }_2){}^t{\gamma }^{*}$$

(31)

where ${\gamma }^{*}$ is the solution to the regularized OT problem (the largest eigenvalue is one, associated with the eigenvector $\mathbf{1}$). We denote the second largest eigenvalue by $1-\eta$; it satisfies $\eta >0$ [36].

Proposition 2.

The SK algorithm converges locally at a linear rate $1-\eta$. For the optimal choice of extrapolation parameter ${\theta }^{*}=2/(1+\sqrt{\eta })$, algorithm SK-SOR converges locally linearly at a rate $(1-\sqrt{\eta })/(1+\sqrt{\eta })$. The local convergence of SK-SOR is guaranteed for $\theta \in ]0,2[$, and the linear rate is given in Figure 6 as a function of $1-\eta$ and θ.

Proof. 

In this proof, we focus on the dual problem, and we recall the relationship ${\gamma }^{\ell }=e^{{\alpha }^{\ell }/\epsilon }{\gamma }^0{e}^{{\beta }^{\ell }/\epsilon }$ between the iterates of the overrelaxed projection algorithm ${\gamma }^{\ell }$ and the iterates $({\alpha }^{\ell },{\beta }^{\ell })$ of the SOR algorithm on the dual problem (21), initialized with $({\alpha }^0,{\beta }^0)=(0,0)$. The dual problem (6) is invariant by translations of the form $(\alpha ,\beta )\mapsto (\alpha -k,\beta +k)$, $k\in \mathbb{R}$, but is strictly convex up to this invariance. In order to deal with this invariance, consider the subspace S of pairs of dual variables $(\alpha ,\beta )$ that satisfy $\sum \alpha =\sum \beta$; let ${\pi }_S$ be the orthogonal projection on S of kernel $(\mathbf{1},-\mathbf{1})$, and let $({\alpha }^{*},{\beta }^{*})\in S$ be the unique dual maximizer in S.

Since one SK-SOR iteration is a smooth map, the local convergence properties of the SK-SOR algorithm are characterized by the local convergence of its linearization, which here corresponds to the SOR method applied to the maximization of the quadratic Taylor expansion of the dual objective E at $({\alpha }^{*},{\beta }^{*})$. This defines an affine map $M_{\theta }:({\alpha }^{\ell },{\beta }^{\ell })\mapsto ({\alpha }^{\ell +1},{\beta }^{\ell +1})$ whose spectral properties are well known [22,37] (see also [38] (Chapter 4) for the specific case of convex minimization and [39] for the non-strictly convex case). For the case $\theta =1$, this is the matrix $M_1$ defined in Equation (31). The operator norm of ${\pi }_S\circ M_1$ is smaller than $1-\eta$, so the operator ${({\pi }_S\circ M_1)}^{\ell }={\pi }_S\circ M_1^{\ell }$ converges at the linear rate $1-\eta$ towards zero (observe that by construction, ${\pi }_S$ and $M_1$ are co-diagonalizable and thus commute): for any $(\alpha ,\beta )\in {\mathbb{R}}^{n_1}\times {\mathbb{R}}^{n_2}$, it holds $\parallel {\pi }_S\circ M_1^{\ell }(\alpha -{\alpha }^{*},\beta -{\beta }^{*}){\parallel }_2\le {\parallel {\pi }_S(\alpha -{\alpha }^{*},\beta -{\beta }^{*})\parallel }_2{(1-\eta )}^{\ell }$. More generally, the convergence of ${\pi }_S\circ M_{\theta }^{\ell }$ is guaranteed for $\theta \in ]0,2[$, with the linear rate:

$$\begin{array}{c}\hfill f(\theta ,\eta )=\left\{\begin{array}{cc}\theta -1\hfill & \mathrm{if}\theta >{\theta }^{*}\hfill \\ \frac{1}{2}{\theta }^2(1-\eta )-(\theta -1)+\frac{1}{2}\sqrt{(1-\eta ){\theta }^2({\theta }^2(1-\eta )-4(\theta -1))}\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(32)

This function is minimized with ${\theta }^{*}:=2/(1+\sqrt{\eta })$, which satisfies $f({\theta }^{*},\eta )=(1-\sqrt{\eta })/(1+\sqrt{\eta })$. The function f is plotted in Figure 6.

To switch from these dual convergence results to primal convergence results, we remark that ${\gamma }^{\ell }\to {\gamma }^{*}$ implies $KL({\gamma }^{\ell },{\gamma }^0)\to KL({\gamma }^{*},{\gamma }^0)$, which in turn implies $E({\alpha }^{\ell },{\beta }^{\ell })\to maxE$, so by invoking the partial strict concavity of E, we have ${\pi }_S({\alpha }^{\ell },{\beta }^{\ell })\to ({\alpha }^{*},{\beta }^{*})$. The converse implication is direct, so it holds $[{\pi }_S({\alpha }^{\ell },{\beta }^{\ell })\to ({\alpha }^{*},{\beta }^{*})]\iff [{\gamma }^{\ell }\to {\gamma }^{*}]$. We conclude by noting the fact that ${\pi }_S({\alpha }^{\ell },{\beta }^{\ell })$ converges at a linear rate, which implies the same rate on ${\gamma }^{\ell }$, thanks to the relationship between the iterates. □

Corollary 1.

The previous local convergence analysis applies to Algorithm 3 with Θ defined as in Equation (29), and the local convergence rate is given by the function of Equation (32) evaluated at the target extrapolation parameter ${\theta }_0$.

Proof. 

What we need to show is that eventually, one always has $\mathsf{\Theta }\left({\gamma }^{\ell }\right)={\theta }_0$. This can be seen from the quadratic Taylor expansion ${\phi }_{{\theta }_0}(1+z)=z^2({\theta }_0-{{\theta }_0}^2/2)+o\left(z^2\right)$, which shows that for any choice of ${\theta }_0\in ]1,2[$, there is a neighborhood of one on which ${\phi }_{{\theta }_0}(·)$ is nonnegative. □

5. Experimental Results

We compared Algorithm 2 to the SK algorithm on two very different optimal transport settings. In the setting of Figure 7a, we considered the domain $[0,1]$ discretized into 100 samples and the squared Euclidean transport cost on this domain. The marginals were densities made of the sum of a base plateau of height $0.1$ and another plateau of height and boundaries chosen uniformly in $[0,1]$, subsequently normalized. In the setting of Figure 7b, the cost was a $100\times 100$ random matrix with entries uniform in $[0,1]$, and the marginals were uniform.

Given an estimation of $1-\eta$, the local convergence rate of the SK algorithm, we define ${\theta }_0$ as the optimal parameter as given in Proposition 2. For estimating $\eta$, we followed two strategies. For strategy “estimated” (in blue on Figure 7), $\eta$ was measured by looking at the local convergence rate of the SK algorithm run on another random problem of the same setting and for the same value of $\epsilon$. For strategy “measured” (in orange on Figure 7), the parameter was set using the local convergence rate of the SK algorithm run on the same problem. Of course, the latter was an unrealistic strategy, but it was interesting to see in our experiments that the “estimated” strategy performed almost as well as the “measured” one, as shown in Figure 7.

Figure 7 displays the ratio of the number of iterations required to reach a precision of ${10}^{-6}$ on the dual variable $\alpha$ for the SK algorithm and Algorithm 2. It is worth noting that the complexity per iteration of these algorithms is the same modulo negligible terms, so this ratio is also the runtime ratio (our algorithm can also be parallelized on GPUs just as the SK algorithm). In both experimental settings, for low values of the regularization parameter $\epsilon$, the acceleration ratio was above 20 with Algorithm 2.

6. Conclusions and Perspectives

The SK algorithm is widely used to solve entropy-regularized OT. In this paper, we first showed that RNA methods are adapted to the numerical acceleration of the SK algorithm. Nevertheless, the global convergence of such approaches is not guaranteed.

Next, we demonstrated that the use of successive overrelaxed projections is a natural and simple idea to ensure and accelerate the convergence, while keeping many nice properties of the SK algorithm (first order, parallelizable, simple). We proposed an algorithm (SK-SOR) that adaptively chooses the overrelaxation parameter so as to guarantee global convergence. The acceleration of the convergence speed was numerically impressive, in particular in low regularization regimes. It was theoretically supported in the local regime by the standard analysis of SOR iterations. Nevertheless, the SK-SOR algorithm was not as performant as RNA, and no guarantee was given on the global computational complexity of either algorithm.

This idea of overrelaxation can be generalized to solve more general problems such as multi-marginal OT, barycenters, gradient flows, and unbalanced OT [38] (Chapter 4), but there is no systematic way to derive globally convergent algorithms. Our work is a step in the direction of building and understanding the properties of robust first-order algorithms for solving OT. More understanding is needed regarding SOR itself (global convergence speed, choice of ${\theta }_0$), as well as its relation to other acceleration methods [13,20].