Computer Science Mathematics
Computer Science-2021.gz:143270
Computer Science
ee60ff0508499991b94e7e2b261cb70ff835bff8
The Last-Iterate Convergence Rate of Optimistic Mirror Descent in Stochastic Variational Inequalities

In this paper, we analyze the local convergence rate of optimistic mirror descent methods in stochastic variational inequalities, a class of optimization problems with important applications to learning theory and machine learning. Our analysis reveals an intricate relation between the algorithm's rate of convergence and the local geometry induced by the method's underlying Bregman function. We quantify this relation by means of the Legendre exponent, a notion that we introduce to measure the growth rate of the Bregman divergence relative to the ambient norm near a solution. We show that this exponent determines both the optimal step-size policy of the algorithm and the optimal rates attained, explaining in this way the differences observed for some popular Bregman functions (Euclidean projection, negative entropy, fractional power, etc.).


Introduction
Variational inequalities -and, in particular, stochastic variational inequalities -have recently attracted considerable attention in machine learning and learning theory as a flexible paradigm for "optimization beyond minimization" -i.e., for problems where finding an optimal solution does not necessarily involve minimizing a loss function. In this context, our paper examines the local rate of convergence of optimistic mirror descent (OMD), a state-of-the-art algorithmic template for solving variational inequalities (VIs) that incorporates an "optimistic" look-ahead step with a "mirror descent" apparatus relying on a suitably chosen Bregman kernel. Our contributions focus exclusively on the stochastic case, which is of central interest to learning theory. To put them in context, we begin with a general overview below and discuss more specialized references in Section 5.
General overview. Algorithms for solving variational inequalities have a very rich history in optimization that goes back at least to the original proximal point algorithm of Martinet (1970) and Rockafellar (1976); for a survey, see Facchinei & Pang (2003). At a high level, if the vector field defining the problem is bounded and strictly monotone, simple forward-backward schemes are known to converge -and if combined with a Polyak-Ruppert averaging scheme, they achieve an O(1/ √ t) rate of convergence without the strictness caveat (Bruck Jr., 1977, Passty, 1979, Nesterov, 2009, Bauschke & Combettes, 2017. If, in addition, the problem's defining vector field is Lipschitz continuous, the extra-gradient (EG) algorithm of Korpelevich (1976) achieves trajectory convergence without strict monotonicity requirements, while the method's ergodic average converges at a O(1/t) rate (Nemirovski, 2004, Nesterov, 2007. Finally, if the problem is strongly monotone, forward-backward methods achieve an O(1/t) convergence speed; and if the operator is also Lipschitz continuous, classical results in operator theory show that simple forward-backward methods suffice to achieve a linear convergence rate (Facchinei & Pang, 2003, Bauschke & Combettes, 2017. The stochastic version of the problem arises in two closely related and interconnected ways: the defining operator could have itself a stochastic structure, or the optimizer could only be able to access a stochastic estimate thereof. Here, the landscape is considerably different: For VI problems with a monotone operator, the stochastic version of the EG algorithm achieves an ergodic O(1/ √ t) convergence rate (Juditsky et al., 2011, Gidel et al., 2019. On the other hand, if the problem is strongly monotone, this rate improves to O(1/t) and it is achieved by simple forward-backward methods, with or without averaging (Nemirovski et al., 2009).
In this context, the optimistic mirror descent algorithm has been designed to meet two complementing objectives: (i) improve the dependence of the above rates on the problem's dimensionality in cases with a favorable geometry; and (ii) minimize the number of oracle queries per iteration. The first of these objectives is achieved by the "mirror descent" component: by employing a suitable distance-generating function (DGF) -like the negative entropy on the simplex or the log-barrier on the orthant) -mirror descent achieves convergence rates that are (almost) dimension-free in problems with a favorable geometry. This idea dates back to Nemirovski & Yudin (1983), and it is also the main building block of the mirror-prox (MP) algorithm which achieves order-optimal rates with two oracle queries per iteration (Nemirovski, 2004).
The "optimistic" module then clicks on top of the MD/MP template by replacing one of the two operator queries by already observed information. This "information reuse" idea was originally due to Popov (1980), and it has recently resurfaced several times in learning theory, cf. Chiang et al. (2012), Rakhlin & Sridharan (2013b,a), Daskalakis et al. (2018), Gidel et al. (2019), Hsieh et al. (2019Hsieh et al. ( , 2021, and references therein. However, the quantitative behavior of the last-iterate in such mirror methods still remains an open question.
Our contributions. In view of the aove, the aim of our paper is to examine the local convergence rate of OMD in stochastic VI problems. For generality, we focus on nonmonotone VIs, and we investigate the algorithm's convergence to local solutions that satisfy a second-order sufficient condition.
In this regard, our first finding is that the algorithm's rate of convergence depends sharply on the local geometry induced by the method's underlying Bregman function. We formalize this by introducing the notion of the Legendre exponent, which can roughly be described as the logarithmic ratio of the volume of a regular ball w.r.t. the ambient norm centered at the solution under study to that of a Bregman ball of the same radius. For example, the ordinary Euclidean version of OMD has a Legendre exponent of β = 0; on the other hand, the entropic variant has a Legendre exponent of β = 1/2 on boundary points.
We then obtain the following rates as a function of β: for arbitrarily small ε > 0. Interestingly, these guarantees undergo a first-order phase transition between the Euclidean-like phase (0 ≤ β < 1/2) and the entropy-like phase (1/2 ≤ β < 1). This coincides with the dichotomy between steep and non-steep DGFs (i.e., DGFs that are differentiable only on the interior of the problem's domain versus those that are differentiable throughout). Moreover, these rate guarantees are reached for different step-size schedules, respectively γ t = Θ(1/t 1−β ) and Θ(1/t (1−ε)/2 ), also depending on the Legendre exponent.
To the best of our knowledge, the only comparable result in the literature is the recent work of Hsieh et al. (2019) that derives an O(1/t) convergence rate for the Euclidean case (i.e., when β = 0). Along with some other recent works on optimistic gradient methods, we discuss this in detail in Section 5, after the precise statement of our results.

Problem setup and preliminaries
Problem formulation, examples, and blanket assumptions. Throughout what follows, we will focus on solving (Stampacchia) variational inequalities of the general form: In the above, K is a closed convex subset of a d-dimensional real space X with norm · , y, x denotes the canonical pairing between y ∈ Y and x ∈ X , and the defining vector field v : K → Y of the problem is a single-valued operator with values in the dual space Y := X * of X . Letting y * := max{ y, x : x ≤ 1} denote the induced dual norm on Y, we will make the following blanket assumption for v.
For concreteness, we provide below some archetypal examples of VI problems.
Example 2.1. Consider the minimization problem with f : K → R assumed C 1 -smooth. Then, letting v(x) = ∇f (x), the solutions of (VI) are precisely the Karush-Kuhn-Tucker (KKT) points of (Min), cf. Facchinei & Pang (2003). ♦ Example 2.2. A saddle-point -or min-max -problem can be stated, in normal form, as where K 1 ⊆ R d1 and K 2 ⊆ R d2 are convex and closed, and the problem's objective function Φ : K 1 × K 2 → R is again assumed to be smooth. In the game-theoretic interpretation of von Neumann (1928), x 1 is controlled by a player seeking to minimize Φ(·, x 2 ), whereas x 2 is controlled by a player seeking to maximize Φ(x 1 , ·). Accordingly, solving (SP) consists of finding a min-max point ( Min-max points may not exist if Φ is not convex-concave. In this case, one looks instead for first-order stationary points of Φ, i.e., action profiles (x * 1 , x * 2 ) ∈ K 1 × K 2 such that x * 1 is a KKT point of Φ(·, x * 2 ) and x * 2 is a KKT point of −Φ(x * 1 , ·). Letting x = (x 1 , x 2 ) and v = (∇ x1 Φ, −∇ x2 Φ), we see that the solutions of (VI) are precisely the first-order stationary points of Φ.

♦
The above examples show that not all solutions of (VI) are desirable: for example, such a solution could be a local maximum of f in the case of (Min) or a max-min point in the case of (MinMax). For this reason, we will concentrate on solutions of (VI) that satisfy the following second-order condition: Assumption 2 (Second-order sufficiency). For a solution x * of (VI), there exists a neighborhood B and a positive constant In the context of (Min), Assumption 2 implies that f grows (at least) quadratically along every ray emanating from x * ; in particular, we have for all x ∈ B (though this does not mean that f is strongly convex in B). Likewise, in the context of (MinMax), Assumption 2 gives , so x * is a local Nash equilibrium of Φ. In general, Assumption 2 guarantees that x * is the unique solution of (VI) in B. (Indeed, any other solutionx = x * of (VI) would satisfy 0

The optimistic mirror descent algorithm
In this section, we recall the optimistic mirror descent (OMD) method for solving (VI). To streamline the flow of our paper, we first discuss the type of feedback available to the optimizer.
3.1. The oracle model. OMD is a first-order method that only requires access to the problem's defining vector field via a "black-box oracle" that returns a (possibly imperfect) version of v(x) at the selected query point x ∈ K. Concretely, we focus on stochastic first-order oracles of the form where θ is a random variable taking values in some abstract measure space Θ and Err(x; θ) is an umbrella error term capturing all sources of uncertainty in the model. In practice, the oracle is called repeatedly at a sequence of query points X t ∈ K with a different random seed θ t at each time. 1 The sequence of oracle signals V t = V(X t ; θ t ) may then be written as V t = v(X t ) + U t (2) where U t = Err(X t ; θ t ) denotes the error of the oracle model at stage t. To keep track of this sequence of events, we will treat X t as a stochastic process on some complete probability space (Ω, F , P), and we will write F t := F (X 1 , . . . , X t ) ⊆ F for the history of play up to stage t (inclusive). Since the randomness entering the oracle is triggered only after a query point has been selected, we will posit throughout that θ t (and hence U t and V t ) is F t+1 -measurable -but not necessarily F t -measurable. We will also make the blanket assumption below for the oracle.
Assumption 3 (Oracle signal). The error term U t in (2) satisfies the following properties: Both assumptions are standard in the literature on stochastic methods in optimization, cf. Polyak (1987), Hazan (2012), Bubeck (2015), and references therein.
3.2. Optimistic mirror descent. With all this in hand, the optimistic mirror descent algorithm is defined in recursive form as follows: In the above, t = 1, 2, . . . , is the algorithm's iteration counter, X t and X t+1/2 respectively denote the algorithm's base and leading states at stage t, V t+1/2 is an oracle signal obtained by querying (SFO) at X t+1/2 , and P x (y) denotes the method's so-called "proxmapping". In terms of initialization, we also take X 1 = X 1/2 ∈ ri K for simplicity. All these elements are defined in detail below. At a high level, (OMD) seeks to leverage past feedback to anticipate the landscape of v and perform more informed steps in subsequent iterations. In more detail, starting at some base state X t , t = 1, 2, . . . , the algorithm first generates an intermediate, leading state X t+1/2 , and then updates the base state with oracle input from X t+1/2 -that is, V t+1/2 . This is also the main idea behind the extra-gradient algorithm of Korpelevich (1976); the key difference here is that OMD avoids making two oracle queries per iteration by using the oracle input V t−1/2 received at the previous leading state X t−1/2 as a proxy for V t (which is never requested or received by the optimizer).
The second basic component of the method is the prox-mapping P , which, loosely speaking, can be viewed as a generalized, proximal/projection operator adapted to the geometry of the problem's primitives. To define it formally, we will need the notion of a distance-generating function on K: (1) The subdifferential of h admits a continuous selection, i.e., there exists a continuous mapping ∇h : dom ∂h → Y such that ∇h(x) ∈ ∂h(x) for all x ∈ dom ∂h. (2) h is continuous and 1-strongly convex on K, i.e., for all x ∈ dom ∂h, x ′ ∈ dom h, we have: (4) For posterity, the set K h := dom ∂h will be referred to as the prox-domain of h. We also define the Bregman divergence of h as and the induced prox-mapping P : K h × Y → K h as P x (y) = arg min In the rest of this section we take a closer look at some commonly used DGFs and the corresponding prox-mappings; for simplicity, we focus on one-dimensional problems on K = [0, 1].
For a (highly incomplete) list of references on the use of this DGF in learning theory and optimization, cf. Auer et al. (1995), Beck & Teboulle (2003), Shalev-Shwartz (2011), Arora et al. (2012), Coucheney et al. (2015), Bomze et al. (2019), Lattimore & Szepesvári (2020). ♦ Example 3.3 (Tsallis entropy). An alternative to the Gibbs-Shannon negentropy is the (Tsallis, 1988, Alvarez et al., 2004, Mertikopoulos & Sandholm, 2016. Formally, to define h for q ← 1, we will employ the continuity convention t q /(q − 1) ← log t, which yields the entropic DGF of Example 3.2. Instead, for q ← 2, we readily get the Euclidean DGF of Example 3.1. The corresponding prox-mapping does not admit a closed-form expression for all p, but it can always be calculated in logarithmic time with a simple binary search. Finally, we note that the prox-domain K h of h is the entire space K for q > 1; on the other hand, for q ∈ (0, 1], we have K h = (0, 1).

♦
The above examples will play a key role in illustrating the analysis to come and we will use them as running examples throughout.

The geometry of distance-generating functions
The Bregman topology. To proceed with our analysis and the statement of our results for (OMD), we will need to take a closer look at the geometry induced on K by the choice of h. The first observation is that, by the strong convexity requirement for h, we have: Topologically, this means that the convergence topology induced by the Bregman divergence of h on K is at least as fine as the corresponding norm topology: if a sequence x t converges to p ∈ K in the Bregman sense (D(p, x t ) → 0), it also converges in the ordinary norm topology ( x t − p → 0). On the other hand, the converse of this statement is, in general, false: specifically, the norm topology could be strictly coarser than the Bregman topology. To see this, let K be the unit Euclidean ball in R d , i.e., K = {x ∈ R d : x 2 ≤ 1}. For this space, a popular choice of DGF is the Hellinger regularizer h( (Alvarez et al., 2004. Then, for all p on the boundary of K (i.e., p 2 = 1), the Bregman divergence is Shadowing the Euclidean definition, a "Hellinger sphere" of radius r centered at p is defined as the corresponding level set of Hence, by (8), this means that a point x ∈ K h is at "Hellinger distance" r relative to p if Now, if we take a sequence x t converging to p in the ordinary Euclidean sense (i.e., x t − p 2 → 0), both sides of (9) converge to 0. This has several counterintuitive consequences: (i) S r (p) is not closed in the Euclidean topology; and (ii) the "Hellinger center" p of S r (p) actually belongs to the Euclidean closure of S r (p).
As a result, we can have a sequence x t with x t − p 2 → 0, but which remains at constant Hellinger divergence relative to p. From a qualitative standpoint, this discrepancy means that D(p, x) is not a reliable measure of convergence of x to p: if D(p, x) is large, this does not necessarily mean that p and x are topologically "far". For this reason, in the analysis of non-Euclidean instances of Bregman proximal methods, it is common to make the so-called "reciprocity" assumption A straightforward verification shows that the DGFs of Examples 3.1-3.3 all satisfy (10); by contrast, this requirement of course fails for the Hellinger regularizer above. For a detailed discussion of this condition, we refer the reader to Chen & Teboulle (1993), Kiwiel The Legendre exponent. From a quantitative standpoint, even if (10)  (1) In the Euclidean case (Example 3.1): (2) In the entropic case (Example 3.2): Consequently, the rate of convergence of x t to p = 0 is drastically different if it is measured relative to the Euclidean divergence of Example 3.1 or the Kullback-Leibler divergence of Example 3.2. As an example, if we take the sequence x t = 1/t, we get D(0, x t ) = Θ(1/t 2 ) in terms of the Euclidean divergence, but D(0, x t ) = Θ(1/t) in terms of the Kullback-Leibler divergence.
This disparity is due to the fact that it might be impossible to invert the inequality (7), even to leading order. To account for this, we introduce below the notion of the Legendre exponent of h.
Definition 2. Let h be a DGF on K. Then the Legendre exponent of h at p ∈ K is defined to be the smallest β ∈ [0, 1) such that there exists a neighborhood V of p and some κ ≥ 0 such that Heuristically, the Legendre exponent measures the difference in relative size between ordinary "norm neighborhoods" in K and the corresponding "Bregman neighborhoods" induced by the sublevel sets of the Bregman divergence. It is straightforward to see that the requirement β ≥ 0 is imposed by the strong convexity of h: since D(p, x) = Ω( p − x 2 ) for x close to p, we cannot also have D(p, x) = O( p − x 2+ε ) for any ε > 0. Note also that the notion of a Legendre exponent is only relevant if h satisfies the reciprocity requirement (10): otherwise, we could have sup x∈U D(p, x) = ∞ for any neighborhood U of p in K h , in which case we say that the Legendre exponent is 1 by convention.
To see the Legendre exponents in action, we compute them for our running examples: (1) Euclidean projection (Example 3.1): D(p, x) = (p − x) 2 /2 for all p, x ∈ K, so β = 0 for all p.
(2) Negative entropy (Example 3.2): If p ∈ K h = (0, 1), a Taylor expansion with Lagrange remainder gives D(p, (3) Tsallis entropy (Example 3.3): Assume that q = 1 (the case q = 1 has been treated above). We then have two cases: if p ∈ K h = (0, 1), h is twice continuously differentiable in a neighborhood of p so that, as previously, by the Taylor formula, D(p, x) On the other hand, if p = 0, so D(0, x) = Θ(x min(q,2) ) when x goes to zero. By symmetry, the situation is the same at p = 1, so the Legendre exponent of the Tsallis entropy is β = max(0, 1 − q/2). (4) Hellinger regularizer: As we saw at the beginning of section, the Hellinger divergence (8) does not satisfy (10), so β = 1 by convention.
Of course, beyond our running examples, the same rationale applies when K is a subset of a d-dimensional Euclidean space: If p ∈ K h , then the associated Legendre exponent is typically 0 (see Lemma A.7 for a precise statement), while if p lies at the boundary of K, the exponent may be positive; we illustrate this below for the important case of the simplex with the entropic DGF.
The associated divergence (also called Kullback-Leibler divergence) can be decomposed as The first sum is simply i:pi=0 x i and the terms of the second are bounded by O( If

Analysis and results
5.1. Main result and discussion. We are finally in a position to state and prove our main result for the convergence rate of (OMD). Since we make no global monotonicity assumptions for the defining vector field v, this result is de facto of a local nature. We provide below the complete statement, followed by a series of remarks and an outline of the main steps of the proof.
Theorem 1. Let x * be a solution of (VI) and fix some confidence level δ > 0. Suppose further that Assumptions 1-3 hold and (OMD) is run with a step-size of the form γ t = γ/(t + t 0 ) η with η ∈ (1/2, 1] and γ, t 0 > 0. Then, (OMD) enjoys the following properties: Stochastic stability: If γ/t 0 is small enough, x * is stochastically stable: for every neighborhood U of x * , there exists a neighborhood U 1 of x * such that the event occurs with probability at least 1 − δ if X 1 ∈ U 1 . Convergence: X t converges to x * at the following rate: Case I: If the Legendre exponent of h at x * is β ∈ (0, 1) and min{γ, 1/t 0 } is small enough, the iterates X t of (OMD) enjoy the convergence rate provided that U is small enough and X 1 ∈ U 1 . In particular, this rate undergoes a phase transition at β = 1/2: • Euclidean-like phase: If β ∈ (0, 1/2) and η = 1 − β, we obtain the optimized rate • Entropy-like phase: If β ∈ [1/2, 1) and η = (1 + ǫ)/2 for any small ε > 0, we get Case II: If β = 0 and γ and t 0 are sufficiently large, the iterates X t of (OMD) enjoy the rate provided that U is small enough and X 1 ∈ U 1 . In particular, for η = 1 we get Related work. The rate guarantees of Theorem 1 should be compared and contrasted to a series of recent results on the convergence speed of OMD and its variants. Most of these results concern the deterministic case, i.e., when U t = 0 for all t. In this regime, and focusing on the last iterate of the method, Golowich et al. (2020a,b) showed that the convergence speed of the (optimistic) extra-gradient algorithm in unconstrained, smooth, If, in addition, v is strongly monotone, it is well known that this rate becomes linear; for a recent take on different (Euclidean) variants of (OMD), see Malitsky (2015) Unsurprisingly, the stochastic regime is fundamentally different. In the merely monotone case (i.e., neither strongly nor strictly monotone), we are not aware of any result on the method's last-iterate convergence speed: existing results in the literature either cover the algorithm's ergodic average (Juditsky et al., 2011, Gidel et al., 2019, Cui & Shanbhag, 2016 or are asymptotic in nature (Iusem et al., 2017). Because of the lack of deterministic correlation between gradients at the base and leading states in stochastic problems, smoothness does not help to achieve better rates in this case -except for mollifying the requirement for bounded gradient signals (Juditsky et al., 2011). By contrast, strong monotonicity -local or global -does help: as was shown by Hsieh et al. (2019), the iterates of the Euclidean version of (OMD) run with a 1/t step-size schedule converge locally to solutions satisfying (SOS) at a O(1/t) rate in terms of the mean square   Fig. 1a indicates that the rate prediction is only valid up to an arbitrarily small exponent or suitable logarithmic factor. The experiments in Fig. 1b were conducted for a linear vector field v(x) = x on K = [0, +∞) with the regularizer setups of Examples 3.1-3.3 (for more details, see Appendix D). We run OMD with the stepsizes prescribed by Theorem 1 and plot D(x * , Xt) × t ν where ν is the predicted rate (ν = 1, 0.5, 0.167, 0.75 respectively for Euclidean, entropy, Tsallis entropy with q = 0.5 and q = 1.5). A horizontal line corresponds to a match between the theoretical and observed rates.
distance to the solution. This result is a special instance of Case II of Theorem 1: in particular, for β = 0, η = 1, (19) indicates that the 1/t step-size schedule examined by Hsieh et al. (2019) is optimal for the Euclidean case; however, as we explain below, it is not optimal for non-Euclidean DGFs.
Step-size tuning. A key take-away from Theorem 1 is that the step-size schedule that optimizes the rate guarantees in Eqs. (17) and (18) is either γ t = Θ(1/t 1−β ) when 0 ≤ β < 1/2, or as close as possible to Θ(t 1/2 ) when 1/2 ≤ β < 1; in both cases, the optimal tuning of the algorithm's step-size depends on the Legendre exponent of x * . We feel that this is an important parameter to keep in mind when deploying non-Euclidean versions of (OMD) in a stochastic setting -e.g., as in the optimistic multiplicative weights update (OMWU) algorithm of Daskalakis & Panageas (2019). Theorem 1 suggests that the "best practices" of the Euclidean setting do not carry over to non-Euclidean ones: Euclidean step-size policy 1/t would lead to a catastrophic drop from 1/ √ t to 1/ log t in the rate guarantees for (OMD).
Geometry-dependent rates. One final observation for Theorem 1 is that the convergence speed of (OMD) is measured in terms of the Bregman divergence D(x * , X t ). Since X t − x * = O( D(x * , X t )), this rate immediately translates to a rate on X t − x * . However, this rate is not tight, again because of the disparity between the Bregman and norm geometries discussed in Section 4. Specifically, in many cases, a non-zero Legendre exponent implies a sharper lower bound on the Bregman divergence of the form D(x * , X t ) = Ω( X t − x * 2−2β ). Thus, for β ∈ [0, 1/2), Jensen's inequality applied to the optimized bound (17) . Importantly, this rate guarantee is the same as in the Euclidean case (β = 0), though it requires a completely different step-size policy to achieve it. Determining when the Legendre exponent also provides a lower bound is a very interesting direction for future research, but one that lies beyond the scope of this work.

5.2.
Main ideas of the proof. The heavy lifting in the proof of Theorem 1 is provided by Proposition 1 below, which also makes explicit some of the hidden constants in the statement of the theorem. Proposition 1. With notation and assumptions as in Theorem 1, fix some r > 0 and let is bounded according to the following table and conditions:

Legendre exponent
Rate We outline here the main ideas of the proof, deferring all details to Appendix C.
Sketch of proof. Our proof strategy is to show that, under the stated conditions for r and the algorithm's initialization, the iterates of (OMD) remain within a neighborhood of x * where (SOS) holds with high probability. Then, conditioning on this event, we use the Bregman divergence as a stochastic potential function, and we derive the algorithm's rate of convergence using a series of lemmas on (random) numerical sequences. The step-by-step process is as follows: (1) The cornerstone of the proof is a descent inequality for the iterates of OMD, which relates D(x * , X t+1 ) to D(x * , X t ). As stated in Lemma B.2, we have 2 for all t ≥ 2, and φ 1 = 0.
(2) The next step is to show that, with high probability and if initialized close to x * , the iterates of OMD remain in U. This is a technical argument relying on a use of the Doob-Kolmogorov maximal inequality for submartingales (in the spirit of Hsieh et al. (2020)) and the condition η > 1/2 (which in turn guarantees that the step-size is square summable, so the corresponding submartingale error terms are bounded in L 1 ). The formal statement is Lemma B.4. (3) To proceed, choose r > 0 small enough so that the second order sufficiency condition (SOS) holds and the Legendre exponent estimate (11) both hold. More precisely, we take r small enough so that U = B(x * , r) ∩ K is included in the domains of validity of both properties (B and V respectively), and we will work on the event E t = {X s+1/2 ∈ U for all s = 1, 2, . . . , t}. Then, conditioning on E t , Assumption 2 gives and hence, by Lemma A.6, we get However, the descent inequality only involves D(x * , X t ) and not X t − x * 2 , so we need to relate the latter to the former. This is where the Legendre exponent comes into play. More precisely, Definition 2 with 1 + (4) Thus, going back to the descent inequality (22), employing the bound derived above, and taking expectations ultimately yields (after some technical calculations) the iterative bound where we used Assumption 3 to bound the expectations of U t+1/2 2 * and U t−1/2 2 * . (5) We are now in a position to reduce the problem under study to the behavior of the sequence a t = E (D(x * , X t ) + φ t )1 Et−1 , t ≥ 1. Indeed, taking a closer look at E t shows that . Consequently, the inequality (23) above can be rewritten as follows: (6) The behavior of this last sequence hinges heavily on whether α > 0 or not. In detail, we have: The long-run behavior of this recursive inequality is described by Lemmas A.8 and A.9: as long as µγ/β > 1 for η = 1 (and only for η = 1), we have a t = O(1/t η ).
• If α > 0, a series of further technical calculations in the same spirit yields The final step of our proof is to upper bound the behavior of a t based on this recursive inequality. The necessary groundwork is provided by Lemmas A.11 and A.12, and hinges on whether η ≥ 1+α 1+2α or not (and is also where the value of c(η, α) comes in).
Putting everything together, we obtain the conditional rates in the statement of the proposition.

Conclusion
In this work, we investigated the rate of convergence of optimistic mirror descent for solving variational inequalities with stochastic oracle feedback. Our results highlight the relationship between the rate of convergence of the last iterate and the local geometry of the distance-generating function at the solution. To capture this local geometry, we introduced the regularity exponent of the divergence relative to the norm that we dubbed the Legendre exponent. This quantity plays a central role in our results: the less regular the Bregman divergence around the solution, the slower the convergence of the last iterate. Furthermore, we show that that the step-size policy that guarantees the best rates depends on the distance-generating function through the associated Legendre exponent.
This work opens the door to various refinements and extensions diving deeper into the geometry of the method's DGF. A key remark is that the method's Legendre exponent seems to depend crucially on which constraints of the problem are active at a given solution. Deriving a precise characterization of the method's convergence rate in these cases is a fruitful direction for future research.

A. Auxiliary lemmas
We gather in this section the basic results used in our proofs.
A.1. Bregman divergences. We recall here classical lemmas on Bregman divergences; see Juditsky et al. (2011. We complement them with two elementary results Lemma A.6 on a specific bounding and Lemma A.7 which shows that the Legendre exponent is typically trivial on K h . Lemma A.3. For p ∈ K, x ∈ dom ∂h, y ∈ Y, y ′ ∈ N K (p) and x + = P x (y), Lemma A.4. For p ∈ K, x ∈ dom ∂h, y, y ′ ∈ Y and x + 1 = P x (y),x + 2 = P x (y ′ ), As a consequence, x ∈ K h , x + ∈ K, y ∈ Y, Moreover, x + = P x (y) implies that x + ∈ K h h.
Lemma A.6. If Assumption 2 holds, x + ∈ B, and x ∈ K, then, However, we are interested in the distance between x and x * , not in the distance between x + and x * . To remedy this, we use Young's inequality, which leads to, Combined with Eq. (A.1), we get, Lemma A.7 (Trivial Legendre exponent on the interior). Assume that ∇h is locally Lipschitz continuous on K h . Then, for any p ∈ K h , the Legendre exponent of h at p is β = 0.
Proof. Take p ∈ K h . As ∇h is assumed to be locally Lipschitz, there exists V a neighborhood of p and κ > 0 such that, for any x ∈ V ∩ K h , Now, as ∇h(p) ∈ ∂h(p), for any x ∈ V ∩ K h , Using the local Lipschitz continuity of ∇h, we get, for any x ∈ V ∩ K h , which reads β = 0.
A.2. Sequences. We recall results on sequences to transform descent-like inequalities into rates; see the classical textbook Polyak (1987). We also need variants of existing results (Lemma A.11 and Lemma A.12) that we state here.
Rearranging and replacing b by its expression gives, Therefore, with c(η, α) = 1 p 1 α 1 p − β > 0 and q 1 α q ′ ≤ c(η, α), we finally get that q ′ ≤ qb p − βb and, for t ≥ 1, Putting this inequality together with the one on the sequence (a t ) t≥1 gives, . Now, let us discuss separately the case when a t > b t and when a t ≤ b t .
More precisely, define T 0 = min{t ≥ 1 : a t ≤ b t } ∈ N * ∪ {+∞}, so that, for any Lemma A.10 with a t ← a t − b t gives, for 1 ≤ T < T 0 , so the statement holds in this case. We now handle the case of T ≥ T 0 . At t = T 0 , we have a T0 ≤ b T0 . We show, by induction, that, for any t ≥ T 0 , a t ≤ b t . The initialization at t = T 0 is trivial. So, now, assume that a t ≤ b t for some t ≥ T 0 . Recall that we have, so we only need to show that, Hence, both a t and b t belong to the interval on which ϕ t is non-decreasing so that ϕ t (a t ) ≤ ϕ t (b t ), which implies that a t+1 ≤ b t+1 . Therefore, we have shown that, for all B. Descent and stability of optimistic mirror descent B.1. Descent inequality. The lemmas we present here link the divergence between two consecutive iterates and the solution.
In particular, the first one only rely on i) the definition of one iteration of optimistic mirror descent; and ii) the core properties of the Bregman regularizer i.e., the 1-strong convexity of h (through Lemma A.1) and the non-expansivity of the proximal mapping (Lemma A.2). It does not require any assumption on the vector field v, on the noise in V , or additional properties of the Bregman divergence (such as a Legendre exponent).
Proof. First, apply Lemma A.4 with (x, y 1 , y 2 ) ← (X t , −γ t V t−1/2 , −γ t V t+1/2 ) and p ← x * , (B.1) In order to upper-bound and the definition of ∆ t+1 to bound where we used Young's inequality to get the last inequality line. Putting the two equations together, we get that For t ≥ 2, we use the non-expansivity of the proximal mapping (Lemma A.2) to bound where we used the definition of φ t and the assumption on the stepsizes. Combining this result with Eq. (B.1) gives the first assertion of the lemma for t ≥ 2. Now, if t = 1, Eq. (B.2) can be rewritten as, which yields the assertion of the lemma for t = 1 as 4γ 2 1 L 2 ≤ 1 − γ 1 µ. We now use the properties of V and v to refine the bound of Lemma B.1 into a descent inequality.
B.2. Stochastic local stability. In order to use the local properties of the Bregman divergence around the solution, we need a local stability result that gives us that sufficiently close iterates will remain in a neighborhood of the solution with high probability. To show such a result, we will use the following lemma Hsieh et al. (2020, Lem. F.1) for bounding a recursive stochastic process. Lemma B.3 (Hsieh et al. (2020, Lem. F.1)). Consider a filtration (F t ) t and four (F t ) tadapted processes (D t ) t , (ζ t ) t , (χ t ) t , (ξ t ) t such that χ t is non-negative and the following recursive inequality is satisfied for all t ≥ 1 For C > 0, we define the events A t by A 1 : We consider also the decreasing sequence of events I t defined by I t := 1≤s≤t A s . If the following three assumptions hold true Lemma B.4. Let Assumptions 1 and 2 hold. Consider the iterates of optimistic mirror descent with: i) an unbiased stochastic oracle V t = v(X t ) + U t with finite variance σ 2 (see (2) and (3)); ii) step-sizes 0 < γ t ≤ 1 4L with ∞ t=1 γ 2 t = Γ < +∞.

C. Convergence of optimistic mirror descent
Proof of Proposition 1.
(2) The first part of the result comes the stochastic stability lemma Lemma B.4 (which relies on Hsieh et al. (2020)), where we use that the stepsize choice is square-summable with +∞ t=1 (3) We now focus on proving the rates of convergence. Take r > 0 small enough so that U := B(x * , r) ∩ K is included in both V and B.
As a consequence, and using the assumption on the noise, The final step will be to study the sequence (a t ) defined by a t = E (D(x * , X t ) + φ t )1 Et−1 for all t ≥ 1. Indeed, a bound on a t = E (D(x * , X t ) + φ t )1 Et−1 implies a bound on the quantity of interest, E (D(x * , X t ) + φ t | E U ) .
Then, for η = 1, Lemma A.8 can be applied with the additional assumption that γ > κ µ to get that for any T ≥ 1, a T = O(1/T ). For 1 2 < η < 1, Lemma A.9 can be directly apply to obtain a T = O(1/T η ). • If α > 0, a little more work has to be done before on (C.4) before concluding.
We will now apply one of Lemmas A.11 and A.12 to the sequence (a t ). To make this step clear, let us introduce the same notations as these lemmas. Define q := µγ 2 α max(1, Φ α )κ and q ′ := 8σ 2 γ 2 .
-If η ≥ 1+α 1+2α , we apply Lemma A.11 to get that, for any T ≥ 1, α . Using asymptotic notations, this simply means that, The final bound on a T now comes from the fact that, if η = 1 .

D. Numerical experiments
We describe here the setup used in the illustration of Theorem 1 and provide additional plots.