Forward-backward-forward methods with variance reduction for stochastic variational inequalities

Introduction

We call $S(T,\mathcal{X})\equiv\mathcal{X}_{\ast}$ the set of (Stampacchia) solutions of $\operatorname{VI}(T,\mathcal{X})$. The variational inequality problem (1.1) arises in many interesting applications in economics, game theory and engineering , and includes as a special case first-order optimality conditions for nonlinear optimization, by choosing $T=\nabla f$ for some smooth function $f$. If $\mathcal{X}$ is unbounded, it can also be used to formulate complementarity problems, systems of equations, saddle point problems and many equilibrium problems. We refer the reader to for an extensive review of applications in engineering and economics.

Find $x^{\ast}\in\mathcal{X}$ satisfying (1.1).

This definition is known as the expected value formulation of the stochastic variational inequality problem. The expected value formulation goes back to the seminal work of . By its very definition, if the operator $T$ defined in (1.2) would be known, then the expected value formulation can be solved by any standard solution technique for deterministic variational inequalities. However, in practice, the operator $T$ is usually not directly accessible, either due to excessive computations involved in performing the integral, or because $T$ itself is the solution of an embedded subproblem. Hence, in most situations of interest, the solution of $\operatorname{SVI}$ relies on random samples of the operator $F(x,\xi)$. In this context, there are two current methodologies available; the sample average approximation (SAA) approach replaces the expected value formulation with an empirical estimator of the form

and use the resulting deterministic map $T^{N}$ as the input in one existing algorithm of choice. We refer to for this solution approach in connection with Monte Carlo simulation. We note that this approach is the standard choice in expected residual minimization problems, when ${\mathsf{P}}$ is unknown but accessible via a Monte Carlo approach.

A different methodology is the stochastic approximation (SA) approach, where samples are obtained in an online fashion, namely, the decision maker chooses one deterministic algorithm to solve the expected value formulation, and draws a fresh random variable whenever needed. The mechanism to draw a fresh sample from ${\mathsf{P}}$ is usually named a stochastic oracle (SO), which report generates a stochastic error $F(x,\xi)-T(x)$.

Until very recently, the SA approach has only been used for the expected value formulation under very restrictive assumptions. To the best of our knowledge, the first formulation of an SA approach for a stochastic VI problem was made by , under the assumption of strong monotonicity and continuity of the operator $T$. There, a proximal point algorithm of the form

is considered, where $\Pi_{\mathcal{X}}$ denotes the Euclidean projection onto $\mathcal{X}$, $(\xi_{n})_{n\geq 0}$ is a sample of ${\mathsf{P}}$, and $(\alpha_{n})_{n\geq 0}$ is a sequence of positive step sizes. Almost sure convergence of the iterates is proven for small step sizes, assuming $T$ is Lipschitz continuous and strongly monotone, and the stochastic error is uniformly bounded. Relaxing strong monotonicity to plain monotonicity, the recent paper incorporated a Tikhonov regularization scheme into the stochastic approximation algorithm (1.3) and proved almost sure convergence of the generated stochastic process. The only established method guaranteeing almost sure convergence under the significantly weaker assumption of pseudo-monotonicty of the mean operator is the extragradient approach of . The original Korpelevich extragradient scheme of consists of two projection steps using two evaluations of the deterministic map $T$ at generated test points $y_{n}$ and $x_{n}$. Extending this to the stochastic oracle case, we arrive at the stochastic extra-gradient (SEG) method

We briefly summarize the main contributions of this work. The most costly part of SEG are the two separate projection steps performed at each single iteration of the method. We show in this paper that a stochastic version of Tseng’s forward-backward-forward , which we call the stochastic forward-backward-forward (SFBF) algorithm, preserves the strong trajectory-based convergence results, while the saving of one projection step allows us to beat SEG significantly in terms of computational overheads and runtimes. In terms of convergence properties the SFBF algorithm developed in this paper has the same good properties as SEG. However, SFBF is potentially more efficient than SEG in each iteration since it relies only on a single euclidean projection step. The price to pay for this is that we obtain an infeasible method (as is typical for primal-dual schemes) with a lower computational complexity count at the positive side. Additionally, the theoretically allowed range for step sizes is by the constant factor $\sqrt{3}$ times larger than the theoretically allowed largest step size in SEG. This constant factor gain results in significant improvements in terms of the convergence speed. This will be illustrated with extensive numerical evidences reported in Section 6.

Preliminaries

The following properties of the euclidean projection on a closed and convex set are well known.

$\Pi_{K}(x)$ is the unique point of $K$ satisfying $\langle x-\Pi_{K}(x),y-\Pi_{K}(x)\rangle\leq 0$ for all $y\in K$;

In the literature on variational inequalities, there exists an alternative solution concept known as weak, or Minty, solutions. In this paper we are only interested in strong, or Stampacchia, solutions of $\operatorname{VI}(T,K)$, defined by inequality (1.1).

Another useful fact we use in this paper is the following elementary identity.

2. Probabilistic Tools

For the convergence analysis we will make use of the following classical lemma (see e.g. [26, Lemma 11, page 50]).

Finally, we need the celebrated Burkholder-Davis-Gundy inequality (see e.g. ).

When combined with Minkowski inequality, we obtain for all $p\geq 2$ a constant $C_{p}>0$ such that for every $N\geq 1$

The stochastic forward-backward-forward algorithm

In this paper we study a forward-backward-forward algorithm of Tseng type under weak monotonicity assumptions. The blanket hypotheses we consider throughout our analysis are summarized here:

The solution set $\mathcal{X}_{\ast}\equiv S(T,\mathcal{X})$ is nonemtpy.

At each iteration, the decision maker has access to a stochastic oracle, reporting an approximation of $T(x)$ of the form

Given the current position $X_{n}$, Algorithm SFBF queries the SO once, to obtain the estimator $A_{n+1}\triangleq\hat{T}_{n+1}(X_{n},\xi_{n+1})$, and then constructs the random variable $Y_{n}=\Pi_{\mathcal{X}}(X_{n}-\alpha_{n}A_{n+1})$. Next, a second query to SO is made to obtain the estimator $B_{n+1}\triangleq\hat{T}_{n+1}(Y_{n},\eta_{n+1})$, followed by the update $X_{n+1}=Y_{n}+\alpha_{n}(A_{n+1}-B_{n+1})$. The pseudocode for SFBF is given in Algorithm 1.

Observe that Algorithm SFBF is an infeasible method: the iterates $(X_{n})_{n\geq 0}$ are not necessarily elements of the admissible set $\mathcal{X}$, but the process $(Y_{n})_{n\geq 0}$ is by construction so. In the stochastic optimization case, i.e. for instances where $A_{n+1}$ is an unbiased estimator of the gradient of a real-valued function, the process $(Y_{n})_{n\geq 0}$ is seen to be a projected gradient step, where $A_{n+1}$ acts as an unbiased estimator for the stochastic gradient. This gradient step is used in an extrapolation step to generate the iterate $X_{n+1}$. We just mention that related popular primal-dual splitting schemes like ADMM are infeasible by nature as well.

The step-size sequence $(\alpha_{n})_{n\geq 0}$ in Algorithm SFBF satisfies

For $n\geq 0$, we introduce the approximation error

and the sub-sigma algebras $(\mathcal{F}_{n})_{n\geq 0},(\hat{\mathcal{F}}_{n})_{n\geq 0}$, defined by $\mathcal{F}_{0}\triangleq\sigma(X_{0})$, and

The batch size sequence $(m_{n})_{n\geq 1}$ satisfies $\sum_{n=1}^{\infty}\frac{1}{m_{n}}<\infty$.

A sufficient condition on the sequence $(m_{n})_{n\geq 1}$ is that for some constant $\mathtt{c}>0$ and integer $n_{0}>0$, we have

for $a>0$ and $b\geq-1$, or $a=0$ and $b>0$. The next assumption is essentially the same as the variance control assumption in .

Before we proceed with the convergence analysis, we want to make some clarifying remarks on this assumption. The most frequently used assumption on the SO’s approximation error, which dates back to the seminal work of Robbins and Monro (see for a textbook reference), asks for uniformly bounded variance (UBV), i.e.

However, assuming a global variance bound is not realistic in cases where the variance of the stochastic oracle depends on the position $x$ (see e.g. Example 1 in ). 7 is much weaker than UBV, as it exploits the local variance of the stochastic oracle rather than, potentially hard to estimate, global mean square variance bounds. The recent papers make similar assumptions on the variance of the stochastic oracle. It is shown there that 7 is most natural in cases where the feasible set $\mathcal{X}$ is unbounded, and it is always satisfied when the Carathéodory functions $F(\cdot,\xi)$ are random Lipschitz (see Example 3.1 below). Since Algorithm 1 is an infeasible method, we are forced to analyze the behavior of the stochastic process $(X_{n},Y_{n})_{n\geq 0}$ on an unbounded domain, which makes 7 the only realistic and convenient choice for us. Example 3.1 illustrates an important instance where 7 holds.

Convergence Analysis

We consider the quadratic residual function defined by

The reader familiar with the literature on finite-dimensional variational inequalities will recognize this immediately as the energy defined by the natural map $F^{\text{nat}}_{a}(x)\triangleq x-\Pi_{\mathcal{X}}(x-aT(x))$ [11, chapter 10]. It is well known that $r_{a}(x)$ is a merit function for $\operatorname{VI}(T,\mathcal{X})$. Moreover, $\{r_{a}(x);a>0\}$ is a family of equivalent merit functions for $\operatorname{VI}(T,\mathcal{X})$, in the sense that $r_{b}(x)\geq r_{a}(x)$ for all $b>a>0$ [11, Proposition 10.3.6]. Denote

We define recursively the process $(V_{n})_{n\geq 0}$ by $V_{0}=0$ and, for all $n\geq 1$,

For all $x^{\ast}\in\mathcal{X}_{\ast}$ and all $n\geq 0$ we have

This recursive relation follows via several simple algebraic steps. Let be $x^{\ast}\in\mathcal{X}_{\ast}$ and $n\geq 0$ fixed.

Using that $\alpha_{n}>0$ as well as the pseudo-monotonicity of $T$, we see

Using the Doob decomposition in equation (3.2), we can rewrite this inequality as

Since $Y_{n}=\Pi_{\mathcal{X}}(X_{n}-\alpha_{n}A_{n+1})$, from Lemma 2.1(i) we conclude that

Step 2

where we have used the definition of $X_{n+1}$ in the last equality. The Pythagoras identity in Lemma 2.2 gives us

Step 3.

Using again the definition of $X_{n+1}$, we see

The first inequality is the Cauchy-Schwarz inequality. The second inequality follows from the $L$-Lipschitz continuity of the averaged operator $T$ (3), and again the Cauchy-Schwarz inequality. Combining this with the last inequality obtained in Step 2, we see that

Step 4

By the definition of the squared residual function, the definition of $Y_{n}$ and Lemma 2.1(iii), we have

Step 5

Combining (4.8) with the last inequality from Step 3 and recalling 5, we conclude

The definitions of the increments associated with the martingales $(U_{n}(x^{\ast}))_{n\geq 0}$ and $(V_{n})_{n\geq 0}$ give the claimed result. $\blacksquare$

One can notice that in the above proof the pseudo-monotonicity of $T$ is used only in Step 1 of the above proof, if order to obtain relation (4.5). Thus, as happened in , the pseudo-monotonicity of $T$ can actually be replaced by the following weaker assumption

In the following, we let $p\geq 2$ be the exponent as specified in 7. Taking conditional expectations in equation (4.4) and using the martingale property (4.3), we see for all $n\geq 0$ that

In order to prove convergence of the process $(X_{n})_{n\geq 0}$, we aim to deduce a stochastic quasi-Fejér relation. For that we need to understand the properties of the conditional expectation

The next lemma provides the required bounds for these expressions, and also highlights the implicit variance reduction of our method.

Let be $p^{\prime}\geq 2$ and $n\geq 0$. We have

Hence, combining this with (4.10) for $p^{\prime}\in[2,p]$ as in Lemma 4.2, we see that

Plugging this inequality into (4.11), after rearranging the terms we see that

such that, for all $n\geq 0$ and $p^{\prime}\in[2,p]$ we obtain the expressions

Let 7 be fulfilled with $p\geq 2$. For $p^{\prime}\in[2,p]$, $q=\frac{p^{\prime}}{2}\geq 1$ and all $n\geq 0$ we have

If (UBV) holds with variance bound $\hat{\sigma}$, then these upper bounds simplify to

Let be $n\geq 0$. For $q\geq 1$, we know that

Using (4.13) and (4.14), and rearranging terms, we obtain (4.3). On the other hand, we have by definition

After applying (4.15) and rearranging terms we arrive at the expression (4.3).

In case UBV holds with uniform variance bound $\hat{\sigma}$, the upper bound for $\lvert\Delta V_{n+1}\rvert^{q}$ follows immediately from the defining expression (4.2) by using the uniform bounds $\frac{C_{p^{\prime}}\hat{\sigma}}{\sqrt{m_{n+1}}}=G_{n,p^{\prime}}\hat{\sigma}$ for the quadratic error terms $\lVert W_{n+1}\rVert^{2}$ and $\lVert Z_{n+1}\rVert^{2}$. The corresponding bound for $\lvert\Delta U_{n}(x^{\ast})\rvert^{q}$ is obtained from (4.3) by setting $\sigma_{0}=0$ and replacing $\sigma(x^{\ast})$ by its uniform upper bound $\hat{\sigma}$. $\blacksquare$

Based on the previous estimates, we can now derive the announced stochastic quasi-Fejér inequality for the sequence $\left(\lVert X_{n}-x^{\ast}\rVert^{2}\right)_{n\geq 0}$.

For all $x^{\ast}\in\mathcal{X}_{\ast}$ and all $n\geq 0$, we have

If (UBV) holds with uniform variance bound $\hat{\sigma}$, then

where now $\kappa_{n}=\alpha^{2}_{n}C_{2}^{2}(8+\rho_{n})$.

Let be $x^{\ast}\in\mathcal{X}_{\ast}$ and $n\geq 0$. Our point of departure is (4.9), together with (4.3). From here we derive that

In the last equality, we have used that $2(4+\rho_{n})+8(1+\alpha_{n}L+\alpha_{n}\sigma_{0}G_{n,2})^{2}\leq 2(4+\rho_{n})+16(1+\alpha_{n}L+\alpha_{n}\sigma_{0}G_{n,2})^{2}$, and that $2(4+\rho_{n})+16+16\alpha_{n}^{2}\sigma_{0}^{2}G_{n,2}^{2}\leq 2(4+\rho_{n})+16(1+\alpha_{n}L+\alpha_{n}\sigma_{0}G_{n,2})^{2}$. Recalling that $G_{n,2}=C_{2}/\sqrt{m_{n+1}}$, the proof is complete.

In the case where (UBV) holds, we just have to combine (4.9) with (4.18) to obtain the claimed result. $\blacksquare$

The scaling factor $\kappa_{n}$ only depends on the step size $\alpha_{n}$, the Lipschitz constant $L$, and the variance bound on the stochastic oracle. Let $\bar{\alpha}\triangleq\sup_{n\geq 0}\alpha_{n}$ and $\underline{\alpha}\triangleq\inf_{n\geq 0}\alpha_{n}$ (both finite and positive according to 5). Using the definition of $\rho_{n}$ in (4.1), we can bound

where $\mathtt{c}_{1}>1$ is a constant. Combined with the batch size condition (3.3), we obtain the existence of constants $\mathtt{c}_{0}$ and $\mathtt{c}_{1}$ such that

for all $n\gg n_{0}$. Such non-asymptotic bounds will be used in the estimation of the rate of convergence of the algorithm.

Next we will prove that the process $(X_{n})_{n\geq 0}$ converges a.s. to a random variable $X$ with values in the set $\mathcal{X}_{\ast}$. This will be obtained as a consequence of the classical Robbins-Siegmund Lemma 2.3, and recent results on the convergence of stochastic quasi-Féjer monotone sequences (Proposition 2.3 in ).

We fix an element $x^{\ast}\in\mathcal{X}_{\ast}$. Let $\delta_{n}(x^{\ast})\triangleq\lVert X_{n}-x^{\ast}\rVert^{2},u_{n}\triangleq\frac{\rho_{n}}{2}r_{\alpha_{n}}(X_{n})^{2},\theta_{n}\!\triangleq\frac{\kappa_{n}\sigma_{0}^{2}}{m_{n+1}}$, and $\beta_{n}=\frac{\kappa_{n}\sigma(x^{\ast})^{2}}{m_{n+1}}$, so that (4.20) can be rewritten for all $n\geq 0$ as

We next show that for all $\omega\in\Omega$ all limit points of $(X_{n}(\omega))_{n\geq 0}$ are points in $\mathcal{X}_{\ast}$, and then apply Proposition 2.3(iii) to conclude that $(X_{n})_{n}$ converges almost surely to a random variable $X$ with values in $\mathcal{X}_{\ast}$. Let $\omega\in\Omega$ be such that $X_{n}(\omega)$ is bounded. Since $(\alpha_{n})_{n\geq 0}$ is bounded as well, we can construct subsequences $(\alpha_{n_{j}})_{j\geq 0}$ and $(X_{n_{j}}(\omega))_{j\geq 0}$ such that $\lim_{j\to\infty}\alpha_{n_{j}}=\alpha\in[\underline{\alpha},\bar{\alpha}]$ and $\lim_{j\to\infty}X_{n_{j}}(\omega)=\chi(\omega)$. Additionally, we have $\lim_{j\to\infty}r_{\alpha_{n_{j}}}(X_{n_{j}}(\omega))=0$, so that

To prove that $r_{\alpha_{n}}(X_{n})$ converges to in mean square as $n\rightarrow\infty$, observe first that

Theorem 4.5 considerably strengthens similar results obtained via different splitting techniques. For SEG, asymptotic convergence of the iterates in the sense of Theorem 4.5 is established in Theorem 3 of . However, different to SFBF, SEG requires two costly projection steps, with the same number of oracle calls. This makes Algorithm SFBF a potentially more efficient tool, and we will demonstrate that this is actually the case empirically, as well as theoretically. Under strong monotonicity assumptions, a version of Theorem 4.5 has been recently established for a stochastic version of the classical forward-backward splitting technique in , assuming a similar variance structure on the stochastic oracle as we do. Theorem th:converge shows convergence of SFBF under the much weaker assumption of pseudo-monotonicity of the mean operator $T$.

We close this section by reporting an improved stochastic quasi-Fejér property in terms of the distance to the solution set $\mathcal{X}_{\ast}$.

Suppose that Assumptions 1-7 hold. For $x^{\ast}\in\mathcal{X}_{\ast}$ set $\hat{\sigma}(x^{\ast})\triangleq\max\{\sigma(x^{\ast}),\sigma_{0}\}$, and define $\operatorname{dist}(x,\mathcal{X}_{\ast})\triangleq\inf_{y\in\mathcal{X}_{\ast}}\lVert y-x\rVert=\lVert\Pi_{\mathcal{X}_{\ast}}(x)-x\rVert$. For all $n\geq 0$ it holds

If (UBV) holds, then we get for all $n\geq 0$ the uniform bound

with $\kappa_{n}=\alpha^{2}_{n}C^{2}_{2}(8+\rho_{n})$.

where the second inequality uses Proposition 4.4. $\blacksquare$

Complexity analysis and rates

The next two propositions provide explicit norm bounds on the iterates $(X_{n})_{n\geq 0}$. These bounds are going to be crucial to assess the convergence rate and the per-iteration complexity of the proposed algorithm. To be sure, the formal appearance of the complexity estimates derived in this section is naturally similar to the corresponding bounds derived in . However, the key observation we would like to emphasize here is that an explicit comparison between the constants involved in the upper bounds obtained for Algorithm SFBF with those appearing in SEG shows that the constants are consistently smaller. This indicates that SFBF should empirically outperform SEG. This fact is consistently observed in all our numerical experiments, and, as we show in Section 6, actually this promised gain can be quite significant.

Suppose that Assumptions 1-7 hold. For all $x^{\ast}\in\mathcal{X}_{\ast}$ let

Using this bound, for all $n\geq n_{0}+1$ the previous display telescopes to

Rearranging, and using $\mathtt{c}_{1}>1$ as well as (5.4), gives

Since $p>\psi_{n_{0}}(x^{\ast})$ has been chosen arbitrarily, we can let $p\rightarrow\infty$ and obtain a contradiction. Therefore, there exists $\hat{p}>\psi_{n_{0}}(x^{\ast})$ such that $\bar{p}\triangleq\sup_{n\geq n_{0}+1}\psi_{n}(x^{\ast})\leq\hat{p}<\infty$. From here we get for all $n\geq n_{0}+1$

Taking the supremum over $n\geq n_{0}+1$, and shifting back to the original expressions of the involved data, we get

which further leads to (5.5). $\blacksquare$

In case where the local variance of the SO is uniformly bounded over the solution set $\mathcal{X}_{\ast}$, we obtain much sharper results, allowing us to bound the distance of the iterates away from the solution set.

Suppose that Assumptions 1-7 hold. Suppose the variance over the solution set $\mathcal{X}_{\ast}$ is bounded: $\hat{\sigma}(x^{\ast})\triangleq\max\{\sigma(x^{\ast}),\sigma_{0}\}\leq\hat{\sigma}$ for all $x^{\ast}\in\mathcal{X}_{\ast}$. Define

Let $\phi\in(0,\frac{\sqrt{5}-1}{2})$ and choose $n_{0}\geq 1$ such that $\sum_{i\geq n_{0}}\frac{1}{m_{i+1}}\leq\frac{\phi}{\mathtt{a}}$. Then

so that $\hat{\sigma}^{2}\kappa_{n}\leq\mathtt{a}(1+\frac{\mathtt{a}}{m_{n+1}\mathtt{c}_{1}})$ for all $n\geq 0$. Hence, for all $n\geq n_{0}+1$

From here proceed, mutatis mutandis, as in the proof of Proposition 5.1. $\blacksquare$

For all $n\geq 0,x^{\ast}\in\mathcal{X}_{\ast}$ and $\phi\in\left(0,\frac{\sqrt{5}-1}{2}\right)$, define

Suppose that Assumptions 1-7 hold. Let $x^{\ast}\in\mathcal{X}_{\ast}$ be arbitrarily chosen, and consider Algorithm SFBF with constant step size $\alpha\in\left(0,\frac{1}{\sqrt{2}L}\right)$. Choose $\phi\in\left(0,\frac{\sqrt{5}-1}{2}\right)$ and $n_{0}\triangleq n_{0}(x^{\ast})$ to be the first integer such that

where $\mathtt{a}(x^{\ast})$ is defined in (5.2). Let

For all $\varepsilon>0$ define the stopping time

Let $\gamma=\frac{\phi}{\mathtt{a}(x^{\ast})}$, with the constant $\mathtt{a}(x^{\ast})$ defined in (5.2), and $n_{0}=n_{0}(x^{\ast})$ as required in the statement of the theorem. From Proposition 5.1, we deduce the bound

Taking expectations in equation (4.20), we get

Using the variance bound $\hat{\sigma}(x^{\ast})=\max\{\sigma(x^{\ast}),\sigma_{0}\}$, which is well defined given the local boundedness of the variance, we get first from Remark 4.2 the bound

Second, recalling that $\mathtt{a}(x^{\ast})=\alpha^{2}\hat{\sigma}(x^{\ast})^{2}C_{2}^{2}\mathtt{c}_{1},$ it yields for all $n\geq 0$

The two cases above can be compactly summarized to statement (5.10). $\blacksquare$

We next turn to the case where the local variance is uniformly bounded over the solution set. In the previous theorem, given $x^{\ast}\in\mathcal{X}_{\ast}$, the constant $\Lambda_{\infty}(x^{\ast},n_{0}(x^{\ast}),\phi)$ in the convergence rate depends on the variance and on the distance of the $n_{0}(x^{\ast})$ initial iterates to $x^{\ast}$, where $n_{0}(x^{\ast})$ and $\phi$ are chosen such that (5.8) holds. Assuming a uniformly bound on the variance of SO over the solution set $\mathcal{X}_{\ast}$, we can obtain much stronger convergence rate estimates, holding uniformly over the solution set.

Assume that $\sup_{x^{\ast}\in\mathcal{X}_{\ast}}\hat{\sigma}(x^{\ast})\leq\hat{\sigma}$, where the function $\hat{\sigma}(\cdot)$ is defined in (5.1). Let $x^{\ast}\in\mathcal{X}_{\ast}$ be arbitrarily chosen, and consider Algorithm SFBF with constant step size $\alpha\in\left(0,\frac{1}{\sqrt{2}L}\right)$. Choose $\phi\in\left(0,\frac{\sqrt{5}-1}{2}\right)$ and $n_{0}\triangleq n_{0}(\hat{\sigma})$ to be the first integer such that

where $\mathtt{a}=\hat{\sigma}^{2}\alpha^{2}C^{2}_{2}\mathtt{c}_{1}$. Let

For all $\varepsilon>0$ consider the stopping time defined in (5.9). Then, either $N_{\varepsilon}=0$, or

The proof is almost identical to the proof of Theorem 5.3, but now we will use the estimates from Proposition 4.6 and Proposition 5.2 . We first remark that the upper variance bound $\hat{\sigma}$ is the only parameter in this statement; hence, the threshold index $n_{0}=n_{0}(\hat{\sigma})$ depends on this parameter only. Once we made this choice, we can repeat all the steps involved in the proof of Theorem 5.3 verbatim, but by using Proposition 4.6 instead of Proposition 4.4, to conclude that

From here, we conclude just as in the proof of Theorem 5.3 that

Choose $\varepsilon>0$ arbitrary, and consider the stopping time (5.9). Then, either $N_{\varepsilon}=0$, or else $N_{\varepsilon}\geq 1$. Focussing on the latter case, we argue just as in the proof of Theorem 5.3, that

Hence, if $N_{\varepsilon}$ not zero, we must have

We now turn to the estimate of the oracle complexity. By this we mean the overall size of the data set needed to be processed in order to make the natural residual function smaller than a given tolerance level $\varepsilon>0$, in mean square. Hence, using the stopping time (5.9), we would like to estimate the number $\sum_{i=0}^{N_{\varepsilon}}2m_{i+1}$.

For simplicity, we will assume that the local variance function $\sigma(x^{\ast})$ is uniformly bounded over the solution set $\mathcal{X}_{\ast}$. That is, we assume that there exists $\hat{\sigma}\in(0,\infty)$ such that $\sup_{x\in\mathcal{X}_{\ast}}\hat{\sigma}(x)\leq\hat{\sigma}$. A more complete argument, without making this strong assumption can be given similar to Proposition 3.23 in . We refrain doing so, since our main aim in this paper is to illustrate the improvement in the convergence rate when using Algorithm SFBF instead of SEG, and the simplest setting is enough for this purpose. We organize the derivation of an oracle complexity estimate in two parts. First, we will show that a specific (though admissible) choice of the sample rate, allows us to give an explicit bound on the number of preliminary iterates $n_{0}\triangleq n_{0}(\hat{\sigma})$ needed to apply the general bounds reported in Proposition 5.4. Building on this insight, we directly estimate the oracle complexity.

As announced, we first establish a bound on the number of iterations we need to meet condition (5.12).

Let $\mathtt{a}$ be the constant defined in (5.6), and $\phi\in(0,\frac{\sqrt{5}-1}{2})$. We choose the sample rate

for $i\geq 1,\theta>0,\mu>1$ and $b>0$. Then, if $n_{0}$ is an integer satisfying

we have $\sum_{i\geq n_{0}}\frac{1}{m_{i+1}}\leq\frac{\phi}{\mathtt{a}}$.

Therefore, if $\frac{1}{\theta b\ln(n_{0}-1+\mu)^{b}}\leq\frac{\phi}{\mathtt{a}}$, we obtain the desired bound. Solving the latter inequality for $n_{0}$ gives the claimed result. $\blacksquare$

Using the sample rate (5.14), we will now bound the constant $\bar{\Lambda}(\hat{\sigma},\phi)$, and the stopping time $N_{\varepsilon}$. Define the constants

This yields the following refined uniform bound on the squared residual function.

For all $\varepsilon>0$, the stopping time $N_{\varepsilon}$ defined in (5.9) is either zero, or

We now turn to the estimation of the oracle complexity. To this end, we have to bound the total number of data points involved in the $N_{\varepsilon}$ batches needed to execute Algorithm SFBF, i.e. we want to upper bound the sum $2\sum_{i=0}^{N_{\varepsilon}}m_{i}$. Given the definition of the sample rate in (5.14), we can perform the following computation:

Let $\varepsilon\in(0,1)$ be arbitrarily chosen, and $\mu\in(1,1/\varepsilon)$. Define

If the sample rate $(m_{i})_{i\geq 1}$ is given by (5.14), then we can bound the oracle complexity by

The proof is patterned after . Using $N_{\varepsilon}<\bar{\Lambda}_{\infty}(\phi,\hat{\sigma})/\varepsilon$, we continue from (5.15) to obtain the bound

Computational Experiments

We provide four examples to verify our theoretical results and compare our methods with the SEG proposed in . All experiments, beside 2, were generated with Matlab R2017a on a Linux OS with a 2.39 Ghz processor and 16 GB of memory. 2 was generated with Mathematica 11 on a MacBook Pro with a 2.9 Ghz processor and 16 GB memory.

Due to its widespread use and applications, fractional programming is instrumental to operations research and engineering, ranging from network science to signal processing, wireless communications and many other related fields . The standard form of a stochastic fractional program is as follows:

where $G$ and $h$ are positive and convex in $x$ for all $\xi$. It is well known that such problems are pseudo-convex , so they fall within the general framework of this paper. In particular, one of the cases most commonly encountered in practice is when $h$ is linear in $x$ and deterministic, i.e.,

for vectors $a$ and $b$ of suitable dimension. Solving this problem directly involves the pseudo-monotone operator $T(x)=\nabla f(x)$. Indeed, $x^{\ast}\in\mathcal{X}$ solves problem (6.1) if and only if $x^{\ast}$ solves $\operatorname{VI}(T,\mathcal{X})$.

In our first experiment, we consider functions $G$ of the form

where $M$ is a random matrix of size $d\times d$ and $\operatorname{Id}$ is the $d\times d$ identity matrix. Finally, the vectors $a$ and $c$ are drawn uniformly at random from $(0,2)^{d}$, $q$ is a random number in $(1,2)$, and $b=1+4d$.

At each sample of the methods, we generate a sample matrix as

where $V(\xi)$ is a $d\times d$ random matrix with iid entries drawn from a normal distribution with zero mean and standard derivation $\sigma=0.1$. Similarly,

where $c_{1}(\xi)$ and $q(\xi)$ are a random vector and a random number with zero mean and normal distribution with derivation $\sigma=0.1$, respectively. Also, for the problem’s feasible region, we consider box constraints of the form

where the lower bound $a_{i}$ is a random vector in $(0,1)^{d}$ and the upper bound $b_{i}=a_{i}+10$. We have implemented SEG and SFBF for this problem, using the random operator $F(x,\xi)=\nabla_{x}\left(\frac{G(x,\xi)}{h(x)}\right)$. The starting point $x_{0}$ is randomly chosen in $(1,10)^{d}$. Both algorithms are run with a constant step-size policy. We fix the stepsize of SFBF and SEG as $\alpha_{FBF}=10/d$ and $\alpha_{EG}=\alpha_{FBF}/\sqrt{3}$. The step-size $\alpha_{EG}$ is the largest one compatible with the theory developed in . We choose the batch size sequence $m_{n+1}=\left[\frac{(n+1)^{1.5}}{d}\right]$, so that Assumption 6 is satisfied. We stop the algorithms when the residual is below a given tolerance $\varepsilon$. Specifically, our stopping criterion is

Our numerical experiments involve dimension $d\in\{200,500,1000,2000\}$, and for each value of $d$ we perform $10$ runs and compare the average number of iterations and CPU time. The results are displayed in Table 1, Figs. 1 and 2. It can be seen that SFBF is constantly about $1.5$ faster than SEG in both computational time and number of iterations. An interesting observation is that the number of iterations seems not to depend on the problem dimension.

Energy efficiency is one of the most important requirements for mobile systems, and it plays a crucial role in preserving battery life and reducing the carbon footprint of multi-antenna devices (i.e., wireless devices equipped with several antennas to multiplex and demultiplex received or transmitted signals).

Following , the problem can be formulated as follows: consider $K$ wireless devices (e.g., mobile phones), each equipped with $M$ transmit antennas and seeking to connect to a common base-station with $N$ receiver antennas. In this case, the users’ achievable throughput (received bits/sec) is given by the familiar Shannon–Telatar capacity formula :

$X_{k}$ is the $M\times M$ Hermitian input signal covariance matrix of user $k$ and $X=(X_{1},\dotsc,X_{K})$ denotes their aggregate covariance profile. As a covariance matrix, each $X_{k}$ is Hermitian positive semi-definite.

$H_{k}$ is the $N\times M$ channel matrix of user $k$, representing the quality of the wireless medium between user $k$ and the receiver.

$\operatorname{Id}$ is the $N\times N$ identity matrix.

In practice, because of fading and other signal attenuation factors, the channel matrices $H_{k}$ are random variables, so the users’ achievable throughput is given by

where the expectation is taken over the (often unknown) law of $H$. The system’s energy efficiency (EE) is then defined as the ratio of the users’ achievable throughput per the unit of power consumed to achieved, i.e.,

$P^{t}_{k}$ is the transmit power of the $k$-th device; by elementary signal processing considerations, it is given by $P^{t}_{k}=\operatorname{tr}(X_{k})$.

$P^{c}_{k}>0$ is a constant representing the total power dissipated in all circuit components of the $k$-th device (mixer, frequency synthesizer, digital-to-analog converter, etc.), except for transmission. For concision, we will also write $P^{c}=\sum_{k}P^{c}_{k}$ for the total circuit power dissipitated by the system.

The users’ transmit power is further constrained by the maximum output of the transmitting device, corresponding to a trace constraint of the form

Hence, putting all this together, we obtain the stochastic fractional problem:

Note that the overall problem dimension is $d=KM^{2}$. The energy efficiency objective of this problem (which, formally, has units of bits/Joule) has been widely studied in the literature and it captures the fundamental trade-off between higher spectral efficiency and increased battery life. Importantly, switching from maximization to minimization, we also see that (6.8) is of the general form (6.1), so it can be solved by applying the SFBF algorithm: in fact, given the costly projection step to the problem’s feasible region, SFBF seems ideally suited to the task.

We do so in a series of numerical experiments reported in Fig. 3. Specifically, we consider a network consisting of $K=16$ users, each with $M=4$ transmit antennas, and a common receiver with $N=128$ receive antennas. To simulate realistic network conditions, the users’ channel matrices are drawn at each update cycle from a COST Hata radio propagation model with Rayleigh fading ; to establish a baseline, we also ran an experiment with static, deterministic channels. For comparison purposes, we ran both SFBF and SEG with the same variance reduction schedule, the same number of iterations, and step-sizes chosen as in 1; also, to reduce statistical error, we performed $S=100$ sample runs for each algorithm. As in the case of 1, the SFBF algorithm performs consistently better than SEG, converging to a given target value between $1.5$ and $3$ times faster.

2. Matrix Games

As numerical illustration we investigate the performance of the algorithm to compute Nash equilibria in random matrix games. To be specific, we revisit in this experiment the problem of computing one Nash equilibrium in random two-player bimatrix games. A bimatrix game presented in its mixed extension consists of a tuple $\mathcal{G}=\left(\{I,II\},(u_{I},u_{II}),(S_{I},S_{II})\right)$, defined by

real valued utility functions $u_{I}(p,y)\triangleq p^{\top}U_{I}q,u_{II}(p,q)\triangleq q^{\top}U_{II}^{\top}y$, defined by the matrices $(U_{I},U_{II})$, both of which are real matrices of dimension $n_{I}\times n_{II}$.

Recall that a pair of mixed actions $(p^{\ast},q^{\ast})$ is called a Nash equilibrium of the bimatrix game $(U_{I},U_{II})$, if

The bimatrix game $\mathcal{G}$ is symmetric if $n_{I}=n_{II}$ and $U_{I}=U_{II}$. In symmetric games, it is natural to focus on symmetric Nash equilibria, which is a Nash equilibrium $(p^{\ast},q^{\ast})$ with $p^{\ast}=q^{\ast}$.

It is a classical fact that a Nash equilibrium $(p^{\ast},q^{\ast})$ can be computed by finding a pair $(x_{1},x_{2})\neq(\mathbf{0}_{n_{I}},\mathbf{0}_{n_{II}})\in\mathcal{X}$ such that

The payoffs of the players in equilibrium can be recovered by looking at $v=\frac{1}{\sum_{j=1}^{n_{I}}x_{1,j}},u=\frac{1}{\sum_{i=1}^{n_{II}}x_{2,i}}$, and the mixed actions defining equilibrium play are recovered by $p=x_{1}\cdot v,q=x_{2}\cdot u$. It is clear that $(\mathbf{0}_{n_{I}},\mathbf{0}_{n_{II}})$ is always a solution to the linear complementarity problem

we can reformulate the conditions (6.11) compactly as

To turn this into a stochastic complementarity problem, we consider a stochastic Nash game , where the player set and the set of mixed actions if fixed, but the payoff functions are realizations of random matrices

In our experiments, $M$ is defined as in (6.9) and $d=n_{I}+n_{II}$. Each element of the matrices $U_{I},U_{II}$ is generated randomly with uniform distribution in $(0,1)$. To setup the experiments, we generate random matrices $M(\xi):=M+V(\xi)$, where $V(\xi)$ is a $d\times d$ random matrix with zero mean and normal distribution with derivation $\sigma=0.1$. Since the operator $T$ is Lipschitz continuous with modulus $L=\|M\|$, we run SEG and SFBF with constant stepsizes $\alpha_{FBF}=\frac{0.99}{\sqrt{2}L}$, and $\alpha_{EG}=\frac{0.99}{\sqrt{6}L}$, respectively. We choose the batch size sequence $m_{n+1}=\left[\frac{(n+1)^{1.5}}{d}\right]$ so that 6 is satisfied. The same stopping criterion as in the previous experiments of Section 6.1 is used.

From the numerical experiments, we observe that the SFBF outperforms the SEG, being on average 1.7 times faster in computational time and 1.5 times faster in number of iterations. The difference becomes larger as the problem dimension increases. There are two reasons for results: firstly, SEG requires two projections per iteration while SFBF only requires one and more importantly, the stepsize of SFBF is $\sqrt{3}$ times larger than that of SEG.

We compare the performance SFBF and SEG for zero sum game, i.e., $U_{I}=-U_{II}^{T}$. The results are displayed in Table 2 and Fig. 4 showing the advantage of SFBF over SEG. On average, SFBF is 1.7 times faster in computational time and 3.4 times faster in number of iterations than SEG.

We compare the performance SFBF and SEG for symmetric game, i.e., $U_{I},U_{II}$ are symmetric and $U_{I}=U_{II}^{T}$. We choose $n_{I}=n_{II}\in\left\{50,100,150,\ldots,500\right\}$ and $d=n_{I}+n_{II}$. The results are displayed in Table 3 and Fig. 4 showing the advantage of SFBF over SEG.

We compare the performance SFBF and SEG for asymmetric game. We choose $n_{I}\in\left\{100,200,\ldots,1000\right\}$ and $n_{II}=2n_{I}$. The results are displayed in Table 4 and Fig. 5 and Fig. 6 showing the advantage of SFBF over SEG.

Conclusion

In this paper we have developed a stochastic version of Tseng’s forward-backward-forward algorithm for solving stochastic variational inequality problems over nonempty closed and convex sets. As in , the current analysis can be generalized to Cartesian $\operatorname{VI}$ problems, though have not done this explicitly. We show that the known theoretical convergence guarantees of SEG carry over to this setting, but our method consistently outperforms SEG in terms of convergence rate and complexity. We therefore believe that SFBF is a serious competitor to SEG in typical primal-dual settings, where feasibility is a minor issue. Interesting directions for the future are to test the performance of the method in other instances where variance reduction is of importance, such as in composite optimization involving a large but finite sum of functions. Another possible extenstion would be to develop an infinite-dimensional Hilbert space version of the algorithm, and modify the basic SFBF scheme to induce strong convergence of the iterates. We will investige these, and other issues, in the future.

Appendix A Auxiliary Results

Setting $\mathcal{G}_{i}\triangleq\sigma(\xi^{(1)},\ldots,\xi^{(i)}),1\leq i\leq N$, we see that the process $\{(M^{N}_{i}(x),\mathcal{G}_{i}),1\leq i\leq N\}$ is a martingale starting at zero.

Using this, together with Lemma 2.4, we get

Observe that $M^{m_{n+1}}_{m_{n+1}}(X_{n})=W_{n+1}$ and $M^{m_{n+1}}_{m_{n+1}}(Y_{n})=Z_{n+1}$. Hence, we immediately obtain from Lemma A.1 that

To prove (4.11), we notice that Lemma A.1 implies that

The tower property of conditional expectations (recall that $\mathcal{F}_{n}\subseteq\hat{\mathcal{F}}_{n}$) gives

Finally, by the Minkowski inequality, we get

and our proof is complete. $\blacksquare$

References