# Overrelaxed Sinkhorn-Knopp Algorithm for Regularized Optimal Transport
## 3 Regularized Nonlinear Acceleration of the Sinkhorn-Knopp Algorithm
### 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 \(\varepsilon\). 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 \(\mathcal{C}_{k},\,k=1,2\), we took as the convergence criterion the error realized on the first marginal of the transport matrix \(\gamma=\operatorname{diag}(\exp(\alpha/\varepsilon))\exp(-c/\varepsilon) \operatorname{diag}(\exp(\beta/\varepsilon))\), that is \(\sum_{i}|\sum_{j}\gamma_{ij}^{\ell}-(\mu_{1})_{i}|<\eta=1e-9\). Note that the variable \(\alpha\) is introduced into the algorithm only for computing the convergence criterion. ```
0:\(\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/\varepsilon\), \(\omega\in\mathbb{R}\), \(N>0\) and \(\eta>0\)  Set \(x_{i}=-\max_{j}(c_{i,j}-y^{0}_{j})\) and \(\alpha_{i}=-\max_{j}(c_{i,j}-\beta^{0}_{j})\), \(i=1\cdots n_{1}\) while\(||\exp(\alpha/\varepsilon)\otimes(\gamma^{0}\exp(\beta^{\ell}/\varepsilon))-\mu_{1}||>\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 }}/(^{t}RR+\lambda\mathrm{Id}_{\tilde{N}})^{-1}\mathbf{1}_{\tilde{N}^{*}} \mathbf{1}_{\tilde{N}})\) \(y^{\ell+1}=\sum_{l=0}^{N-1}w_{l}((1-\omega)y^{\ell-1}+\omega\beta^{\ell-l})\) \(\tilde{x}_{i}=-\max_{j}(c_{i,j}-y^{\ell+1}_{j})\) \(i=1\cdots n_{1}\) \(x_{i}=\tilde{x}_{i}-\varepsilon\log(\sum_{j}\exp\left((-c_{i,j}+\tilde{x}_{i} +y^{\ell+1}_{j})/\varepsilon-\log(\mu_{1})_{i}\right))\), \(i=1\cdots n_{1}\) \(\beta^{\ell+1}_{j}=y^{\ell+1}-\varepsilon\log(\sum_{i}\exp\left((-c_{i,j}+x_{i }+y^{\ell+1}_{j})/\varepsilon-\log(\mu_{2})_{j}\right))\), \(j=1\cdots n_{2}\) \(\alpha_{i}=-\max_{j}(c_{i,j}-\beta^{\ell+1}_{j})\), \(i=1\cdots n_{1}\) \(\ell\leftarrow\ell+1\) endwhile return\(\gamma_{i,j}=\operatorname{diag}(\exp(\alpha/\varepsilon))\gamma^{0} \operatorname{diag}(\exp(\beta/\varepsilon))\)
```  
**Algorithm 1** RNA SK Algorithm in the Log Domain. 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 \(\varepsilon\in\{0.003,0.01,0.03,0.1\}\). Figure 1 first illustrates that the convergence is improved with higher RNA orders \(N\). Figure 1: Convergence of RNA schemes for a relaxation parameter \(\omega=1\) and different orders \(N\in\{1,2,4,8\}\). All approaches lead to similar convergence as the original SK algorithm (order \(N=1\), with blue dots) for high values of the entropic parameter such as \(\varepsilon=0.1\). When facing more challenging regularized optimal transport problems, high-order RNA schemes (\(N>1\)) realize important accelerations. This behavior is highlighted in the bottom row for \(\varepsilon=0.01\) and \(\varepsilon=0.003\). In these settings, with respect to SK, RNA of order \(N=8\) (plain red curves) reduces by a factor 100 the number of iterations required to reach the same accuracy.