Ioan Radu

PANAYOTIS MERTIKOPOULOS • Mathias Boţ

T U Staudigl
email: m.staudigl@maastrichtuniversity.nl

Vuong

MINI-BATCH FORWARD-BACKWARD-FORWARD METHODS FOR SOLVING STOCHASTIC VARIATIONAL INEQUALITIES

Keywords: 2010 Mathematics Subject Classification. 65K15, 62L20, 90C15, 90C33 variational inequalities, stochastic approximation, forward-backwardforward algorithm, mini-batch

We develop a new stochastic algorithm for solving pseudo-monotone stochastic variational inequalities. Our method builds on Tseng's forward-backwardforward (FBF) algorithm, which is known in the deterministic literature to be a valuable alternative to Korpelevich's extragradient method when solving variational inequalities over a convex and closed set governed by pseudo-monotone, Lipschitz continuous operators. The main computational advantage of Tseng's algorithm is that it relies only on a single projection step and two independent queries of a stochastic oracle. Our algorithm incorporates a mini-batch sampling mechanism and leads to almost sure (a.s.) convergence to an optimal solution. To the best of our knowledge, this is the first stochastic look-ahead algorithm achieving this by using only a single projection at each iteration.

Introduction

In this paper we consider the following variational inequality problem, denoted as VI(T, X ), or simply VI: given a nonempty closed and convex set X ⊆ R d and a single valued map T : R d → R d , find x * ∈ X such that T (x * ), x -x * ≥ 0 for all x ∈ X .

(1.1)

We call S(T, X ) ≡ X * the set of (Stampacchia) solutions of VI(T, X ). The variational inequality problem (1.1) arises in many interesting applications in economics, game theory and engineering [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF], [START_REF] Kannan | Distributed computation of equilibria in monotone nash games via iterative regularization techniques[END_REF], [START_REF] Mertikopoulos | Stochastic mirror descent dynamics and their convergence in monotone variational inequalities[END_REF], [START_REF] Ravat | On the characterization of solution sets of smooth and nonsmooth convex stochastic nash games[END_REF], [START_REF] Scutari | Convex optimization, game theory, and variational inequality theory[END_REF], and includes as a special case first-order optimality conditions for nonlinear optimization, by choosing T = ∇f for some smooth function f . If 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems -Volume I and[END_REF] for an extensive review of applications in engineering and economics. In many instances the problem VI arises as the expected value of an underlying stochastic optimization problem whose primitives are defined on a probability space (Ω, F, P) carrying a random variable ξ : (Ω, F) → (Ξ, A) taking values in a measurable space (Ξ, A) and inducing a law P = P • ξ -1 . Given the random element ξ, consider the measurable mapping F : X × Ξ → R d , defining an integrable random vector F (x, ξ) : Ω → R d via the composition F (x, ξ)(ω) = F (x, ξ(ω)). The stochastic variational inequality problem on which we will focus in this paper is denoted by SVI and defined as follows:

Definition 1.1. Let the operator T : R d → R d be defined by

T (x) := E ξ [F (x, ξ)] := Ω F (x, ξ(ω)) dP(ω) = Ξ F (x, z) dP(z).
(1.2) Find x * ∈ 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 [START_REF] King | Asymptotic theory for solutions in statistical estimation and stochastic programming[END_REF]. 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 SVI relies on random samples of the operator F (x, ξ). 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

T N (x) = 1 N N j=1 F (x, ξ j ),
and use the resulting deterministic map T N as the input in one existing algorithm of choice. We refer to [START_REF] Shapiro | Lectures on stochastic programming: modeling and theory[END_REF] 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 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, and key terms in a deterministic algorithm, such as gradients, are replaced by unbiased estimators by drawing a fresh random variable whenever needed.The mechanism to draw a fresh sample from P is usually named a stochastic oracle (SO), which report generates a stochastic error F (x, ξ) -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 [START_REF] Jiang | Stochastic approximation approaches to the stochastic variational inequality problem[END_REF], under the assumption of strong monotonicity and continuity of the operator T . There, a proximal point algorithm of the form

X n+1 = Π X [X n + α n F (X n , ξ n )]
(1.3) is considered, where Π X denotes the Euclidean projection onto X , (ξ n ) n≥0 is a sample of P, and (α n ) n≥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 in mean square. Relaxing strong monotonicity to plain monotonicity, the recent paper [START_REF] Yousefian | On smoothing, regularization, and averaging in stochastic approximation methods for stochastic variational inequality problems[END_REF] 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 [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF] and [START_REF] Kannan | Optimal stochastic extragradient schemes for pseudomonotone stochastic variational inequality problems and their variants[END_REF]. The original Korpelevich extragradient scheme of [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF] 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

Y n = Π X [X n -α n A n+1 ] X n+1 = Π X [X n -α n B n+1 ] (SEG)
where (A n ) n≥1 , (B n ) n≥1 are stochastic estimators of T (X n ), and T (Y n ), respectively. The paper [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF] constructs these estimators by relying on a dynamic sampling strategy, where noise reduction of the estimators is achieved via a mini-batch sampling of the stochastic operators F (X n , ξ) and F (Y n , ξ). Within this mini-batch formulation, almost sure convergence of the stochastic process (X n ) n∈N to the solution set can be proven even with constant step size implementations of SEG. In addition, optimal convergence rates of O(1/N ) in terms of the mean squared residual of the VI are obtained. 1

1.1. Our Contribution. 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 method [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF], which we call the stochastic forwardbackward-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 overhead and runtime. 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. With infeasibility of a method we mean that parts of the algorithm's outputs may not satisfy state-space constraints present in the underlying optimization problem. While feasibility is a big concern in many applications (in particular in engineering and economics, where such constraints may represent technological constraints), it is not really a big problem for our method. Our numerical scheme will always provide one sequence respecting state-space constraints, and we will show that this feasible "shadow sequence" is an equally good proposal for an approximate solution.

We will make this somewhat loose statement precise in the paper. Additionally, the 1 Iusem et al. (2017), as well as our working paper, called this sampling process a variance reduction strategy. We follow the suggestion of the Associate Editor, and do not use this potentially confusing terminology anymore, and simply use the term mini-batch instead.

theoretically allowed range for step sizes is by the constant factor √ 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

2.1. Notation. For x, y ∈ R d , we denote by x, y the standard inner product, and by

x ≡ x 2 := x, x 1 2 the corresponding norm. For p ∈ [1, ∞], the p norm on R d is defined for x = (x 1 , •, x p ) as x p := ( n i=1 |x i | p ) 1 p . For a nonempty, closed and convex set E ⊆ R d , the Euclidean projector is defined as Π E (x) := argmin y∈E y -x for x ∈ R d .
All random elements are defined on a given probability space (Ω, F, P). An E-valued random variable is a (F, E)-measurable mapping f : Ω → E; we write f ∈ L 0 (Ω, F, P; E).

For every p ∈ [1, ∞], define the equivalence class of random variables f ∈ L 0 (Ω, F, P; E) with E( f p ) 1/p < ∞ as L p (Ω, F, P; E). If G ⊆ F, the conditional expectation of the random variable f ∈ L p (Ω, F, P; E) is denoted by E[f |G]. For f 1 , . . . , f k ∈ L p (Ω, F, P; E)
we denote the sigma-algebra generated by these random variables by σ(f 1 , . . . , f k ), this is the smallest sigma-algebra measuring the random variables f 1 , . . . , f k . Let (Ω, F, F = (F n ) n≥0 , P) be a complete stochastic basis. We denote by 0 (F) the set of random sequences

(ξ n ) n≥1 such for each n ∈ N, ξ n ∈ L 0 (Ω, F n , P; R). For p ∈ [1, ∞], we set p (F) {(ξ n ) n≥1 ∈ 0 (F)| n≥1 |ξ n | p < ∞ P-a.s.}.
The following properties of the Euclidean projection onto a closed convex set are well known.

Lemma 2.1. Let K ⊆ R d be a nonempty, closed and convex set. Then:

(i) Π K (x) is the unique point of K satisfying x -Π K (x), y -Π K (x) ≤ 0 for all y ∈ K; (ii) for all x ∈ R d and y ∈ K, we have Π K (x) -y 2 + Π K (x) -x 2 ≤ x -y 2 ; (iii) for all x, y ∈ R d , Π K (x) -Π K (y) ≤ x -y ; (iv) given α > 0 and T : K → R d , the set of solutions of the variational problem VI(T, K) can be expressed as S(T, K) = {x ∈ R d |x = Π K (x -αT (x))}.
Remark 2.1. 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 VI(T, K), defined by inequality (1.1). For the problems of interest in this paper, Minty and Stampacchia solutions coincide [START_REF] Cottle | Pseudo-monotone complementarity problems in hilbert space[END_REF].

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

Lemma 2.2 (Pythagorean identity). For all x, x n , x n+1 ∈ R d we have

x n+1 -x 2 + x n+1 -x n 2 -x n -x 2 = 2 x n+1 -x n , x n+1 -x .
2.2. Probabilistic Tools. We recall the Minkowski inequality: for f, g ∈ L p (Ω, F, P; E), G ⊆ F and p ∈ [1, ∞], we have

E[ f + g p |G] 1/p ≤ E[ f p |G] 1/p + E[ g p |G] 1/p . (2.1)
For the convergence analysis we will make use of the following classical lemma (Polyak 1987, Lemma 11, page 50).

Lemma 2.3 (Robbins-Siegmund). Let (Ω, F, F = (F n ) n≥0 , P) be a discrete stochastic basis. Let (v n ) n≥1 , (u n ) n≥1 ∈ 0 + (F) and (θ n ) n≥1 , (β n ) n≥1 ∈ 1 + (F) be such that for all n ≥ 0 E[v n+1 |F n ] ≤ (1 + θ n )v n -u n + β n P -a.s. .
Then (v n ) n≥0 converges a.s. to a random variable v, and (u n ) n≥1 ∈ 1 + (F). Finally, we need the celebrated Burkholder-Davis-Gundy inequality (see e.g. [START_REF] Stroock | Probability Theory: An Analytic View[END_REF]).

Lemma 2.4. Let (Ω, F, (F n ) n≥0 , P) be a discrete stochastic basis and (U n ) n≥0 a vectorvalued martingale relative to this basis. Then, for all p ∈ [1, ∞), there exists a universal constant C p > 0 such that for every N ≥ 1

E sup 0≤i≤N U i p 1/p ≤ C p E   N i=1 U i -U i-1 2 p/2   1 p .
When combined with Minkowski inequality, we obtain for all p ≥ 2 a constant C p > 0 such that for every N ≥ 1

E sup 0≤i≤N U i p 1/p ≤ C p N i=1 E ( U i -U i-1 p ) 2/p

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:

Assumption 1 (Consistency). The solution set X * ≡ S(T, X ) is non-emtpy.

Assumption 2 (Stochastic Model). The set X ⊆ R d is nonempty, closed and convex, (Ξ, A) is a measurable space and

F : R d × Ξ → R d is a Carathéodory map. 2 Assumption 3 (Lipschitz continuity). The averaged operator T (•) = E ξ [F (•, ξ)] : R d → R d is Lipschitz continuous with modulus L > 0. Assumption 4 (Pseudo-Monotonicity). The averaged operator T (•) = E ξ [F (•, ξ)] is pseudo- monotone on R d , which means ∀x, y ∈ R d : T (x), y -x ≥ 0 ⇒ T (y), y -x ≥ 0.
At each iteration, the decision maker has access to a stochastic oracle, reporting an approximation of T (x) of the form

Tn+1 (x, ξ n+1 ) 1 m n+1 mn+1 i=1 F (x, ξ (i) n+1 ) for x ∈ R d .
(3.1)

2 The mapping x → F (x, ξ) is continuous for a. e. ξ ∈ Ξ, and ξ → F (x, ξ) is measurable for all

x ∈ R d ; ξ is a random variable with values in Ξ, defined on a probability space (Ω, F , P).

The sequence (m n ) n≥1 ⊆ N determines the batch size of the stochastic oracle. The random sequence ξ n = (ξ

(1) n , . . . , ξ (mn) n
) is an i.i.d draw from P. Approximations of the form (3.1) are very common in Monte-Carlo simulation approaches, machine learning and computational statistics (see e.g. [START_REF] Atchadé | On perturbed proximal gradient algorithms[END_REF], [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF], and references therein); they are easy to obtain in case we are able to sample from the measure P. The forward-backward-forward algorithm requires two queries from the stochastic oracle in which mini-batch estimators of the averaged map T are revealed. This dynamic sampling strategy requires a sequence of integers (m n ) n≥1 (the batch size) determining the size of the data set to be processed at each iteration. The random sample on each mini-batch consists of two independent stochastic processes ξ n and η n drawn from the law P, and explicitly given by

ξ n (ξ (1) n , . . . , ξ (mn) n ) and η n (η (1) n , . . . , η (mn) n ) ∀n ≥ 1.
Given the current position X n , Algorithm SFBF queries the SO once, to obtain the estimator A n+1 Tn+1 (X n , ξ n+1 ), and then constructs the random variable Y n = Π X (X n -α n A n+1 ). Next, a second query to SO is made to obtain the estimator B n+1

Tn+1 (Y n , η n+1 ), followed by the update

X n+1 = Y n + α n (A n+1 -B n+1 ).
The pseudocode for SFBF is given in Algorithm 1.

Algorithm 1 Stochastic forward-backward-forward (SFBF) Require: step-size sequence α n ; batch size sequence m n 1: Initialize X # initialization 2: for n = 1, 2, . . . do

3:

Draw samples ξ i and η i from P (i = 1, . . . , m n ) 4:

Oracle returns A ← 1 m n mn i=1 F (X, ξ i ) # first oracle query 5: Set Y ← Π X (X -α n A) # forward-backward step 6: Oracle returns B ← 1 m n mn i=1 F (Y, η i ) # second oracle query 7: Set X ← Y + α n (A -B)
# second forward step 8: end for Observe that Algorithm SFBF is an infeasible method: the iterates (X n ) n≥0 are not necessarily elements of the admissible set X , but the "shadow sequence" (Y n ) n≥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≥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 [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Chen | Stochastic accelerated alternating direction method of multipliers with importance sampling[END_REF]) are infeasible by nature as well. In concrete applications, the infeasibility of Algorithm SFBF is not really a big problem. First, if feasibility is a strict requirement, we can always propose the shadow sequence (Y n ) n≥0 as an approximate solution. This is justified by Proposition 4.7. Moreover, Theorem 4.5 shows that the random process (X n ) n≥0 will converge to a solution almost surely. Hence, SFBF is for sure "asymptotically feasible", and always contains feasible approximate solutions in case of early stopping.

Assumption 5 (Step-size choice). The step-size sequence

(α n ) n≥0 in Algorithm SFBF satisfies 0 < α inf n≥0 α n ≤ ᾱ sup n≥1 α n < 1 √ 2L .
For n ≥ 0, we introduce the approximation error

W n+1 A n+1 -T (X n ), and Z n+1 B n+1 -T (Y n ), (3.2)
and the sub-sigma algebras (F n ) n≥0 , ( Fn ) n≥0 , defined by F 0 σ(X 0 ), and F n σ(X 0 , ξ 1 , ξ 2 , . . . , ξ n , η 1 , . . . , η n ) for all n ≥ 1, and Fn σ(X 0 , ξ 1 , . . . , ξ n , ξ n+1 , η 1 , . . . , η n ) for all n ≥ 0, respectively. Observe that F n ⊆ Fn for all n ≥ 0. We also define the filtrations F (F n ) n≥0 and F ( Fn ) n≥0 . The introduction of these two different subsigma algebras is important for many reasons. First, observe that they embody the information the learner has about the optimization problem. Indeed, the sub-sigma algebra (F n ) n≥0 corresponds to the information the decision maker has at the beginning the n-th iteration, whereas ( Fn ) n≥0 is the information the decision maker has after the first (projection)-step of the iteration. Therefore, (Y n ) n≥0 is measurable with respect to the sub-sigma algebra ( Fn ) n≥0 and (X n ) n≥0 is measurable with respect to the subsigma algebra (F n ) n≥0 . Second, we see that the process (W n ) n≥1 is F-adapted, whereas the process (Z n ) n≥1 is F-adapted, unbiased approximations relative to the respective information structures are provided:

E[W n+1 |F n ] = 0 and E[Z n+1 | Fn ] = 0 ∀n ≥ 0.
Assumption 6 (Batch Size). The batch size sequence (m n ) n≥1 satisfies ∞ n=1 1 mn < ∞. A sufficient condition on the sequence (m n ) n≥1 is that for some constant c > 0 and integer n 0 > 0, we have

m n = c • (n + n 0 ) 1+a ln(n + n 0 ) 1+b (3.3)
for a > 0 and b ≥ -1, or a = 0 and b > 0.

The next assumption is essentially the same as the variance control assumption in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF].

Assumption 7 (Variance Control). For all x ∈ R d and p ≥ 1, let

s p (x) E ξ [ F (x, ξ) -T (x) p ] 1/p .
There exist p ≥ 2, σ 0 ≥ 0, and a measurable locally bounded function σ :

X * → R + such that for all x ∈ R d and all x * ∈ X * s p (x) ≤ σ(x * ) + σ 0 x -x * .
(3.4)

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 (Duflo 1996, Kushner and[START_REF] Kushner | Stochastic Approximation Algorithms and Applications[END_REF], asks for uniformly bounded variance (UBV), i.e.

sup x∈X s 2 (x) ≤ σ. (UBV)
UBV is covered by Assumption 7 when σ 0 = 0 and sup x∈X * σ(x * ) ≤ σ. UBV is for instance valid when additive noise with finite p-th moment is assumed, that is, for some random variable ξ with E[ ξ p ] 1/p ≤ σ < ∞, we have F (x, ξ) = T (x) + ξ P-a.s. 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 [START_REF] Jofré | On variance reduction for stochastic smooth convex optimization with multiplicative noise[END_REF]). Assumption 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 [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF], [START_REF] Jofré | On variance reduction for stochastic smooth convex optimization with multiplicative noise[END_REF] make similar assumptions on the variance of the stochastic oracle. It is shown there that Assumption ( 7) is most natural in cases where the feasible set X is unbounded, and it is always satisfied when the Carathéodory functions F (•, ξ) are random Lipschitz, as illustrated with the example below.

Example 3.1. Assume for the Carathéodory map

F : R d × Ξ → R d that there exists L ∈ L p (Ω, F, P; R + ) with F (x, ξ) -F (y, ξ) ≤ L(ξ) x -y ∀x, y ∈ R d . Call L the Lipschitz constant of the map x → T (x) = E ξ [F (x, ξ)].
Then, a repeated application of the Minkowski inequality, shows that for all x ∈ R d and all x * ∈ X * we have

s p (x) ≤ E ξ [ F (x, ξ) -F (x * , ξ) p ] 1/p + s p (x * ) + T (x) -T (x * ) ≤ (E ξ [L(ξ) p ] 1/p + L) x -x * + s p (x * ).
Let σ(x * ) denote a bound on s p (x * ) and set σ 0 L + E ξ [L(ξ) p ] 1/p , to get a variance bound as required in Assumption 7.

Convergence Analysis

We consider the quadratic residual function defined by

r a (x) 2 x -Π X (x -aT (x)) 2 ∀x ∈ R d .
The reader familiar with the literature on finite-dimensional variational inequalities will recognize this immediately as the energy defined by the natural map F nat a (x)

x -Π X (x -aT (x)) (Facchinei and Pang 2003, chapter 10). It is well known that r a (x) is a merit function for VI(T, X ). Moreover, {r a (x); a > 0} is a family of equivalent merit functions for VI(T, X ), in the sense that r b (x) ≥ r a (x) for all b > a > 0 (Facchinei and Pang 2003, Proposition 10.3.6). Denote

ρ n 1 -2L 2 α 2 n ∀n ≥ 0. (4.1)
We define recursively the process (V n ) n≥0 by V 0 0 and, for all n ≥ 1,

V n+1 V n + (4 + ρ n )α 2 n W n+1 2 + 4α 2 n Z n+1 2 , so that ∆V n V n+1 -V n = (4 + ρ n )α 2 n W n+1 2 + 4α 2 n Z n+1 2 ∀n ≥ 0. (4.2)
Additionally, we define for all x ∈ R d the process (U n (x)) n≥0 given by U 0 (x) 0, and

U n+1 (x) U n (x) + 2α n Z n+1 , x -Y n ∀n ≥ 1 with corresponding increment ∆U n (x) 2α n Z n+1 , x -Y n ∀n ≥ 0.
For any reference point x ∈ R d we see that E[∆U n (x)| Fn ] = 0 for all n ≥ 0, i.e. the process (U n (x)) n≥0 is a martingale with respect to the filtration F. Since F n ⊆ Fn , the tower property implies that

E[∆U n (x)|F n ] = 0 ∀x ∈ R d ∀n ≥ 0, (4.3)
showing that it is also a F-martingale. (V n ) n≥0 is an increasing process, with increments ∆V n whose expected value is determined by the variance of the approximation error of the stochastic oracle feedback. In terms of these increment processes, we establish the following fundamental recursion.

Lemma 4.1. For all x * ∈ X * and all n ≥ 0 we have

X n+1 -x * 2 ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + ∆U n (x * ) + ∆V n P -a.s.. (4.4)
Proof. Proof. This recursive relation follows via several simple algebraic steps. Let x * ∈ X * and n ≥ 0 be fixed.

Step 1.. We have

T (x * ), y -x * ≥ 0 ∀y ∈ X .
Using that α n > 0 as well as the pseudo-monotonicity of T , we see

α n T (Y n ), Y n -x * ≥ 0.
Using the Doob decomposition in equation ( 3.2), we can rewrite this inequality as

α n B n+1 , Y n -x * ≥ α n Z n+1 , Y n -x * . (4.5) Since Y n = Π X (X n -α n A n+1 ), from Lemma 2.1(i) we conclude that x * -Y n , Y n -X n + α n A n+1 ≥ 0. (4.6)
Adding (4.5) and (4.6) gives

α n (A n+1 -B n+1 ) -X n + Y n , x * -Y n ≥ α n Z n+1 , Y n -x * , which is equivalent to x * -Y n , X n+1 -X n ≥ α n Z n+1 , Y n -x * . (4.7)
Step 2.. Using (4.7), we get

X n+1 -X n , X n+1 -x * = X n+1 -X n , Y n -x * + X n+1 -X n , X n+1 -Y n ≤ α n Z n+1 , x * -Y n + X n+1 -X n 2 + X n+1 -X n , X n -Y n = α n Z n+1 , x * -Y n + X n+1 -X n 2 -X n -Y n 2 + α n A n+1 -B n+1 , X n -Y n
where we have used the definition of X n+1 in the last equality. The Pythagorean identity (Lemma 2.2) gives us

X n+1 -x * 2 = X n -x * 2 -X n+1 -X n 2 + 2 X n+1 -X n , X n+1 -x * ≤ X n -x * 2 + X n+1 -X n 2 -2 X n -Y n 2 + 2 α n Z n+1 , x * -Y n + 2α n A n+1 -B n+1 , X n -Y n .
Step 3. Using again the definition of X n+1 , we see

X n+1 -X n 2 = Y n + α n (A n+1 -B n+1 ) -X n 2 = X n -Y n 2 + α 2 n A n+1 -B n+1 2 + 2α n A n+1 -B n+1 , Y n -X n ≤ X n -Y n 2 + 2α 2 n T (X n ) -T (Y n ) 2 + 2α 2 n W n+1 -Z n+1 2 + 2α n A n+1 -B n+1 , Y n -X n ≤ X n -Y n 2 + 2L 2 α 2 n X n -Y n 2 + 4α 2 n W n+1 2 + 4α 2 n Z n+1 2 + 2α n A n+1 -B n+1 , Y n -X n .
The first inequality is the Cauchy-Schwarz inequality. The second inequality follows from the L-Lipschitz continuity of the averaged operator T (Assumption 3), and again the Cauchy-Schwarz inequality. Combining this with the last inequality obtained in Step 2, we see that

X n+1 -x * 2 ≤ X n -x * 2 -(1 -2L 2 α 2 n ) X n -Y n 2 + 4α 2 n W n+1 2 + 4α 2 n Z n+1 2 + 2 α n Z n+1 , x * -Y n .
Step 4.. By the definition of the squared residual function, the definition of Y n and Lemma 2.1(iii), we have

r αn (X n ) 2 = X n -Π X (X n -α n T (X n )) 2 ≤ 2 X n -Y n 2 + 2 Y n -Π X (X n -α n T (X n )) 2 = 2 X n -Y n 2 + 2 Π X (X n -α n A n+1 ) -Π X (X n -α n T (X n ) 2 ≤ 2 X n -Y n 2 + 2 α n W n+1 2 . Hence, -2 X n -Y n 2 ≤ 2α 2 n W n+1 2 -r αn (X n ) 2 .
(4.8) Step 5.. Combining (4.8) with the last inequality from Step 3 and recalling Assumption 5, we conclude

X n+1 -x * 2 ≤ X n -x * 2 - 1 2 (1 -2L 2 α 2 n )r αn (X n ) 2 + (1 -2L 2 α 2 n )α 2 n W n+1 2 + 4α 2 n W n+1 2 + 4α n Z n+1 2 + 2 α n Z n+1 , x * -Y n = X n -x * 2 - ρ n 2 r αn (X n ) 2 + (4 + ρ n )(α n ) 2 W n+1 2 + 4α 2 n Z n+1 2 + 2 α n Z n+1 , x * -Y n .
The definitions of the increments associated with the martingale (U n (x * )) n≥0 and the non-decreasing process (V n ) n≥0 give the claimed result.

Remark 4.1. One can notice that in the above proof the pseudo-monotonicity of T is used only in Step 1, in order to obtain relation (4.5). Thus, as happened in [START_REF] Dang | On the convergence properties of non-euclidean extragradient methods for variational inequalities with generalized monotone operators[END_REF], [START_REF] Solodov | A new projection method for variational inequality problems[END_REF], the pseudo-monotonicity of T can actually be replaced by the weaker assumption

T (x), x -x * ≥ 0 ∀x ∈ X , x * ∈ X * .
See also [START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF] for a similar condition.

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

E[ X n+1 -x * 2 |F n ] ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + E[∆V n |F n ]. (4.9)
In order to prove convergence of the process (X n ) n≥0 , we aim to deduce a stochastic quasi-Fejér relation. For that, we need to understand the properties of the conditional expectation

E[∆V n |F n ] = (4 + ρ n )α 2 n E[ W n+1 2 |F n ] + 4α 2 n E[ Z n+1 2 |F n ] ∀n ≥ 0. Let be q ∈ [1, ∞]. The monotonicity of L q (P) L q (Ω, F, P; R) norms gives E[∆V n |F n ] ≤ E[|∆V n | q |F n ]
1 q for all n ≥ 0. By Minkowski inequality,

E[|∆V n | q |F n ] 1 q ≤ (4 + ρ n )α 2 n E[ W n+1 2q |F n ] 1/q + 4α 2 n E[ Z n+1 2q |F n ] 1/q ∀n ≥ 0.
The next lemma provides the required bounds for these expressions, and also highlights the implicit variance reduction of our method.

Lemma 4.2. Let be p ∈ [2, p]. For all n ≥ 0 we have P-a.s.

E[ W n+1 p |F n ] 1 p ≤ C p (σ(x * ) + σ 0 X n -x * ) √ m n+1 (4.10)
and

E[ Z n+1 p |F n ] 1 p ≤ C p √ m n+1 σ(x * ) + σ 0 E[ Y n -x * p |F n ] 1 p . (4.11)
In particular, in case of (UBV) with σ 0 = 0 and sup x∈X * σ(x * ) ≤ σ, both approximation errors are bounded in L p (P) by the common factor

C p σ √ mn+1 .
Proof. Proof. See Appendix A of the Appendix.

Let be p ≥ 2 and n ≥ 0. We have

E[ Y n -x * p |F n ] 1 p ≤ (1 + α n L) X n -x * + α n E[ W n+1 p |F n ] 1/p .
Hence, combining this with (4.10) for p ∈ [2, p] as in Lemma 4.2, we see that

E[ Y n -x * p |F n ] 1 p ≤ (1 + α n L) X n -x * + α n C p (σ(x * ) + σ 0 X n -x * ) √ m n+1 = 1 + α n L + α n C p σ 0 √ m n+1 X n -x * + α n C p σ(x * ) √ m n+1 . (4.12)
Plugging this inequality into (4.11), after rearranging the terms we see that

E[ Z n+1 p |F n ] 1 p ≤ C p σ(x * ) √ m n+1 1 + α n σ 0 C p √ m n+1 + X n -x * C p σ 0 √ m n+1 1 + α n L + α n C p σ 0 √ m n+1 .
We denote

G n,p C p √ m n+1 , (4.13)
such that, for all n ≥ 0 and p ∈ [2, p] we obtain the expressions

E[ W n+1 p |F n ] 1 p ≤G n,p (σ(x * ) + σ 0 X n -x * ) , (4.14) E[ Z n+1 p |F n ] 1 p ≤σ(x * )G n,p (1 + α n σ 0 G n,p ) (4.15) + σ 0 G n,p X n -x * (1 + α n L + α n σ 0 G n,p ), E[ Y n -x * p |F n ] 1 p ≤(1 + α n L + α n σ 0 G n,p ) X n -x * + α n σ(x * )G n,p . (4.16)
In case of a (UBV), we obtain from the above estimates simple upper bounds, by setting σ 0 = 0, and replacing σ(x * ) with the uniform upper bound σ. We next use these derived expressions to obtain L q (P) bounds for the error increments (∆U n (x * )) n≥1 and

(∆V n ) n≥1 , when q ∈ [1, p/2].
Lemma 4.3. Let Assumption 7 be fulfilled with p ≥ 2. For p ∈ [2, p], q = p 2 ≥ 1 and all n ≥ 0 we have

E[|∆V n | q |F n ] 1 q ≤ α 2 n G 2 n,p σ(x * ) 2 [2(4 + ρ n ) + 16 + 16α 2 n σ 2 0 G 2 n,p ] + α 2 n G 2 n,p σ 2 0 X n -x * 2 × [2(4 + ρ n ) + 8(1 + α n L + α n σ 0 G n,p ) 2 ]
(4.17)

and

E[|∆U n (x * )| q |F n ] 1 q ≤ 2α 2 n G 2 n,p σ(x * ) 2 (1 + α n G n,p σ 0 ) + 2α n G n,p σ(x * ) X n -x * [1 + α n L + α n σ 0 G n,p (3 + 2α n L) + 2α 2 n σ 2 0 G 2 n,p ] + 2α n G n,p σ 0 X n -x * 2 (1 + α n L + α n σ 0 G n,p ) 2 .
(4.18)

If (UBV) holds with variance bound σ, then these upper bounds simplify to

E[|∆V n | q |F n ] 1 q ≤ α 2 n σ2 G 2 n,p (8 + ρ n ) (4.19)
and, respectively,

E[|∆U n (x * )| q |F n ] 1 q ≤ 2α n σG n,p (1 + Lα n ) X n -x * + 2α 2 n σ2 G 2 n,p . (4.20)
Proof. Proof. Let be n ≥ 0. For q ≥ 1, we know that

E[|∆V n | q |F n ] 1 q ≤ (4 + ρ n )α 2 n E[ W n+1 p |F n ] 2 p + 4α 2 n E[ Z n+1 p |F n ] 2 p .
Using (4.14) and (4.15), and rearranging terms, we obtain (4.17). On the other hand, we have by definition

E[|∆U n (x * )| q | Fn ] 1 q ≤ 2α n Y n -x * • E[ Z n+1 q | Fn ] 1 q ≤ 2α n Y n -x * • E[ Z n+1 p | Fn ] 1 p ≤ 2α n Y n -x * G n,p σ(x * ) + 2α n G n,p σ 0 Y n -x * 2 ,
where the first estimate follows from the Cauchy-Schwarz inequality, the second one uses the monotonicity of L q (P) norms, and the third one uses eq. (A.4). Applying the operator E[•|F n ] on both sides, and using again the monotonicity of L q (P) norms, we obtain

E[|∆U n (x * )| q |F n ] 1 q ≤ 2α n G n,p σ(x * )E[ Y n -x * q |F n ] 1 q + 2α n G n,p σ 0 E[ Y n -x * p |F n ] 2 p ≤ 2α n G n,p σ(x * )E[ Y n -x * p |F n ] 1 p + 2α n G n,p σ 0 E[ Y n -x * p |F n ] 2 p .
After applying (4.16) and rearranging terms we arrive at the expression (4.18).

In case UBV holds with uniform variance bound σ, the upper bound for |∆V n+1 | q follows immediately from the defining expression (4.2) by using the uniform bounds

C p σ √ mn+1 = G n,p σ for the quadratic error terms W n+1
2 and Z n+1 2 . The corresponding bound for |∆U n (x * )| q is obtained from (4.18) by setting σ 0 = 0 and replacing σ(x * ) by its uniform upper bound σ.

Based on the previous estimates, we can now derive the announced stochastic quasi-Fejér inequality for the sequence X n -x * 2 n≥0 . Proposition 4.4. For all x * ∈ X * and all n ≥ 0, we have

E[ X n+1 -x * 2 |F n ] ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + κ n m n+1 σ 2 0 X n -x * 2 + σ(x * ) 2 , (4.21) where κ n α 2 n C 2 2 [2(4 + ρ n ) + 16(1 + α n L + α n σ 0 G n,2 ) 2 ]. If (UBV) holds with uniform variance bound σ, then E[ X n+1 -x * 2 |F n ] ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + κ n σ2 m n+1 , (4.22)
where now κ n = α 2 n C 2 2 (8 + ρ n ). Proof. Proof. Let be x * ∈ X * and n ≥ 0. Our point of departure is (4.9), together with (4.17). From here we derive that

E[ X n+1 -x * 2 |F n ] ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + α 2 n G 2 n,2 σ(x * ) 2 [2(4 + ρ n ) + 16 + 16α 2 n σ 2 0 G 2 n,2 ] + α 2 n G 2 n,2 σ 2 0 X n -x * 2 [2(4 + ρ n ) + 8(1 + α n L + α n σ 0 G n,p2 ) 2 ] ≤ X n -x * 2 - ρ n 2 r αn (X n ) 2 + σ 2 0 X n -x * 2 + σ(x * ) 2 2(4 + ρ n ) + 16(1 + α n L + α n σ 0 G n,2 ) 2 α 2 n G 2 n,2 .
In the last equality, we have used that 2(4

+ ρ n ) + 8(1 + α n L + α n σ 0 G n,2 ) 2 ≤ 2(4 + ρ n ) + 16(1 + α n L + α n σ 0 G n,2 ) 2 , and that 2(4 + ρ n ) + 16 + 16α 2 n σ 2 0 G 2 n,2 ≤ 2(4 + ρ n ) + 16(1 + α n L + α n σ 0 G n,2 ) 2 . Recalling that G n,2 = C 2 /
√ m n+1 , the proof is complete.

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

Remark 4.2. The scaling factor κ n only depends on the step size α n , the Lipschitz constant L, and the variance bound on the stochastic oracle. Let ᾱ sup n≥0 α n and α inf n≥0 α n (both finite and positive according to Assumption 5). Using the definition of ρ n in (4.1), we can bound

κ n = α 2 n C 2 2 2(4 + ρ n ) + 16(1 + α n L + α n σ 0 C 2 √ m n+1 ) 2 ≤ α 2 n C 2 2 10 + 32(1 + α n L) 2 + 32α 2 n σ 2 0 C 2 2 m n+1 ≤ ᾱ2 C 2 2 c 1 1 + ᾱ2 σ 2 0 C 2 2 m n+1 ∀n ≥ 0,
where c 1 > 1 is a constant. Combined with the batch size condition (3.3), we obtain the existence of constants c 0 and c 1 such that

κ n ≤ c 1 1 + ᾱ2 σ 2 0 C 2 2 c 0 (n + n 0 ) 1+a ln(n + n 0 ) 1+b
for all n 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≥0 converges a.s. to a random variable X with values in the set X * . 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 [START_REF] Combettes | Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping[END_REF]).

Given a stochastic process (f n ) n≥0 ⊆ L 0 (Ω, F, P; R d ), we define the (random) set of cluster points

Lim(f )(ω) {x ∈ R d |(∃(n j ) ↑ ∞) : lim nj →∞ f nj (ω) = x}
Theorem 4.5. Consider the stochastic process (X n , Y n ) n≥0 generated by Algorithm SFBF under Assumptions (1)-( 7). Then, (X n ) n≥0 converges as n → ∞ almost surely to a limit random variable X with values in X * , and

lim n→∞ E[r αn (X n ) 2 ] = 0. Proof. Proof. We fix an element x * ∈ X * . Let δ n (x * ) X n -x * 2 , u n ρn 2 r αn (X n ) 2 , θ n κnσ 2 0 mn+1
, and β n = κnσ(x * ) 2 mn+1 , so that (4.21) can be rewritten for all n ≥ 0 as

E[δ n+1 (x * )|F n ] ≤ (1 + θ n )δ n (x * ) -u n + β n P -a.s. .
Hence, by Lemma 2.3, there exists a random variable δ(

x * ) ∈ [0, ∞) such that (δ n (x * )) n≥1 → δ(x * ) a.s. as n → ∞, and P n≥0 u n < ∞ = 1. In particular, (X n ) n≥0 is bounded for almost every ω ∈ Ω. Since n≥0 u n = n≥0 ρ n r αn (X n ) 2 ≥ ρ n≥0 r αn (X n ) 2 , where ρ = 1 -2 ᾱ2 L 2 > 0, it follows that lim n→∞ r αn (X n ) = 0 P-a.s.
We next show that for all ω ∈ Ω all limit points of (X n (ω)) n≥0 are points in X * , and then apply Proposition 2.3(iii) in [START_REF] Combettes | Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping[END_REF] to conclude that (X n ) n converges almost surely to a random variable X with values in X * . Let ω ∈ Ω be such that X n (ω) is bounded. Since (α n ) n≥0 is bounded as well, we can construct subsequences (α nj ) j≥0 and (X nj (ω)) j≥0 such that lim j→∞ α nj = α ∈ [α, ᾱ] and lim j→∞ X nj (ω) = χ(ω). Additionally, we have lim j→∞ r αn j (X nj (ω)) = 0, so that

lim j→∞ X nj (ω) = lim j→∞ Π X (X nj (ω) -α nj T (X nj (ω))).
Therefore, by continuity of the projection operator and of the averaged map T , Lemma 2.1(iv) allows us to conclude that χ(ω) ∈ X * . Since the subsequence is arbitrary, it follows that Lim((X n ) n≥0 )(ω) ⊆ X * for P-almost all ω ∈ Ω. Now apply Proposition 2.3(iv) of [START_REF] Combettes | Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping[END_REF] to conclude that X n → X ∈ L 0 (Ω, F, P; X * ) P-a.s.

To prove that r αn (X n ) converges to 0 in mean square as n → ∞, observe first that

E[δ n+1 (x * )] ≤ E[δ n (x * )] - ρ n 2 E[r αn (X n ) 2 ] + κ n m n+1 σ 2 0 E[δ n (x * )] + σ(x * ) 2 ∀n ≥ 0. Let z n = E[δ n (x * )], u n = ρn 2 E[r αn (X n ) 2
] and θ n and β n be defined as in the previous paragraph. The deterministic version of the Robbins-Siegmund Lemma (Lemma 2.3) gives (u n ) n≥1 ∈ 1 + (N). Hence, lim n→∞ E[r αn (X n ) 2 ] = 0. 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 [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. 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 [START_REF] Rosasco | Stochastic forward-backward splitting for monotone inclusions[END_REF], assuming a similar variance structure on the stochastic oracle as we do. Theorem 4.5 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 X * .

Proposition 4.6. Suppose that Assumptions (1)-( 7) hold. For x * ∈ X * , set σ(x * ) max{σ(x * ), σ 0 }, and define dist(x, X * ) inf y∈X * y -x = Π X * (x) -x . For all n ≥ 0, it holds

E[dist(X n+1 , X * ) 2 |F n ] ≤ dist(X n , X * ) 2 - ρ n 2 r αn (X n ) 2 + κ n σ(Π X * (X n )) 2 m n+1 [1 + dist(X n , X * ) 2 ].
If (UBV) holds, then we get for all n ≥ 0 the uniform bound

E[dist(X n+1 , X * ) 2 |F n ] ≤ dist(X n , X * ) 2 - ρ n 2 r αn (X n ) 2 + κ n σ2 m n+1 , with κ n = α 2 n C 2 2 (8 + ρ n ). Proof. Proof. Let be π n (ω) = Π X * (X n (ω)
) for all n ≥ 0 and all ω ∈ Ω. Since the projection operator onto the closed and convex set X * is nonexpansive, we have (π n ) n≥0 ∈ 0 (F). For all n ≥ 0, we have

E[dist(X n+1 , X * ) 2 |F n ] ≤ E[ X n+1 -π n 2 |F n ] ≤ X n -π n 2 - ρ n 2 r αn (X n ) + κ n m n+1 [σ 2 0 X n -π n 2 + σ(π n ) 2 ] ≤ dist(X n , X * ) 2 - ρ n 2 r αn (X n ) + κ n σ2 (π n ) m n+1 [dist(X n , X * ) 2 + 1],
where the second inequality uses Proposition 4.4.

We are now in the position to make our claim precise that we can always provide the current iterate of the shadow sequence (Y n ) n≥0 as the approximate solution of our SVI.

Proposition 4.7. Consider the stochastic process (X n , Y n ) n≥0 generated by Algorithm SFBF under Assumptions (1)-( 7). Let (π n ) n≥0 be the X * -valued F-adapted stochastic process defined by π n = Π X * (X n ) for all n ≥ 0. Define the sequences (γ n ) n≥0 , (β n ) n≥0 by γ n = 1 + α n L + α n C2σ0 √ mn+1 , and

β n = α n C2σ(πn) √ mn+1 . If β n → 0 almost surely, then (Y n ) n converges almost surely to a X * valued random variable Y .
Proof. Proof. Let d(x) dist(x, X * ), and consider the stochastic process

(c n ) n≥0 defined as c n E[d(Y n ) 2 |F n ] 1/2
. Furthermore, we define the X * -valued random process by π n Π X * (X n ). From (4.12), we see

c n ≤ E[ Y n -π n 2 |F n ] 1 2 ≤ 1 + α n L + α n C 2 σ 0 √ m n+1 d(X n ) + α n C 2 σ(π n ) √ m n+1 = γ n d(X n ) + β n .
Taking expectations on both sides,

E[c n ] ≤ γ n E[d(X n )] + E[β n ] By Theorem 4.5, X n → X almost surely, an X * -valued random variable. Therefore, we know that lim n→∞ E[d(X n )] = 0. By hypothesis β n → 0 almost surely, so that lim sup n→∞ E[c n ] ≤ 0. Now, c 2 n = E[d(Y n ) 2 |F n ],
and by Jensen's inequality

E[c 2 n ] = E[d(Y n ) 2 ] ≤ E[c n ] 2 . Hence, lim sup n→∞ E[d(Y n ) 2 ]
≤ 0, and consequently d(Y n ) → 0 almost surely. The convergence to an X * -valued limit random variable then follows from Proposition 2.3 in [START_REF] Combettes | Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping[END_REF].

We remark that the assumption β n → 0 a.s. is rather mild. A sufficient condition is that (β n ) n≥0 ∈ 1 + (F). It trivially holds under the (UBV) assumption on the stochastic oracle's variance.

Complexity analysis and rates

The next two propositions provide explicit norm bounds on the iterates (X n ) n≥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 [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. 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.

Proposition 5.1. Suppose that Assumptions (1)-( 7) hold. For all x * ∈ X * let σ(x * ) max{σ(x * ), σ 0 }, (5.1)

a(x * ) σ2 (x * )ᾱ 2 C 2 2 c 1 .
(5.2)

Choose n 0 ∈ N and γ > 0 such that n≥n0 1 m n+1 ≤ γ (5.3) and β(x * ) γa(x * ) + γ 2 a(x * ) 2 ∈ (0, 1).

(5.4)

Then sup n≥n0+1 E[ X n -x * 2 ] ≤ E[ X n0 -x * 2 ] + 1 1 -β(x * ) .
(5.5)

Proof. Proof. Thanks to Assumption 6, for every γ > 0, we can find an index n 0 ∈ N such that (5.3) holds. Consequently, we fix n 0 ∈ N to be the smallest positive integer so that (5.3) holds. For all n ≥ 0 we denote

ψ n (x * ) E[ X n -x * 2 ].
From Proposition 4.4, we obtain

ψ n+1 (x * ) ≤ ψ n (x * ) - ρ n 2 E[r αn (X n ) 2 ] + κ n m n+1 [σ 2 0 ψ n (x * ) + σ(x * ) 2 ] ∀n ≥ 0.
Recall from Remark (4.2) that

κ n ≤ ᾱ2 C 2 2 c 1 1 + ᾱ2 σ 2 0 C 2 2 m n+1 ≤ ᾱ2 C 2 2 c 1 1 + a(x * ) c 1 m n+1 .
Using this bound, for all n ≥ n 0 + 1 the previous display telescopes to

ψ n (x * ) ≤ ψ n0 (x * ) + n-1 k=n0 (1 + ψ k (x * )) a(x * ) m k+1 + n-1 k=n0 (1 + ψ k (x * )) a(x * ) 2 c 1 m 2 k+1 .
For p > ψ n0 (x * ), define τ p (x * ) inf{n ≥ n 0 + 1|ψ n (x * ) ≥ p} ∈ N ∪ {+∞}. We claim that there exists p > ψ n0 (x * ) such that τ p(x * ) = ∞. Assuming that this is not the case, then we must have that τ p (x * ) < ∞ for all p > ψ n0 (x * ). Therefore, by definition of τ p (x * ) and (5.3), we get

p ≤ ψ τp(x * ) (x * ) ≤ ψ n0 (x * ) + τp(x * )-1 k=n0 (1 + ψ k (x * )) a(x * ) m k+1 + τp(x * )-1 k=n0 (1 + ψ k (x * )) 1 c 1 a(x * ) m k+1 2 ≤ ψ n0 (x * ) + (1 + p)γa(x * ) + (1 + p) γ 2 a(x * ) 2 c 1 .
Rearranging, and using c 1 > 1 as well as (5.4), gives

p ≤ ψ n0 (x * ) + 1 1 -γa(x * ) -γ 2 c1 a(x * ) 2 ≤ ψ n0 (x * ) + 1 1 -γa(x * ) -γ 2 a(x * ) 2 .
Since p > ψ n0 (x * ) has been chosen arbitrarily, we can let p → ∞ and obtain a contradiction. Therefore, there exists p > ψ n0 (x * ) such that p sup n≥n0+1 ψ n (x * ) ≤ p < ∞.

From here we get for all n ≥ n 0 + 1

ψ n (x * ) ≤ ψ n0 (x * ) + n-1 k=n0 (1 + ψ k (x * )) a(x * ) m k+1 + n-1 k=n0 (1 + ψ k (x * )) 1 c 1 a(x * ) m k+1 2 ≤ ψ n0 (x * ) + (1 + p)γa(x * ) + (1 + p) γ 2 a(x * ) 2 c 1 .
Taking the supremum over n ≥ n 0 + 1, and shifting back to the original expressions of the involved data, we get

p = sup n≥n0+1 E[ X n -x * 2 ] ≤ E[ X n0 -x * 2 ] + 1 1 -β(x * ) ,
which further leads to (5.5).

In case where the local variance of the SO is uniformly bounded over the solution set X * , we obtain much sharper results, allowing us to bound the distance of the iterates away from the solution set.

Proposition 5.2. Suppose that Assumptions (1)-( 7) hold. Suppose the variance over the solution set X * is bounded: σ(x * ) max{σ(x * ), σ 0 } ≤ σ for all x * ∈ X * . Define

a ᾱ2 σ2 C 2 2 c 1 . (5.6) Let φ ∈ (0, √ 5-1 2 ) and choose n 0 ≥ 1 such that i≥n0 1 mi+1 ≤ φ a . Then sup n≥n0+1 E[dist(X n , X * ) 2 ] ≤ 1 + E[dist(X n0 , X * ) 2 ] 1 -φ -φ 2 .
(5.7)

Proof. Proof. We denote by d(x) dist(x, X * ) : R d → R + the distance function of the solution set X * . Since X * is a nonempty, closed and convex subset of R d , the function d(X n ) : Ω → R + given by ω → d(X n )(ω) dist(X n (ω), X * ) is F n -measurable for all n ≥ 0. Indeed, letting π n (ω) Π X * (X n (ω)) for all n ≥ 0, then first, (π n ) n≥0 ∈ 0 + (F), and second d(X n )(ω) = X n (ω) -π n (ω) is a well-defined random process in 0 + (F), being a composition of continuous and measurable functions. Therefore, for all n ≥ 0,

E[d(X n+1 ) 2 |F n ] ≤ E[ X n+1 -π n 2 |F n ] ≤ X n -π n 2 - ρ n 2 E[r αn (X n ) 2 ] + κ n m n+1 σ 2 0 d(X n ) 2 + σ(π n ) 2 . Call ψ n E[d(X n ) 2
] for all n ≥ 0. Taking expectations in the previous display, and using the assumed uniform bound of the variance, we arrive at

ψ 2 n+1 ≤ ψ 2 n - ρ n 2 E[r αn (X n ) 2 ] + σ2 κ n m n+1 1 + ψ 2 n ∀n ≥ 0.
From Remark (4.2), we know that

κ n ≤ ᾱ2 C 2 2 c 1 1 + ᾱ2 σ2 C 2 2 m n+1 , so that σ2 κ n ≤ a(1 + a mn+1c1
) for all n ≥ 0. Hence, for all n ≥ n 0 + 1

ψ 2 n ≤ ψ 2 n0 + n-1 k=n0 (1 + ψ 2 k ) a m k+1 + n-1 k=n0 (1 + ψ 2 k ) a 2 c 1 m 2 k+1 .
From here proceed, mutatis mutandis, as in the proof of Proposition 5.1.

We next give explicit estimates of the rate of convergence and the oracle complexity of SFBF. The reported results are very similar to the extragradient method, with the important remark that all numerical constants can be improved under our forwardbackward-forward scheme. For that purpose, it is sufficient to consider Algorithm SFBF with a constant step size

α n = α ∈ 0, 1 √ 2L
for all n ≥ 0. 3 As in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF],

we can provide non-asymptotic convergence rates for the sequence

(E[r α (X n ) 2 ]) n≥0 .
For all n ≥ 0, x * ∈ X * and φ ∈ 0,

√ 5-1 2
, define

Γ n n i=0 1 m i+1 , Γ 2 n n i=0 1 m 2 i+1 , ρ = 1 -2α 2 L 2 , δ n (x * ) X n -x * 2 ,
and H(x * , n, φ)

1 + max 0≤i≤n E[δ i (x * )] 1 -φ -φ 2 .
Theorem 5.3. Suppose that Assumptions (1)-( 7) hold. Let x * ∈ X * be arbitrarily chosen, and consider Algorithm SFBF with constant step size α ∈ 0,

1 √ 2L . Choose φ ∈ 0, √ 5-1 2
and n 0 n 0 (x * ) to be the first integer such that

i≥n0 1 m i+1 ≤ φ a(x * ) , (5.8)
where a(x * ) is defined in (5.2). Let

Λ n (x * , φ) 2 ρ E[δ 0 (x * )] + (1 + H(x * , n 0 , φ)) a(x * )Γ n + a(x * ) 2 Γ 2 n , Λ ∞ (x * , φ) sup n≥0 Λ n (x * , φ).
For all ε > 0, define the stopping time

N ε inf{n ≥ 0|E[r α (X n ) 2 ] ≤ ε}.
(5.9)

Then, either N ε = 0, or E[r α (X Nε ) 2 ] ≤ ε < Λ ∞ (x * , φ) N ε .
(5.10)

Proof. Proof. Let γ = φ a(x * )
, with the constant a(x * ) defined in (5.2), and n 0 = n 0 (x * ) as required in the statement of the theorem. From Proposition 5.1, we deduce the bound

sup n≥n0+1 E[δ n (x * )] ≤ 1 + E[δ n0 (x * )] 1 -φ -φ 2 ≤ H(x * , n 0 , φ). Since 1 -φ -φ 2 ∈ (0, 1), sup 0≤i≤n0 E[δ i (x * )] ≤ H(x * , n 0 , φ). Therefore, sup n≥0 E[δ n (x * )] ≤ H(x * , n 0 , φ).
(5.11)

3 The reason for this is that {ra(x); a > 0} is a family of equivalent merit functions of VI(T, X ) (see Proposition 10.3.6 in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems -Volume I and[END_REF], and the opening to Section 4). Hence, as long as the step size policy (αn) n≥0 obeys Assumption 5, we obtain the same rate estimates.

Taking expectations in equation ( 4.21), we get

ρ 2 E[r α (X n ) 2 ] ≤ E[δ n (x * )] -E[δ n+1 (x * )] + κ n m n+1 σ(x * ) 2 + σ 2 0 E[δ n (x * )] ∀n ≥ 0.
Therefore, for all n ≥ 0,

ρ 2 n i=0 E[r α (X i ) 2 ] ≤ E[δ 0 (x * )] + n i=0 κ i m i+1 σ(x * ) 2 + σ 2 0 E[δ i (x * )] .
Using the variance bound σ(x * ) = max{σ(x * ), σ 0 }, which is well defined given the local boundedness of the variance, we get first from Remark (4.2) the bound

κ i ≤ α 2 C 2 2 c 1 1 + α 2 C 2 2 σ(x * ) 2 m i+1 ∀i ≥ 0.
Second, recalling that a(x

* ) = α 2 σ(x * ) 2 C 2 2 c 1 , it yields for all n ≥ 0 ρ 2 n i=0 E[r α (X i ) 2 ] ≤ E[δ 0 (x * )] + n i=0 a(x * ) m i+1 (1 + E[δ i (x * )]) + n i=0 1 c 1 a(x * ) m i+1 2 (1 + E[δ i (x * )]) ≤ E[δ 0 (x * )] + 1 + max 0≤i≤n E[δ i (x * )] a(x * )Γ n + a(x * ) 2 Γ 2 n .
From (5.11), we conclude

ρ 2 n i=0 E[r α (X i ) 2 ] ≤ E[δ 0 (x * )] + (1 + H(x * , n 0 , φ)) a(x * )Γ n + a(x * ) 2 Γ 2 n = ρ 2 Λ n (x * , φ) ∀n ≥ 0.
In conclusion,

n i=0 E[r α (X i ) 2 ] ≤ Λ n (x * , φ) ∀n ≥ 0.
From Theorem 4.5, we know that for all ε > 0 there exists

M ε ∈ N such that E[r α (X n ) 2 ] ≤ ε for all n ≥ M ε .
Hence, the (deterministic) stopping time N ε defined in (5.9) is either 0, or an integer bounded from above. Focussing on the latter case N ε ≥ 1, then for every

0 ≤ k ≤ N ε -1, we have ε < E[r α (X i ) 2 ].
From here, it follows

εN ε < Nε-1 i=0 E[r α (X i ) 2 ] ≤ Λ Nε-1 (x * , φ).
Hence,

E[r α (X Nε ) 2 ] ≤ ε < Λ ∞ (x * , φ) N ε .
The two cases above can be compactly summarized to statement (5.10).

We next turn to the case where the local variance is uniformly bounded over the solution set. In the previous theorem, given x * ∈ X * , the constant Λ ∞ (x * , n 0 (x * ), φ) in the convergence rate depends on the variance and on the distance of the n 0 (x * ) initial iterates to x * , where n 0 (x * ) and φ are chosen such that (5.8) holds. Assuming a uniformly bound on the variance of SO over the solution set X * , we can obtain much stronger convergence rate estimates, holding uniformly over the solution set.

Proposition 5.4. Assume that sup x * ∈X * σ(x * ) ≤ σ, where the function σ(•) is defined in (5.1). Let x * ∈ X * be arbitrarily chosen, and consider Algorithm SFBF with constant

step size α ∈ 0, 1 √ 2L . Choose φ ∈ 0, √ 5-1 2
and n 0 n 0 (σ) to be the first integer such that

i≥n0 1 m i+1 ≤ φ a ,
(5.12)

where a = σ2 α 2 C 2 2 c 1 . Let Λn (σ, φ) 2 ρ E[dist(X 0 , X * ) 2 ] + (1 + H(σ, n 0 , φ))(aΓ n + a 2 Γ 2 n ) , Λ∞ (σ, φ) = sup n≥0
Λn (φ, σ), and

H(σ, n 0 , φ)) 1 + max 0≤i≤n0(σ) E[dist(X i , X * )] 1 -φ -φ 2 .
For all ε > 0 consider the stopping time defined in (5.9). Then, either

N ε = 0, or E[r α (X Nε ) 2 ] ≤ ε < Λ∞ (σ, φ) N ε , (5.13)
Proof. Proof. 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 σ is the only parameter in this statement; hence, the threshold index n 0 = n 0 (σ) 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

n i=0 E[r α (X i ) 2 ] ≤ E[dist(X 0 , X * ) 2 ] + a n i=0 1 + E[dist(X i , X * ) 2 ] m i+1 + a 2 n i=0 1 + E[dist(X i , X * ) 2 ] m 2 i+1 ∀n ≥ 0. Proposition 5.2 gives us sup n≥n0+1 E[dist(X n , X * ) 2 ] ≤ 1 + E[δ(X n0 , X * ) 2 ] 1 -φ -φ 2 ≤ H(σ, n 0 , φ), from which it follows sup n≥0 E[dist(X n , X * ) 2 ] ≤ H(σ, n 0 , φ).
From here, we conclude just as in the proof of Theorem 5.3 that

n i=0 E[r α (X i ) 2 ] ≤ Λn (σ, φ) ≤ Λ∞ (σ, φ) ∀n ≥ 0.
Choose ε > 0 arbitrary, and consider the stopping time (5.9). Then, either N ε = 0, or else N ε ≥ 1. Focussing on the latter case, we argue just as in the proof of Theorem 5.3, that

εN ε < Nε-1 i=0 E[r α (X i ) 2 ] ≤ Λ Nε-1 (x * , φ).
Hence, if N ε not zero, we must have

E[r α (X Nε ) 2 ] ≤ ε < Λ ∞ (σ, φ) N ε .
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 ε > 0, in mean square. Hence, using the stopping time (5.9), we would like to estimate the number Nε i=0 2m i+1 . For simplicity, we will assume that the local variance function σ(x * ) is uniformly bounded over the solution set X * . That is, we assume that there exists σ ∈ (0, ∞) such that sup x∈X * σ(x) ≤ σ. A more complete argument, without making this strong assumption can be given similar to Proposition 3.23 in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. 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 n 0 (σ) 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).

Lemma 5.5. Let a be the constant defined in (5.6), and φ ∈ (0, √ 5-1

2 ). We choose the sample rate

m i = θ(µ -1 + i) ln(µ + i -1) 1+b , (5.14) for i ≥ 1, θ > 0, µ > 1 and b > 0. Then, if n 0 is an integer satisfying n 0 ≥ 1 -µ + e ( a φθb ) 1/b , we have i≥n0 1 mi+1 ≤ φ a . Proof. Proof. For n 0 ≥ 1, we compute i≥n0 1 m i+1 ≤ 1 θ i≥n0 1 (i + µ) ln(i + µ) 1+b ≤ 1 θ ∞ n0-1 1 (t + µ) ln(t + µ) 1+b dt = 1 θb ln(n 0 -1 + µ) b . Therefore, if 1 θb ln(n0-1+µ) b ≤ φ
a , we obtain the desired bound. Solving the latter inequality for n 0 gives the claimed result.

Using the sample rate (5.14), we will now bound the constant Λ(σ, φ), and the stopping time N ε . Define the constants

A µ,b α 2 C 2 2 c 1 b ln(µ -1) b , B µ,b α 4 C 4 2 c 2 1 (1 + 2b)(µ -1) ln(µ -1) 1+2b . Since, Γ ∞ ≤ 1 θb 1 ln(µ -1) b , and Γ 2 ∞ ≤ 1 θ 2 1 (2b + 1)(µ -1) ln(µ -1) 1+2b , we conclude aΓ ∞ + a 2 Γ 2 ∞ ≤ max{1, θ -2 }(A µ,b σ2 + B µ,b σ4 ).
Therefore,

Λ(σ, φ) ≤ max{1, θ -2 } 2 ρ E[dist(X 0 , X * ) 2 ] + 2 ρ (1 + H(σ, n 0 , φ)) A µ,b σ2 + B µ,b σ4 max{1, θ -2 }Q(φ, σ).
This yields the following refined uniform bound on the squared residual function.

Corollary 5.6. For all ε > 0, the stopping time N ε defined in (5.9) is either zero, or

E[r α (X Nε ) 2 ] ≤ ε < max{1, θ -2 }Q(φ, σ) N ε .
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 ε batches needed to execute Algorithm SFBF, i.e. we want to upper bound the sum 2 Nε i=0 m i . Given the definition of the sample rate in (5.14), we can perform the following computation:

Nε+1 i=1 m i ≤ max{1, θ} ln(N ε + µ + 1) 1+b Nε+1 i=1 (i -1 + µ) + (N ε + 1) ≤ max{1, θ} ln(N ε + 1 + µ) 1+b (N ε + 1) 2 (N ε + 2µ) + (N ε + 1) .
Hence,

2 Nε i=1 m i ≤ max{1, θ}(N ε + 1)(N ε + 2µ) ln(N ε + 1 + µ) 1+b + 2 N ε + 2µ
.

(5.15)

Proposition 5.7. Let ε ∈ (0, 1) be arbitrarily chosen, and µ ∈ (1, 1/ε). Define

I(σ, φ) 3 2 ρ E[dist(X 0 , X * ) 2 ] + 2 2 + 12 ρ 2 (1 + H(σ, n 0 , φ)) 2 A 2 µ,b σ4 + 12 ρ 2 (1 + H(σ, n 0 , φ)) 2 B 2 µ,b σ8 , J (σ, φ) Λ∞ (σ, n 0 , φ) + 2.
If the sample rate (m i ) i≥1 is given by (5.14), then we can bound the oracle complexity by

2 Nε+1 i=1 m i ≤ 2 max{1, θ} max{1, θ -4 }I(σ, φ) ln(J (σ, φ)/ε) 1+b + µ -1 ε 2 .
Proof. Proof. The proof is patterned after [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. Using N ε < Λ∞ (φ, σ)/ε, we continue from (5.15) to obtain the bound

2 Nε+1 i=1 m i ≤ max{1, θ} ( Λ∞ ( σ , φ) + 1)( Λ∞ (σ, φ) + 2) ε 2 ln Λ∞ (σ, φ) + 2 ε 1+b + µ -1 ≤ max{1, θ} ( Λ∞ (σ, φ) + 2) 2 ε 2 ln(ε -1 J (σ, φ)) 1+b + µ -1 .
Since,

(Λ ∞ (σ, φ) + 2) 2 ≤ max{1, θ -4 } 2 ρ E[dist(X 0 , X * ) 2 ] + 2 ρ (1 + H(σ, n 0 , φ)) A µ,b σ2 + B µ,b σ4 + 2 2 ≤ max{1, θ -4 }3 2 ρ E[dist(X 0 , X * ) 2 ] + 2 2 + 12 ρ 2 (1 + H(σ, n 0 , φ)) 2 A 2 µ,b σ4 + 12 ρ 2 max{1, θ -4 }(1 + H(σ, n 0 , φ)) 2 B 2 µ,b σ8 = max{1, θ -4 }I(σ, φ),
the result follows.

Computational Experiments

We provide three examples to verify our theoretical results and compare our methods with the SEG proposed in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. All experiments, beside Experiment 2, were generated with Matlab R2017a on a Linux OS with a 2.39 Ghz processor and 16 GB of memory. Experiment 2 was generated with Mathematica 11 on a MacBook Pro with a 2.9 Ghz processor and 16 GB memory.

6.1. Fractional programming and applications to communication networks. 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 [START_REF] Shen | Fractional programming for communication systems-part i: Power control and beamforming[END_REF]. The standard form of a stochastic fractional program is as follows:

minimize f (x) = E G(x; ξ) h(x; ξ) , subject to x ∈ X (6.1)
where G and h are positive and convex in x for all ξ. It is well known that such problems are pseudo-convex [START_REF] Boyd | Convex Optimization[END_REF], 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., h(x; ξ) h(x) = a x + b for vectors a and b of suitable dimension. Solving this problem directly involves the pseudo-monotone operator T (x) = ∇f (x). Indeed, x * ∈ X solves problem (6.1) if and only if x * solves VI(T, X ).

Experiment 1 (Quadratic fractional programming). In our first experiment, we consider functions G of the form

G(x, ξ) = 1 2 x Q(ξ)x + c(ξ) x + q(ξ),
where Q(ξ) ∈ R d×d , c(ξ) ∈ R d and q(ξ) ∈ R are randomly generated, and Q is further assumed to be positive semi-definite. More specifically, the problem data for Q is randomly generated as follows:

Q = M M + Id,
where M is a random matrix of size d × d and Id is the d × 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

Q(ξ) = Q + 1 2 V (ξ) + V (ξ) T ,
where V (ξ) is a d × d random matrix with iid entries drawn from a normal distribution with zero mean and standard derivation σ = 0.1. Similarly,

c(ξ) := c + c 1 (ξ), q(ξ) = q + q 1 (ξ), (6.2)
where c 1 (ξ) and q(ξ) are a random vector and a random number with zero mean and normal distribution with derivation σ = 0.1, respectively. Also, for the problem's feasible region, we consider box constraints of the form ..., d}, (6.3) where the lower bound a i is a random vector in (0, 1) d and the upper bound b i = a i + 10.

X = {x ∈ R d : a i ≤ x i ≤ b i i = 1,
We have implemented SEG and SFBF for this problem, using the random operator

F (x, ξ) = ∇ x G(x,ξ) h(x)
. 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 α F BF = 10/d and α EG = α F BF / √ 3. The step-size α EG is the largest one compatible with the theory developed in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF]. We choose the batch size sequence m n+1 = (n+1) 1.5 d , so that Assumption 6 is satisfied. We stop the algorithms when the residual is below a given tolerance ε. Specifically, our stopping criterion is

r n x n -Π X (x n -T (x n )) ≤ ε = 10 -3 .
Our numerical experiments involve dimension d ∈ {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 and Fig. 1. 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.

Experiment 2 (Energy efficiency in multi-antenna communications). 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 [START_REF] Feng | A survey of energy-efficient wireless communications[END_REF], [START_REF] Isheden | Framework for link-level energy efficiency optimization with informed transmitter[END_REF], [START_REF] Mertikopoulos | Learning to be green: Robust energy efficiency maximization in dynamic MIMO-OFDM systems[END_REF], 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 [START_REF] Telatar | Capacity of multi-antenna Gaussian channels[END_REF]: (6.4) where:

R(X; H) = log det Id + K k=1 H k X k H † k ,
(1) X k is the M × M Hermitian input signal covariance matrix of user k and X = (X 1 , . . . , X K ) denotes their aggregate covariance profile. As a covariance matrix, each X k is Hermitian positive semi-definite.

(2) H k is the N ×M channel matrix of user k, representing the quality of the wireless medium between user k and the receiver.

(3) Id is the N × 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

R(X) = E H [R(X; H)], (6.5)
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.,

EE(X) = R(X) K k=1 [P c k + P t k ] , (6.6)
where

(1) P t k is the transmit power of the k-th device; by elementary signal processing considerations, it is given by P t k = tr(X k ).

(2) 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 = 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 tr(X k ) ≤ P max ∀k = 1, . . . , K.

(6.7)

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

maximize EE(X) = E H [R(X; H)] P c + K k=1 tr(X k ) subject to X k 0,
tr(X k ) ≤ P max ∀k = 1, 2, . . . , K.

(6.8)

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 [START_REF] Cui | Energy-efficiency of MIMO and cooperative MIMO techniques in sensor networks[END_REF], [START_REF] Isheden | Framework for link-level energy efficiency optimization with informed transmitter[END_REF] and it captures the fundamental tradeoff 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. 2. 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 [START_REF] Hata | Empirical formula for propagation loss in land mobile radio services[END_REF]; 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 Experiment 1; also, to reduce statistical error, we performed S = 100 sample runs for each algorithm. As in the case of Experiment 1, the SFBF algorithm performs consistently better than SEG, converging to a given target value between 1.5 and 3 times faster. 6.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

○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○ ○○○○○○○○○○ ○○○○ ○ ○○ ○○ ○○○ ○○ ○○ ○○ ○○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○ □□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□ □□□□□□ □□ □ □ □□ □□ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ○ �����-�������� □ ��� � ���� ���� ���� ���� ���� �� -� ����� ����� ����� � �� ��� ���������� �������� ����� ����������� ���� ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ ◇◇◇◇◇◇◇◇◇ ◇◇◇◇◇ ◇◇◇ ◇◇◇ ◇◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇◇ ◇◇ ◇◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ ◇◇◇◇◇◇◇◇◇◇◇◇◇ ◇◇◇◇◇◇◇ ◇◇◇◇◇ ◇◇◇ ◇◇◇ ◇◇ ◇◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇

����� ����������� ���� (�=��� �������)

Figure 2. Comparison of the extra-gradient and forward-backward-forward (FBF) methods in the energy efficiency maximization problem (6.8). On the left, we considered static channels, and we ran SFBF and SEG with the same initialization. On the right, we considered ergodic channels following a Rayleigh fading model and we performed S = 100 sample runs for each algorithm; we then plotted a sample run, the sample mean, and the best and worst values at each iteration for each algorithm. In all cases, SFBF exhibits significant performance gains over SEG. bimatrix games. A bimatrix game presented in its mixed extension consists of a tuple G = ({I, II}, (u I , u II ), (S I , S II )), defined by

• the set of players {I, II};

• strategy sets

S I {p ∈ R n I + | n I i=1 p i = 1}, S II {q ∈ R n II + | n II
i=1 q i = 1}; • real valued utility functions u I (p, y) p U I q, u II (p, q) p U II q, defined by the matrices (U I , U II ), both of which are real matrices of dimension n I × n II . Recall that a pair of mixed actions (p * , q * ) is called a Nash equilibrium of the bimatrix game (U I , U II ), if

p * i > 0 ⇒ (U I q) i = max 1≤j≤n I (U I q) j and q * i > 0 ⇒ (U II p) i = max 1≤j≤n II (U II p) j .
The bimatrix game 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 * , q * ) with p * = q * . Let d n I + n II , and note that

R d ∼ = R n I × R n II , via the usual embedding of a pair (p, q) to a stacked vector in R d . Define the d × d matrix M 0 n I ,n I -U I -U II 0 n II ,n II , ( 6
.9) and consider the set

X {(x 1 , x 2 ) ∈ R n I + × R n II + |U I x 2 ≤ 1 n I and U II x 1 ≤ 1 n II }. (6.10)
It is a classical fact that a Nash equilibrium (p * , q * ) can be computed by finding a pair

(x 1 , x 2 ) = (0 n I , 0 n II ) ∈ X such that x 1 (1 n I -U I x 2 ) = 0, and x 2 (1 n II -U II x 1 ) = 0.
The payoffs of the players in equilibrium can be recovered by looking at v = (6.11) This the so-called artificial equilibrium of the game, and serves as the initial point in the most used algorithm for computing Nash equilibria in bimatrix games, the Lemke-Howson algorithm, as masterly surveyed in Von [START_REF] Stengel | Computing equilibria for two-person games[END_REF]. Defining the mapping

= x 1 •v, q = x 2 • u. It is clear that (0 n I , 0 n II ) is always a solution to the linear complementarity problem x 1 (1 n I -U I x 2 ) = 0, 1 n I -U I x 2 ≥ 0 n I , x 2 (1 n II -U II x 1 ) = 0, 1 n II -U II x 1 ≥ 0 n II .
T : R d ∼ = R n I × R n II → R d ∼ = R n I × R n II , by T (x) 1 n I 1 n II + M x (6.12)
we can reformulate the conditions (6.11) compactly as

x * ≥ 0 n and T (x * ) ≥ 0 n , x * , T (x * ) = 0. (6.13)

To turn this into a stochastic complementarity problem, we consider a stochastic Nash game [START_REF] Duvocelle | Learning in time-varying games[END_REF], [START_REF] Kannan | Distributed computation of equilibria in monotone nash games via iterative regularization techniques[END_REF], where the player set and the set of mixed actions if fixed, but the payoff functions are realizations of random matrices

U n I = U I (ξ n ), U n II = U II (ξ n ) and (ξ n
) is a random process in some set Ξ, defined on a probability space (Ω, F, P). For each n ≥ 1, we look at that random operator (6.14) and run Algorithm SFBF.

F (x, ξ n ) 1 n I 1 n II + M (ξ n )x,
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 (ξ) := M + V (ξ), where V (ξ) is a d × d random matrix with zero mean and normal distribution with derivation σ = 0.1. Since the operator T is Lipschitz continuous with modulus L = M , we run SEG and SFBF with constant stepsizes α F BF = 0.99 √ 2L , and α EG = 0.99 √ 6L , respectively. We choose the batch size sequence m n+1 = (n+1) 1.5 d so that Assumption 6 is satisfied. The same stopping criterion as in the previous 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 √ 3 times larger than that of SEG.

Experiment 3 (Zero-Sum games). We compare the performance SFBF and SEG for zero sum game, i.e., U I = -U T II . The results are displayed in Table 2, 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. Experiment 4 (Symmetric game). We compare the performance SFBF and SEG for symmetric game, i.e., U I , U II are symmetric and U I = U T II . We choose n I = n II ∈ {50, 100, 150, . . . , 500} and d = n I + n II . The results are displayed in Table 3 showing the advantage of SFBF over SEG. On average, sFBF is 1.4 times faster in computational time and 1.8 times faster in number of iterations.

Experiment 5 (Bimatrix Games). We compare the performance SFBF and SEG for asymmetric game. We choose n I ∈ {100, 200, . . . , 1000} and n II = 2n I . The results are displayed in Table 4 and Fig. 3, showing the advantage of SFBF over SEG.

Conclusion

In this paper we have developed a stochastic version of Tseng's forward-backwardforward algorithm for solving stochastic variational inequality problems over nonempty closed and convex sets. As in [START_REF] Iusem | Extragradient method with variance reduction for stochastic variational inequalities[END_REF], the current analysis can be generalized to Cartesian 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 extension would be to develop an infinite-dimensional Hilbert space version of the algorithm, and modify the basic SFBF scheme to allow for inertial and relaxation effects. We will investigate these, and other issues, in the future.

Appendix A. Proof of Lemma 4.2

We start with a general result. Let N ∈ N and ξ (1) , . . . , ξ (N ) be an i.i.d sample from the measure P. Define the process M N i (x) N i=0 by M 0 (x) 0, and for 1 ≤ i ≤ N , by

M N i (x) 1 N i n=1 F (x, ξ (n) ) -T (x) ∀x ∈ R d . (A.1)
Setting G i σ(ξ (1) , . . . , ξ (i) ), 1 ≤ i ≤ N , we see that the process {(M N i (x), G i ), 1 ≤ i ≤ N } is a martingale starting at zero. Lemma A.1. Let p ≥ 2 be as specified in Assumption 7. For all 1 ≤ q ≤ p, N ∈ N and x ∈ R d , we have

E M N N (x) q 1 q ≤ C q √ N (σ(x * ) + σ 0 x -x * ). (A.2)
Proof. Proof. For i ∈ {1, 2, . . . , N }, the monotonicity of L p (P) norms implies that

E ∆M N i-1 (x) q 1 q = 1 N E F (x, ξ (i) ) -T (x) q 1 q ≤ 1 N E F (x, ξ (i) ) -T (x) p 1 p ≤ σ(x * ) + σ 0 x -x * N .
Using this, together with Lemma 2.4, we get

E M N N (x) q 1/q ≤ C q N k=1 E F (x, ξ (k) ) -T (x) N q 2/q ≤ C q N -2 N k=1
E F (x, ξ (k) ) -T (x) q 2/q ≤ C q (σ(x * ) + σ 0 x -x * ) √ N .

Proof. Proof of Lemma (4.2). Observe that M The tower property of conditional expectations (recall that F n ⊆ Fn ) gives

E Z n+1 p |F n = E E[ Z n+1 p | Fn ]|F n ≤ C p √ m n+1 p E (σ(x * ) + σ 0 Y n -x * ) p |F n .
Finally, by the Minkowski inequality, we get

E Z n+1 p |F n 1/p ≤ C p √ m n+1 σ(x * ) + σ 0 E Y n -x * p |F n 1/p
, and our proof is complete.

Figure 1 .

 1 Figure 1. Comparison between SFBF and SEG for solving the fractional programming. We represent the residual vs. running time (left) and number of iterations (right) for one random example n = 5000.

  i , and the mixed actions defining equilibrium play are recovered by p

Figure 3 .

 3 Figure 3. Comparison between SFBF and SEG for solving the asymmetric game. We represent the Residual vs. running time (left) and number of iterations (right) for one random example nI = 1000, nII = 2000.

  mn+1 mn+1 (X n ) = W n+1 and M mn+1 mn+1 (Y n ) = Z n+1 . Hence, we immediately obtain from Lemma A.1 that E W n+1 p |F n 1/p ≤ C p (σ(x * ) + σ 0 X n -x * ) σ(x * ) + σ 0 Y n -x * ) √ m n+1 . (A.4)

Table 1 .

 1 Averaged over 100 runs for fractional problems of different size

		SFBF		SEG
	d	number of iterations time(sec.) number of iterations time(sec.)
	200	29.88	0.0473	43.96	0.0835
	500	29.84	0.2647	44.49	0.3793
	1000	30.14	1.1650	44.99	1.7017
	2000	30.54	8.0487	45.68	11.4803

Table 2 .

 2 Averaged over 100 runs for zero sum game of different size

	Dimension	SFBF		SEG
	d = n I + n II	Iterations time(sec.) Iterations time(sec.)
	n I = n II = 100	84.38	0.4421	172.42	1.4768
	n I = n II = 250	214.09	9.2088	372.80	32.4321
	n I = n II = 500	430.18	73.9068	749.65	270.5911
	n I = n II = 1000	865.67	672.0806	1508.50	2535.50

Table 3 .

 3 Averaged over 100 runs for symmetric game of different size

	Dimension	SFBF		SEG
	d = n I + n II	Iterations time(sec.) Iterations time(sec.)
	n I = n II = 100	52.00	0.3882	68.68	0.6293
	n I = n II = 250	97.96	2.589	142.55	5.1276
	n I = n II = 500	173.30	10.5297	247.30	21.0797
	n I = n II = 1000	319.92	92.0417	455.48	191.6854

Table 4 .

 4 Averaged over 100 runs for asymmetric game of different size

		SFBF		SEG
	d = n I + n II	Iterations time(sec.) Iterations time(sec.)
	n I = 100, n II = 200	100.28	1.9553	155.28	4.8202
	n I = 300, n II = 600	293.36	32.3010	466.01	90.2339
	n I = 500, n II = 1000	492.21	136.7019	779.86	394.7606
	n I = 1000, n II = 2000	992.64	1597.7266 1564.12 4655.9213

Acknowledgments R. I. Boţ and P. T. Vuong acknowledge support by the Austrian Science Fund (FWF) within the project I2419-N32 ("Employing Recent Outcomes in Proximal Theory Outside the Comfort Zone"). P. Mertikopoulos has received financial support from the French National Research Agency (ANR) under grant ORACLESS (ANR-16-CE33-0004-01). M. Staudigl and P. Mertikopoulos have been sponsored by the COST Action CA16228 "European Network for Game Theory".