New First-Order Algorithms for Stochastic Variational Inequalities

Given a constraint set $\mathcal{Z}\subseteq\mathbb{R}^{n}$ and a mapping $F:\mathcal{Z}\rightarrow\mathbb{R}^{n}$, the classical variational inequality (VI) problem is to find $z^{*}\in\mathcal{Z}$ such that

$\displaystyle F(z^{*})^{\top}(z-z^{*})\geq 0,\quad\forall z\in\mathcal{Z}.$

(1)

For an introduction to VI and its applications, we refer the readers to Facchinei and Pang [4] and the references therein.

In this paper, we consider a stochastic version of problem (1), where the exact evaluation of the mapping $F(\cdot)$ is inaccessible. Instead, only a stochastic oracle is available. The stochasticity in question may stem from, e.g., the non-deterministic nature of mixed strategies of the players in a game-setting, or simply because of the difficulty in evaluating the mapping itself. The latter has become more pronounced in the literature, due to its recent-found application as a training/learning subproblem in machine learning and/or statistical learning. The so-called stochastic oracle is a noisy estimation of the mapping $F(\cdot)$, and an iterative scheme that incorporates such oracle is known as stochastic approximation (SA). As far as we know, the first proposal to use such approach for stochastic optimization can be traced back to the seminal work of Robbins and Monro [31]. In 2008, Jiang and Xu [10] followed the SA approach to solve VI models. Since then, efforts have been made to extend existing deterministic methods to the stochastic VI models; see e.g. [11, 41, 12, 15, 9, 7, 8].

Let us start our discussion by introducing the assumptions made throughout the paper. We consider VI model (1) where $\mathcal{Z}$ is a closed convex set. Moreover, the following two conditions are assumed:

$$\left(F(z)-F(z^{\prime})\right)^{\top}(z-z^{\prime})\geq\mu\|z-z^{\prime}\|^{2},\quad\forall z,z^{\prime}\in\mathcal{Z},$$

(2)

for some $\mu>0$, and

$$\|F(z)-F(z^{\prime})\|\leq L\|z-z^{\prime}\|,\quad\forall z,z^{\prime}\in\mathcal{Z},$$

(3)

for some $L\geq\mu>0$. Condition (2) is known as the strong monotonicity of $F$, while Condition (3) is known as the Lipschitz continuity of $F$. If Condition 2 is met with $\mu=0$ then $F$ is known as monotone. VI problems satisfying (2) with positive $\mu$ can be easily shown to have a unique solution $z^{*}$. Let us denote $\kappa:=\frac{L}{\mu}\geq 1$. Parameter $\kappa$ is usually known as the condition number of the VI model (1). We also assume

$\displaystyle\max\limits_{z,z^{\prime}\in\mathcal{Z}}\|z-z^{\prime}\|\leq D,$

(4)

namely the constraint set is bounded. Remark that this assumption can actually be removed without affecting the results, but then the analysis becomes lengthier and tedious without conceivable conceptual benefit, and so we shall not pursue that generality in this paper.

The stochastic oracle of the mapping, denote by $\hat{F}(z,\xi)$, takes a random sample $\xi\in\Xi$ from some sample space $\Xi$. The oracle is required to satisfy:

	$\displaystyle\mathbb{E}\left[\|\hat{F}(z,\xi)-F(z)\|\right]$	$\displaystyle\leq$	$\displaystyle\delta$		(5)
	$\displaystyle\mathbb{E}\left[\|\hat{F}(z,\xi)-F(z)\|^{2}\right]$	$\displaystyle\leq$	$\displaystyle\sigma^{2}$		(6)

for all $z\in\mathcal{Z}$, where $\delta,\sigma\geq 0$ are some constants. In other words, we assume both the bias and the deviation are uniformly upper-bounded.

In this paper, we propose two stochastic first-order schemes: the stochastic extra-point scheme and the stochastic extra-momentum scheme. The first scheme maintains two sequences of iterates featuring several well-known first-order search directions such as: the extra-gradient [13, 36], the heavy-ball [29], Nesterov’s extrapolation [23, 25], and the optimism direction [30, 20]. The second scheme, on the other hand, specifically combines the heavy-ball momentum and the optimism momentum in its updating formula, and maintains only one sequence throughout the iterations, therefore requiring only one projection per iteration. These two approaches require different types of analysis. Both schemes render a wider range of search directions than the existing first-order methods, and the parameters associated with each search direction could and should be tuned differently from problem-class to problem-class in order to secure good practical performances. The deterministic counterpart of these methods can be found in our previous work [6]. In the stochastic context, we show that as long as the variance can be reduced throughout the iterations, they yield the optimal iteration complexity $\mathcal{O}(\kappa\ln(1/\epsilon))$ (cf. [42]) to reach $\epsilon$-solution: $\|z^{k}-z^{*}\|^{2}\leq\epsilon$, with an additional biased term depending on $\delta$. In a later section, we demonstrate an application to the stochastic black-box minimax saddle-point problem where only noisy function values $f(x,y)$ are accessible. This application is particularly relevant, given its applications in machine learning, where the training data set may be very large and evaluating exact gradient/function value is usually impractical. Through a smoothing technique, we propose a stochastic zeroth-order gradient as our update directions in either the stochastic extra-point scheme or the stochastic extra-momentum scheme. We show that both approaches yield an iteration complexity of $\mathcal{O}(\kappa\ln(1/\epsilon))$ and a sample-complexity of $\mathcal{O}\left(\frac{\kappa^{2}}{\epsilon}\ln\left(\frac{1}{\epsilon}\right)\right)$.

The rest of the paper is organized as follows. In Section 2, we survey the relevant literature with a focus on stochastic VI. In Section 3, we present the main results in this paper, i.e. the convergence results of the two proposed methods, while the technical proofs are relegated to the appendices. In Section 4, we introduce a stochastic black-box saddle-point problem and present the sample complexity results of our methods. We present some promising preliminary numerical results in Section 5 and conclude the paper in Section 6.

The first-order algorithms for deterministic VI (1) serve as a basis for the developments of their stochastic counterparts. These algorithms include the projection method [4], the proximal method [18, 32, 36], the extra-gradient method [13, 36], the optimistic gradient descent ascent (OGDA) method [30, 20, 21], the mirror-prox method [22], the extrapolation method [26, 24], and the extra-point method [6].

In this section, we shall focus on the developments of algorithms for stochastic VI, starting with a paper of Jiang and Xu [10], where the authors propose a stochastic projection method for solving strongly monotone and Lipschitz continuous VI problems and present an almost-sure convergence result. Koshal et al. [14] propose iterative Tikhonov regularization method and iterative proximal point method and show almost-sure convergence with monotone and Lipschitz continuous VI problems. Both methods solve a strongly monotone VI subproblem at each iteration. Yousefian et al. [39] further introduce local smoothing technique to the above-mentioned regularized methods to account for non-Lipschitz mappings and show almost-sure convergence. A survey on these methods, as well as applications and the theory behind stochastic VI can be found in Shanbhag [35].

Juditsky et al. [11] are among the first to show an iteration complexity bound for stochastic VI algorithms. They extend the mirror-prox method [22] to stochastic settings and prove an optimal iteration complexity bound for monotone VI: $\mathcal{O}(\frac{1}{\epsilon^{2}})$ , or $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ when the variance can be controlled small enough. Yousefian et al. [41] further extend the stochastic mirror-prox method with a more general step size choice and show an $\mathcal{O}\left(\frac{1}{\epsilon^{2}}\right)$ iteration complexity, where they also show an $\mathcal{O}(\frac{1}{\epsilon})$ complexity for the stochastic extra-gradient method for solving strongly monotone VI problems. Yousefian et al. [40] use randomized smoothing technique for non-Lipschitz mapping and show an $\mathcal{O}\left(\frac{1}{\epsilon^{6}}\right)$ iteration complexity. Chen et al. [3] consider a specific class of VI model: a mapping that consists of a Lipschitz continuous and monotone operator, a Lipschitz continuous gradient mapping of a convex function, and a subgradient mapping of a simple convex function. They propose a method that combines Nesterov’s acceleration [25] with the stochastic mirror-prox method to exploit this special structure, resulting in an optimal iteration complexity for such class of problem: $\mathcal{O}\left(\frac{1}{\epsilon^{2}}\right)$, or $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ when the variance can be controlled small enough, or $\mathcal{O}\left(\sqrt{\frac{1}{\epsilon}}\right)$ when the operator consists only of gradient/subgradient mappings from some convex function. Kannan and Shanbhag [12] analyze a general variant of extra-gradient method (which uses general distance-generating functions) and show that under a slightly weaker assumptions than the strongly monotonicity, the optimal $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ iteration bound still hold. Kotsalis et al. [15] extend the OGDA method to strongly monotone stochastic VI with iteration complexity $\mathcal{O}\left(\frac{1}{\epsilon}\right)$, or $\mathcal{O}\left(\kappa\ln\left(\frac{1}{\epsilon}\right)\right)$ when the variance can be controlled small enough.

We shall note that the detailed implementation of variance-reduction is in general not considered in the above-mentioned methods (although some do present additional complexity term when the variance is small, such as in [11, 3]). Therefore, the optimal iteration complexity bound is $\mathcal{O}\left(\frac{1}{\epsilon^{2}}\right)$ for monotone VI and $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ for strongly monotone VI, as compared to $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ and $\mathcal{O}(\kappa\ln(1/\epsilon))$ for their deterministic counterpart. By increasing the sample size (aka mini-batch) in each iteration, the variance can be reduced as the algorithm progresses, therefore attaining the same optimal iteration complexity bound as the deterministic problems.

There have been developments for variance-reduction-based methods in recent years. Jalilzadeh and Shanbhag [9] extend the method [26] for deterministic strongly monotone VI to stochastic VI and show that with variance reduction the optimal iteration complexity $\mathcal{O}(\kappa\ln(1/\epsilon))$ can be achieved, together with a total sample complexity of $\mathcal{O}\left(\frac{1}{\epsilon^{\beta}}\right)$ for some constant $\beta>1$. With this method as a subroutine, they also propose a variance-reduced proximal point method with iteration complexity $\mathcal{O}\left(\frac{1}{\epsilon}\ln\left(\frac{1}{\epsilon}\right)\right)$ and sample complexity $\mathcal{O}(\frac{1}{\epsilon^{1+2\alpha\beta}})$ for some constants $\alpha,\beta>1$. Iusem et al. [7] propose a variance-reduced extra-gradient-based method for monotone VI and show $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ iteration complexity and $\mathcal{O}\left(\frac{1}{\epsilon}\right)$ sample complexity. They further extend the method [8] by incorporating line-search for unknown Lipschitz constant, while preserving similar bounds. Palaniappan and Bach [28] propose variance-reduced stochastic forward-backward methods based on (accelerated) stochastic gradient descent methods in optimization and show $\mathcal{O}(\kappa\ln(1/\epsilon))$ iteration complexity. For another line of work, which includes the concept of differential privacy in the stochastic VI, we refer the readers to a recent paper [2] and the references therein. The stochastic oracle maybe man-made. For instance, the technique of randomized smoothing has been applied in the so-called zeroth-order methods (i.e. derivative-free methods), refer to [27, 34] or the survey [16] in the context of optimization and [37, 38, 17, 33, 19] in the context of minimax saddle-point problems.

Let us start this section by introducing the notations to facilitate our analysis. We shall denote the stochastic oracle as $\hat{F}(\cdot)$, suppressing the notation $\xi$ whenever it is clear from the context. For example, $\hat{F}(z^{k})$ is associated with the random sample $\xi^{k}\in\Xi$. In addition, we denote $P_{\mathcal{Z}}(\cdot)$ as the projection operator onto the feasible set $\mathcal{Z}$.

In this section, we shall apply the proposed stochastic extra-point/extra-momentum scheme to solve the following saddle-point problem without needing to compute the gradients of $f$:

$\displaystyle\min\limits_{x\in\mathcal{X}}\max\limits_{y\in\mathcal{Y}}f(x,y),$

(30)

where $\mathcal{X}\subseteq\mathbb{R}^{n}$, $\mathcal{Y}\subseteq\mathbb{R}^{m}$ are convex sets, $f(x,y)$ is strongly convex (with fixed $y$) and strongly concave (with fixed $x$) with modulus $\mu$, and the partial gradients $\nabla_{x}f(x,y)$/$\nabla_{y}f(x,y)$ are Lipschitz continuous with constant $L_{x}/L_{y}$ for fixed $y/x$, and with constant $L_{xy}$ with fixed $x/y$. We let $L=2\cdot\max(L_{x},L_{y},L_{xy})$. We further assume that the function $f(x,y)$ is Lipschitz continuous for either fixed $x$ or $y$ with constant $M$. This implies that the norms of the partial gradients are bounded by $M$:

$$\|\nabla_{x}f(x,y)\|\leq M,\quad\|\nabla_{y}f(x,y)\|\leq M.$$

In particular, we consider the settings when the partial gradients $\nabla_{x}f(x,y)$ and $\nabla_{y}f(x,y)$ (and any higher-order information) are not available. Furthermore, the exact evaluation of the function value itself is also not available; instead, we can only access a stochastic oracle $\hat{f}(x,y,\xi)$, which satisfies the following assumption:

$\displaystyle\left\{\begin{array}[]{l}\mathbb{E}\left[\hat{f}(x,y,\xi)\right]=f(x,y),\\ \\ \mathbb{E}\left[\nabla_{x}\hat{f}(x,y,\xi)\right]=\nabla_{x}f(x,y),\\ \\ \mathbb{E}\left[\nabla_{y}\hat{f}(x,y,\xi)\right]=\nabla_{y}f(x,y),\\ \\ \mathbb{E}\left[\|\nabla_{x}\hat{f}(x,y,\xi)-\nabla_{x}f(x,y)\|^{2}\right]\leq\sigma^{2},\\ \\ \mathbb{E}\left[\|\nabla_{y}\hat{f}(x,y,\xi)-\nabla_{y}f(x,y)\|^{2}\right]\leq\sigma^{2}.\end{array}\right.$

(40)

Now, we shall use the so-called smoothing technique to approximate the first-order information, which then enables us to apply the proposed stochastic methods for VI, which includes the saddle-point model as a special case. In particular, we use a randomized smoothing scheme using uniform distributions $U_{b}$/$V_{b}$ over the unit Euclidean ball $B$ in the $\mathbb{R}^{n}$/$\mathbb{R}^{m}$ space, respectively. The smoothing functions with parameters $\rho_{x},\rho_{y}>0$ are defined as follows:

	$\displaystyle f_{\rho_{x}}(x,y)$	$\displaystyle:=$	$\displaystyle\mathbb{E}_{u\sim U_{b}}\left[f(x+\rho_{x}u,y)\right]=\frac{1}{\alpha(n)}\int_{B}f(x+\rho_{x}u,y)du,$
	$\displaystyle f_{\rho_{y}}(x,y)$	$\displaystyle:=$	$\displaystyle\mathbb{E}_{v\sim V_{b}}\left[f(x,y+\rho_{y}v)\right]=\frac{1}{\alpha(m)}\int_{B}f(x,y+\rho_{y}v)dv,$

where $\alpha(n)$/$\alpha(m)$ is the volume of the unit ball in $\mathbb{R}^{n}$/$\mathbb{R}^{m}$.

Let us summarize main properties of the smoothing functions $f_{\rho_{x}},\,f_{\rho_{y}}$ below:

We are now ready to define the stochastic zeroth-order gradient as follows:

$\displaystyle\left\{\begin{array}[]{l}F_{\rho_{x}}(x,y,\xi,u):=\frac{n}{\rho_{x}}\left(\hat{f}(x+\rho_{x}u,y,\xi)-\hat{f}(x,y,\xi)\right)u,\\ F_{\rho_{y}}(x,y,\xi,v):=\frac{m}{\rho_{y}}\left(\hat{f}(x,y+\rho_{y}v,\xi)-\hat{f}(x,y,\xi)\right)v,\end{array}\right.$

(51)

where $u$ and $v$ are the uniformly distributed random vectors over the unit spheres in $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$ respectively.

The next lemma shows that such stochastic zeroth-order gradients are unbiased with respect to the gradients of the smoothing functions and have uniformly bounded variance.

Before applying the stochastic extra-point/extra-momentum scheme to solve (30), let us first introduce the connections between these two models. As we regard the saddle-point model as a special case of VI, we shall treat the variables $x,y$ in the saddle-point problem as one variable and denote $\mathcal{Z}=\mathcal{X}\times\mathcal{Y},z=(x,y)$. Additionally, we define:

$\displaystyle F(z):=\begin{pmatrix}\nabla_{x}f(x,y)\\ -\nabla_{y}f(x,y)\end{pmatrix},\quad F_{\rho}(z):=\begin{pmatrix}\nabla_{x}f_{\rho_{x}}(x,y)\\ -\nabla_{y}f_{\rho_{y}}(x,y)\end{pmatrix},\quad\hat{F}_{\rho}(z,u,v,\xi):=\begin{pmatrix}F_{\rho_{x}}(x,y,\xi,u)\\ -F_{\rho_{y}}(x,y,\xi,v)\end{pmatrix}.$

These terms correspond to the gradient of $f(x,y)$, the gradient of the smoothing functions $f_{\rho_{x}}(x,y)$ and $f_{\rho_{y}}(x,y)$, and the stochastic zeroth-order gradient, respectively. Note that we have flipped the sign on partial gradient correspond to $y$ to account for the concavity of $f$ with respect to $y$.

Finally, as we shall use a sample size of $t_{k}\in\mathbb{N}$ (a natural number) at iteration $k$, we reserve the subscripts for the random vectors $\xi,u,v$ for the sample index $(i),\,i=1,...,t_{k}$, and denote:

$$\hat{F}^{k}_{\rho}(z^{k})=\frac{1}{t_{k}}\sum\limits_{i=1}^{t_{k}}\hat{F}_{\rho}(z^{k},u^{k}_{(i)},v^{k}_{(i)},\xi^{k}_{(i)}).$$

In the above definition we suppress the notation of the random vectors $u,v,\xi$ on the LHS for cleaner presentation. Note that by the law of large numbers, together with (55)-(59), we have

	$\displaystyle\mathbb{E}\left[\hat{F}^{k}_{\rho}(z^{k})\right]=F_{\rho}(z^{k}),$		(60)
	$\displaystyle\mathbb{E}\left[\|\hat{F}^{k}_{\rho}(z^{k})-F_{\rho}(z^{k})\|^{2}\right]\leq\frac{2\tilde{\sigma}^{2}}{t_{k}}.$		(61)