Gradient Methods for Problems with Inexact Model of the Objective

TL;DR

This paper develops gradient methods for convex optimization problems with inexact objective models, providing convergence analysis and applying the framework to clustering, optimal transport, and barycenter problems.

Contribution

It introduces a general inexact model framework encompassing various inexact oracle conditions and analyzes gradient methods within this setting.

Findings

Convergence rates are established for convex and strongly convex cases.

Proximal Sinkhorn algorithm effectively approximates optimal transport distance.

Proximal Iterative Bregman Projections algorithm efficiently computes optimal transport barycenters.

Abstract

We consider optimization methods for convex minimization problems under inexact information on the objective function. We introduce inexact model of the objective, which as a particular cases includes $(\delta,L)$ inexact oracle and relative smoothness condition. We analyze gradient method which uses this inexact model and obtain convergence rates for convex and strongly convex problems. To show potential applications of our general framework we consider three particular problems. The first one is clustering by electorial model introduced in [Nesterov, 2018]. The second one is approximating optimal transport distance, for which we propose a Proximal Sinkhorn algorithm. The third one is devoted to approximating optimal transport barycenter and we propose a Proximal Iterative Bregman Projections algorithm. We also illustrate the practical performance of our algorithms by numerical experiments.

Figures

Tables

Table 1. Table 1: Results for Example 1 .

	Non-adaptive		Adaptive
k	Time	Estimate	Time	Estimate
160	0:01:19	0.19827	0:05:25	0.02110
180	0:01:27	0.16220	0:05:55	0.01258
200	0:01:36	0.13264	0:07:11	0.00750
220	0:01:55	0.10849	0:07:19	0.00474
240	0:01:57	0.08873	0:07:56	0.00282

Table 2. Table 2: Results for Example 2 .

	Adaptive
k	Time	Estimate
50	0:07:37	0.71273
100	0:14:27	0.51241
150	0:23:00	0.372301
200	0:28:07	0.27334
250	0:34:32	0.19699
300	0:43:10	0.14456

Equations

$$\displaystyle f(x)\rightarrow\min_{x\in Q}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}$$

$$0\le f(x)-(f(y)+{\psi }_{\delta }(x,y))\le LV[y](x)+\delta$$

$$\mu V[y](x)\le f(x)-(f(y)+{\psi }_{\delta }(x,y))\le LV[y](x)+\delta$$

$$0\le f(x)-f(y)-⟨g_{\delta }(y),x-y⟩\le \frac{L}{2}{\Vert x-y\Vert }^2+\delta$$

$$\frac{\mu }{2}{\Vert x-y\Vert }^2\le f(x)-f(y)-⟨g_{\delta ,L,\mu }(y),x-y⟩\le \frac{L}{2}{\Vert x-y\Vert }^2+\delta \forall x\in Q,$$

$$\varphi (x)\to \underset{x\in Q\subseteq {\mathbb{R}}^n}{\min }.$$

$${\varphi }_{k+1}(x):={\psi }_{\delta }(x,x_k)+LV[x_k](x),x_{k+1}:={\arg \underset{x\in Q}{\min }}^{\tilde{\delta }}{\varphi }_{k+1}(x).$$

$$f({\bar{x}}_N)-f(x_{\ast })\le \frac{L{R}^2}{N}+\tilde{\delta }+\delta ,$$

$$f(x_{k+1})\le f(x_k)+{\psi }_{\delta }(x_{k+1},x_k)+L_{k+1}V[x_k](x_{k+1})+\delta ,$$

$${\varphi }_{k+1}(x):={\psi }_{\delta }(x,x_k)+L_{k+1}V[x_k](x),x_{k+1}:={\arg \underset{x\in Q}{\min }}^{\tilde{\delta }}{\varphi }_{k+1}(x).$$

$$1-\frac{\mu }{\hat{L}}=\sqrt[k+1]{\left(1-\frac{\mu }{L^{k+1}}\right)\left(1-\frac{\mu }{L_k}\right)\dots \left(1-\frac{\mu }{L_1}\right)}.$$

$$\min\limits_{1\leq i\leq k+1}L_{i}\leq\hat{L}\leq\max\limits_{1\leq i\leq k+1}L_{i}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\leq 2L}.$$

$$V[x^{k+1}](x_{\ast })\le \frac{2L(\delta +\tilde{\delta })}{{\mu }^2}\left(1-{\left(1-\frac{\mu }{2L}\right)}^{k+1}\right)+{\left(1-\frac{\mu }{\hat{L}}\right)}^{k+1}V[x^0](x_{\ast }),$$

$$f(x^{k+1})-f(x_{\ast })\le \frac{4L^2(\delta +\tilde{\delta })}{{\mu }^2}\left(1-{\left(1-\frac{\mu }{2L}\right)}^{k+1}\right)+2L{\left(1-\frac{\mu }{\hat{L}}\right)}^{k+1}V[x^0](x_{\ast }).$$

$$f_{{\mu }_1,{\mu }_2}(x=(z,p))=g(x)+{\mu }_1\sum _{k=1}^n{z}_k\ln z_k+\frac{{\mu }_2}{2}\Vert p{\Vert }_2^2\to \underset{z\in S_n(1),p\in {\mathbb{R}}_{+}^m}{\min },$$

$$\Vert \nabla g(x)-\nabla g(y){\Vert }_{\ast }\le L_g\Vert x-y\Vert \forall x,y\in S_n(1)\times {\mathbb{R}}_{+}^m,$$

$$\Vert x+y\Vert =\sqrt{\Vert z_x+z_y{\Vert }_1^2+\Vert p_x+p_y{\Vert }_2^2}\le \sqrt{(\Vert z_x{\Vert }_1+\Vert z_y{\Vert }_1{)}^2+(\Vert p_x{\Vert }_2+\Vert p_y{\Vert }_2{)}^2}\le$$

$$\le \sqrt{\Vert z_x{\Vert }_1^2+\Vert p_x{\Vert }_2^2}+\sqrt{\Vert z_y{\Vert }_1^2+\Vert p_y{\Vert }_2^2}=\Vert x\Vert +\Vert y\Vert ,$$

$${\psi }_{\delta }(x,y)=⟨\nabla g(y),x-y⟩-L_g\cdot KL(z_x|z_y)-\frac{L_g}{2}\Vert p_x-p_y{\Vert }_2^2+$$

$$+{\mu }_1(KL(z_x|\text{1})-KL(z_y|\text{1}))+\frac{{\mu }_2}{2}\left(\Vert p_x{\Vert }_2^2-\Vert p_y{\Vert }_2^2\right)$$

$$V[y](x)=KL(z_x|z_y)+\frac{1}{2}\Vert p_x-p_y{\Vert }_2^2.$$

$$|g(x)-g(y)-⟨\nabla g(y),x-y⟩|\le \frac{L_g}{2}\Vert x-y{\Vert }^2.$$

$$=g(x)-g(y)-⟨\nabla g(y),x-y⟩+{\mu }_1\cdot KL(z_x|\text{1})-{\mu }_1\cdot KL(z_y|\text{1})+$$

$$+\frac{{\mu }_2}{2}\Vert p_x{\Vert }_2^2-\frac{{\mu }_2}{2}\Vert p_y{\Vert }_2^2-{\mu }_1\cdot KL(z_x|\text{1})+{\mu }_1\cdot KL(z_y|\text{1})+L_g\cdot KL(z_x|z_y)+$$

$$+\frac{L_g}{2}\Vert p_x-p_y{\Vert }_2^2-\frac{{\mu }_2}{2}\Vert p_x{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert p_y{\Vert }_2^2=$$

$$=g(x)-g(y)-⟨\nabla g(y),x-y⟩+L_g\cdot KL(z_x|z_y)+\frac{L_g}{2}\Vert p_x-p_y{\Vert }_2^2.$$

$$f_{{\mu }_1,{\mu }_2}(x)-f_{{\mu }_1,{\mu }_2}(y)-{\psi }_{\delta }(x,y)\le \frac{L_g}{2}\Vert x-y{\Vert }^2+L_g\cdot KL(z_x|z_y)+\frac{L_g}{2}\Vert p_x-p_y{\Vert }_2^2,$$

$$f_{{\mu }_1,{\mu }_2}(x)-f_{{\mu }_1,{\mu }_2}(y)-{\psi }_{\delta }(x,y)\ge -\frac{L_g}{2}\Vert x-y{\Vert }^2+L_g\cdot KL(z_x|z_y)+\frac{L_g}{2}\Vert p_x-p_y{\Vert }_2^2.$$

$$0\le f_{{\mu }_1,{\mu }_2}(x)-f_{{\mu }_1,{\mu }_2}(y)-{\psi }_{\delta }(x,y)\le 2L_g\cdot KL(z_x|z_y)+L_g\Vert p_x-p_y{\Vert }_2^2=2L_{g}V[y](x),$$

$${\psi }_{\delta }(x,y)={\psi }_{\delta }^{lin}(x,y)+({\mu }_1-L_g)\cdot KL(z_x|z_y)+\frac{{\mu }_2-L_g}{2}\Vert p_x-p_y{\Vert }_2^2\ge$$

$$\ge (\min \{{\mu }_1,{\mu }_2\}-L_g)\cdot V[y](x),$$

$${\psi }_{\delta }^{lin}(x,y)=⟨\nabla g(y),x-y⟩+{\mu }_1⟨\nabla KL(z_y|1),z_x-z_y⟩+{\mu }_2⟨p_y,p_x-p_y⟩$$

Full text

1 Introduction

In this paper we consider optimization methods for convex problems under inexact information on the objective function. This information is given by an object, which we call inexact model. Inexact model generalizes the inexact oracle introduced in [19], where inexactness is assumed to be present in the objective value and its gradient. The authors show that, based on these two objects, it is possible to construct a linear function, which is a lower approximation and, up to a quadratic term, an upper approximation of the objective, and these two approximations are enough to obtain convergence rates for gradient method and accelerated gradient method. We go beyond and assume that the approximations of the objective are given through some function, which is not necessarily linear.

This allows us to construct general gradient-type method which is applicable in for different problem classes and allows to obtain convergence rates in these situations as a corollary of our general theorem. Besides convex problems we focus also on strongly convex objectives and illustrate the application of our general theory by two examples. The first example is data clustering by electoral model [49]. The second example relates to Wasserstein distance and barycenter, which are widely used in data analysis [15, 16].

Many optimization methods use some model of the objective function to define a step by minimization of this model. Usually the model is constructed using exact first-order [46, 21, 52], second-order [51], or higher-order information [11, 48] information on the objective. The influence of inexactness on the convergence of gradient-type methods have being studied at least since [55]. Accelerated first-order methods with inexact oracle are studied in [17, 44, 19, 24, 14]. Some recent works study also non-convex problems in this context [10, 22]. Randomized methods with inexact oracle are also studied in the literature, e.g. coordinate descent in [61, 32], random gradient-free methods and random directional derivative methods in [27, 26]. A method with inexact oracle for variational inequalities can be found in [31].

The contributions of this paper can be summarized as follows.

• $\square$

We introduce an inexact model of the objective function for convex optimization problems and strongly convex optimization problems.

• $\square$

We introduce and theoretically analyze a gradient-type method for convex and strongly convex problems with an inexact model of the objective function. For the latter case we prove linear rate of convergence.

• $\square$

We apply our method to, generally speaking, non-convex optimization problem which arises in clustering model introduced in [49]. To do this we construct an inexact model and apply our general algorithms and convergence theorems.

• $\square$

We apply our general framework for Wasserstein distance and barycenter problems and show that it allows to construct a proximal á la [12] version of the Sinkhorn’s algorithm [58] and Iterative Bregman Projection algorithm [7].

**Notation. We define $\textbf{1}=(1,...,1)^{T}\in{\mathbb{R}}^{n}$, $KL(z|t)$ to be the Kullback-Leibler divergence: $KL(z|t)=\sum\limits_{k=1}^{n}z_{k}\ln(z_{k}/t_{k})$, $\forall z,t\in S_{n}(1)$, where $S_{n}(1)$ is the standard simplex in ${\mathbb{R}}^{n}$. We also denote by $\odot$ the entrywise product of two matrices. **

2 Gradient Methods with Inexact Model of the Objective

Consider the convex optimization problem

[TABLE]

where function $f$ is convex and $Q\subseteq{\mathbb{R}}^{n}$ is a simple convex compact set. Moreover, assume that $\min_{x\in Q}f(x)=f(x_{*})$ for some $x_{*}\in Q$.

To solve this problem, we introduce a norm $\|\cdot\|$ on ${\mathbb{R}}^{n}$ and a prox-function $d(x)$ which is continuous and convex. We underline that, unlike most of the literature, we do not require $d$ to be strongly convex. Without loss of generality, we assume that $\min\limits_{x\in{\mathbb{R}}^{n}}d(x)=0$. Further, we define Bregman divergence $V[y](x):=d(x)-d(y)-\langle\nabla d(y),x-y\rangle$. Next we define the inexact model of the objective function, which generalizes the inexact oracle of [19] (see also [24, 10, 28, 35, 60, 62]).

Definition 1

Let function $\psi_{\delta}(x,y)$ be convex in $x\in Q$ and satisfy $\psi_{\delta}(x,x)=0$ for all $x\in Q$.

i) We say that $\psi_{\delta}(x,y)$ is a $(\delta,L)$-model of the function $f$ at a given point $y$ with respect to $V[y](x)$ iff, for all $x\in Q$, the inequality

[TABLE]

holds for some $L,\delta>0$.

ii) We say that $\psi_{\delta}(x,y)$ is a $(\delta,L,\mu)$-model of the function $f$ at a given point $y$ with respect to $V[y](x)$ iff, for all $x\in Q$, the inequality

[TABLE]

Note that we allow $L$ to depend on $\delta$. We refer to the case i) as convex case and to the case ii) as strongly convex case.

Remark 1

In the particular case of function $f$ possessing $(\delta,L)$-oracle [19] at a given point $y$, one has

[TABLE]

and $\psi_{\delta}(x,y)=\langle g_{\delta}(y),x-y\rangle$. In the same way, if function $f$ is equipped with $(\delta,L,\mu)$-oracle [20], i.e.,

[TABLE]

we have $\psi_{\delta}(x,y)=\langle g_{\delta,L,\mu}(y),x-y\rangle$.

The algorithms we develop are based on solving auxiliary simple problems on each iteration. We assume that these problems can be solved inexactly and, following [6] introduce a definition of inexact solution of a problem.

Definition 2

Consider a convex minimization problem

[TABLE]

If $\phi$ is smooth, we say that we solve it with $\widetilde{\delta}$-‘precision’ ($\widetilde{\delta}\geq 0$) if we find $\tilde{x}$ s.t. $\max_{x\in Q}\langle\nabla\phi(\tilde{x}),\tilde{x}-x\rangle=\widetilde{\delta}$. If $\phi$ is general convex, we say that we solve this problem with $\widetilde{\delta}$-‘precision’ if we find $\tilde{x}$ s.t. \exists h\in\partial{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}(\widetilde{x}),\,\,\,\langle h,x_{*}-\widetilde{x}\rangle\geq-\widetilde{\delta}. In both cases we denote this $\tilde{x}$ as \operatorname*{argmin}_{x\in Q}^{\widetilde{\delta}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}(x).

We notice that the case $\widetilde{\delta}=0$ corresponds to the case when $\tilde{x}$ is an exact solution of convex optimization problem (4) [6, 46]. The connection of Definition 2 with standard definitions of inexact solution, e.g. in terms of the objective residual, can be found in Appendix 0.G.

2.1 Convex Case

In this subsection we describe a gradient-type method for problems with $(\delta,L)$-model of the objective. This algorithm is a natural extension of gradient method, see [35, 60, 62].

Theorem 2.1

Let $V[x_{0}](x_{*})\leq R^{2}$, where $x_{0}$  is the starting point, and $x_{*}$  is the nearest minimum point to the point $x_{0}$ in the sense of Bregman divergence $V[y](x)$. Then, for the sequence, generated by Algorithm 1 the following inequality holds:

[TABLE]

In appendix 0.A we prove this theorem and provide an adaptive version of Algorithm 2.1, which does not require knowledge of the constant $L$.

2.2 Strongly Convex Case

In this subsection we consider problem (1) with $(\delta,L,\mu)$-model of the objective function satisfying (3). This more strong assumption allows us to obtain linear rate of convergence of the proposed algorithm. Our algorithm is listed as Algorithm 2 and it is a version of Algorithm 1, which is adaptive to possibly unknown constant $L$.

Let’s introduce average parameter $\hat{L}$:

[TABLE]

Note that by $L_{i}\geq\mu~{}(i=1,2,\ldots)$

[TABLE]

The following result holds.

Theorem 2.2

Let $\psi_{\delta}(x,y)$ is a $(\delta,L,\mu)$-model for $f$ w.r.t. $V[y](x)$. Then, after $k$ iterations of Algorithm 2, we have

[TABLE]

[TABLE]

The details of proof can be found in Appendix 0.B. Note that Algorithm 1 also has linear convergence rate for the strongly convex case. The details can be found in Appendix 0.C. The benefit of Algorithm 1 is that there is no need to know the strong convexity parameter $\mu$ for the algorithm to work. On the other hand, this parameter is needed for assessing the quality of the solution returned by the algorithm. The benefit of the adaptive version is that it does not require to know the value of the parameter $L$ and adapts to it. Moreover, the parameter $L$ can be different for the model at different points and the algorithm adapts also for the local value of this parameter.

3 Clustering by Electorial Model

In this section we consider clustering model introduced in [49]. In this model voters (data points) choose a party (cluster) in an iterative manner by alternative minimization of the following function.

[TABLE]

where ${\mathbb{R}}_{+}^{m}$ is a non-negative orthant and $S_{n}(1)$ is the standard $n$-dimensional simplex in ${\mathbb{R}}^{n}$. The vector $z$ contains probabilities with which voters choose the considered party, and vector $p$ describes the position of the party in the space of voter opinions. The minimized potential is the result of combining two optimization problems into one: voters choose the party whose position is closest to their personal opinion and the party adjusts its position minimizing dispersion and trying not to go too far from its initial position. Yu. Nesterov in [49] used sequential elections process to show that under some natural assumptions the process convergence and gives the clustering of the data-points. This was done for a particular choice of the function $g$ which has limited interpretability. We show, how our framework of inexact model of the objective allows to construct a gradient-type method for the case of general function $g$, which is not necessarily convex.

Assume that $g(x)$ (generally, non-convex) is an function with $L_{g}$-Lipschitz continuous gradient:

[TABLE]

and, following [49], the numbers $\mu_{1},\mu_{2}$ are chosen such that $L_{g}\leq\mu_{1}$ and $L_{g}\leq\mu_{2}$.

The norm $\|\cdot\|$ in $S_{n}(1)\times{\mathbb{R}}_{+}^{m}$ is defined as $\|(z,p)\|^{2}=\|z\|^{2}_{1}+\|p\|^{2}_{2}$, where $\|z\|_{1}=\sum\limits_{k=1}^{n}z_{k}$ and $\|p\|_{2}=\sqrt{\sum\limits_{k=1}^{m}p_{k}^{2}}$. This is indeed a norm since, for $x=(z_{x},p_{x})$ and $y=(z_{y},p_{y})$ we have:

[TABLE]

[TABLE]

because $\sqrt{(a+b)^{2}+(c+d)^{2}}\leq\sqrt{a^{2}+c^{2}}+\sqrt{b^{2}+d^{2}}$ for each $a,b,c,d\geq 0$.

Let us show that

[TABLE]

[TABLE]

is a $(0,2L_{g})$-model of $f_{\mu_{1},\mu_{2}}(x)$ in $x$ with respect to the following Bregman divergence

[TABLE]

It is easy to see that $\psi_{\delta}(x,x)=0$. Let us show, that inequality (2) holds for $\psi_{\delta}(x,y)$. For the function $g(x)$ satisfying (12) we have:

[TABLE]

It means that $f_{\mu_{1},\mu_{2}}(x)-f_{\mu_{1},\mu_{2}}(y)-\psi_{\delta}(x,y)=$

[TABLE]

[TABLE]

Along with (13) and $KL(z_{x}|z_{y})\geq\frac{\|z_{x}-z_{y}\|_{1}^{2}}{2}$ it leads to

[TABLE]

Finally, by definition of the norm $\|\cdot\|$, we have

[TABLE]

i.e. $\psi_{\delta}(x,y)$ is a $(0,2L_{g})$-model of the function $f_{\mu_{1},\mu_{2}}$.

Further, for the case $\min\{\mu_{1},\mu_{2}\}>L_{g}$ $\psi_{\delta}(x,y)$ is a strongly convex w.r.t. $V[y](x)$:

[TABLE]

[TABLE]

where

[TABLE]

is linear in $y$. Indeed,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus, $\psi_{\delta}^{lin}(x,y)$ is a $(0,\max\{\mu_{1},\mu_{2}\}+L_{g},\min\{\mu_{1},\mu_{2}\}-L_{g})$-model of the function $f_{\mu_{1},\mu_{2}}$:

[TABLE]

and

[TABLE]

So, we can apply our Algorithms 1 and 2 to the problem (11).

4 Proximal Sinkhorn Algorithm for Optimal Transport

In this section we consider the problem of approximating an optimal transport (OT) distance. Recently optimal transport distances has gained a lot of interest in machine learning and statistical applications [4, 8, 18, 33, 40, 54, 59]. To state the OT problem, assume that we are given two discrete probability measures $p,q\in S_{n}(1)$ and ground cost matrix $C\in{\mathbb{R}}_{+}^{n\times n}$, then the optimal transport problem is

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes Frobenius dot product of matrices, $\pi$ is a transportation plan. The above optimal transport problem is the Kantorovich [37] linear program (LP) formulation of the problem, which goes back to the Monge’s problem [45]. The best known theoretical complexity for this linear program is 111Here and below for all (large) $n$: $\widetilde{O}(g(n))\leq\tilde{C}\cdot(\ln n)^{r}g(n)$ with some constants $\tilde{C}>0$ and $r\geq 0$. Typically, $r=1$, but not in this particular case. If $r=0$, then $\widetilde{O}(\cdot)=O(\cdot)$. $\widetilde{O}(n^{2.5})$, see [42]. However, there is no known practical implementation of this algorithm. In practice, the simplex method gives complexity $O(n^{3}\ln n)$ [53]. We follow the alternative approach based on entropic regularization of the OT problem [15]. We show how our general framework of inexact model of the objective allows to construct Proximal Sinkhorn algorithm with better computational stability in comparison with the standard Sinkhorn algorithm.

For any optimization problem (1), $\psi_{\delta}(x,y)=f(x)-f(y)$ satisfies Definition 1 with any $L\geq 0$. In this case, our Algorithm 1 becomes inexact Bregman proximal gradient method

[TABLE]

Our idea is to apply this proximal method for the OT problem and approximately find the next iterate $x^{k+1}$ by Sinkhorn’s algorithm [58, 15, 2, 29]. The latter is made possible by the choice of $V$ as KL divergence, which makes the problem of finding the point $x^{k+1}$ to be an entropy-regularized OT problem, which, in turn, is efficiently solvable by the Sinkhorn algorithm.

Consider the iterates

[TABLE]

which we call outer iterations. On each outer iteration we use Sinkhorn’s algorithm 3, which solves the minimization problem in (4) with accuracy $\tilde{\varepsilon}$ in terms of its objective residual. Notice that unlike [29] we provide a slightly refined theoretical bounds for the Sinkhorn’s algorithm not depending on vectors $p$, $q$.

Theorem 4.1

Let $\bar{\pi}^{N}=\frac{1}{N}\sum_{k=1}^{N}\pi^{k}$, where $\pi^{k}$ are the iterates of (4). Then, after $N=\frac{4L\ln n}{\varepsilon}$ iterations, it holds that $\langle C,\bar{\pi}^{N}\rangle\leq\min_{\pi\in\mathcal{U}(p,q)}\langle C,\pi\rangle+\varepsilon$. Moreover, the accuracy $\tilde{\varepsilon}$ for the solution of (4) is sufficient to be set as \widetilde{O}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varepsilon^{4}/(Ln^{4})}) and the complexity of Sinkhorn’s Algorithm on $k$-th iteration is bounded as

[TABLE]

*where222This bound is rough and typically $\bar{c}_{k}$ is smaller in practice. By proper rounding of $\pi^{k}$ one can guarantee (without loss of generality) that \pi^{k}_{ij}\geq\varepsilon/(2n^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left\lVert C\right\rVert_{\infty}}), which gives

{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{\bar{c}_{k}}{L}=\frac{\left\lVert C\right\rVert_{\infty}}{L}+\ln\left(\frac{2n^{2}\left\lVert C\right\rVert_{\infty}}{\varepsilon}\right).}

But, in practice there often is no need to make ‘rounding’ after each outer iteration. *

[TABLE]

Proof

The estimate for the number of iterations $N$ follows from Theorem 2.1 since $V[\pi_{0}](\pi_{*})\leq\ln n^{2}$ as $\pi\in S_{n^{2}}(1)$. The first component of (24) is proved in [34], and the second component basically follows from [5, 29]. Proofs of the second component and bound on $\bar{c}_{k}$ (25) are provided in Appendix 0.E (Theorem 0.E.1). Let us show that it is sufficient to solve minimization problem (4) on each iteration with accuracy \tilde{\varepsilon}=\widetilde{O}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varepsilon^{4}/(Ln^{4})}) in terms of the objective residual to guarantee $\tilde{\delta}=\varepsilon/2$ accuracy in terms of Definition 2.

To prove this fact, we use relation (49) in Theorem 0.G.1 of Appendix 0.G with $\|\cdot\|=\|\cdot\|_{1}$, $\tilde{R}=2$, $\mu=L$. To bound $\Delta=\tilde{L}\widetilde{R}+\|\nabla\phi(\tilde{x}^{*})\|_{*}$ (in notations of Theorem 0.G.1) we modify $\mathcal{U}(p,q)$ by adding constraints: \pi_{ij}\geq\varepsilon/(4{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}n^{2}}), $i,j=1,...,n$. The solution of the changed problem is still an $O(\varepsilon)$-solution of the original problem. For the modified problem \Delta=5L{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}n^{2}}\tilde{R}/\varepsilon. According to (49) one should solve auxiliary problem with accuracy by function value $\widetilde{\varepsilon}$, which is chosen such that \varepsilon/2=\tilde{\delta}=(5L{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}n^{2}}/\varepsilon)\tilde{R}\sqrt{2\widetilde{\varepsilon}/L}. The only problem is that now we cannot directly apply Sinkhorn’s algorithm. This problem can be solved by trivial affine transformation of $\pi$-space. This transformation reduces modified polyhedral to the standard one and we can use Sinkhorn’s Algorithm. Such a transformation doesn’t change (in terms of $O(~{})$) the requirements to the accuracy.

Figure 1: Adaptive choice of L

Remark 2

The standard Sinkhorn’s method can be seen as a particular case of our algorithm (4) with only one step. To obtain an $\varepsilon$-approximate solution of (21), the regularization parameter $L$ needs to be chosen $O\left(\varepsilon/\ln n\right)$ [2, 29, 36]. This can lead to instability of the Sinkhorn’s algorithm [57]. On the opposite, our Proximal Sinkhorn algorithm allows to run Sinkhorn’s algorithm with larger regularization parameter. This parameter can be chosen by minimization of the theoretical bound (24), which gives $L=\widetilde{O}(\|C\|_{\infty})$. In practice one can choose this constant adaptively since we have a $(\delta,L)$-model for any $L$ and can vary $L$ from iteration to iteration. First, the inner problem (4) is solved with overestimated $L$. Then, we set $L:=L/2$ and the problem is solved with the updated value of the parameter and so on until a significant increase (e.g. 10 times) in the complexity of the auxiliary entropy-linear programming problem in comparison with the initial complexity is detected, see Figure 1, where $N(L)$ is a number of required iterations of Sinkhorn algorithm to solve the inner problem with accuracy $\varepsilon$.

From the Theorem 4.1 and Remark 2 one can roughly estimate the total complexity of Proximal Sinkhorn algorithm as333Our experiments on MNIST data set show (see Figures 2, 3, 7) that in practice the bound is better. {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{O}(n^{4}/\varepsilon^{2})}. We also mention several recent complexity bounds444Strictly speaking for the moment we can not verify all the details of the proof of estimate $\tilde{O}(n^{2}/\varepsilon)$. Also the proposed in [9, 56] methods are mainly theoretical, like Lee–Sidford’s method for OT problem with the complexity $\tilde{O}(n^{2.5})$ [42]. For the moment it is hardly possible to implement these methods such that theirs practical efficiencies correspond to the theoretical ones. for the OT problem $\tilde{O}(n^{2}/\varepsilon^{3})$ [2], $\tilde{O}(n^{2}/\varepsilon^{2})$ and $\tilde{O}(n^{2.5}/\varepsilon)$ [29], $\tilde{O}(n^{2}/\varepsilon)$ [9, 56], $\tilde{O}(n/\varepsilon^{3+d})$, $d\geq 1$ [1].

4.1 Numerical Illustration

In this subsection we provide numerical illustration of the Proximal Sinkhorn algorithm.555The code is available at https://github.com/dmivilensky/Proximal-Sinkhorn-algorithm In the experiments we use a standard MNIST dataset with images scaled to a size $10\times 10$. The vectors $p$ and $q$ contain the pixel intensities of the first and second images respectively. The value of $c_{ij}$ is equal to the Euclidean distance between the $i$-th pixel from the vector $p$ and the $j$-th pixel from the vector $q$ on the image pixel grid. For experiments with varying number of pixels $n$ the images are resized to be images of $10\cdot m\times 10\cdot m$ pixels, where $m\in\mathbb{N}$. We replace all the zero elements in $p$ and $q$ with $10^{-3}$ and, then, normalize these vectors.

Fig. 2 shows that the growth rate of the iteration number with increasing accuracy or size of the problem for the Sinkhorn’s algorithm is greater than for the Proximal Sinkhorn’s method. At the same time, with a higher value of $L$ in proximal method, the iteration number is greater, and the growth rates with some precision are equal. The same type of dependence on the accuracy and the size of the problem can be seen for the working time (fig. 3):

5 Proximal IBP Algorithm for Wasserstein Barycenter

In this section we consider a more complicated problem of approximating an OT barycenter. OT barycenter is a natural definition of a mean in a space endowed with an OT distance. Such barycenters are used in the analysis of data with geometric structure, e.g. images, and other machine learning applications [16, 7, 38, 54, 39]. For a set of probability measures $\{p_{1},\dots,p_{m}\}$, cost matrices $C_{1},\dots,C_{m}\in{\mathbb{R}}_{+}^{n\times n}$, and $w\in S_{n}(1)$, the weighted barycenter of these measures is defined as a solution of the following convex optimization problem

[TABLE]

[TABLE]

The idea is similar to the one in Sect. 4, namely, we use our framework to define a Proximal Iterative Bregman Projections algorithm. The algorithm starts from the point $\boldsymbol{\pi}$ s.t. \pi^{0}_{l}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{n}p_{l}\mathbbm{1}^{T}\in\mathcal{U}(p_{l},\mathbbm{1}/n)}, $l=1,...,m$ and iterates

[TABLE]

These iterations are called outer iterations and on each such iteration, the Iterative Bregman Projections algorithm [7] listed as Algorithm 4 below is used to solve the auxiliary minimization problem.

Theorem 5.1

Let $\bar{\boldsymbol{\pi}}^{N}=\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{\pi}^{k}$, where $\boldsymbol{\pi}^{k}$ are the iterates of (5). Then, after $N=\frac{4Lm\ln n}{\varepsilon}$ iterations, it holds that

[TABLE]

Moreover, the accuracy $\tilde{\varepsilon}$ for the solution of (5) is sufficient to be set as $\tilde{\varepsilon}=\widetilde{O}(\varepsilon^{2}/(mn^{3}))$ and the complexity of IBP on $k$-th iteration is bounded as

[TABLE]

The proof of Theorem 5.1 is based on Theorem 2.1 and [38]. All the remarks from Section 4 for Proximal Sinkhorn algorithm also hold for Proximal IBP. In [38] it was shown that complexity of IBP is $\widetilde{O}\left(n^{2}/\varepsilon^{2}\right)$. Despite the theoretical complexity of Proximal IBP is worse than this bound, we show in the next section that in practice Proximal IBP beats the standard IBP algorithm. As an alternative to the IBP algorithm we mention primal-dual accelerated gradient descent [23, 63].

5.1 Numerical Illustration

In this section, we present preliminary computational results for the numerical performance analysis of the Proximal Iterative Bregman Projection (ProxIBP) method discussed above asthe iterates (5).

Initially, we show the results for the computation of a non-regularized Wasserstein barycenter of a set of $10$ truncated Gaussian distributions with finite support. For the finite support $x=[-5,-4.9,-4.8,\ldots,-0.1,0,0.1,\dots,4.8,4.9,5]$, we set the finite distribution $p_{l}$ such that $p_{l}(i)=\mathcal{N}(x_{i};\mu_{i},\sigma_{i})$, that is, the value at coordinate $i$ of the distribution $p_{l}$, for $1\leq l\leq m$, is the value of the Normal distribution with mean $\mu_{i}$ and standard deviation $\sigma_{i}$. The values $\{\mu_{i}\}\sim\text{Uniform}[-5,5]$, are uniformly chosen in the line segment $[-5,5]$, and the values are selected as $\{\mu_{i}\}\sim\text{Uniform}[0.25,1.25]$. For simplicity of exposition, we select uniform weighting for all distributions, i.e., $w_{l}=1/m$.

Figure 4 shows the numerical results for a number of comparative scenarios between the Iterative Bregman Projection (IBP) algorithm proposed in [7] and its Proximal variant in (5). For both algorithms, we show the function values achieved by the generated iterates, and the final approximated barycenter. The results for the IBP algorithm are shown in Figure 4(a) and Figure 4(b). Figure 4(a) shows the weighted distance between the generated barycenter and the original distributions for three different desired accuracy values. It is clear that a bigger $\varepsilon$ generates a faster convergence, but the final cost is slightly higher than in other cases. Figure 4(b) shows the resulting barycenter for the three values of the accuracy parameter. For higher accuracy, the effects of the regularization constant are smaller and thus we obtain a “spikier” barycenter. Figure 4(c) and Figure 4(d) shows a similar analysis for the proposed Proximal IBP in (5), in Figure 4(c) we observe the function value of the generated barycenter, for a fixed number of inner loop iterations, and changing values of $L$, note that here $L$ is not a regularization parameter but the weight on the Bregman function. For larger values of $L$, the inner loop problem is easier to solve, requires less iterations to achieve certain accuracy, with the price in a larger number of iterations in the outer loop. For the particular problem studied, $200$ iterations in the outer loop are sufficient to achieve good performance even with relavively smaller values of $L$. Figure 4(c) shows the generated barycenters for the Proximal IBP algorithm. Finally, Figure 4(e) and Figure 4(f) show the results, for the analogous adaptive stopping condition described in Line $11$ of Algorithm 3 with $\varepsilon=1\cdot 10^{-10}$. We test two different values of the parameter $L$, namely $1$ and $0.1$. Additionally, we explore the suggested adaptive search procedure, where one decreases the value of the parameter $L$ at each iteration, until the inner problem has become particularly hard to solve. This last approach is shown a fast convergence as it reaches a comparable value in around $10$ iterations. Figure 4(f) shows the resulting barycenters.

Figure 5 shows the result of applying the proximal IBP algorithm to the computation of the barycenter of $20$ images of the number $7$ from the MNIST dataset [41]. As shown in Figure 5(b), the Euclidean mean among the images does not preserve the geometric properties of the images, making the optimal transport distances suitable for this particular problem. Figure 5(a) shows, from left to right, and top to bottom, the sequence of generated barycenters for some of the initial iterations of the algorithm. It is evident, that the generated barycenter maintains the geometric structure of the images, as the barycenter is itself a prototypical image of the number $7$. Figure 5(c) shows the average distance of the generated barycenter for the first $400$ iterations of the algorithm.

Figure 6 shows the influence of the weights in the computation of the Wasserstain barycenter of four images via the Prox IBP algorithm. We have used four images corresponding to the digits [math], $1$, $2$, and $3$, and positioned each one of them at the corner of a square. Each of the images shown inside the square corresponds to the obtained barycenter whith weights proportional to the distance to each of the corners. For example, the upper edge of the square assumes zero weights to the digits $2$ and $3$ and only shows the changes in the weights between [math] and $1$. The center image corresponds to equal weights to all images. The other images are generated similarly.

Acknowledgments. The work in sections 4 and 5 was funded by Russian Science Foundation (project 18-71-10108). The work in section 3 was supported by RFBR 18-31-20005 mol$\_$a$\_$ved. The work of D. Dvinskikh and D. Pasechnyk partially conducted in Sochi-Sirius, July, 2018.

Appendix 0.A Adaptive gradient method and proof of Theorem 2.1

Before we prove the theorem we should note that that we can reduce Algorithm 5 to Algorithm 1. Indeed, let us always choose $L_{k+1}=L$ instead of $L_{k+1}=2^{i_{k}-1}L_{k}$ in Algorithm 5. In this case it is guaranteed that we exit the inner loop in Algorithm 5 after the first step due to $(\delta,L)$-model definition. Moreover, with this choice of $L_{k+1}$ Algorithm 5 generates the same sequences as in Algorithm 1. We prove Thus, Theorem 2.1 is a corollary of Theorem 0.A.1.

Theorem 0.A.1

Let $V[x_{0}](x_{*})\leq R^{2}$, where $x_{0}$ is the starting point, and $x_{*}$ is the nearest minimum point to the point $x_{0}$ in the sense of Bregman divergence $V[y](x)$. Then, for the sequence, generated by Algorithm 5 the following inequality holds:

[TABLE]

Moreover, for Algorithm 5 the total number of attempts to solve (32) is bounded by $2N+\log_{2}\frac{L}{L_{0}}$.

First, we need to prove two lemmas in order to obtain the final result. Let us prove Lemma 1.

Lemma 1

Let $\psi(x)$ be a convex function and

[TABLE]

where $\beta\geq 0$. Then

[TABLE]

Proof

Using Definition 2, we have:

[TABLE]

Then inequality

[TABLE]

and equality

[TABLE]

complete the proof.

Lemma 2

We have the following inequality:

[TABLE]

Proof

From the stopping criterion (31), we have:

[TABLE]

Using Lemma 1 with $\psi(x)=\psi_{\delta}(x,x_{k})$ and $\beta=L_{k+1}$, we obtain:

[TABLE]

In view of the model definition (2), we have:

[TABLE]

Remark 3

Let us show that $L_{k}\leq 2L\quad\forall k\geq 0$. For $k=0$ this is true from the fact that $L_{0}\leq L$. For $k\geq 1$ this follows from the fact that we leave the inner cycle earlier than $L_{k}$ will be greater than $2L$. The exit from the cycle is guaranteed by the condition that there is an $(\delta,L)$-model for $f(x)$ at any point $x\in Q$.

Finally, we prove the theorem.

Proof

Let us sum up the inequality from Lemma 2 from [math] to $N-1$:

[TABLE]

Since $V[x_{N}](x_{*})\geq 0$ and $V[x_{0}](x_{*})\leq R^{2}$, we obtain inequality

[TABLE]

Let us divide both parts by $S_{N}$.

[TABLE]

Using the convexity of $f(x)$ we can show that

[TABLE]

Remains only to prove that

[TABLE]

As it follows from Definition 1 and Remark 3 for all $k\geq 0$ $L_{k}\leq 2L$. Thus, we have that

[TABLE]

and

[TABLE]

The total number of attempts to solve (32) is bounded in the same way as in [47].

In the same way as it is done in [3, 13, 25, 30, 50], one can show that the proposed method is primal-dual.

Appendix 0.B Analysis of Algorithm 2 in the case of $(\delta,L,\mu)$-model

Now we consider a proof of Theorem 2.2. To analyze Algorithm 2, assume that it works for $k$ iterations. By Lemma 1, for each $x\in Q$:

[TABLE]

It means that

[TABLE]

[TABLE]

Further, $\psi_{\delta}(x,y)$ is a $(\delta,L)$-model w.r.t. $V[y](x)$ and from

[TABLE]

we get

[TABLE]

Now (34) means

[TABLE]

Since $\psi_{\delta}(x,y)$ is a $(\delta,L,\mu)$-model for $f$, we have:

[TABLE]

Considering (35), we obtain:

[TABLE]

Set $x=x_{*}$. Since $L_{0}\leq 2L$, we have $L_{k+1}\leq 2L$ for each $k\geq 0$. We also assume in Algorithm 2 that $L_{k+1}\geq\mu$. Thus, we have

[TABLE]

Then we have $\forall i\in\mathbb{N}:i<k$

[TABLE]

Therefore, we obtain:

[TABLE]

and

[TABLE]

For further reasoning we introduce average parameter $\hat{L}$:

[TABLE]

Note that by $L_{i}\geq\mu~{}(i=1,2,\ldots)$

[TABLE]

Now, taking into account (37), we have:

[TABLE]

[TABLE]

Finally, we have

[TABLE]

and by (38) – (39) and $L^{k+1}\leq 2L$ means:

[TABLE]

Appendix 0.C Analysis of Algorithm 1 in the case of $(\delta,L,\mu)$-model

Theorem 0.C.1

Let $\psi_{\delta}(x,y)$ be a $(\delta,L,\mu)$-model for $f$ w.r.t. $V[y](x)$ and $y_{k}=\operatorname*{argmin}_{i=1,...,k}(f(x_{i}))$. Then, after $k$ iterations of Algorithm 1, we have

[TABLE]

and

[TABLE]

Proof

Clearly, $f(x_{*})\leq f(x^{k+1})$ and

[TABLE]

i.e.

[TABLE]

Further,

[TABLE]

Therefore, taking into account the following fact

[TABLE]

we obtain

[TABLE]

Now we consider the question on convergence by function:

[TABLE]

Therefore, we have

[TABLE]

Denote by $y_{k}=\operatorname*{argmin}_{i=1,...,k}(f(x_{i}))$. Then, taking into account

[TABLE]

we obtain

[TABLE]

[TABLE]

Appendix 0.D Some Numerical Tests for Algorithms 1 and 2

We consider two numerical examples for Algorithms 1 and 2 for minimizing $\mu$-strongly convex objective function of $N$ variables on a unit ball $B_{1}(0)$ with center at zero with respect to the standard Euclidean norm. It is clear that such functions admit $(\delta,L,\mu)$-model of the standard form $\psi_{\delta}(x,y)=\langle\nabla f(y),x-y\rangle$ for the case of Lipschitz-continuous gradient $\nabla f$. In the first of the considered examples, it is easy to estimate $L$ and $\mu$, and the ratio $\frac{\mu}{L}$ is not very small, which ensures a completely acceptable rate of convergence of the non-adaptive method (see Table 1 below). In the second example, the objective is ill-conditioned meaning that the ratio $\frac{\mu}{L}$ so small that the computer considers the value $1-\frac{\mu}{L}$ to be equal to 1 and Theorem 0.C.1 for the non-adaptive algorithm does not allow to estimate the rate of convergence at all. In this case, the use of adaptive Algorithm 2 leads to noticeable results (see the Table 2 below).

Example 1

Consider a function

[TABLE]

where $N=100$ and input data

[TABLE]

$\mu=2$, $L^{0}=2\mu$, $L=2N$.

The results of the comparison of the work of algorithms 1 and 3 are presented in the comparative Table 1, where $k$ is the number of iterations of these algorithms.

As we can see from the Table 1, in the previous example the non-adaptive method converges no worse than the adaptive one. However, it is possible that $\frac{\mu}{L}$ is too small, which leads to $1-\frac{\mu}{L}\approx 1$. In this case, Theorem 0.C.1 cannot estimate the rate of convergence of the method. We give another example.

Example 2

Consider the target functional

[TABLE]

It is easy to verify that for such a function one can choose $\mu=2+\frac{1}{e}$ and $L=2N+N^{2}e$ and the program calculates the value of $1-\frac{\mu}{L}$ equal to 1. However, applying Algorithm 2 with adaptive tuning to the constant $L$ and Theorem 2.2 we obtain meaningful results, which we present in Table 2.

Experiments were performed using CPython 3.7 software on a computer with a 3-core AMD Athlon II X3 450 processor with a clock frequency of 803.5 MHz per core. The computer’s RAM was 8 GB.

Appendix 0.E Complexity Analysis of Sinkhorn’s Algorithm

Let us consider regularized optimal transport problem

[TABLE]

Recall that the dual problem to (46) is equivalent to

[TABLE]

where $B(u,v):=\operatorname{diag}(e^{u})e^{-C/\gamma}\operatorname{diag}(e^{v})$ [29]. Below we present a slightly refined complexity analysis of Algorithm 3 based on the same approach as in [29]. First, we prove that Sinkhorn’s iterations are contractant for $e^{u^{t}-u^{*}}$ and $e^{v^{t}-v^{*}}$ in Hilbert’s projective metric (cf. [34]).

Lemma 3

Let us define

[TABLE]

where $(u^{*},v^{*})$ is the solution of problem (47). Then for any $t\geq 0$ it holds $R_{t+1}\leq R_{t}$.

Proof

W.l.o.g. consider even $t$. Let us denote $\pi^{*}=B(u^{*},v^{*})\in\mathcal{U}(p,q)$. Then for any $i$

[TABLE]

and since $\sum_{j}\frac{\pi^{*}_{ij}}{p_{i}}=1$ one obtains

[TABLE]

therefore,

[TABLE]

Now repeating the proof of Theorem 1 from [29] we obtain the following complexity bound.

Theorem 0.E.1

The inner cycle of Algorithm 3 stops in number of iterations

[TABLE]

and

[TABLE]

Nptice that now we require an approximated solution of regularized problem (46), thus the choice of $\varepsilon^{\prime}$ in Algorithm 3 differs from the one from [2, 29].

Theorem 0.E.2

Algorithm 3 returns $\hat{\pi}\in\mathcal{U}(p,q)$ s.t.

[TABLE]

where $\pi^{*}$ is the solution of problem (46).

Proof

Notice that for any $\pi\in\mathcal{U}(B(u,v)\mathbbm{1},B(u,v)^{T}\mathbbm{1})$ it holds

[TABLE]

It is easy to see that for any pair $\pi,\tilde{\pi}\in S_{n\times n}(1)$

[TABLE]

thus

[TABLE]

Now, for any $\pi\in S_{n\times n}(1)$ and $r,c\in S_{n}(1)$ there exists $\tilde{\pi}\in\mathcal{U}(r,c)$ given by Algorithm 2 from [2] s.t. $\left\lVert\pi-\tilde{\pi}\right\rVert_{1}\leq\left\lVert\pi\mathbbm{1}-r\right\rVert_{1}+\left\lVert\pi^{T}\mathbbm{1}-c\right\rVert_{1}$. Combining all these facts together we obtain for $\hat{\pi}$ defined in Algorithm 3 the following estimate

[TABLE]

Substituting

[TABLE]

we obtain

[TABLE]

Appendix 0.F Additional Experiments for Prox-Sinkhorn Algorithm

Figure 7 shows the dependence of the mean inner method iteration number upon accuracy and size of the vector $\pi$. With the growth of $L$, there is a decrease in the mean inner method iteration number. However, the type of dependence from the accuracy or size of the problem is the same.

Consider a graph of change of the auxiliary problem solution complexity with increasing external method iteration number (fig. 8). Note that at the first interval there is an increase in the inner method iteration number with two peaks at different levels. On subsequent iterations of the external method the complexity of the solution of the auxiliary problem decreases, approaching a constant.

Appendix 0.G On Inexact Solution of Auxiliary Subproblems

Our goal is to provide a relation between the accuracy of the solution of an optimization problem in terms of the objective residual and $\widetilde{\delta}$-‘precision’ in the sense of Definition 2.

Theorem 0.G.1

Assume that we find a point $\tilde{x}$ such that $\phi(\tilde{x})-\phi(\tilde{x}^{*})\leq\widetilde{\varepsilon}$, where $\tilde{x}^{*}$ is an exact solution of problem (4). Assume also that $\phi$ has $\tilde{L}$-Lipschitz continuous gradient in $Q$.

If $\nabla\phi(\widetilde{x}^{*})=0$, then \tilde{x}=\operatorname*{argmin}_{x\in Q}^{\widetilde{\delta}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}(x) with $\widetilde{\delta}=\widetilde{R}\sqrt{2\tilde{L}\tilde{\varepsilon}}$, where $\widetilde{R}=\max_{x,y\in Q}\|y-x\|$.

*If $\phi$ is $\mu$-strongly convex on $Q$, then \tilde{x}=\operatorname*{argmin}_{x\in Q}^{\widetilde{\delta}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}(x) with *

[TABLE]

Proof

1. Assume that $\nabla\phi(\tilde{x}^{*})=0$. Then

[TABLE]

[TABLE]

2. Let us now assume that $\nabla\phi(\tilde{x}^{*})\neq 0$. For strongly convex function $\phi(x)$ we have

[TABLE]

Hence,

[TABLE]

Using this and Lipschitz gradient condition, we obtain

[TABLE]

Hence,

[TABLE]