Golden Ratio Algorithms for Variational Inequalities

Introduction

We are interested in the variational inequality (VI) problem:

where $V$ is a finite dimensional vector space and we assume that

$g\colon V\to(-\infty,+\infty]$ is a proper convex lower semicontinuous (lsc) function;

$F\colon\operatorname{dom}g\to V$ is monotone: $\langle F(u)-F(v),u-v\rangle\geq 0$ $\forall u,v\in\operatorname{dom}g$.

The function $g$ can be nonsmooth, and it is very common to consider VI with $g=\delta_{C}$, the indicator function of a closed convex set $C$. In this case (1) reduces to

which is a more widely-studied problem. It is clear that one can rewrite (1) as a monotone inclusion: $0\in(F+\partial g)(x^{*})$. Henceforth, we implicitly assume that we can (relatively simply) compute the resolvent of $\partial g$ (the proximal operator of $g$), that is $(\operatorname{Id}+\partial g)^{-1}$, but cannot do this for $F$, in other words computing the resolvent $(\operatorname{Id}+F)^{-1}$ is prohibitively expensive.

VI is a useful way to reduce many different problems that arise in optimization, PDE, control theory, games theory to a common problem (1). We recommend as excellent references for a broader familiarity with the subject.

As a motivation from the optimization point of view, we present two sources where VI naturally arise. The first example is a convex-concave saddle point problem:

Saddle point problems are ubiquitous in optimization as this is a very convenient way to represent many nonsmooth problems, and this in turn often allows to improve the complexity rates from $O(1/\sqrt{k})$ to $O(1/k)$. Even in the simplest case when $K$ is bilinear form, the saddle point problem is a typical example where the two simplest iterative methods, the forward-backward method and the Arrow-Hurwicz method (see ), will not work. Korpelevich in and Popov in resolved this issue by presenting two-step methods that converge for a general monotone $F$. In turn, these two papers gave birth to various improvements and extensions, see .

Another important source of VI is a simpler problem of composite minimization

For a general VI, even when $F$ is Lipschitz continuous but nonlinear, computing its Lipschitz constant is not an easy task. Moreover, the curvature of $F$ can be quite different, so the stepsizes governed by the global Lipschitz constant will be very conservative. Thus, most practical methods for VI are required to use linesearch — an auxiliary iterative procedure which runs in each iteration of the algorithm until some criterion is satisfied, and it seems this is the only option for the case when $F$ is not Lipschitz. To this end, most known methods for VI with a fixed stepsize have their analogues with linesearch. This is still an active area of research rich in diverse ideas, see . The linesearch can be quite costly in general as it requires computing additional values of $F$ or $\operatorname{prox}_{g}$, or even both in every linesearch iteration. Moreover, the complexity estimates become not so informative, as they only say how many outer iterations one needs to reach the desired accuracy in which the number of linesearch iterations is of course not included.

Contributions. In this paper, our aim is to propose an adaptive algorithm for solving problem (1) with $F$ locally Lipschitz continuous. By adaptive we mean that the method does not require a linesearch to be run, and its stepsizes are computed using current information about the iterates. These stepsizes approximate an inverse local Lipschitz constant of $F$, thus they are separated from zero. Each iteration of the method needs only one evaluation of the proximal operator and one value of $F$. Moreover, the stepsizes are allowed to increase from iteration to iteration. To our knowledge, it is the first adaptive method with such properties. The method is easy to implement and it satisfies all standard rates for monotone VI: ergodic $O(1/k)$ and $R$-linear if the error bound condition holds. In particular, one of possible instances of the proposed algorithm can be given just in a few lines

where $z^{0},z^{1}\in V$, $\bar{z}^{0}=z^{1}$, $\lambda_{0}=\lambda_{-1}>0$, $\bar{\lambda}>0$.

Our approach is to start from the simple case when $F$ is $L$-Lipschitz continuous. For this case, we present the Golden Ratio Algorithm with a fixed stepsize, which is interesting on his own right and gives us an intuition for the more difficult case with dynamic steps. Sect. 2 collects these results. In Sect. 3 we show how one can derive new algorithms for fixed point problems based on the proposed framework. In particular, instead of working with the standard class of nonexpansive operators, we consider a more general class of demi-contractive operators. Sect. 4 collects two extensions of the adaptive Golden Ratio Algorithm. The first proposes an extension of the adaptive algorithm enhanced by two auxiliary metrics. Although it is simple theoretically, it is nevertheless still very important in applications, where it is preferable to use different weights for different coordinates. For the second extension we do not need monotonicity of $F$, but instead we require that the Minty variational inequality associated with (1) has a solution. In Sect. 5 we illustrate the performance of the method for several problems including the aforementioned nonmonotone case. Finally, Sect. 6 concludes the paper by presenting several directions for further research.

Preliminaries. Let $V$ be a finite-dimensional vector space equipped with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|=\sqrt{\langle\cdot,\cdot\rangle}$. For a lsc function $g\colon V\to(-\infty,+\infty]$ by $\operatorname{dom}g$ we denote the domain of $g$, i.e., the set $\{x\colon g(x)<+\infty\}$. Given a closed convex set $C$, $P_{C}$ stands for the metric projection onto $C$, $\delta_{C}$ denotes the indicator function of $C$ and $\operatorname{dist}(x,C)$ the distance from $x$ to $C$, that is $\operatorname{dist}(x,C)=\|P_{C}x-x\|$. The proximal operator $\operatorname{prox}_{g}$ for a proper lsc convex function $g\colon V\to(-\infty,+\infty]$ is defined as $\operatorname{prox}_{g}(z)=\operatorname{argmin}_{x}\{g(x)+\frac{1}{2}\|x-z\|^{2}\}$. The following characteristic property (prox-inequality) will be frequently used:

A simple identity important in our analysis is

The following important lemma will simplify the proofs of the main theorems.

Let $(z^{k})\subset V$ be a bounded sequence and suppose $\lim_{k\to\infty}\|z^{k}-z\|$ exists whenever $z$ is a cluster point of $(z^{k})$. Then $(z^{k})$ is convergent.

Assume that $u,v$ are two arbitrary cluster points of $(z^{k})$. From

we see that there exists $\lim_{k\to\infty}\langle z^{k}-u,u-v\rangle$. Assume now that $z^{k_{i}}\to u$ and $z^{k_{j}}\to v$. Then one can observe that

and hence, $u=v$. Since $u,v$ were arbitrary, we can conclude that $(z^{k})$ converges to some element in $V$. ∎

Golden Ratio Algorithms

Let $\varphi=\frac{\sqrt{5}+1}{2}$ be the golden ratio, that is $\varphi^{2}=1+\varphi$. The proposed Golden RAtio ALgorithm (GRAAL for short) reads as a simple recursion:

Suppose that $F\colon\operatorname{dom}g\to V$ is $L$–Lipschitz and conditions (C1)–(C3) are satisfied. Let $z^{1},\bar{z}^{0}\in V$ be arbitrary and $\lambda\in(0,\frac{\varphi}{2L}]$. Then $(z^{k})$, $(\bar{z}^{k})$, generated by (10), converge to a solution of (1).

Note that $z^{k}-\bar{z}^{k-1}=\frac{1+\varphi}{\varphi}(z^{k}-\bar{z}^{k})=\varphi(z^{k}-\bar{z}^{k})$. Hence, we can rewrite (12) as

Expressing the first two terms in (14) through norms, we arrive at

Choose $z=z^{*}\in S$. By (C3), the rightmost term in (2) is nonnegative:

By (7) and $z^{k+1}=(1+\varphi)\bar{z}^{k+1}-\varphi\bar{z}^{k}$ we have

Combining (2) and (16) with (2), we deduce

From $\lambda\leq\frac{\varphi}{2L}$ it follows that

From (20) we have that \bigl{(}(1+\varphi)\|\bar{z}^{k+1}-z^{*}\|^{2}+\frac{\varphi}{2}\|z^{k+1}-z^{k}\|^{2}\bigr{)} is bounded and so is $(\bar{z}^{k})$, and $\lim_{k\to\infty}\|z^{k}-\bar{z}^{k}\|=0$. Hence, $(z^{k})$ has at least one cluster point. By definition of $\bar{z}^{k}$, $z^{k+1}-\bar{z}^{k}=\varphi(z^{k+1}-\bar{z}^{k+1})\to 0$ and, hence, $z^{k+1}-z^{k}\to 0$ as well. Taking the limit in (11) (going to the subsequences if needed) and using (C2), we prove that all cluster points of $(z^{k})$ (and thus, of $(\bar{z}^{k})$) belong to $S$. From (20) one can see that the sequence \bigl{(}(1+\varphi)\|\bar{z}^{k}-z^{*}\|^{2}+\frac{\varphi}{2}\|z^{k}-z^{k-1}\|\bigr{)} is non-increasing and, hence, it is convergent. As $\lim_{k\to\infty}\|z^{k}-z^{k-1}\|=0$, there must exist $\lim_{k\to\infty}\|\bar{z}^{k}-z^{*}\|$. Since $z^{*}$ is an arbitrary point in $S$, from Lemma 1 it follows that $(\bar{z}^{k})$ and $(z^{k})$ converge to a point in $S$. ∎

Notice that the constant $\varphi$ is chosen not arbitrary, but as the largest constant $c$ that satisfies $\frac{1}{c}\geq c-1$ in order to get rid of the term $\|z^{k+1}-\bar{z}^{k}\|^{2}$ in (2). It is interesting to compare the proposed GRAAL with the reflected projected (proximal) gradient method . At a first glance, they are quite similar: both need one $F$ and one $\operatorname{prox}_{g}$ per iteration. The advantage of the former, however, is that $F$ is computed at $z^{k}$, which is always feasible $z^{k}\in\operatorname{dom}g$, due to the properties of the proximal operator. In the reflected projected (proximal) gradient method $F$ is computed at $2z^{k}-z^{k-1}$ which might be infeasible. Sometimes, as it will be illustrated in Sect. 5, this can be important.

In this section we introduce our fully adaptive algorithm. The algorithm preserves the same computational cost per iteration as (10) (i.e., no linesearch) and the stepsizes approximate an inverse local Lipschitz constant of $F$. Furthermore, for our purposes, the locally Lipschitz continuity of $F$ is sufficient. The algorithm, which we call the adaptive Golden Ratio Algorithm (aGRAAL), is presented as Alg. 1. For simplicity, we adopt the convention $\frac{0}{0}=+\infty$.

Notice that in the $k$-th iteration we need only to compute $F(z^{k})$ once and reuse the already computed value $F(z^{k-1})$. The constant $\bar{\lambda}$ in (21) is given only to ensure that $(\lambda_{k})$ is bounded. Hence, it makes sense to choose $\bar{\lambda}$ quite large. Although $\lambda_{0}$ can be arbitrary from the theoretical point of view, from (21) it is clear that $\lambda_{0}$ will influence further steps, so in practice we do not want to take it too small or too large. The simplest way is to choose $z^{0}$ as a small perturbation of $z^{1}$ and take $\lambda_{0}=\frac{\|z^{1}-z^{0}\|}{\|F(z^{1})-F(z^{0})\|}$. This gives us an approximationWe assume that $F(z^{1})\neq F(z^{0})$, otherwise choose another $z^{0}$. of the local inverse Lipschitz constant of $F$ at $z^{1}$.

As one can see from Alg. 1, for $\phi<\varphi$ one has $\rho>1$ and hence, $\lambda_{k}$ can be larger than $\lambda_{k-1}$. This is probably the most important feature of the proposed algorithm. When $F$ is Lipschitz, there are some other methods without linesearch that do not require to know the Lipschitz constant, e.g.,. However, such methods use nonincreasing stepsizes, which is rather restrictive.

Condition (21) leads to two key estimations. First, one has $\lambda_{k}\leq\lambda_{k-1}(\frac{1}{\phi}+\frac{1}{\phi^{2}})$, which in turn implies $\theta_{k}-1-\frac{1}{\phi}\leq 0$. Second, from $\frac{\phi\theta_{k-1}}{\lambda_{k-1}}=\frac{\theta_{k}\theta_{k-1}}{\lambda_{k}}$ one can derive

Alg. 1 is quite generic, and in order to simplify it we can predefine some of the constants. Let $\phi=\frac{3}{2}$. Then $\rho=\frac{10}{9}$ and as $\theta_{k-1}=\frac{\lambda_{k-1}}{\lambda_{k-2}}\phi$, we have $\frac{\phi\theta_{k-1}}{4\lambda_{k-1}}=\frac{9}{16\lambda_{k-2}}$. It is easy to see that with such choice of constants, Alg. 1 reduces to the one we presented in the introduction.

Suppose that $F\colon\operatorname{dom}g\to V$ is locally Lipschitz continuous. If the sequence $(z^{k})$ generated by Alg. 1 is bounded, then both $(\lambda_{k})$ and $(\theta_{k})$ are bounded and separated from .

It is obvious that $(\lambda_{k})$ is bounded. Let us prove by induction that it is separated from . As $(z^{k})$ is bounded there exists some $L>0$ such that $\|F(z^{k})-F(z^{k-1})\|\leq L\|z^{k}-z^{k-1}\|$. Moreover, we can take $L$ large enough to ensure that $\lambda_{i}\geq\frac{\phi^{2}}{4L^{2}\bar{\lambda}}$ for $i=0,1$. Now suppose that for all $i=0,\dots k-1$, $\lambda_{i}\geq\frac{\phi^{2}}{4L^{2}\bar{\lambda}}$. Then we have either $\lambda_{k}=\rho\lambda_{k-1}\geq\lambda_{k-1}\geq\frac{\phi^{2}}{4L^{2}\bar{\lambda}}$ or

Hence, in both cases $\lambda_{k}\geq\frac{\phi^{2}}{4L^{2}\bar{\lambda}}$. The claim that $(\theta_{k})$ is bounded and separated from now follows immediately. ∎

Define the bifunction $\Psi(u,v):=\langle F(u),v-u\rangle+g(v)-g(u)$. It is clear that (1) is equivalent to the following equilibrium problem: find $z^{*}\in V$ such that $\Psi(z^{*},z)\geq 0$ $\forall z\in V$. Notice that for any fixed $z$, the function $\Psi(z,\cdot)$ is convex.

Suppose that $F\colon\operatorname{dom}g\to V$ is locally Lipschitz and conditions (C1)–(C3) are satisfied. Then $(z^{k})$ and $(\bar{z}^{k})$, generated by Alg. 1, converge to a solution of (1).

Let $z\in V$ be arbitrary. By the prox-inequality (6) we have

Multiplying (27) by $\frac{\lambda_{k}}{\lambda_{k-1}}\geq 0$ and using that $\frac{\lambda_{k}}{\lambda_{k-1}}(z^{k}-\bar{z}^{k-1})=\theta_{k}(z^{k}-\bar{z}^{k})$, we obtain

Expressing the first two terms in (2.1) through norms, we derive

Recall that $\theta_{k}\leq 1+\frac{1}{\phi}$. Using (24), the rightmost term in (32) can be estimated as

Applying the obtained estimation to (32), we deduce

Iterating the above inequality, we derive

For a variational inequality in the form (2), condition (C3) can be relaxed to the following

$\langle F(z),z-\bar{z}\rangle\geq 0\quad\forall z\in C\quad\forall\bar{z}\in S$.

This condition, used for example in , is weaker than the standard monotonicity assumption (C3) or pseudomonotonicity assumption:

It is straightforward to see that Theorems 1, 2 hold under (C4). In fact, in the proof of the theorems, we choose $z=z^{*}\in S$ only to ensure that $\Psi(z^{*},z^{k})\geq 0$. However, it is sufficient to show that the first left-hand side in (2.1) is nonnegative. But for $z=z^{*}\in S$ this is true, since by (C4), $\langle F(z^{k}),z^{k}-z^{*}\rangle\geq 0$.

2 Ergodic convergence

It is known that many algorithms for monotone VI (or for more specific convex-concave saddle point problems) exhibit an $O(1/k)$ rate of convergence, see, for instance, , where such an ergodic rate was established. Moreover, Nemirovski has shown in that this rate is optimal. In this section we prove the same result for our algorithm. When the set $\operatorname{dom}g$ is bounded establishing such a rate is a simple task for most methods, including aGRAAL, however the case when $\operatorname{dom}g$ is unbounded has to be examined more carefully. To deal with it, we use the notion of the restricted merit function, first proposed in .

The proof is almost identical to Lemma 1 in . The only difference is that we have to consider VI (1) with a general $g$ instead of $\delta_{C}$. ∎

Now we can obtain something meaningful. Since $F$ is continuous and $g$ is lsc, there exist some constant $M>0$ that majorizes the right-hand side of (35) for all $z\in U$. From this follows that $\sum_{i=1}^{k}\lambda_{i}\Psi(z,z^{i})\leq M$ for all $z\in U$ (we ignore the constant $2$ before the sum). Let $Z^{k}$ be the ergodic sequence: $Z^{k}=\frac{\sum_{i=1}^{k}\lambda_{i}z_{i}}{\sum_{i=1}^{k}\lambda_{i}}$. Then using convexity of $\Psi(z,\cdot)$, we obtain

Taking into account that $(\lambda_{k})$ is separated from zero, we obtain the $O(1/k)$ convergence rate for the ergodic sequence $(Z^{k})$.

For the case of the composite minimization problem (5), instead of using the merit function it is simpler to use the energy residual: $J(z^{k})-J(z^{*})$. For this we need to use in (2.1) that

In this way, we may proceed analogously to obtain

3 Linear convergence

For many VI methods it is possible to derive a linear convergence rate under some additional assumptions. The most general tool for that is the use of error bounds. For a survey of error bounds, we refer the reader to and for their applications to the VI algorithms to .

Let us briefly recall the terminology associated with linear convergence. Suppose that $(u^{k})\subset V$ is a sequence that converges to $u\in V$. We say that convergence is $Q$-linear, if there is $q\in(0,1)$ such that $\|u^{k+1}-u\|<q\|u^{k}-u\|$ for all $k$ large enough. We also say that convergence is $R$-linear, if for all $k$ large enough, $\|u^{k}-u\|\leq\gamma_{k}$ and $(\gamma_{k})$ converges $Q$-linearly to zero.

Let us fix some $\lambda>0$ and define the natural residual $r(z,\lambda):=z-\operatorname{prox}_{\lambda g}(z-\lambda F(z))$. Evidently, $z\in S$ if and only if $r(z,\lambda)=0$. We say that problem (1) satisfies an error bound condition if there exist positive constants $\mu$ and $\eta$ such that

The function $\lambda\mapsto\|r(z,\lambda)\|$ is nondecreasing and $\lambda\mapsto\frac{\|r(z,\lambda)\|}{\lambda}$is nonincreasing (see [21, Proposition 10.3.6]), and thus all natural residuals $r(\cdot,\lambda)$ are equivalent. Hence the choice of $\lambda$ in the above definition is not essential. No doubt, it is not an easy task to decide whether (42) holds for a particular problem. Several examples are known, see for instance and it is still an important and active area of research.

In the analysis below we are not interested in sharp constants, but rather in showing the linear convergence for $(z^{k})$. This will allow us to keep the presentation simpler. For the same reason we assume that $F$ is $L$-Lipschitz continuous.

Choose any $\lambda>0$ such that $\lambda_{k}\geq\lambda$ for all $k$. Without loss of generality, we assume that $\lambda$ is the same as in (42). As $\lambda\mapsto\|r(\cdot,\lambda)\|$ is nondecreasing and $\operatorname{prox}_{\lambda_{k}g}$ is nonexpansive, using the triangle inequality, we obtain

Let $\beta=\frac{\phi}{\phi-1}$. If $(\theta_{k})$ is separated from zero, then the above inequality ensures that for any $z^{k},\bar{z}^{k}$ and $\varepsilon\in(0,1)$ there exists $m\in(0,1)$ such that

The presence of so many constants in (45) will be clear later. In order to proceed, we have to modify Alg. 1. Now instead of (21), we choose the stepsize by

This modification basically means that we slightly bound the stepsize. However, this is not crucial for the steps, as we can choose $\delta$ arbitrary close to one. An argument completely analogous to that in the proof of Lemma 2 (up to the factor $\delta$) shows that both $(\lambda_{k})$ and $(\theta_{k})$ are bounded and separated from zero. This confirms correctness of our arguments about $(\lambda_{k})$ and $(\theta_{k})$ in (2.3) and (44). It should be also obvious that Alg. 1 with (46) instead of (21) has the same convergence properties, since (33) — the only place where this modification plays some role — will be still valid. For any $\delta\in(0,1)$ there exist $\varepsilon\in(0,1)$ and $m\in(0,1)$ such that $\delta=(1-\varepsilon)(1-m)$ and (45) is fulfilled for any $z^{k},\bar{z}^{k}$. Now using (46), one can derive a refined version of (33):

With this inequality, instead of (34), for $z=z^{*}\in S$ we have

where in the last inequality we have used (45). As $(\bar{z}^{k})$ converges to a solution, $r(\bar{z}^{k},\lambda)$ goes to , and hence $\|r(\bar{z}^{k},\lambda)\|\leq\eta$ for all $k\geq k_{0}$. Setting $z^{*}=P_{S}(\bar{z}^{k})$ in (2.3) and using (42) and that $\Psi(z^{*},z^{k})\geq 0$, we obtain

From this the $Q$-linear rate of convergence for the sequence \bigl{(}\operatorname{dist}(\bar{z}^{k},S)^{2}+\frac{\theta_{k-1}}{2}\|z^{k}-z^{k-1}\|^{2}\bigr{)} follows. Since $(\theta_{k})$ is separated from zero, we conclude that $\|z^{k}-z^{k-1}\|$ converges $R$-linearly and this immediately implies that the sequence $(z^{k})$ converges $R$-linearly. We summarize the obtained result in the following statement.

Suppose that conditions (C1)–(C3) are satisfied, $F\colon\operatorname{dom}g\to V$ is $L$–Lipschitz and the error bound (42) holds. Then $(z^{k})$, generated by Alg. 1 with (46) instead of (21), converges to a solution of (1) at least $R$-linearly.

Fixed point algorithms

Although in general it is a very standard way to formulate a VI (1) as a fixed point equation $x=\operatorname{prox}_{g}(x-F(x))$, sometimes other way around might also be beneficial. In this section we show how one can apply general algorithms for VI to find a fixed point of some operator $T\colon V\to V$. Clearly, any fixed point equation $x=Tx$ is equivalent to the equation $F(x)=0$ with $F=\operatorname{Id}-T$. The latter problem is of course a particular instance of (1) with $g\equiv 0$. Hence, we can work under the assumptions of Remark 1.

By $\operatorname{Fix}T$ we denote the fixed point set of the operator $T$. Although, with a slight abuse of notation, we will not use brackets for the argument of $T$ (this is common in the fixed point literature), but we continue doing that for the argument of $F$.

We are interested in the following classes of operators:

Therefore, (e) is the most general class of the aforementioned operators. Sometimes in the literature the authors consider operators that satisfy a more restrictive condition

and call them $\rho$–demi-contractive for $\rho\in[0,1)$ and hemi-contractive for $\rho=1$. We consider only the most general case with $\rho=1$, but still for simplicity will call them as demi-contractive. It is also tempting to call the class (e) as quasi-pseudocontractive by keeping the analogy between (b) and (c), but this tends to be a bit confusing due to “quasi-pseudo”. Notice that for the case $\rho<1$ in (51), one can consider a relaxed operator $S=\rho\operatorname{Id}+(1-\rho)T$. It is easy to see that $S$ belongs to the class (c) and $\operatorname{Fix}S=\operatorname{Fix}T$. However, with $\rho=1$, the situation becomes much more difficult.

When $T$ belongs to (a) or (b), the standard way to find a fixed point of $T$ is by applying the celebrated Krasnoselskii-Mann scheme (see e.g. ): $x^{k+1}=\alpha x^{k}+(1-\alpha)Tx^{k}$, where $\alpha\in(0,1)$. The same method can be applied when $T$ is in class (c), but to the averaged operator $S=\beta\operatorname{Id}+(1-\beta)T$, $\beta\in(0,1)$, instead of $T$. However, things become more difficult when we consider broader classes of operators. In particular, Ishikawa in proposed an iterative algorithm when $T$ is Lipschitz continuous and pseudo-contractive. However, its convergence requires a compactness assumption, which is already too restrictive. Moreover, the convergence is rather slow, since the scheme uses some auxiliary slowly vanishing sequence (as in the subgradient method), also the scheme uses two evaluations of $T$ per iteration. Later, this scheme was extended in to the case when $T$ is Lipschitz continuous and demi-contractive but with the same assumptions as above.

Obviously, one can rewrite condition (e) as

which means that the angle between vectors $Tx-x$ and $\bar{x}-x$ must always be nonobtuse.

We know that $T$ is pseudo-contractive if and only if $F$ is monotone [10, Example 20.8]. It is not more difficult to check that $T$ is demi-contractive if and only if $F$ satisfies (C4). In fact, in this case $S=\operatorname{Fix}T$, thus, (52) and (C4) are equivalent: $\langle F(x),x-\bar{x}\rangle=\langle x-Tx,x-\bar{x}\rangle\geq 0$. The latter observation allows one to obtain a simple way to find a fixed point of a demi-contractive operator $T$. In particular, in order to find a fixed point of $T$, one can apply the aGRAAL. Moreover, since in our case $g\equiv 0$ and $\operatorname{prox}_{g}=\operatorname{Id}$, (23) simplifies to

Let $T\colon V\to V$ be locally Lipschitz and demi-contractive operator. Define $F=\operatorname{Id}-T$. Then the sequence $(x^{k})$ defined by aGRAAL with (3) instead of (23) converges to a fixed point of $T$.

Since $F=\operatorname{Id}-T$ is obviously locally Lipschitz, the proof is an immediate application of Theorem 2. ∎

The obtained theorem is interesting not only for the very general class of demi-contractive operators, but for a more limited class of non-expansive operators. The scheme (3) requires roughly the same amount of computations as the Krasnoselskii-Mann algorithm, but the former method defines $\lambda_{k}$ from the local properties of $T$, and hence, depending on the problem, it can be much larger than $1$. Recently, there appeared some papers on speeding-up the Krasnoselskii-Mann (KM) scheme for nonexpansive operators . The first paper proposes a simple linesearch in order to accelerate KM. However, as it is common for all linesearch procedures, each inner iteration requires an evaluation of the operator $T$, which in the general case can eliminate all advantages of it. The second paper considers a more general framework, inspired by Newton methods. However, its convergence guarantees are more restrictive.

The demi-contractive property can be useful when we want to analyze the convergence of the iterative algorithm $x^{k+1}=Tx^{k}$ for some operator $T$. It might happen that we cannot prove that $(x^{k})$ is Fejér monotone w.r.t. $\operatorname{Fix}T$, but instead we can only show that $\|x^{k+1}-\bar{x}\|^{2}\leq\|x^{k}-\bar{x}\|^{2}+\|x^{k+1}-x^{k}\|^{2}$ for all $\bar{x}\in\operatorname{Fix}T$. This estimation guarantees that $T$ is demi-contractive and hence one can apply Theorem 4 to obtain a sequence that converges to $\operatorname{Fix}T$.

Finally, we can relax condition in (e) to the following

which in turn is equivalent to $\langle F(x),x-P_{S}x\rangle\geq 0$ for all $x\in V$. For instance, (54) might arise when we know that $F=\operatorname{Id}-T$ satisfies a global error bound condition

where $\varepsilon>0$ is some constant. One can easily show that (54) follows from (55) and (56), whenever $\varepsilon<\frac{1}{\mu^{2}}$. The proof of convergence for aGRAAL with this new condition will be almost the same: in the $k$-th iteration instead of using arbitrary $x^{*}\in S$, we have to set $x^{*}=P_{S}x^{k}$. Of course, the global error bound condition (55) is too restrictive and is difficult to verify [21, Chapter 6]. On the other hand, property (56) is very attractive: it is often the case that for some iterative scheme one can only show (56), which is not enough to proceed further with standard arguments like the Krasnoselskii-Mann theorem or Banach contraction principle. We believe it might be an interesting avenue for further research to consider the above settings with $\mu$ and $\varepsilon$ dependent on $x$ in order to eliminate the limitation of (55). Condition (56) also relates to the recent work where the generalization of nonexpansive operators for multivalued case was considered.

Generalizations

This section deals with two generalizations. First, we present aGRAAL in a general metric settings. Although this extension should be straightforward for most readers, we present it to make this work self-contained. Next, by revisiting the proof of convergence for aGRAAL we relax the monotonicity condition to a new one and discuss its consequences.

Throughout section 2 we have been working in standard Euclidean metric $\langle\cdot,\cdot\rangle$ and assumed that $F$ is monotone with respect to it. There are at least two possible generalizations of how one can incorporate some metric into Alg. 1. Firstly, the given operator $F$ may not be monotone in metric $\langle\cdot,\cdot\rangle$, but it is so in $\langle\cdot,\cdot\rangle_{P}$, induced by some symmetric positive definite operator $P$. For example, the generalized proximal method $z^{k+1}=(\operatorname{Id}+P^{-1}G)^{-1}z^{k}$ for some monotone operator $G$ gives us the operator $T=(\operatorname{Id}+P^{-1}G)^{-1}$ which is nonexpansive in metric $\langle\cdot,\cdot\rangle_{P}$, and hence $F=\operatorname{Id}-T$ is monotone in that metric but is not so in $\langle\cdot,\cdot\rangle$. This is an important example, as it incorporates many popular methods: ADMM, Douglas-Rachford, PDHG, etc. Secondly, it is often desirable to consider an auxiliary metric $\langle\cdot,\cdot\rangle_{M}$, induced by a symmetric positive definite operator $M$, that will enforce faster convergence. For instance, for saddle point problems, one may give different weights for primal and dual variables, in this case one will consider some diagonal scaling matrix $M$. The standard analysis of aGRAAL (and this is common for other known methods) does not take into account that $F$ is derived from a saddle point problem; it treats the operator $F$ as a black box, and just use one stepsize $\lambda$ for both primal and dual variables. For example, the primal-dual hybrid gradient algorithm , which is a very popular method for solving saddle point problems with a linear operator, uses different steps $\tau$ and $\sigma$ for primal and dual variables; and the performance of this method drastically depends on the choice of these constants. This should be kept in mind when one applies aGRAAL for such problems.

Another possibility is if we already work in metric $\langle\cdot,\cdot\rangle_{P}$, then a good choice of the matrix $M$ can eliminate some undesirable computations, like computing the proximal operator in metric $\langle\cdot,\cdot\rangle_{P}$. Of course, the idea to incorporate a specific metric to the VI algorithm for a faster convergence is not new and was considered, for example, in . Our goal is to show that the framework of GRAAL can easily adjust to these new settings.

Let $M,P\colon V\to V$ be symmetric positive definite operators. We consider two norms induces by $M$ and $P$ respectively

For a symmetric positive definite operator $W$, we define the generalized proximal operator as $\operatorname{prox}_{g}^{W}=(\operatorname{Id}+W^{-1}\partial g)^{-1}$. Since now we work with the Euclidean metric $\langle\cdot,\cdot\rangle_{P}$, we have to consider a more general form of VI:

the solution set $S$ of (58) is nonempty.

$g\colon V\to(-\infty,+\infty]$ is a convex lsc function;

$F\colon\operatorname{dom}g\to V$ is $P$–monotone: $\langle F(u)-F(v),u-v\rangle_{P}\geq 0$ $\forall u,v\in\operatorname{dom}g$.

The modification of Alg. 1 presented below assumes that the matrices $M,P$ are already given. In this note we do not discuss how to choose the matrix $M$, as it depends crucially on the problem instance.

Suppose that $F\colon\operatorname{dom}g\to V$ is locally Lipschitz continuous and conditions (D1)–(D3) are satisfied. Then $(z^{k})$ and $(\bar{z}^{k})$, generated by Alg. 2, converge to a solution of (58).

Fix any $z^{*}\in S$. By the prox-inequality (6) we have

Multiplying (63) by $\frac{\lambda_{k}}{\lambda_{k-1}}\geq 0$ and using that $\frac{\lambda_{k}}{\lambda_{k-1}}(z^{k}-\bar{z}^{k-1})=\theta_{k}(z^{k}-\bar{z}^{k})$, we obtain

Expressing the first two terms in (4.1) through norms, we derive

Notice that $\theta_{k}\leq 1+\frac{1}{\phi}$. Using (24), the rightmost term in (4.1) can be estimated as

Applying the obtained estimation to (4.1), we obtain

It is obvious to see that the statement of Lemma 2 is still valid for Alg. 2, hence one can finish the proof by the same arguments as in the end of Theorem 2. ∎

2 Beyond monotonicity

A more careful examination of the proof of Theorem 2 can help us to relax assumption (C3) (or (C4)) even more. In particular, one can impose

The above problem is known as Minty variational inequality (MVI) associated with VI (1). Let $S_{MVI},S_{VI}$ denote the solution sets of MVI (71) and VI (1) respectively. Essentially, condition (71) asks for $S_{MVI}\neq\varnothing$. It is a standard fact that when $F$ is monotone, both problems are equivalent, see [29, Lemma 1.5]. In general case ($F$ is continuous), one can only claim that $S_{MVI}\subset S_{VI}$.

Suppose that $F\colon\operatorname{dom}g\to V$ is locally Lipschitz continuous, $g\colon V\to(-\infty,+\infty]$ is convex lsc, and $S_{MVI}\neq\varnothing$. Then all cluster points of $(z^{k})$, generated by Alg. 1, are solutions of (1).

Fix any $\bar{z}\in S_{\text{MVI}}$. Then in the proof of Theorem 2 instead of taking arbitrary $z$, choose $z=\bar{z}$. Then instead of (2.1), we obtain

where the last inequality holds because of (71). Proceeding as in (2.1)–(35), we can deduce

From this we obtain that $(\bar{z}^{k})$ is bounded and $\theta_{k}\|z^{k}-\bar{z}^{k}\|\to 0$. Then in the same way as in Theorem 2, one can show that all cluster points of $(z^{k})$ are elements of $S_{VI}$. ∎

Of course, we obtain slightly weaker convergence guarantees than in Theorem 2, but on the other hand our assumptions are much more general.

Numerical experiments

This section collects several numerical experimentsAll codes can be found on https://gitlab.gwdg.de/malitskyi/graal.git. to confirm our findings. Computations were performed using Python 3.6 on a standard laptop running 64-bit Debian GNU/Linux. In all experiments we take $\phi=1.5$ for Alg. 1.

Here we study a Nash–Cournot oligopolistic equilibrium model. We give only a short description, for more details we refer to . There are $n$ firms, each of them supplies a homogeneous product in a non-cooperative fashion. Let $q_{i}\geq 0$ denote the $i$th firm’s supply at cost $f_{i}(q_{i})$ and $Q=\sum_{i=1}^{n}q_{i}$ be the total supply in the market. Let $p(Q)$ denote the inverse demand curve. A variational inequality that corresponds to the equilibrium is

where $F(q^{*})=(F_{1}(q^{*}),\dots,F_{n}(q^{*}))$ and

As a particular example, we assume that the inverse demand function $p$ and the cost function $f_{i}$ take the form:

with some constants that will be defined later. This is a classical example of the Nash-Cournot equilibrium first proposed in for $n=5$ players and later reformulated as monotone VI in . To make the problem even more challenging, we set $n=1000$ and generate our data randomly by two scenarios. Each entry of $\beta$, $c$ and $L$ are drawn independently from the uniform distributions with the following parameters:

$\gamma=1.1$, $\beta_{i}\sim\mathcal{U}(0.5,2)$, $c_{i}\sim\mathcal{U}(1,100)$, $L_{i}\sim\mathcal{U}(0.5,5)$;

$\gamma=1.5$, $\beta_{i}\sim\mathcal{U}(0.3,4)$ and $c_{i}$, $L_{i}$ as above.

2 Convex feasibility problem

Given a number of closed convex sets $C_{i}\subset V$, $i=1,\dots,m$, the convex feasibility problem (CFP) aims to find a point in their intersection: $x\in\cap_{i=1}^{m}C_{i}$. The problem is very general and allows one to represent many practical problems in this form. Projection methods are a standard tool to solve such problem (we refer to for a more in-depth overview). In this note we study a simultaneous projection method: $x^{k+1}=Tx^{k}$, where $T=\frac{1}{m}(P_{C_{1}}+\dots+P_{C_{m}})$. Its main advantage is that it can be easily implemented on a parallel computing architecture. Although it might be slower than the cyclic projection method $x^{k+1}=P_{C_{1}}\dots P_{C_{m}}x^{k}$ in terms of iterations, for large-scale problems it is often much faster in practice due to parallelization and more efficient ways of computing $Tx$.

One can look at the iteration $x^{k+1}=Tx^{k}$ as an application of the Krasnoselskii-Mann scheme for the firmly-nonexpansive operator $T$. By that $(x^{k})$ converges to a fixed point of $T$, which is either a solution of CFP (consistent case) or a solution of the problem $\min_{x}\sum_{i=1}^{m}\operatorname{dist}(x,C_{i})^{2}$ (inconsistent case when the intersection is empty).

To illustrate Remark 3, we show how in many cases aGRAAL with $F=\operatorname{Id}-T$ can accelerate convergence of the simultaneous projection algorithm. We believe, this is quite interesting, especially if one takes into account that our framework works as a black box: it does not require any tuning or a priori information about the initial problem.

Tomography reconstruction. The goal of the tomography reconstruction problem is to obtain a slice image of an object from a set of projections (sinogram). Mathematically speaking, this is an instance of a linear inverse problem

3 Sparse logistic regression

In this section we demonstrate that even for nice problems as convex composite minimization (5) with Lipschitz $\nabla f$, aGRAAL can be faster than the proximal gradient method or accelerated proximal gradient method. This is interesting, especially since the latter method has a better theoretical convergence rate.

Our problem of interest is a sparse logistic regression

We compare our method with the proximal gradient method (PGM) and FISTA with a fixed stepsize. We have not included into comparison the extensions of these methods with linesearch, as we are interested in methods that require one evaluation of $\nabla f$ and one of $\operatorname{prox}_{g}$ per iteration. For both methods we set the stepsize as $\lambda=\frac{1}{L_{\nabla f}}=\frac{4}{\|K^{T}K\|}$. We take two popular datasets from LIBSVM : real-sim with $m=72309$, $n=20958$ and kdda-2010 with $m=510302$, $n=2014669$. For both problems we set $\gamma=0.005\|A^{T}b\|_{\infty}$. We are aware of that neither PGM nor FISTA can be considered as the state of the art for (78), stochastic methods seem to be more competitive as the size of the problem is quite large. Overall our motivation is not to propose the best method for (78) but to demonstrate the performance of aGRAAL on some real-world problems.

We run all methods for sufficiently many iterations and compute the energy $J(x^{k})$ in each iteration. If after $k$ iterations the residual was small enough: $\|x^{k}-\operatorname{prox}_{g}(x^{k}-\nabla f(x^{k}))\|\leq 10^{-6}$, we choose the smallest energy value among all methods and set it to $J_{*}$. In Fig. 4 we show how the energy residual $J(x^{k})-J_{*}$ is changing w.r.t. the iterations. Since the dimensions in both problems are quite large, the CPU time for all methods is approximately the same.

We have presented the results of only two test problems merely for compactness, in fact similar results were observed for other datasets that we tested: rcv1, a9a, ijcnn1, covtype. An explanation for such a good performance of aGRAAL is of course that for this problem the global Lipschitz constant of $\nabla f$ is too conservative. Notice also that our algorithm did not take into account that in this case $F=\nabla f$ is a potential operator. It would be interesting to see how we can enhance aGRAAL with this information.

4 Beyond monotonicity

Finally, we illustrate numerically that aGRAAL can work even for nonmonotone problems.

For each $n\in\{100,500,1000,5000\}$ we generate $100$ random problems. We run aGRAAL for $10000$ iterations and stopped it whenever the accuracy $\|F(z^{k})\|\leq 10^{-6}$ reached. Table 1 shows the success rate of solving this problem, i.e., we counted only those problems where $\|z^{k}\|$ was large enough to make sure that this is not a trivial solution. We also report the average number of iterations (among all successful instances) the method needs to find a non-trivial solution. The entries of $A$, $B$ are drawn independently from the normal distribution $\mathcal{N}(0,1)$. The starting point is always $z^{1}=(1,\dots,1)$. It is clear that $F$ is not monotone; moreover, the transcendental functions $\sin$, $\exp$ make this problem highly nonlinear.

Conclusions and further directions

We conclude our work by presenting some possible directions for future research.

Fixed point iteration. It is interesting to represent scheme (10) as a fixed point iteration. To this end, let

Then it is not difficult to see that we can rewrite (10) as

If we now set $\mathbf{u}^{k}=(\bar{z}^{k-1},z^{k})$, the above equation simplifies to $\mathbf{u}^{k+1}=GR\mathbf{u}^{k}$. However, so far it not clear how to derive convergence of (10) from the fixed point perspective. Although, $G$ is a firmly-nonexpansive operator, $R$ is definitely not; thus, it is difficult to say something meaningful about $G\circ R$.

Inertial extensions. Starting from the paper , it is observed that using some inertia for the optimization algorithm often accelerates the latter. Later, papers extend this idea to a more general case with monotone operators. In our case it will be in particular interesting to do so, since our scheme

uses $\bar{z}^{k}$ as a convex combination of all previous iterates $z^{1}\dots,z^{k}$. This is completely opposite to the inertial methods, where one uses $z^{k}+\alpha(z^{k}-z^{k-1})$ for some $\alpha>0$.

Bregman distance. For many VI methods it is possible to derive their analogues for the Bregman distances, as this is done for the the extragradient method by extensions . It is possible to do so for GRAAL? This extension is not trivial since, for example, in (2) we have used the identity (7), where the linear structure was explicitly used.

Stochastic settings. For large-scale problems it is often the case that even computing $F(z^{k})$ becomes prohibitively expensive. For this reason, the stochastic VI methods that compute $F(z^{k})$ approximately can be advantageous over their deterministic counterparts, as it was shown in . It is interesting to derive similar extensions for GRAAL. The same concerns the coordinate extensions of GRAAL.

Continuous dynamic. Many discrete optimization/VI methods can be studied from the continuous perspective which may shed light on the discrete method. Originated from 60-s, this line of research was popularized by many authors, see and references therein. This research brings new and often deeper understanding of the respective iterative schemes, as, for instance, the case with the Nesterov acceleration in . For aGRAAL it is a challenging task to derive a continuous scheme, since the function $\lambda(t)$ cannot be defined in advance.

Y. Maltsky was supported by German Research Foundation grant SFB755-A4. The author would like to thank Anna–Lena Martins, Panayotis Mertikopoulos, Matthew Tam, associate editor and anonymous referees for their useful comments that have significantly improved the quality of the paper.

References