# Overrelaxed Sinkhorn-Knopp Algorithm for Regularized Optimal Transport
## 3 Regularized Nonlinear Acceleration of the Sinkhorn-Knopp Algorithm
### Regularized Nonlinear Acceleration  
To introduce RNA, we first rewrite the SK algorithm (7) and (8) as:  
\[\beta^{\ell+1}=\operatorname{SK}(\beta^{\ell}):=\operatorname*{argmax}_{ \beta}E\biggl{(}\operatorname*{argmax}_{\alpha}E(\alpha,\beta^{\ell}),\beta \biggr{)}. \tag{11}\]  
The goal of this algorithm is to build a sequence \((\beta^{\ell})_{l\geq 1}\) converging to a fixed point of SK, i.e., to a point \(\beta^{*}\) satisfying:  
\[\beta^{*}=SK(\beta^{*}). \tag{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 \((x^{n})_{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 \((x^{k})_{k}\) and \((Tx^{k})_{k}\) to provide a value of \(x^{n}\) such that \(\|x^{n}-Tx^{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\
abla 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 \((x^{n})_{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:  
\[(s_{l}^{\ell+1})_{l=0}^{N-1} =\operatorname*{argmin}_{\mathbf{s};\,\sum_{l=0}^{N-1}s_{l}=1} \mathcal{P}_{\lambda}\Big{(}\mathbf{s},(\beta^{\ell-l}-y^{\ell-l})_{l=0}^{N-1} \Big{)} \tag{13}\] \[y^{\ell+1} =\sum_{l=0}^{N-1}s_{l}^{\ell+1}(y^{\ell-l}+\omega(\beta^{\ell-l}- y^{\ell-l}))\] (14) \[\beta^{\ell+1} =\operatorname{SK}(y^{\ell+1}), \tag{15}\]  
where extrapolation weights \(\mathbf{s}=(s_{l})_{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\geqslant 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 [PERSON]'s method [33], relaxation leads to numerical convergence improvements in gradient descent schemes [34]. [PERSON] suggested underrelaxing with \(\omega\in(0;1]\), while the authors of [21] proposed to take \(\omega=1\).  
**Extrapolation:** Let us define the residual \(r(y)=\operatorname{SK}(y)-y\). As the objective is to estimate the fixed point of SK, the extrapolation step builds a vector \(y\) such that \(||r(y)||\) 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 \((s_{l}^{\ell+1})_{l=0}^{N-1}\) as the unique solution of the regularized problem:  
\[(s_{l}^{\ell+1})_{l=0}^{N-1} =\operatorname*{argmin}_{\mathbf{s};\,\sum_{l=0}^{N-1}s_{l}=1} \mathcal{P}_{\lambda}(\mathbf{s},R):=||R\mathbf{s}||^{2}+\lambda||\mathbf{s}|| ^{2} \tag{16}\] \[=\frac{(t^{\prime}RR+\lambda\mathrm{Id}_{N})^{-1}\mathbf{1}_{N}} {((t^{\prime}RR+\lambda\mathrm{Id}_{N})^{-1}\mathbf{1}_{N},\mathbf{1}_{N})}. \tag{17}\]  
where the columns of \(R:=[r(y^{\ell}),\cdots,r(y^{\ell+1-N})]\) 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.