New Version of Mirror Prox for Variational Inequalities with Adaptation to Inexactness

TL;DR

This paper introduces an adaptive Mirror Prox method for variational inequalities that adjusts to both Lipschitz constant and oracle error, achieving near-optimal complexity and demonstrating effectiveness through experiments.

Contribution

It presents a novel adaptive Mirror Prox algorithm that accounts for inexact oracle information and adapts to problem parameters, improving convergence guarantees.

Findings

Achieves near $O(1/\varepsilon \log(1/\varepsilon))$ complexity for certain monotone operators.

Demonstrates superior performance in experiments compared to existing methods.

Shows effectiveness in matrix games with non-Euclidean setups.

Abstract

Some adaptive analogue of the Mirror Prox method for variational inequalities is proposed. In this work we consider the adaptation not only to the value of the Lipschitz constant, but also to the magnitude of the oracle error. This approach, in particular, allows us to prove a complexity near $O\left(\frac{1}{\varepsilon}\log_2\frac{1}{\varepsilon}\right)$ for variational inequalities for a special class of monotone bounded operators. This estimate is optimal for variational inequalities with monotone Lipschitz-continuous operators. However, there exists some error, which may be insignificant. The results of experiments on the comparison of the proposed approach with some known analogues are presented. Also, we discuss the results of the experiments for matrix games in the case of using non-Euclidean proximal setup.

Tables

Table 1. Table 1 : The results of Algorithm 1 for objective function ( 27 ).

Iterations	General estimate	Time (sec.)
17	0.1051	0.264
19	0.0527	0.291
21	0.0266	0.315
22	0.0212	0.342
23	0.0177	0.354
24	0.0133	0.364
25	0.0106	0.380
26	0.0082	0.427
27	0.0063	0.467
28	0.0048	0.423
29	0.0044	0.443

Table 2. Table 2 : The results of Algorithm 1 for objective function ( 29 ).

$n=600,m=400,N=25$			$n=1000,m=500,N=50$
Iteration	General estimate	Time (sec.)	Iteration	General estimate	Time (sec.)
22	0.122	67.955	19	0.1343	252.151
23	0.061	70.587	20	0.0672	252.673
24	0.0305	75.107	21	0.0336	266.636
25	0.0153	72.917	22	0.0168	279.883
26	0.0076	76.686	23	0.0084	293.866

Table 3. Table 3 : The results of Algorithm 1 for objective function ( 29 ), with points A k subscript 𝐴 𝑘 A_{k} in the unit ball.

$n=100,m=50,N=25$			$n=200,m=100,N=50$
Iteration	General estimate	Time (sec.)	Iteration	General estimate	Time (sec.)
318	0.2539	35.805	684	0.2522	441.744
468	0.1702	52.809	1016	0.1688	645.082
618	0.1279	71.484	1350	0.1267	857.185
768	0.1026	87.851	1682	0.1015	1049.885
918	0.0857	103.686	2016	0.0846	1305.006
1068	0.0736	123.056	2349	0.0727	1534.333
1218	0.0645	141.044	2683	0.0637	1753.489
1368	0.0575	153.877	3014	0.0567	2026.402
1518	0.0518	173.538	3348	0.0511	2210.362
1668	0.0472	190.714	3680	0.0465	2301.327
1818	0.0434	208.966	4014	0.0427	2611.942
1968	0.0401	228.198	4346	0.0394	2970.324
2118	0.0373	243.970	4678	0.0367	3153.905
2268	0.0349	253.112	5012	0.0343	3387.451
2426	0.0323	266.583	5346	0.0322	3619.831

Equations

$$\Vert g(x)-g(y){\Vert }_{\ast }⩽L_{\nu }\Vert x-y{\Vert }^{\nu }\forall x,y\in Q\text{ for some }\nu \in [0;1],$$

$$⟨g(y)-g(x),y-z⟩⩽LV(y,x)+LV(y,z)+\delta \Vert y-z\Vert \forall x,y,z\in Q,$$

$$||y|{|}_{\ast }=\underset{x}{\max }\{⟨y,x⟩,||x||\le 1\},$$

$$V(x,y)=d(x)-d(y)-⟨\nabla d(y),x-y⟩.$$

$$⟨g(x_{\ast }),x_{\ast }-x⟩⩽0\forall x\in Q.$$

$$⟨g(x)-g(y),x-y⟩\ge 0\forall x,y\in Q,$$

$$⟨g(x),x_{\ast }-x⟩⩽0\forall x\in Q.$$

$$\Vert g(x)-g(y){\Vert }_{\ast }⩽L\Vert x-y\Vert \forall x,y\in Q.$$

$$\Vert \tilde{g}(x)-g(x){\Vert }_{\ast }⩽\frac{\delta }{2}\forall x\in Q.$$

$$|⟨\tilde{g}(y)-\tilde{g}(x),y-z⟩-⟨g(y)-g(x),y-z⟩|=|⟨\tilde{g}(y)-g(y),y-z⟩+⟨g(x)-\tilde{g}(x),y-z⟩|⩽$$

$$⩽\Vert \tilde{g}(y)-g(y){\Vert }_{\ast }\cdot \Vert y-z\Vert +\Vert \tilde{g}(x)-g(x){\Vert }_{\ast }\cdot \Vert y-z\Vert ⩽\left(\frac{\delta }{2}+\frac{\delta }{2}\right)\Vert y-z\Vert =\delta \Vert y-z\Vert .$$

$$⟨\tilde{g}(y)-\tilde{g}(x),y-z⟩⩽⟨g(y)-g(x),y-z⟩+\delta \Vert y-z\Vert ⩽\Vert g(y)-g(x){\Vert }_{\ast }\cdot \Vert y-z\Vert +\delta \Vert y-z\Vert ⩽$$

$$⩽L\Vert y-x\Vert \cdot \Vert y-z\Vert +\delta \Vert y-z\Vert ⩽\frac{L}{2}\Vert y-x{\Vert }^2+\frac{L}{2}\Vert y-z{\Vert }^2+\delta \Vert y-z\Vert ⩽$$

$$⩽LV(y,x)+LV(y,z)+\delta .$$

$$\{y_{t_k}{\}}_{k=1}^{\infty }:t_1<t_2<t_3<\dots \text{ and}\underset{k\to \infty }{\lim }{t}_k=1;$$

$$diam\partial f(y_{t_k})=:{\delta }_k>0,\text{ where }\sum _{k=1}^{+\infty }{\delta }_k=:\delta <+\infty ,$$

$$\underset{\begin{array}{c}\hat{\partial }f(x)\in \partial f(x),\\ \hat{\partial }f(y)\in \partial f(y)\end{array}}{\min }\Vert \hat{\partial }f(x)-\hat{\partial }f(y){\Vert }_{\ast }⩽L\Vert x-y\Vert .$$

$$\Vert g(y)-g(x){\Vert }_{\ast }⩽L\Vert y-x\Vert +\delta .$$

$$⟨g(y)-g(x),y-z⟩⩽\Vert g(y)-g(x){\Vert }_{\ast }\cdot \Vert y-z\Vert ⩽L\Vert y-x\Vert \cdot \Vert y-z\Vert +\delta \Vert y-z\Vert ⩽$$

$$⩽\frac{L}{2}(\Vert y-x{\Vert }^2+\Vert y-z{\Vert }^2)+\delta \Vert y-z\Vert ⩽LV(y,x)+LV(z,y)+\delta \Vert y-z\Vert .$$

$$f(x)=\{\begin{array}{ll}kx& ;0⩽x⩽\frac{1}{2},\\ \left(k+{\sum }_{i=1}^n\frac{\delta }{2^i}\right)x-{\sum }_{i=1}^n\frac{\delta }{2^i}\left(1-\frac{1}{2^i}\right)& ;1-\frac{1}{2^n}<x⩽1-\frac{1}{2^{n+1}},\\ {\lim }_{x\to +1}f(x)& ;x=1.\end{array}$$

$$y^{N+1}:=\arg \underset{x\in Q}{\min }\{⟨g(x^N),x-x^N⟩+L^{N+1}V(x,x^N)\},$$

$$x^{N+1}:=\arg \underset{x\in Q}{\min }\{⟨g(y^{N+1}),x-x^N⟩+L^{N+1}V(x,x^N)\}.$$

$$⟨g(y^{N+1})-g(x^N),y^{N+1}-x^{N+1}⟩\le L^{N+1}V(y^{N+1},x^N)+$$

$$+L^{N+1}V(x^{N+1},y^{N+1})+{\delta }^{N+1}\Vert y^{N+1}-x^{N+1}\Vert ,$$

$$\sum _{k=0}^{N-1}\frac{1}{L^{k+1}}⟨g(y^{k+1}),y^{k+1}-x⟩⩽V(x,x^0)-V(x,x^N)+\sum _{k=0}^{N-1}\frac{{\delta }^{k+1}}{L^{k+1}}\Vert y^{k+1}-x^{k+1}\Vert .$$

$$\left\langle\nabla_{x}V(x,x^{k})\big{|}_{x=x^{k+1}},\,x-x^{k+1}\right\rangle=V(x,x^{k})-V(x,x^{k+1})-V(x^{k+1},x^{k}),$$

$$\left\langle\nabla_{x}V(x,x^{k})\big{|}_{x=y^{k+1}},\,x-y^{k+1}\right\rangle=V(x,x^{k})-V(x,y^{k+1})-V(y^{k+1},x^{k}).$$

$$\left\langle\nabla_{x}\left(\left\langle g(x^{k}),\,x-x^{k}\right\rangle+L^{k+1}V(x,x^{k})\right)\big{|}_{x=y^{k+1}},\ x-y^{k+1}\right\rangle\geqslant 0,$$

$$\left\langle\nabla_{x}\left(\left\langle g(y^{k+1}),\,x-x^{k}\right\rangle+L^{k+1}V(x,x^{k})\right)\big{|}_{x=x^{k+1}},\ x-x^{k+1}\right\rangle\geqslant 0.$$

Full text

1 Introduction

Variational inequalities (VI) and saddle point problems often arise in a variety of important applications [1]. For solving such problems a lot of algorithmic schemes are known (see e.g. [1, 2, 3, 4]). The Mirror Prox method proposed by A. S. Nemirovski [4] is currently one of the most popular of such methods. This method goes back to the well-known extragradient method proposed by G. M. Korpelevich in [5]. At the same time, unlike the standard extragradient method, Mirror Prox allows to effectively solve problems with non-Euclidean norms, as well as with the Holder-continuous operators:

[TABLE]

all notations are explained in Sect. 2 below.

Recently, a universal analogue of the A. S. Nemirovski method [4] was proposed in [6, 7]. The universality is understood as an adaptive adjustment to the optimal smoothness level $\nu$ in (1), as well as the constant value $L_{\nu}>0$. Note that universal gradient method for unconstrained convex optimization problems was proposed by Yu.E. Nesterov [8] (see also Sect. 5 of the textbook [9]). And it is possible to observe the convergence rate of the proposed method, which is typical for the smooth case $\nu=1$ (Lipschitz-continuous operators), for some problems with bounded operators ($L_{0}<+\infty$ and $L_{\nu}=+\infty$ for all $\nu>0$).

This paper is devoted to the modification of Mirror Prox method [6, 7] for the following analogue of the Lipschitz condition for the operator $g$ with constant $L>0$

[TABLE]

where $\delta>0$ is a fixed value and $V(y,x)$ is the Bregman divergence (see Sect. 2 below).

We propose an analogue of аdaptive Mirror Prox method [6, 7]. At the same time, we consider adaptive tuning both for the value of the parameters $L$ and $\delta$. One of the features of the proposed method which are important for applications is the possibility for the value of $\delta$ to reflect the inexactness of operator $g$. In addition, the value of $\delta$ can indicate the degree of discontinuity of the operator $g$. Adaptive tuning to its value can approximate the convergence rate for variational inequalities with bounded operators ($\nu=0$) to the convergence rate for variational inequalities with Lipshitz-continuous operators ($\nu=1$). Effects of this approach can be observed for the universal method but without a theoretical justification for the convergence rate $O\left(\frac{1}{\varepsilon}\right)$ for non-smooth operators [6]. This means that the proposed approach in this paper is an alternative to the universal method.

The contribution of the paper can be summarized as follows:

• An analogue of the Mirror Prox method for variational inequalities with a monotone Lipschitz-continuous operator, which allows for adaptive tuning to the value of the Lipschitz constant $L$, as well as the limit value of the error $\delta$ of the specifying operator $g$, is proposed.

• The applicability of the proposed method to a certain class of variational inequalities with bounded operators ($\nu=0$) is discussed. The rate of convergence $O\left(\frac{1}{\varepsilon}\right)$ of this method is proved with some finite error associated with the non-smoothness of the operator. Thus, some alternative to the universal method has been proposed, but with a clearer theoretical rationale for acceleration.

• The results of numerical experiments for finding the equilibrium in a bilinear matrix game (or VI with Lipschitz-continuous operator) with a bounded error in the definition of the operator are given. A comparison of the quality of the calculated solution is given depending on the number of iterations for the adaptive Mirror Prox method from [7] and the method proposed in this paper with an adaptive setting for the magnitude of the error.

• The results of numerical experiments for a variational inequality with a bounded ($\nu=0$) operator (related to the Fermat-Torricelli-Steiner problem) are presented. These results show that the method due to the proposed adaptation of the non-smoothness error can converge much faster than the optimal lower estimate $O\left(\frac{1}{\varepsilon^{2}}\right)$ for the corresponding class of problems.

2 Problem Statement and Some Examples

Let $(E,\,\|\cdot\|)$ be a normed finite-dimensional vector space and $E^{*}$ be the conjugate space of $E$ with the norm:

[TABLE]

where $\langle y,\,x\rangle$ is the value of the continuous linear functional $y$ at $x\in E$.

Let $Q\subset E$ be a (simple) closed convex set, $d:Q\rightarrow\mathbb{R}$ be a distance generating function (d.g.f.), which is continuously differentiable and $1$-strongly convex with respect to the norm $\|\cdot\|$ and assume that $\min\limits_{x\in Q}d(x)=d(0).$

For all $x,\,y\in Q\subset E$ we consider the corresponding Bregman divergence

[TABLE]

Let $g:Q\to E^{*}$ be a continuous operator. We consider the problem of finding a solution to a variational inequality of the form

[TABLE]

Under the assumption of the monotony of the operator $g$, i.e.

[TABLE]

the inequality (4) is equivalent to the following weak variational inequality

[TABLE]

Assume that the operator $g$ satisfies the condition (2). In this section, we show some examples of problems for which a condition of the form (2) naturally arises. First of all, this is due to the inexactness of the oracle for the operator of a variational inequality. But also the value of $\delta\|y-z\|$ can describe the degree of discontinuity of the operator $g$ (i.e. using considered approach, one can propose an approach to the solution of some VI’s with bounded operators).

Example 1

Let $g:Q\rightarrow\mathbb{R}^{n}$ be a Lipschitz-continuous operator with constant $L>0$, i.e.

[TABLE]

However, suppose that the exact value of the operator $g$ is not available, and only an approximate value of $g(x)$, i.e. $\tilde{g}(x)$, is known:

[TABLE]

Then for each $x,\,y,\,z\in Q$ we have:

[TABLE]

[TABLE]

Therefore,

[TABLE]

[TABLE]

[TABLE]

Example 2

Note that the term $\delta\|y-z\|$ in (2) can describe non-smoothness for the operator $g$ along any fixed vector segment $\{ty+(1-t)x\}_{0\leq t\leq 1}$. In general (if you combine all possible vector segments), on the domain of non-smoothness points there can be an infinite number.

For example, assume that for some subset $Q_{0}\subset Q$ the function $f$ is differentiable at all points of $Q\setminus Q_{0}$ and that for an arbitrary $x\in Q_{0}$ there exists a finite subdifferential $\partial f(x)$ in the sense of convex analysis.

For fixed $x,y\in Q$ with $t\in[0;1]$ we denote $y_{t}:=(1-t)x+ty$.

Definition 1

([10]) Fix $\delta>0$ and $L>0$. We say that the convex function $f:Q\rightarrow\mathbb{R}\;(Q\subset\mathbb{R}^{n})$ has $(\delta,L)$-Lipschitz subgradient ($f\in C_{L,\delta}^{1,1}(Q)$), if:

• (i)

for arbitrary $x,y\in Q$ the function $f$ is differentiable at all points of the set $\{y_{t}\}_{0\leqslant t\leqslant 1}$, with the exception of the sequence (possibly finite)

[TABLE]

• (ii)

for a sequence of points from (6) there exist finite subdifferentials $\{\partial f(y_{t_{k}})\}_{k=1}^{\infty}$ and

[TABLE]

and   $diam\;\partial f(x)=\max\{\|y-z\|_{*}\,\mid\,y,z\in\partial f(x)\};$

• (iii)

if for $x,y\in Q$ the function $f$ is differentiable at each point $y_{t}$, $t\in(0;1)$, then the following inequality holds:

[TABLE]

Indeed, the property of $(\delta,L)$-Lipschitzness for each subgradient $g(x)=\hat{\partial}f(x)$ means that

[TABLE]

To prove (9), it suffices to split the segment $\{y_{t}\}_{0\leqslant t\leqslant 1}$ into the intervals of smoothness and take into account the boundedness of diameters of the subdifferentials at non-smoothness points of $f$.

The inequality (9) means that

[TABLE]

[TABLE]

Let us give a concrete example of a non-smooth functional with a ($\delta$, L)-Lipschitz subgradient with an arbitrarily large Lipschitz constant.

Example 3

We fix some $k>0$, the value $\delta>0$ and consider the piecewise linear function $f:[0;1]\rightarrow\mathbb{R}$ (here $Q=[0;1]\subset\mathbb{R}$) defined as follows

[TABLE]

In this case, $Q_{0}=\{1-\frac{1}{2^{n}}\}_{n=1}^{\infty}$, $\partial f(q_{n})=\left[k+\sum\limits_{i=1}^{n-1}\frac{\delta}{2^{i}};k+\sum\limits_{i=1}^{n}\frac{\delta}{2^{i}}\right]$ with $n>1$, $\partial f(q_{1})=\left[k;k+\frac{\delta}{2}\right]$ (here $q_{n}=1-\frac{1}{2^{n}}$ with $n=1,2,3,\ldots$). It is clear that $\partial f(q_{n})=\frac{\delta}{2^{n}}$, which is true for the entered value $\delta>0$. Moreover, on the intervals $(0;q_{1}),\;(q_{n};q_{n+1})$ the function $f$ has a Lipschitz-continuous gradient with the constant $L=0$. Therefore, for the function $f$ from (10), we find that $f\in C_{0,\delta}^{1,1}(Q)$.

Any functional with a finite set of non-smooth points along an arbitrary segment will satisfy the proposed Lipschitz condition for the subgradient. Obviously, this condition holds for each objective function with finite points of non-smoothnes on each vector segment $[x;y]$. Thus, it is possible to apply this technique to problems of minimization for sum distances to several balls in Hilbert spaces [11]. Such an objective function, obviously, will not be differentiable in the usual sense at the points of the boundaries of the balls of which there are infinitely many. Note that among points of each vector segment $[x;y]$ such an objective function have finite points of non-smoothness. However, the considered Lipschitz condition for a special choice of subgradient holds for some functions with infinitely many points of non-smoothness (e.g. for maximum of linear functions).

3 Adaptive Method for Variational Inequalities with Adaptation to Inexactness

In this section, we introduce a new version of the Mirror Prox method for variational inequalities (see [7]), which we call Mirror Prox with Adaptation to Inexactness (MPAI). In this version, which is listed as Algorithm 1 below, we consider the adaptation not only to the level of operator smoothness, but also to the magnitude of the oracle error, which may allow to receive complexity near $O\left(\frac{1}{\varepsilon}\right)$ for VI with bounded operators, i.e. the optimal complexity for VI with Lipshitz-continuous operators.

We evaluate the solution quality of the problem (4), produced by Algorithm 1, by using the Bregman divergence (3).

Theorem 3.1

After $N$ iterations of Algorithm 1, the following estimate holds:

[TABLE]

Proof

One can directly check the following inequalities:

[TABLE]

[TABLE]

Further, for each $x\in Q$ and $k=\overline{0,\ N-1}$:

[TABLE]

[TABLE]

Thus,

[TABLE]

and

[TABLE]

Taking into account (13), we have for each $k=\overline{0,\ N-1}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

i.e.

[TABLE]

After summing (16) by $k=\overline{0,\ N-1}$, we have

[TABLE]

$\blacksquare$

Let us denote

[TABLE]

Theorem 3.2

For monotone operator $g$ after $N$ iterations of Algorithm 1, the following estimate holds:

[TABLE]

Assume that for fixed $\varepsilon$

[TABLE]

Then the following inequality holds:

[TABLE]

If $L^{0}\leqslant 2L$, then inequality (18) holds at no more than

[TABLE]

iterations of Algorithm 1.

Proof

By monotony of $g$ we have for each $k=0,1,...$:

[TABLE]

so the inequality

[TABLE]

can be replaced by

[TABLE]

Remark 1

Due to adaptive choice of parameters $L^{k+1}$ and $\delta^{k+1}$ at each iteration of Algorithm 1 the expression

[TABLE]

in (21) may be small enough even in the case of $L=+\infty$ or $\delta=+\infty$ in (2).

Remark 2

Clearly, for each $k$, we have $\delta_{k}\leqslant C_{L}\delta$ ($C_{L}=\max\left\{1,\frac{2L}{L^{0}}\right\}$) and:

[TABLE]

This means that the value associated with the error in the specifying operator $g$ is bounded on the set $Q$ of a finite diameter.

Remark 3

If $g\not\equiv 0$, then the condition $L^{0}\leqslant 2L$ can be satisfied by choosing

[TABLE]

Remark 4

Note that the estimate of the number of iterations $N=\left\lceil\frac{2LR^{2}}{\varepsilon}\right\rceil$ with accuracy to a numerical factor is optimal for variational inequalities with a Lipschitz-continuous operator [12]. Note that the evaluation of the inexactness of the value of the operator, as we see from the previous remark, is bounded and does not accumulate.

Note that similarly Remark 4 in [13] the total number of attempts to solve (11) and (12) is bounded by $4N+\max\left\{\log_{2}\frac{2L}{L^{0}},~{}\log_{2}\frac{2\delta}{\delta^{0}}\right\}$.

Remark 5

It was pointed earlier that $\delta$ can describe the degree of non-smoothness of the operator $g$. Let us show how Algorithm 1 can be slightly modified in order to provably obtain an approximate complexity estimate

[TABLE]

to achieve the quality of the solution $\widetilde{y}$:

[TABLE]

in the case of a bounded operator $g$ (generally speaking, non-smooth).

Suppose that at each iteration $L\leqslant L_{k+1}\leqslant 2L$ (this can always be achieved with a constant number of calculations in item 2 of the listing of Algorithm 1).

We propose such a procedure for some positive integer $p$: repeat the operations item 2 of Algorithm 1, $p$ times, each time

[TABLE]

(for each iteration of the algorithm 1, $k=0,\,1,\,2,\,\ldots,\,N-1$). A procedure of type (22) will be stopped if one of the following two inequalities holds:

[TABLE]

or

[TABLE]

Note that the procedure (22) involves updating $x^{k+1}$ and $y^{k+1}$ during repetitions. Let us estimate $p$ (the number of times required for executing (23) or (24) to repeat a procedure of type (22)).

It is clear that $\forall\,x^{k},x^{k+1},y^{k+1}\in Q$ it is true that

[TABLE]

Moreover, $\delta_{k+1}\leqslant 2\delta$. If (23) is not true, then $\left\|y^{k+1}-x^{k+1}\right\|>\displaystyle\frac{\varepsilon}{4\delta}$ and the inequality (24) is obviously satisfied at

[TABLE]

since in that case

[TABLE]

So, after repeating the $p$ procedures of type (22) at each of the $N$ iterations of the method, the following inequality holds:

[TABLE]

Further,

[TABLE]

Therefore, $\displaystyle\frac{R^{2}}{S_{N}}\leqslant\frac{2^{p+1}LR^{2}}{N}\leqslant\frac{\varepsilon}{2}$ at $N\geqslant\displaystyle\frac{2^{p+2}LR^{2}}{\varepsilon}$, whence, taking into account (25)

[TABLE]

Generally speaking, we need $O\left(\log\displaystyle\frac{1}{\varepsilon}\right)$ additional steps of item 2 of listing 1 of Algorithm 1 at each iteration. So, the final estimate of complexity to achieve the quality of $\displaystyle\max_{x\in Q}\langle g(x),\,\widetilde{y}-x\rangle\leqslant\varepsilon$ will be $O\left(\displaystyle\frac{1}{\varepsilon^{2}}\log_{2}\displaystyle\frac{1}{\varepsilon}\right)$. It is well-known that this estimate is optimal up to a logarithmic factor.

However, for small enough $\delta$ in (26) we have complexity near $O\left(\displaystyle\frac{1}{\varepsilon}\log_{2}\displaystyle\frac{1}{\varepsilon}\right)$.

4 Numerical Experiments for Non-smooth Problem: Variational Inequality for Some Analogue of Fermat-Torricelli-Steiner Problem

In this section, to show the advantages of the proposed Algorithm 1, we consider some numerical experiments for the saddle point problem (and the corresponding VI), which corresponds to the convex programming problem for some analogues of the Fermat-Torricelli-Steiner problem with functional constraints. Note that the objective functions are non-smooth and the corresponding operators of the variational inequality of the problem under consideration are bounded ($\nu=0$). However, experimentally, due to adaptation, we can observe an estimate of the complexity inherent in the case of Lipshitz-continuous operators of VI.

All experiments in this section were implemented in Python 3.4, on a computer equipped with Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz, 1992 Mhz, 4 Core(s), 8 Logical Processor(s). RAM of the computer is 8GB.

For a given set of $N$ points $\{A_{k}=(a_{1k},a_{2k},\ldots,a_{nk});\,k=\overline{1,N}\}$, that represent the centers of the balls $\omega_{k}$ with radii $r_{k}$, in the $n$-dimensional Euclidean space $\mathbb{R}^{n}$, we need to find such a point $X=(x_{1},x_{2},\ldots,x_{n})$ of the objective function [11]

[TABLE]

where

[TABLE]

would take the minimal value on the set $Q$, which is given by several functional constraints:

[TABLE]

where the coefficients $\alpha_{11},\alpha_{12},\ldots,\alpha_{mn}$ are represented by the matrix

[TABLE]

in which one element of each row is an integer belonging to the interval $(1;10)$, and the remaining elements of the row are equal to 1.

To solve such a problem, we can consider a saddle point problem $\min\limits_{x}\max\limits_{\lambda}L(x,\lambda)$, where

[TABLE]

Consider the corresponding variational inequality:

[TABLE]

where

[TABLE]

[TABLE]

We give an example for $n=100$, $m=20$, $N=5$, initial approximation

[TABLE]

and $\delta_{0}=\frac{1}{20}$. The coordinates of the points $A_{k}$ are chosen in such a way that $\left\|A_{k}\right\|\in[1;\,2]$. We choose the standard Euclidean proximal setup as a prox-function.

Note that the centers of the balls were chosen with the norm in the interval $[1;2]$, and the radii of the balls are 1. Therefore, in a single ball with the center at zero there will be points of the boundary of the balls in which the objective function (27) and the operator of the corresponding variational inequality will be bounded. At the same time, it can be shown that the diameter of the subdifferential at such points will be at least 1. It means that theoretically $\delta$ can be at least 1. However, experimentally, we can see significantly better solution quality (see the column "General estimate" in Table 1).

The results of the work of Algorithm 1, for objective function (27) are represented in Table 1 below.

It is known ([12, 14]) that for variational inequalities with bounded operators, the theoretical estimate of the complexity (the convergence rate) $O\left(\frac{1}{\varepsilon^{2}}\right)$ is theoretically optimal. However, experimentally we see from Table 1 that, for example, an accuracy of $0.1051$ is achieved in $17$ iterations, and a 10-fold greater accuracy of $0.0106$ is achieved in $25$ iterations and approximately in the same time. If the method worked without adaptation and strictly according to optimal lower bounds for the specified class of problems, then the increase in costs could be approximately $100$ times. Thus, due to the adaptability of the proposed method, we observe a convergence rate close to $O\left(\ln\left(\frac{1}{\varepsilon}\right)\right)$.

Now for a given set of $N$ points $\{A_{k}=(a_{1k},a_{2k},\ldots,a_{nk});\,k=\overline{1,N}\}$ in $n$-dimensional Euclidean space $\mathbb{R}^{n}$ we need to find such a point $x=(x_{1},x_{2},\ldots,x_{n})$, that the objective function

[TABLE]

would take the minimal value on the set $Q$, which is defined by the previous constraints (28). The coordinates of the points $A_{k}$ for $k=\overline{1,N}$, are chosen as the rows of the matrix $A\in\mathbb{R}^{N\times n}$. The entries of the matrix $A$ are random integers in the closed interval $[-10;10]$, which are drawn from the discrete uniform distribution.

The results of the work of Algorithm 1, for objective function (29) and for some different values of $n,m$ and $N$, are presented in Table 2 below. These results demonstrate the number of iterations produced by Algorithm 1 to reach the solution of the problem, the quality of the solution "General estimate", which is the right side of inequality (17), and the running time of the algorithm in seconds.

As we see from Table 2 we also observe a convergence rate close to $O\left(\ln\left(\frac{1}{\varepsilon}\right)\right)$, due to the adaptability of the proposed method.

Remark 6

Now we take all previous parameters but with points $A_{k}(k=\overline{1,N})$ in the unit ball. The results of Algorithm 1, for objective function (29) and for some different values of $n$, $m$ and $N$, with constraints (28), are presented in Table 3 below.

In this case, since the points $A_{k}$ are chosen in the unit ball, the operator of the variational inequality is bounded. The results of experiments in Table 3 show that the rate of convergence of the proposed method is close to $O\left(\frac{1}{\varepsilon}\right)$, which is significantly better than the optimal one, which is $O\left(\frac{1}{\varepsilon^{2}}\right)$, for non-smooth convex optimization problems and bounded variational inequalities [14].

5 Numerical Experiments for Matrix Games with Inexactness

We continue our experiments with computing a Nash equilibrium of a matrix game. For that purpose one should solve the saddle point problem

[TABLE]

where $x=(x_{1},\,x_{2},\,\ldots,\,x_{n})\,\in\,\mathbb{R}^{n}$, $y=(y_{1},\,y_{2},\,\ldots,\,y_{n})\,\in\,\mathbb{R}^{m}$, $\Delta_{n}$ is a standard simplex in $\mathbb{R}^{n}$, i.e. $\Delta_{n}=\{x\,\in\,\mathbb{R}^{n}\,|\,x\,\geq\,0,\,\,\sum\limits_{i=1}^{n}x_{i}=1\}$, $\Delta_{m}$ is a standard simplex in $\mathbb{R}^{m}$, $A$ is the payoff matrix for the $y$ player. In all experiments we use payoff distributions centered at zero. Consider the following operator

[TABLE]

The operator $g(u)$ from (31) is monotone on $Q$, and with this operator the VI (5) has the same solution as the saddle point problem (30). So, Mirror Prox methods could be used for solving it.

In all experiments with matrix games we use the entropy prox-function $d(x)=\sum\limits_{i=1}^{n}x_{i}\ln x_{i}$ in Bregman divergence (3). Entropy prox for matrix games on simplex is the best option (see [15]).

First of all, we calculate experimental numbers of iterations for adaptive proximal method for VI [7] and compare these numbers with theoretical estimation $N=\frac{C}{\varepsilon}$, for some $C>0$ [4]. For that kind of experiments we run simulations on two classes of random matrix games: $10\times 10$ and $100\times 100$ normally distributed payoff matrices. For the first setting we create 50 games at random and calculate average experimental number of iterations over all games. Fig. 1 shows the results for different $\varepsilon$. The experimental results are better than theoretical estimation in all cases.

In the second part of experiments with matrix games we compare the proposed Algorithm MPAI with adaptive Mirror Prox method from [7]. In this part of experiments we modified problem (30) by adding inexactness (a bounded by $\delta$) random noise) to the operator $g$ of VI (5). Fig. 2-3 show the results of these experiments.

For comparison we calculate the specific values that determine the degree of influence of the inexactness on the final estimate of the decisions’ accuracy. So, not the whole error estimations are compared, but only improved by our approach part of the error estimations.

For Algorithm 1 this specific value is equal to (see (19)):

[TABLE]

and for adaptive Mirror Prox [7] method we can estimate the analogous value in the following way:

[TABLE]

According to the scheme of the proofs of Theorems 3.1 and 3.2, we can obtain the analogous value for non-adaptive Mirror Prox with constant step $1/L$:

[TABLE]

It is worth mentioning that estimate (32) should be less than (33) because of adaptive reduction $\delta^{k+1}<\delta$.

We may see that the accumulated error for proposed MPAI method is smaller than the error for adaptive Mirror Prox method from  [7].

6 Conclusion

The paper proposes an analogue of the Mirror Prox method for variational inequalities with adaptive tuning not only for the constant $L_{\nu}$, but also for the magnitude of the operator’s error. Moreover, such an error can set the degree of discontinuity of the operator. It is proved that the proposed method converges with the optimal estimate of complexity $O\left(\frac{1}{\varepsilon}\right)$ for the Lipschitz-continuous operator and the magnitude of the error is limited. It is important that the result applies to a certain class of variational inequalities with bounded (generally, discontinuous) operators ($\nu=0$). The paper also presents the results of experiments that demonstrate the ability to work with the estimate of the complexity close to $O\left(\ln\left(\frac{1}{\varepsilon}\right)\right)$ even for a problem with bounded operator ($\nu=0$) and the experimental comparison between adaptive Mirror Prox and the proposed algorithm for a special application (matrix games with inexactness). We also note that all result of this paper is applicable for VI with relatively smooth operators (with the exception of Remark 5). More precisely, the prox-function and Bregman divergence in (2) may not be strongly convex (for convex optimization problems this situation was studied e.g. in [16, 17]).