Generalized Mirror Prox Algorithm for Monotone Variational Inequalities: Universality and Inexact Oracle

Abstract

We introduce an inexact oracle model for variational inequalities with monotone operators, propose a numerical method that solves such variational inequalities, and analyze its convergence rate. As a particular case, we consider variational inequalities with Hölder-continuous operator and show that our algorithm is universal. This means that, without knowing the Hölder exponent and Hölder constant, the algorithm has the least possible in the worst-case sense complexity for this class of variational inequalities. We also consider the case of variational inequalities with a strongly monotone operator and generalize the algorithm for variational inequalities with inexact oracle and our universal method for this class of problems. Finally, we show how our method can be applied to convex–concave saddle point problems with Hölder-continuous partial subgradients.

Appendices

Appendix A

Proof of Lemma 3.1

Proof

Let us fix some $\nu \in [0,1]$. Then, for any $x \in [0,1 ]$, we have that $x^{2\nu }\le 1$. On the other hand, for any $x \ge 1$, we have that $\ x^{2\nu }\le x^{2}$. Thus, for any $x \ge 0$, we have that $x^{2\nu }\le x^2+ 1$. Hence, for any $\alpha , \beta \ge 0$,

$$\begin{aligned} \alpha ^{\nu }\beta \le \frac{\alpha ^{2\nu }}{2}+\frac{\beta ^{2}}{2}\le \frac{\alpha ^{2}}{2}+\frac{\beta ^{2}}{2}+\frac{1}{2}. \end{aligned}$$

Substituting $\alpha =\frac{ba^{\frac{1}{1+\nu }}}{\delta ^{\frac{1}{1+\nu }}}$ and $\beta =\frac{ca^{\frac{1}{1+\nu }}}{\delta ^{\frac{1}{1+\nu }}}$, we obtain

$$\begin{aligned} \frac{b^{\nu }a^{\frac{\nu }{1+\nu }}}{\delta ^{\frac{\nu }{1+\nu }}}\frac{ca^{\frac{1}{1+\nu }}}{\delta ^{\frac{1}{1+\nu }}}\le \frac{b^{2}a^{\frac{2}{1+\nu }}}{2\delta ^{\frac{2}{1+\nu }}}+\frac{c^{2}a^{\frac{2}{1+\nu }}}{2\delta ^{\frac{2}{1+\nu }}}+\frac{1}{2} \end{aligned}$$

and

$$\begin{aligned} ab^{\nu }c \le \left( \frac{1}{\delta }\right) ^{\frac{1-\nu }{1+\nu }} \frac{a^{\frac{2}{1+\nu }}}{2} \left( b^2+c^2\right) + \frac{\delta }{2}. \square \end{aligned}$$

Appendix B

To show the practical performance of the proposed Algorithm 1, we performed a series of numerical experiments for the Lagrange saddle point problem induced by the Fermat-Torricelli-Steiner problem.

All experiments were made using Python 3.4, on a computer with Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz, 1992 Mhz, 4 Core(s), 8 Logical Processor(s), and 8 GB RAM.

We consider an example of a variational inequality with a non-smooth, i.e., with $\nu = 0$, and non-strongly monotone operator. For this VI, the proposed universal method, due to its adaptivity to the smoothness level of the problem, works in practice with iteration complexity much smaller than the one predicted by the theory. This example is inspired by the well-known Fermat–Torricelli–Steiner problem, to which we add some non-smooth functional constraints. This problem can be solved by a switching subgradient scheme [8, 55] with the complexity $O(1/\varepsilon ^2)$. But, as we will see, our method allows us to obtain much faster convergence in practice than the one given by this bound.

More precisely, for a given set of N points $A_k \in \mathbb {R}^n, k=1,\ldots , N$, consider the optimization problem

$$\begin{aligned} \min _{x \in Q} \left\{ f(x):= \sum \limits _{k=1}^N \Vert x - A_k\Vert _2 \left| \; \varphi _p(x):= \sum _{i=1}^n \alpha _{pi}|x_i| - 1 \le 0 \right. , \; p=1, \ldots , m \right\} , \end{aligned}$$

where Q is a convex and compact set, $\alpha _{pi}$ are drawn from the standard normal distribution and then truncated to be positive. The corresponding Lagrange saddle point problem is defined as

$$\begin{aligned} \min _{x \in Q} \max _{\lambda = (\lambda _1,\lambda _2,\ldots ,\lambda _m)^T \in \mathbb {R}^m_+} L(x,\lambda ):=f(x)+\sum \limits _{p=1}^m\lambda _p\varphi _p(x). \end{aligned}$$

As it was described in Sect. 6, this problem is equivalent to the variational inequality with the following monotone non-smooth operator:

$$\begin{aligned} G(x,\lambda )= \begin{pmatrix} \nabla f(x)+\sum \limits _{p=1}^m\lambda _p\nabla \varphi _p(x), \\ (-\varphi _1(x),-\varphi _2(x),\ldots ,-\varphi _m(x))^T \end{pmatrix}. \end{aligned}$$

For simplicity, we assume that there exists (potentially very large) bound on the norm of the optimal Lagrange multiplier $\lambda ^*$, which allows us to compactify the feasible set for the pair $(x,\lambda )$ to be a Euclidean ball of some radius. We also believe that the approach of [14, 44] for dealing with unbounded feasible sets can be extended to our setting and we leave this for future work.

Fig. 1

Results of Algorithm 1 for Fermat–Torricelli–Steiner problem with different values of m and n

We run Algorithm 1 with different values of n, m, and N using the standard Euclidean prox-setup and the starting point $(x^0, \lambda ^0) = \frac{1}{\sqrt{m+n}} {\textbf {1}} \in \mathbb {R}^{n+m}$, where ${\textbf {1}}$ is the vector of all ones. The points $A_k$, $k=1,\ldots , N$ are drawn randomly from the standard normal distribution. For each value of the parameters, the random data was drawn 10 times and the results were averaged. The results of the work of Algorithm 1 are represented in Fig. 1. For different values of the accuracy $\varepsilon \in \{ 1/2^i, i=1,2,3,4,5,6\}$, we report the number of iterations and the running time in seconds required by Algorithm 1 to reach an $\varepsilon$-solution of the considered problem.

As it is known [46], for a VI with a non-smooth operator, the theoretical iteration complexity estimate $O\left( \frac{1}{\varepsilon ^2}\right)$ is optimal. However, experimentally we see from the slope of the lines in Fig. 1 that, due to the adaptivity, the proposed Algorithm 1 has iteration complexity $O\left( \frac{1}{\root 4 \of {\varepsilon }}\right)$.

Fig. 2

Results of Algorithm 2 and Modified Projection Method with $n = 10^7$

Appendix C

In this appendix, in order to demonstrate the performance of the Generalized Mirror Prox with restarts (Algorithm 2), we consider the variational inequality with Lipschitz-continuous and strongly monotone operator (see Example 5.2 in [40]):

$$\begin{aligned} g: Q \subset \mathbb {R}^n \rightarrow \mathbb {R}^n, \quad g(x) = x. \end{aligned}$$

(35)

We compare the work of the proposed Algorithm 2 with Modified Projection Method, which was proposed in [40]. We run Algorithm 2 with different values of the accuracy $\varepsilon \in \{10^{-i}, i = 3,4,\ldots , 10 \}$ and for the dimension $n = 10^7$. We take $Q = \{x \in \mathbb {R}^n, \Vert x\Vert _2 \le 2\}$. The results of the comparison are presented in Fig. 2, which shows the norm $\Vert x_{\text {out}} - x_*\Vert _2$ as a function of iteration number, where $x_{\text {out}}$ is the output of each algorithm, and $x_*$ is the solution of the problem (1) with the operator g given in (35). Note that $x_* = {\textbf {0}} \in \mathbb {R}^n$. In the conducted experiments, we first run Algorithm 2 and calculate $\Vert x_{\text {out}} - x_*\Vert _2$ for the different previously mentioned values of $\varepsilon$ and the corresponding number of iterations required by the algorithm. For the calculated numbers of iterations of Algorithm 2, we rum Modified Projection Method and calculate the corresponding values $\Vert x_{\text {out}} - x_*\Vert _2$. From Fig. 2, we can see the higher efficiency of the proposed Algorithm 2 and the significant difference between the performance of the compared algorithms.