Waïss Azizian

AND Panayotis Mertikopoulos

THE RATE OF CONVERGENCE OF BREGMAN PROXIMAL METHODS: LOCAL GEOMETRY VS. REGULARITY VS. SHARPNESS

Keywords: 2020 Mathematics Subject Classification. Primary 65K15, 90C33, secondary 68Q25, 68Q32 Legendre exponent, optimistic mirror descent, variational inequalities

We examine the last-iterate convergence rate of Bregman proximal methods -from mirror descent to mirror-prox and its optimistic variants -as a function of the local geometry induced by the prox-mapping defining the method. For generality, we focus on local solutions of constrained, non-monotone variational inequalities, and we show that the convergence rate of a given method depends sharply on its associated Legendre exponent, a notion that measures the growth rate of the underlying Bregman function (Euclidean, entropic, or other) near a solution.

In particular, we show that boundary solutions exhibit a stark separation of regimes between methods with a zero and non-zero Legendre exponent: the former converge at a linear rate, while the latter converge, in general, sublinearly. This dichotomy becomes even more pronounced in linearly constrained problems where methods with entropic regularization achieve a linear convergence rate along sharp directions, compared to convergence in a finite number of steps under Euclidean regularization.

Introduction

Bregman proximal methods (BPMs) have a long and rich history in optimization, going back at least to the introduction of the influential mirror descent algorithm by Nemirovski & Yudin [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF]. In plain terms, BPMs are first-order (constrained) optimization algorithms that forego Euclidean projections in favor of a more sophisticated "proxmapping" that minimizes a certain distance-like functional known as the Bregman divergence [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF][START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF][START_REF] Kiwiel | Proximal minimization methods with generalized Bregman functions[END_REF][START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF]. When this Bregman divergence is the Euclidean distance squared, one recovers the standard projection-based methods; other than that, depending on the problem's feasible region, different Bregman setups lead to a diverse collection of algorithms, from exponentially weighted gradient descent in the simplex [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF][START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF], to matrix multiplicative weights on the positive-semidefinite cone [START_REF] Kakade | Regularization techniques for learning with matrices[END_REF][START_REF] Tsuda | Matrix exponentiated gradient updates for on-line Bregman projection[END_REF], variants of Karmarkar's affine scaling algorithm for linear programs [START_REF] Vanderbei | A modification of Karmarkar's linear programming algorithm[END_REF], etc.

From an operational standpoint, one of the most appealing features of BPMs is that they achieve almost dimension-free convergence rates in problems with a convex structure and a favorable geometry -such as the L 1 ball, the spectraplex, second-order cones, etc. [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Bubeck | Convex optimization: Algorithms and complexity[END_REF][START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF]. This is owed to a delicate interplay between the algorithms' non-Euclidean update scheme and the global geometry of the problem's domain. However, these (almost) dimension-free guarantees also come with some strings attached: they do not concern the sequence of iterates generated by the method, but only its time average; as a result, the best guarantee that can be achieved after t iterations is O(1/t). In terms of oracle complexity, this is sufficient for problems that are not strongly convex / strongly monotone, but if one targets finer, geometric convergence rates, the lag induced by averaging cannot be compensated. And, on the other extreme, if the problem is not convex / monotone to begin with, iterate averaging does not provide any quantifiable benefits whatsoever; in this case, only the actual iterates generated by the method can be used as candidate solutions.

Our contributions. In this context, our paper seeks to provide a precise characterization of the last-iterate convergence rate of Bregman proximal methods, as a function of the Bregman divergence defining the method and the local geometry that it induces. To treat this question in as general a setting as possible, we focus throughout on variational inequality (VI) problems of the form Find x * ∈ X such that ⟨g(x * ), x -x * ⟩ ≥ 0 for all x ∈ X , (VI)

where X is a closed convex subset of a finite-dimensional normed space V, and g : X → V * is a (possibly non-monotone) single-valued operator on X with values in V * , the dual of V. This problem is a staple of many areas of mathematical programming, game theory and data science, as it provides a template for "optimization beyond minimization" -i.e., for problems where finding an optimal solution does not necessarily involve minimizing a loss function. In particular, in addition to standard minimization problems -which are recovered when g = ∇f for some smooth function f -the general formulation (VI) includes saddle-point problems, games, complementarity problems, etc.; for an introduction, see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] and references therein.

In this broad context, we examine the rate of convergence of a wide class of Bregman proximal methods to local solutions of (VI) that satisfy a second-order sufficient condition. Specifically, the class of algorithms we consider includes as special cases (i) the original mirror descent (MD) algorithm of [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF]; (ii) the mirror-prox (MP) method of Nemirovski [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] -which has the same update structure as the Bregman-based algorithm of [START_REF] Auslender | Interior projection-like methods for monotone variational inequalities[END_REF] and contains as a special case the extra-gradient (EG) algorithm of [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF]; (iii) the so-called optimistic mirror descent (OMD) method of [START_REF] Rakhlin | Optimization, learning, and games with predictable sequences[END_REF] -itself a Bregman analogue of the modified Arrow-Hurwicz algorithm of [START_REF] Popov | A modification of the Arrow-Hurwicz method for search of saddle points[END_REF]; etc.

Our first finding is a crisp characterization of last-iterate convergence rate of BPMs in terms of the local geometry induced by the underlying Bregman function near a given solution of (VI). We make this dependence precise via the notion of the Legendre exponent, a regularity measure for Bregman methods due to [START_REF] Azizian | The last-iterate convergence rate of optimistic mirror descent in stochastic variational inequalities[END_REF], which can roughly be described as the logarithmic ratio of the volume of a Euclidean ball to that of a Bregman ball of the same radius. For example, Euclidean methods have a Legendre exponent of β = 0 and they converge at a linear rate; entropic methods have a Legendre exponent of β = 1/2 at boundary points, and they converge at a rate of O(t -1 ); more generally, as we show in Theorem 1, methods with a Legendre exponent β > 0 converge at a rate of O(t 1-1/β ). The Euclidean regime (β = 0) is perfectly aligned with existing results for the geometric last-iterate convergence rate of the EG algorithm and its variants [START_REF] Gidel | A variational inequality perspective on generative adversarial networks[END_REF][START_REF] Hsieh | On the convergence of single-call stochastic extra-gradient methods[END_REF][START_REF] Malitsky | Projected reflected gradient methods for monotone variational inequalities[END_REF][START_REF] Mokhtari | A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: proximal point approach[END_REF]. By contrast, the Legendre regime (β > 0) indicates a significant drop in the algorithm's last-iterate convergence speed, even though ergodic convergence rates [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] and results for bilinear games [START_REF] Wei | Linear last-iterate convergence in constrained saddle-point optimization[END_REF] might suggest otherwise.

Subsequently, motivated by applications to game theory and linear programming, we take a closer look at the convergence rate of BPMs across the constraints that are active at a solution x * of (VI) depending on the position of g(x * ) relative to said constraints. This analysis reveals that Bregman proximal methods have a particularly fine structure: along sharp directions (i.e., constraints along which g(x * ) is strictly inward-pointing), BPMs converge (i) at a rate of O(1/t 1/(2β-1) ) if 1/2 < β < 1; (ii) at a geometric rate if 0 < β ≤ 1/2 (e.g., for entropic methods); and (iii ) in a finite number of iterations if β = 0 (cf. Theorem 2). Thus, even though the estimates of Theorem 1 are in general tight, the actual rate of convergence of a Bregman method along different coordinates / constraints could be starkly different -and, in fact, dramatically faster if the solution under study is itself sharp.

The closest antecedent of our work is the conference paper [START_REF] Azizian | The last-iterate convergence rate of optimistic mirror descent in stochastic variational inequalities[END_REF] where the Legendre exponent was introduced to analyze the convergence of OMD in stochastic VI problems (without considering sharp directions and/or faster identification rates). The stochastic and deterministic settings are obviously very different, both in the challenges involved as well as the rates obtained, so there is no overlap in our analysis and results. Other than that, we are not aware of any comparable results in the literature concerning the radically different convergence landscape of BPMs along active and inactive constraints.

Problem setup and preliminaries

To fix notation and terminology, in the rest of our paper V will denote an n-dimensional real space with norm ∥•∥ and X will be a closed convex subset thereof. We will also write Y := V * for the dual of V, ⟨y, x⟩ for the canonical pairing between y ∈ Y and x ∈ V, and ∥y∥ * := max{⟨y, x⟩ : ∥x∥ ≤ 1} for the induced dual norm on Y.

2.1. Blanket assumptions. Throughout the sequel, we will make the following assumptions for the defining vector field g : X → Y of (VI) and the solution x * ∈ X under study: Assumption 1 (Lipschitz continuity). The vector field g is L-Lipschitz continuous, i.e.,

∥g(x ′ ) -g(x)∥ * ≤ L∥x ′ -x∥ for all x, x ′ ∈ X . (LC)
Assumption 2 (Second-order sufficiency). There exists a convex neighborhood B of x * in X and a positive constant µ > 0 such that

⟨g(x) -g(x * ), x -x * ⟩ ≥ µ∥x -x * ∥ 2 for all x ∈ B. (SOS)
In general, Assumption 2 guarantees that x * is the unique solution of (VI) in B; we illustrate this in two special cases of interest:

(1) Minimization problems: suppose that g = ∇f for some Lipschitz smooth objective function f on X . Then, Assumption 2 implies that f grows (at least) quadratically along every ray emanating from x * , i.e., f

(x) -f (x * ) ≥ ⟨∇f (x * ), x -x * ⟩ + (µ/2)∥x - x * ∥ 2 = Ω(∥x -x * ∥ 2
) for all x ∈ B, i.e., x * is an isolated local minimizer of f . (2) Min-max problems: suppose that X factorizes as X = X 1 × X 2 for suitable factor sets X 1 , and X 2 , let L : X → R be a smooth function on X , and write g = (∇ x1 L, -∇ x2 L) for the min-max gradient of L (with respect to x 1 ∈ X 1 and x 2 ∈ X 2 respectively). Then, any solution x * = (x * 1 , x * 2 ) of (VI) that satisfies Assumption 2 enjoys the local growth bounds

L(x 1 , x * 2 ) -L(x * 1 , x * 2 ) = Ω(∥x 1 -x * 1 ∥ 2 ) and L(x * 1 , x * 2 ) -L(x * 1 , x 2 ) = Ω(∥x 2 -x * 2 ∥ 2 ), implying in turn that x * is an isolated, hyperbolic saddle-point of L.
More examples satisfying (SOS) include strict Nash equilibria in finite games [START_REF] Fudenberg | Game Theory[END_REF], deterministic Nash policies in (generic) stochastic games [START_REF] Solan | Stochastic games[END_REF], etc. Overall, Assumptions 1 and 2 apply to a very wide range of problems, so we will treat them as blanket assumptions throughout.

2.2. Bregman proximal methods. As we discussed in the introduction, the main algorithmic template that we will examine for solving (VI) is a general class of first-order algorithms known as Bregman proximal methods (BPMs). The defining ingredient of this class is the notion of Bregman regularizer, which is defined below as follows:

Definition 1 (Bregman regularizers and related notions). A proper, lower semi-continuous, strictly convex function

h : V → R ∪ {∞} is a Bregman regularizer on X if (1) h is supported on X , i.e., dom h = X .
(2) 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.

(3) h is 1-strongly convex relative to ∥•∥, i.e., for all x ∈ dom ∂h, x ′ ∈ dom h, we have

h(x ′ ) ≥ h(x) + ⟨∇h(x), x ′ -x⟩ + 1 2 ∥x ′ -x∥ 2 . (2.1)
The set X h := dom ∂h will be referred to as the prox-domain of h. In addition, we also define the Bregman divergence of h as

D(p, x) = h(p) -h(x) -⟨∇h(x), p -x⟩ for all x ∈ X h , p ∈ X (2.2)
and the induced prox-mapping as

P x (y) = arg min x ′ ∈X {⟨y, x -x ′ ⟩ + D(x ′ , x)} for all x ∈ X h , y ∈ Y. (2.3)
Examples of Bregman regularizers are given in Section 3, where we also take an in-depth look at their properties. For now, given a Bregman regularizer on X , the general class of Bregman proximal methods (BPMs) that we will consider is defined via the generic recursion

x t+1/2 = P xt (-γ t g t ) x t+1 = P xt (-γ t g t+1/2 ) (BPM)
where (i ) t = 1, 2, . . . denotes the method's iteration counter; (ii) γ t > 0 is a (non-increasing) step-size sequence; (iii ) g t and g t+1/2 are sequences of "oracle signals" that we discuss in detail below. In terms of vocabulary, the iterates x t , t = 1, 2, . . . , will be referred to as the "base states" of the method, while the "half-iterates" x t+1/2 , t = 1, 2, . . . , will be referred to as the method's "leading states". Finally, in terms of initialization, we will take for convenience

x 1 = x 1/2
g t = ag(x t ) + bg(x t-1/2 ) (2.5) for some a, b ∈ [0, 1] with a + b ≤ 1 and a + b = 1 if b > 0. 1
For concreteness, we illustrate below three archetypal Bregman methods that serve as the backbone of the above framework:

(1) Mirror descent: following [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF][START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF], the mirror descent algorithm proceeds recursively as x + = P x (-γg(x)), so it can be recovered from (BPM) by taking a = 0, b = 0 or, equivalently g t = 0 for all t = 1, 2, . . . (MD)

(2) Mirror-prox: following [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF], the mirror-prox algorithm corresponds to the choice a = 1, b = 0 or, equivalently g t = g(x t ) for all t = 1, 2, . . . (MP)

(3) Optimistic mirror descent: originally due to [START_REF] Popov | A modification of the Arrow-Hurwicz method for search of saddle points[END_REF] (in the Euclidean case) and [START_REF] Chiang | Online optimization with gradual variations[END_REF][START_REF] Rakhlin | Optimization, learning, and games with predictable sequences[END_REF] (for the general case), the optimistic mirror descent algorithm is obtained by setting a = 0, b = 1 or, equivalently g t = g(x t-1/2 ) for all t = 1, 2, . . . (OMD)

These three algorithms are the most widely studied Bregman methods in the literature, so we will use them as running examples throughout.

Motivating examples

We now proceed to take a closer look at some commonly used Bregman regularizers (and the induced prox-mappings) with the goal of determining the rate of convergence of the associated Bregman method. For concreteness, we focus on one-dimensional problems where X is a closed interval of R and g is the affine vector field

g(x) = x -x * , x ∈ R, (3.1)
for different choices of x * ∈ R. To streamline our presentation, we will only examine the mirror descent recursion (MD) with constant step-size schedules γ t ≡ γ for some γ > 0.

In this case, we obtain the general recursive scheme

x t+1 = F (x t ) with F (x) = P x (-γg(x)), (3.2)
and we will examine the convergence speed of x t by analyzing the behavior of F near x * .

Example 3.1 (Euclidean regularization). We begin with the quadratic regularizer h(x) = x 2 /2 for x ∈ X . Concretely, taking X = [0, ∞) and noting that h ′ (x) = x, we have:

a) Prox-domain: X h = X b) Bregman divergence: D(p, x) = (p -x) 2 /2 c) Prox-mapping: P x (y) = [x + y] + (3.3)
Consider now the case x * = 0, i.e., g(x) = x. Then, for γ ∈ (0, 1), the update (3.2) becomes F (x) = x -γx = (1 -γ)x (3.4) i.e., F is contracting. We thus conclude that x t converges to x * = 0 at a geometric rate, viz.

D(x * , x t ) = 1 2 x 2 t = Θ (1 -γ) 2t or, in absolute value, |x t -x * | = Θ((1 -γ) t ). ♦ 1 Note that the requirement "a + b = 1 if b > 0"
is only intended to ease notation and does not lead to a loss in generality: if b > 0, we can always rescale γt by a + b so the condition a + b = 1 is satisfied automatically.

Example 3.2 (Entropic regularization). Another popular choice when X = [0, ∞) is the entropic regularizer h(x) = x log x [START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF][START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF]. In this case, h ′ (x) = 1 + log x, which gives a) Prox-domain:

X h = ri X = (0, ∞) b) Bregman divergence: D(p, x) = p log(p/x) + x -p c) Prox-mapping: P x (y) = x exp(y). (3.5)
Now, taking g(x) = x as in the previous example, the update rule (3.2) becomes

F (x) = x exp(-γx) = x(1 -γx + o(x)) = x -γx 2 + o(x 2 ) as x → 0. (3.6)
In contrast to (3.4), we now have F (x) ∼ x instead of (1 -γ)x, so F is no longer a contraction. Instead, the iterates of (3.6) may be analyzed by means of the following lemma:

Lemma 1. Suppose that f : R + → R + admits the asymptotic expansion

f (x) = x -λx 1+r + o(x 1+r ) as x → 0 (3.7)
for positive constants λ, r > 0. Then, for u 1 > 0 small enough, the sequence u t+1 = f (u t ), t = 1, 2, . . . , converges to 0 at a rate of u t ∼ (λrt) -1/r .

Thanks to this lemma (which we prove in Appendix A), we readily conclude that x t converges to 0 at a rate of D(x * , x t ) = x t = |x t -x * | ∼ 1/(γt). ♦ Example 3.3 (Fractional power). Take X = [0, ∞) and g(x) = x as in Examples 3.1 and 3.2 above. Then, for a given q > 0, q ̸ = 1, the fractional power regularizer -or Tsallis entropy -on X is defined as h(x) = [q(1 -q)] -1 (x -x q ) [1, [START_REF] Mertikopoulos | Learning in games via reinforcement and regularization[END_REF][START_REF] Tsallis | Possible generalization of Boltzmann-Gibbs statistics[END_REF]. For this choice of regularizer, we have h ′ (x) = (1 -qx q-1 )/[q(1 -q)], and a series of direct calculations gives: 2 a) Prox-domain:

X h = (0, ∞) if q ∈ (0, 1) and X h = [0, ∞) if q > 1 b) Bregman divergence: D(p, x) = x q -p q q(1 -q) -x q-1 x -p 1 -q c) Prox-mapping: P x (y) = x q-1 -(1 -q)y 1 q-1
for q ∈ (0, 1).

(3.8)

Now, when applied to g(x) = x, the fractional power variant of (3.2) for q ∈ (0, 1) gives

F (x) = x [1 + γ(1 -q)x 2-q ] 1/(q-1) = x -γx 3-q + o(x 3-q ) as x → 0. (3.9)
Hence, by Lemma 1, we conclude that x t converges to 0 at a rate of

D(x * , x t ) = Θ t -q/(2-q) or, in absolute value, |x t -x * | = Θ t -1/(2-q) .
♦ Example 3.4 (Hellinger distance). Our last example concerns the Hellinger regularizer

h(x) = - √ 1 -x 2 on X = [-1, 1]. Since h ′ (x) = x/ √ 1 -x 2 ,
we readily obtain the 2 Strictly speaking, the expression we provide for Px(y) is only valid when y < x q-1 /(1 -q). The reason for this is that this regularizer is not strongly convex over [0, ∞), so the prox-mapping Px(y) is not well-defined for all values of y. This detail is not important in the calculations that follow, so we disregard it for now. Iterations (t)

|x t -x * | Entropy theory: Θ(t -1 )
Fractional Power (q = 3/4) theory: Θ(t -4/5 ) Hellinger (shifted) theory: following:

Θ(t -2/3 )
a) Prox-domain: X h = ri X = (-1, 1) b) Bregman divergence: D(p, x) = 1 -px -(1 -p 2 )(1 -x 2 ) √ 1 -x 2 c) Prox-mapping: P x (y) = x + y √ 1 -x 2 1 -x 2 + (x + y √ 1 -x 2 ) 2 . ( 3.10)
In this case, taking g(x) = x as per the previous examples, yields

F (x) = x -γx √ 1 -x 2 1 -x 2 + (x -γx √ 1 -x 2 ) 2 ∼ x -γx as x → 0, (3.11)
i.e., x t converges to x * = 0 at a geometric rate, as in Example 3.1. On the other hand, if we consider the shifted operator g(x) = x + 1, a somewhat tedious calculation (which we detail in Appendix B) gives the following Taylor expansion near x * = -1:

F (x) = x * + (x -x * ) -2 √ 2γ(x -x * ) 5/2 + o (x -x * ) 5/2 . (3.12)
Hence, by Lemma 1, we conclude that x t converges to x * = -1 at a rate of

D(x * , x t ) = Θ(t -1/3 ) or, in absolute value, |x t -x * | = Θ(t -2/3
).

♦

Albeit one-dimensional, the above examples provide a representative view of the geometry of Bregman proximal methods near a solution. Specifically, they show that the divergence induced by a given regularizer may exhibit a very different behavior at the boundary of X : when x * is a boundary point, D(x * , x) grows as Θ(∥x -x * ∥ 2 ) in the Euclidean case, as Θ(∥x -x * ∥) for the negative entropy, or, more generally, as Θ(∥x -x * ∥ q ) for the q-th power regularizer. As a result, when used as a measure of convergence, it is important to rescale the Bregman divergence in order to avoid inflating -or deflating -an algorithm's rate of convergence.

Nonetheless, even if we take this rescaling into account, different instances of (MD) may lead to completely different rates of convergence. Specifically, in terms of absolute values, we observe a geometric rate in the Euclidean and shifted Hellinger cases, a rate of Θ(1/t) for the negative entropy, and a rate of Θ(1/t 1/(2-q) ) for the q-th power regularizer (cf. Fig. 1 above). This is due to the different first-order behavior of the iterative update map x ← F (x) that underlies (MD), which is itself intimately related to the growth rate of the Bregman divergence near a solution x * of (VI). We will make this relation precise in the next section.

The Legendre exponent and convergence rate analysis

Our goal in this section is to provide a precise link between the geometry induced by a Bregman regularizer near a solution and the convergence rate of the associated Bregman proximal method. The key notion in this regard is that of the Legendre exponent, which we define and discuss in detail below.

4.1. The Legendre exponent. Our starting point is the observation that the strong convexity requirement for h can be equivalently expressed as

D(p, x) ≥ 1 2 ∥p -x∥ 2 for all p ∈ X , x ∈ X h . (4.1)
Qualitatively, this means that the convergence topology induced by the Bregman divergence of h on X is at least as fine as the ambient norm topology: if a sequence x t ∈ X h , t = 1, 2, . . . , converges to p ∈ X in the Bregman sense (D(p, x t ) → 0), then it also converges in the ambient norm topology (∥x t -p∥ → 0). On the other hand, from a quantitative standpoint, the rate of this convergence could be quite different: as we saw in the previous section, the reverse inequality D(p, x) = O(∥p -x∥ 2 ) may fail to hold, in which case D(p, x t ) and ∥x -x t ∥ would exhibit a different asymptotic behavior.

To quantify this gap, we use the notion of the Legendre exponent introduced in [START_REF] Azizian | The last-iterate convergence rate of optimistic mirror descent in stochastic variational inequalities[END_REF].

Definition 2. Let h be a Bregman regularizer on X . Then the Legendre exponent of h at p ∈ X is defined as

β h (p) := inf β ∈ [0, 1] : lim sup x→p D(p, x) ∥x -p∥ 1-β < ∞ (4.2)
and we say that h is tight at p if the infimum is attained in (4.2), i.e., if

β h (p) is the minimal β ∈ [0, 1] such that D(p, x) = O ∥p -x∥ 2(1-β) for x near p. (4.3)
Informally, the Legendre exponent measures the deficit in relative size between ordinary "norm neighborhoods" in X and the corresponding "Bregman neighborhoods" induced by the sublevel sets of the Bregman divergence. Specifically, (i) the case β h (p) = 0 corresponds to the "norm-like" behavior D(p, x) = Θ(∥p -x∥ 2 ); (ii) any other value β h (p) ∈ (0, 1) indicates a different limiting behavior for D(p, x) as x → p; and, finally, (iii) when β h (p) = 1 we may have lim sup x→p D(p, x) > 0. In this last case, the ambient norm topology is strictly coarser than the Bregman topology in the sense that D(p, x t ) may remain bounded away from zero even if x t → p as t → ∞; for a detailed discussion of this phenomenon, cf. [START_REF] Azizian | The last-iterate convergence rate of optimistic mirror descent in stochastic variational inequalities[END_REF]Sec. 4] and [START_REF] Pauwels | On the nature of bregman functions[END_REF].

For illustration purposes, we compute below the Legendre exponent for each of the running examples of Section 3 (see also Table 1):

(1) Quadratic regularization (Example 3.1): Since D(p, x) = (p -x) 2 /2 for all p, x ∈ X , we have β h (p) = 0 for all p ∈ X .

Domain (X ) Regularizer (h) Legendre Exponent (β h (p)) Convergence Rate Euclidean arbitrary x 2 /2 0 Linear Entropic [0, ∞) x log x 1/2 O(1/t) Tsallis [0, ∞) [q(1 -q)] -1 (x -x q ) max{0, 1 -q/2} O(1/t q/(2-q) ) Hellinger [-1, 1] - √ 1 -x 2 3/4 O(1/t 1/3 )
Table 1: Summary of the Legendre exponents for the 1-dimensional examples of Section 2 at a boundary point p of X ⊂ R, and the associated convergence rates.

(2) Negative entropy (Example 3.2): For p = 0, Eq. (3.5) gives D(0, x) = x, so

β h (0) = 1/2.
Otherwise, for all p ∈ X h = (0, ∞), a Taylor expansion with Lagrange remainder yields D(p, x) = O((p -x) 2 ), so β h (p) = 0 for all p ∈ (0, ∞). (3) Tsallis entropy (Example 3.3): For p = 0, Eq. (3.8) gives D(0, x) = x q /q, so β h (0) = max{0, 1 -q/2}. Otherwise, for all p ∈ X h = (0, ∞), a Taylor expansion yields D(p, x) = O((p -x) 2 ), so β h (p) = 0 in this case. (4) Hellinger regularizer (Example 3.4): For p = ±1, Eq. (3.10) readily gives

D(±1, x) = (1 ∓ x)/(1 ± x) = Θ(|x ∓ 1| 1/2 ), so β h (±1) = 1 -1/4 = 3/4. Instead, if p ∈ (-1, 1) a Taylor expansion again yields D(p, x) = O((p -x) 2 ), so β h (p) = 0 in this case.
A common pattern that emerges above is that β h (p) = 0 whenever p is an interior point. We make this observation precise in Appendix A (Lemma A.5), where we show more generally that β h (p) = 0 whenever ∇h is (locally) Lipschitz continuous in a neighborhood of p in X .

4.2.

Convergence rate analysis. We are now in a position to state our first general result for the convergence rate of (BPM). To do so, we will make the blanket assumption that h is tight at x * with Legendre exponent β * := β h (x * ). In particular, this means that there exists a neighborhood U of x * in X and a positive constant κ > 0 such that

D(x * , x) ≤ κ 2 ∥x -x * ∥ 2(1-β * ) for all x ∈ U. (4.4)
We then have the following result.

Theorem 1. Suppose that Assumptions 1-3 hold and (BPM) is run with a constant step-size γ t ≡ γ, t = 1, 2, . . . , such that

γ ≤ 1 2φL and γ(1 -a -b) 2 ≤ µ 8L 2 (4.5)
where φ = ( √ 5 + 1)/2 is the golden ratio. If x 1 is initialized sufficiently close to x * , then the iterates x t of (BPM) enjoys the bound:

D(x * , x t ) ≤ D(x * , x 1 ) •    1 -µγ 2κ t-1 if β * = 0, [1 + ρµγ(t -1)] 1-1/β * if β * ∈ (0, 1), (4.6) where ρ = β * 1-β * max 2κ 1 1-β * D(x * , x 1 ) -β * 1-β * , 2 β * 1-β * -1 .
Before moving on to the proof of Theorem 1, some remarks and corollaries are in order (see also Table 1 for an explicit illustration of the derived rates for Examples 3.1-3.4): Remark 1. A first point of note is the sharp drop in the convergence rate of (BPM) from geometric, when β * = 0, to a power law when β * > 0. As we saw in Section 3, this drop is unavoidable, even when X is 1-dimensional and g is affine; in fact, the calculations of Section 3 show that the guarantees of Theorem 1 are, in general, unimprovable. ♦ Remark 2. We note also that the guarantees of Theorem 1 concern the Bregman divergence, not the ambient norm. Since D(x * , x t ) = Ω(∥x t -x * ∥ 2 ), these bounds can be restated in terms of ∥x t -x * ∥, but this conversion is not without loss of information: if the bound D(x * , x t ) = Ω(∥x t -x * ∥ 2 ) is not tight, the actual rate in terms of the norm may be significantly different. This phenomenon was already observed in the 1-dimensional examples of Section 3 where

D(x * , x t ) = Θ(∥x t -x * ∥ 2(1-β * ) ), in which case Theorem 1 gives ∥x t -x * ∥ = O t -1/(2β * ) (4.7
) whenever β * > 0 (see also Table 1). In general however, the Bregman divergence may grow at different rates along different rays emanating from x * , so it is not always possible to translate a Bregman-based bound to a norm-based bound (or vice versa). This analysis requires a much closer look at the geometric structure of X , depending on which constraints are active at x * ; we examine this issue at depth in Section 5. ♦ Remark 3. Even though Theorem 1 is stated for a constant step-size, our proof allows for a variable step-size γ t , provided that the step-size conditions (4.5) are satisfied. In this case, the bounds (4.6) becomes

D(x * , x t ) ≤ D(x * , x 1 ) t-1 s=1 1 -µγs 2κ for the Euclidean case (β * = 0), 1 + ρµ t-1 s=1 γ s 1-1/β * for the Legendre case (0 < β * < 1).
The sharpness of these guarantees can be checked as before. ♦ 4.3. Proof of Theorem 1. We now proceed to the proof of Theorem 1, beginning with a series of intermediate results tailored to the update structure of (BPM). The first of these lemmas relates Bregman divergences before and after a prox-step modulo an element of the polar cone PC(p) := {v ∈ Y : ⟨v, x -p⟩ ≤ 0 for all x ∈ X } of X at the reference point p. Lemma 2. Let x + = P x (y) for x ∈ X h , y ∈ Y. Then, for all p ∈ X , v ∈ PC(p), we have:

D(p, x + ) ≤ D(p, x) + ⟨y -v, x + -p⟩ -D(x + , x) (4.8a) ≤ D(p, x) + ⟨y -v, x -p⟩ + 1 2 ∥y -v∥ 2 * . (4.8b)
The next lemma extends Lemma 2 to emulate the two-step structure of (BPM):

Lemma 3. Let x + i = P x (y i ) for some x ∈ X h and y i ∈ Y, i = 1, 2.
Then, for all p ∈ X and all v ∈ PC(p), we have:

D(p, x + 2 ) ≤ D(p, x) + ⟨y 2 -v, x + 1 -p⟩ + 1 2 ∥y 2 -y 1 -v∥ 2 * -1 2 ∥x + 1 -x∥ 2 . (4.9)
Versions of the above inequalities already exist in the literature, see e.g., [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF]. The novelty in Lemmas 2 and 3 is the extra term involving the polar vector v ∈ PC(p); this term plays an important role in the sequel, so we provide complete proofs in Appendix A.

With these preliminaries in hand, we proceed to derive two further inequalities that play a pivotal role in the analysis of (BPM). The first is an immediate corollary of Lemma 3: Corollary 1. Let x * be a solution of (VI). Then, for all c ≥ 0 and all t = 1, 2, . . . , the iterates of (BPM) satisfy the template inequality

D(x * , x t+1 ) ≤ D(x * , x t ) -γ t ⟨g t+1/2 -cg(x * ), x t+1/2 -x * ⟩ + 1 2 γ 2 t ∥g t+1/2 -g t -cg(x * )∥ 2 * -1 2 ∥x t+1/2 -x t ∥ 2 . (4.10)
Proof. Since x * is a solution of (VI), we have g(x * ) ∈ -PC(x * ). Eq. (4.10) then follows from Lemma 3 with

x ← x t , p ← x * , v ← -cγ t g(x * ) ∈ PC(x * ) and (y 1 , y 2 ) ← (-γ t g t , -γ t g t+1/2
).

■

The second inequality that we derive provides an "energy function" for (BPM), namely

E t = D t + f t (4.11)
where

D t = D(x * , x t ) and f t = γ 2 t-1 ∥(a + b)g t-1/2 -g t-1 ∥ 2 *
(by convention, we take f 1 = 0). The lemma below outlines the Lyapunov properties of E t .

Proposition 1. Suppose that Assumptions 1 and 3 hold and (BPM) is run with a step-size such that

λγ t + 4γ 2 t L 2 ≤ 1 for some λ ≥ 0 and all t = 1, 2, . . . (4.12)
Then, for all t = 1, 2, . . . , the iterates x t of (BPM) satisfy the inequality

E t+1 ≤ E t -λγ t f t -γ t ⟨g(x t+1/2 ) -g(x * ), x t+1/2 -x * ⟩ -γ t (a + b)⟨g(x * ), x t+1/2 -x * ⟩ -1 2 ∥x t+1/2 -x t ∥ 2 + γ 2 t (1 -a -b) 2 L 2 ∥x t+1/2 -x * ∥ 2 + 2γ 2 t (a + b) 2 L 2 ∥x t+1/2 -x t ∥ 2 . ( 4
D t+1 ≤ D t -γ t ⟨g t+1/2 -g(x * ), x t+1/2 -x * ⟩ -γ t (a + b)⟨g(x * ), x t+1/2 -x * ⟩ -1 2 ∥x t+1/2 -x t ∥ 2 + γ 2 t 2 ∥g t+1/2 -g t -cg(x * )∥ 2 * . (4.14)
Since g t+1/2 = (a + b)g t+1/2 + cg t+1/2 , the last term above may be bounded as

1 2 γ 2 t ∥g t+1/2 -g t -cg(x * )∥ 2 * ≤ γ 2 t c 2 ∥g t+1/2 -g(x * )∥ 2 * + γ 2 t ∥(a + b)g t+1/2 -g t ∥ 2 * ≤ γ 2 t c 2 L 2 ∥x t+1/2 -x * ∥ 2 * + γ 2 t ∥(a + b)g t+1/2 -g t ∥ 2 * = γ 2 t (1 -a -b) 2 L 2 ∥x t+1/2 -x * ∥ 2 + f t+1 , (4.15)
where we used Assumption 1 in the second line and the definition (4.11) of f t in the last one. Thus, combining Eqs. (4.14) and (4.15) and comparing to (4.13), it suffices to show that

2f t+1 ≤ (1 -λγ t )f t + 4γ 2 t (a + b) 2 L 2 ∥x t+1/2 -x t ∥ 2 for all t = 1, 2, . . . (4.16)
We consider two distinct cases for this below.

Case 1: t = 1. By the definition (4.11) of f t and Eqs. (2.4) and (2.5), we have:

f 2 = γ 2 1 ∥(a + b)g 3/2 -g 1 ∥ 2 * = γ 2 1 (a + b) 2 ∥g(x 3/2 ) -g(x 1 )∥ 2 * ≤ γ 2 1 (a + b) 2 L 2 ∥x 3/2 -x 1 ∥ 2 , (4.17)
where we used the initialization assumption x 1 = x 1/2 in the second equality and the Lipschitz continuity of g in the last one. Since f 1 = 0 by construction, our claim is immediate.

Case 2: t > 1. By Young's inequality and the Lipschitz continuity of g, we readily obtain

f t+1 = γ 2 t ∥(a + b)g t+1/2 -g t ∥ 2 * = γ 2 t (a + b)[g(x t+1/2 ) -g(x t )] + b[g(x t ) -g(x t-1/2 )] 2 * ≤ 2γ 2 t (a + b) 2 ∥g(x t+1/2 ) -g(x t )∥ 2 * + 2γ 2 t b 2 ∥g(x t ) -g(x t-1/2 )∥ 2 * ≤ 2γ 2 t (a + b) 2 L 2 ∥x t+1/2 -x t ∥ 2 + 2γ 2 t b 2 L 2 ∥x t -x t-1/2 ∥ 2 ≤ 2γ 2 t (a + b) 2 L 2 ∥x t+1/2 -x t ∥ 2 + 2γ 2 t b 2 L 2 γ 2 t-1 ∥g t-1/2 -g t-1 ∥ 2 *
where, in the last line, we used Lemma A.4 to bound the difference x t -x t-1/2 as

∥x t -x t-1/2 ∥ = ∥P xt-1 (-γ t-1 g t-1/2 )-P xt-1 (-γ t-1 g t-1 )∥ ≤ γ t-1 ∥g t-1/2 -g t-1 ∥ * . (4.18)
Finally, by Assumption 3, we have c = 0

whenever b > 0, so b 2 γ 2 t-1 ∥g t-1/2 -g t-1 ∥ 2 * = b 2 γ 2 t-1 ∥(1 -c)g t-1/2 -g t-1 ∥ 2 * = b 2 f t for all t > 1.
Hence, putting everything together, we get

f t+1 ≤ 2γ 2 t (a + b) 2 L 2 ∥x t+1/2 -x t ∥ 2 + 2γ 2 t b 2 L 2 f t . (4.19)
Eq. (4.16) then follows by the requirement (4.12), which implies that 2γ 2 t L 2 ≤ (1λγ t )/2.

■

Moving forward, since x * is a solution of (VI), the first line of (4.13) yields a negative O(γ t ) contribution to E t , whereas the third and fourth lines collectively represent a subleading O(γ 2 t ) "error term". This decomposition would suffice for the analysis of (BPM) if the coupling term ⟨g(x * ), x t+1/2 -x * ⟩ did not incur an additional O(γ t ) positive contribution to E t+1 . This error term is difficult to control but if x * satisfies (SOS), we have the following bound. Lemma 4. Suppose that Assumption 2 holds. Then, for all x ∈ X , x ′ ∈ B and all c ∈ [0, 1], we have:

⟨g(x ′ ) -cg(x * ), x ′ -x * ⟩ ≥ 1 2 µ∥x -x * ∥ 2 -µ∥x ′ -x∥ 2 . (4.20)
Proof. Since x * is a solution of (VI) and c ∈ [0, 1], we have (1 -c)⟨g(x * ), x ′ -x * ⟩ ≥ 0 for all x ′ ∈ X . Hence, by Assumption 2, we get

⟨g(x ′ ) -cg(x * ), x ′ -x * ⟩ ≥ ⟨g(x ′ ) -g(x * ), x ′ -x * ⟩ ≥ µ∥x ′ -x * ∥ 2 (4.21)
and our assertion follows from the basic bound ∥x -

x * ∥ 2 ≤ 2∥x -x ′ ∥ 2 + 2∥x ′ -x * ∥ 2 . ■
With this ancillary estimate in hand, we may finally sharpen Proposition 1 to obtain a bona fide energy inequality for solutions satisfying (SOS): Proposition 2. Suppose that Assumptions 1-3 hold and (BPM) is run with γ t such that

2µγ t + 4γ 2 t L 2 ≤ 1 and (1 -a -b) 2 γ t ≤ µ 8L 2 for all t = 1, 2, . . . (4.22)
Then, for all t ≥ 1 such that x t+1/2 ∈ B, we have

E t+1 ≤ E t -µγ t f t -1 4 µγ t ∥x t -x * ∥ 2 . (4.23)
Proof. Assume that x t+1/2 ∈ B and set c = 1 -a -b. Then, invoking Lemma 4 with x ← x t and x ′ ← x t+1/2 , we get

⟨g(x t+1/2 )-g(x * ), x t+1/2 -x * ⟩+(a+b)⟨g(x * ), x t+1/2 -x * ⟩ ≥ 1 2 µ∥x t -x * ∥ 2 -µ∥x t+1/2 -x t ∥ 2 .
(4.24) Thus, taking λ ← µ in Proposition 1 (in terms of step-size conditions, the first part of (4.22) implies (4.12)) and combining with the above, the bound (4.13) becomes

E t+1 ≤ E t -µγ t f t -1 2 µγ t ∥x t -x * ∥ 2 + γ 2 t c 2 L 2 ∥x t+1/2 -x * ∥ 2 -1 2 1 -4γ 2 t (a + b) 2 L 2 -2µγ t ∥x t+1/2 -x t ∥ 2 . (4.25)
Hence, writing ∥x t+1/2 -x * ∥ 2 ≤ 2∥x t+1/2 -x t ∥ 2 + 2∥x t -x * ∥ 2 and rearranging, we obtain

E t+1 ≤ E t -µγ t f t -1 2 µγ t -4γ 2 t c 2 L 2 ∥x t -x * ∥ 2 -1 2 1 -4γ 2 t ((a + b) 2 + c 2 )L 2 -2µγ t ∥x t+1/2 -x t ∥ 2 . (4.26)
Since a, b, c ≥ 0 and a + b + c = 1, we also have (a + b) 2 + c 2 ≤ 1, so the step-size assumption (4.22) guarantees that the last term in (4.26) is nonpositive. Likewise, the second part of (4.22) gives µγ t -4γ 2 t c 2 L 2 ≥ 1 2 µγ t , so the energy inequality (4.23) follows and our proof is complete.

■

We finally have all the required building blocks in place to prove Theorem 1.

Proof of Theorem 1. Our proof strategy consists of the following basic steps: (1) We first show that, if the step-size of (BPM) satisfies (4.5) and x 1 is initialized sufficiently close to x * , the sequence of leading states x t+1/2 , t = 1, 2, . . . , remains within the neighborhood U ∩ B of x * where the Legendre bound (4.4) and (SOS) both hold. (2) By virtue of this stability result, the energy inequality (4.23) and the definition of the Legendre exponent allow us to express

D t = D(x * , x t ) as D t+1 ≤ D t -O D 1/(1-β * ) t
up to an error term that vanishes at a geometric rate. The rates (4.6) are then derived by analyzing this recursive inequality for β * = 0 and β * > 0 respectively. We now proceed to detail the two steps outlined above.

Step 1: Stability. Take r > 0 such that B X r (x * ) := {x ∈ X : ∥x -x * ∥ ≤ r} ⊂ U ∩ B and assume further that

x 1/2 = x 1 ∈ U ∩ B is such that D(x * , x 1/2 ) = D(x * , x 1 ) ≤ (1 -λ)r 2 /4,
where λ ∈ (0, 1) is a constant to be determined later. Letting E 0 := D(x * , x 1 ), we will prove by induction on t that

∥x t-1/2 -x * ∥ ≤ r and E t ≤ E t-1 , (4.27)
which will show in particular that x t+1/2 ∈ U ∩ B for all t ≥ 1. Indeed:

• For the base case (t = 1), we have ∥x 1/2 -x * ∥ ≤ 2D(x * , x 1/2 ) ≤ r and E 1 = E 0 by construction, so there is nothing to show.

• For the induction step, assume (4.27) holds. Then (4.1) yields

1 2 ∥x t -x * ∥ 2 ≤ D t ≤ E t ≤ E 1 = D(x * , x 1 ) ≤ 1 4 (1 -λ)r 2 (4.28)
i.e., x t ∈ B X r (x * ). Now, to show that x t+1/2 ∈ B X r (x * ), Lemma 2 with p ← x * , x ← x t , y ← -γ t g t and v ← -(a + b)γ t g(x * ) gives

D t+1/2 ≤ D t -γ t ⟨g t -(a + b)g(x * ), x t+1/2 -x * ⟩ ≤ D t -aγ t ⟨g(x t ) -g(x * ), x t+1/2 -x * ⟩ -bγ t ⟨g(x t-1/2 ) -g(x * ), x t+1/2 -x * ⟩ (4.29)
and hence, by Young's inequality and (4.1), we get

1 2 ∥x t+1/2 -x * ∥ 2 ≤ D t + γ 2 t a∥g(x t ) -g(x * )∥ 2 * + γ 2 t b∥g(x t-1/2 ) -g(x * )∥ 2 * + 1 4 (a + b)∥x t+1/2 -x * ∥ 2 . (4.30)
Since a + b ≤ 1, using Assumption 1 and rearranging gives

∥x t+1/2 -x * ∥ 2 ≤ 4D t + 4γ 2 t L 2 max{∥x t -x * ∥ 2 , ∥x t-1/2 -x * ∥ 2 } ≤ (1 -λ)r 2 + 4γ 2 t L 2 r 2 (4.31)
where we used the fact that ∥x t-1/2 -x * ∥ 2 ≤ r 2 and ∥x t -x * ∥ 2 ≤ 2D t ≤ 1 2 (1 -λ)r 2 (by the inductive hypothesis and (4.28) respectively). Thus, with 2γ t L ≤ 1/φ < 1 by assumption, choosing λ = 1/φ 2 gives ∥x t+1/2 -x * ∥ 2 ≤ r 2 , which completes the first part of the induction. Finally, for the second part, our step-size assumption gives

2µγ t + 4γ 2 t L 2 ≤ 2γ t L + 4γ 2 t L 2 ≤ 1/φ + 1/φ 2 = 1.
(4.32) Thus, since x t+1/2 ∈ B, Proposition 2 readily gives

E t+1 ≤ E t -µγ t f t -1 4 µγ t ∥x t -x * ∥ 2 ≤ E t , (4.33)
and the induction is complete.

Step 2: Convergence rate analysis. From (4.33) and the Legendre bound (4.4), we get

E t+1 ≤ E t -µγ t f t - µγ t 2 1-α κ 1+α D 1+α t with α = β * /(1 -β * ).
(4.34)

We now distinguish two cases, depending on whether β * = 0 or β * > 0.

Case 1: If β * = 0, we have α = 0 by definition and κ ≥ 1 by (4.1). Eq. (4.34) then gives

E t+1 ≤ E t - µγ t 2κ D t -µγ t f t ≤ 1 - µγ t 2κ E t (4.35)
so the case β = 0 of (4.6) follows immediately by setting γ t ≡ γ for all t. Case 2: If β * > 0, then α > 0 too, so we will proceed by rewriting all terms in Eq. (4.34) in terms of E t . To that end, we have:

E t+1 ≤ E t -µγ t f t - µγ t 2 1-α κ 1+α D 1+α t ≤ E t - µγ t D(x * , x 1 ) α f 1+α t - µγ t 2 1-α κ 1+α D 1+α t ≤ E t - µγ t max(2 1-α κ 1+α , D(x * , x 1 ) α ) D(x * , x t ) 1+α + f 1+α t ≤ E t - µγ t max(2κ 1+α , 2 α D(x * , x 1 ) α ) E 1+α t (4.36)
where, in the second line, we used (4.27) to get f t ≤ D(x * , x t )+f t ≤ D(x * , x 1 ), and, in the last line, we used the convexity of x 1+α . The case β ∈ (0, 1) of (4.6) then follows from Lemma 6 of [34, p. 46] (recreated in the appendix as Lemma A.1). This concludes our proof.

Finer results for linearly constrained problems

Our goal in this last section is to take a closer look at the convergence rate of (BPM) for different solution configurations that arise in linearly constrained problems. To that end, we begin by revisiting the examples of Section 3.

Motivating examples, redux.

A common feature of Examples 3.1-3.4 is that the problem's defining vector field vanishes at the solution point under study. In the series of examples below, we examine the rate of convergence achieved when this is not the case.

Example 5.1 (Euclidean regularization). Consider again the quadratic regularizer of Example 3.1 over X = [0, ∞), but with g(x) = x + 1. The solution of (VI) is still x * = 0 but the update (3.2) now becomes

F (x) = [x -γ(x + 1)] + = [(1 -γ)x -γ] + .
(5.1)

Since F (x) = 0 for all sufficiently small x > 0, we readily conclude that x t converges to x * in a finite number of iterations.

F (x) = x exp(-γ(x + 1)) = xe -γ + o(x) ∼ xe -γ (5.2)
i.e., F is a contraction for small x > 0. Hence, in contrast to Example 3.2, x t converges to 0 at a geometric rate, even though the problem's solution lies on the boundary of X . ♦ Example 5.3 (Fractional power). Finally, consider the fractional power regularizer of Example 3.3, again with g(x) = x + 1. Then, for q ∈ (0, 1), the update rule (3.2) gives

F (x) = [x q-1 + γ(1 -q)(x + 1)] 1/(q-1) = x -γx 2-q + o(x 2-q ) (5.3)
for small x > 0. Thus, by Lemma 1, we get that x t converges to 0 as |x t -x * | = Θ t -1/(1-q) and D(x * , x t ) = Θ t -q/(1-q) , which is again faster than the rate given by Theorem 1.

♦

Examples 5.1-5.3 show that the convergence rate of (BPM) can change drastically depending on whether g(x * ) is zero or not. In the example below, we examine in more detail the behavior of the individual coordinates of x t as a function of the position of g(x * ) relative to X .

Example 5.4 (Higher-dimensional simplices). Consider the canonical two

-dimensional simplex X = {(x 1 , x 2 , x 3 ) ∈ R 3 + : x 1 + x 2 + x 3 = 1} of R 3
equipped with the entropic regularizer h(x) = 3 i=1 x i log x i . Consider also the vector field g(x) = x -p with p = (-ν 1 , -ν 2 , 1) for some ν 1 , ν 2 ≥ 0, so the solution of (VI) is x * = (0, 0, 1), an extreme point of X .

Since the Legendre exponent of h at x * is easily seen to be β h (x * ) = 1/2, Theorem 1 would indicate a rate of convergence of D(x * , x t ) = O(1/t) or, in terms of norms, ∥x t -x * ∥ = O(1/t). However, this rate can be very pessimistic if, for example, ν 1 > 0. Indeed, in this case, since x t converges to x * = (0, 0, 1), the relevant coordinates of g(x t ) will evolve as g 1 (x t ) = x 1,t + ν 1 = ν 1 + o(1) and g 3 (x t ) = x 3,t -1 = o(1). Accordingly, since entropic regularization on the simplex leads to the exponential weights update [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF] x i,t+1 ∝ x i,t exp(-γg i (x t )) for all t ≥ 1, i = 1, 2, 3,

(5.4) the fact that lim t→∞ x 3,t = 1 readily yields

x 1,t+1 ∼ x 1,t+1 x 3,t+1 = x 1,t x 3,t exp(-γg 1 (x t ) + γg 3 (x t )) = x 1,t x 3,t exp(-γν 1 + o(1)) (5.5)
i.e., x 1,t converges to 0 at a geometric rate whenever ν 1 > 0.

By symmetry, the argument above yields the same rate for x 2,t if ν 2 > 0. However, as we show in Appendix B, if ν 2 = 0, we would have x 2,t = Θ(1/t) no matter the value of ν 1 (and likewise for the rate of x 1,t if ν 1 = 0). In other words, the rate provided by Theorem 1 is tight for the coordinate i ∈ {1, 2} with a vanishing drift coefficient ν i , but not otherwise; we will devote the rest of this section to deriving a formal statement (and proof) of the general principle underlying this observation. ♦ 5.2. Linearly constrained problems. For concreteness, we focus below on linearly constrained problems: the structural configurations outlined in the previous examples are there more prominent. To simplify the presentation, we identify V with R n endowed with the Euclidean scalar product ⟨•, •⟩, and we do not distinguish between primal and dual vectors (meaning in particular that the distinction between normal and polar cones are likewise blurred).

Formally, we consider polyhedral domains written in normal form as

X = {x ∈ R n + : Ax = b} (5.6)
for some matrix A ∈ R m×n and b ∈ R n . 3 Moreover, to avoid trivialities, we further assume that X admits a Slater point, i.e., there exists some x ∈ X such that x i > 0 for all i = 1, . . . , n. This setup is particularly flexible, as it allows us to identify the active constraints at x ∈ X with the zero components of x.

Elaborating further on this, since ⟨g(x * ), x -x * ⟩ ≥ 0 for all x ∈ X and any solution x * of (VI), we directly infer that -g(x * ) is an element of the normal cone NC(x * ) to X at x * . In our polyhedral setting, NC(x * ) admits an especially simple representation as

NC(x * ) = row(A) -{(ν 1 , . . . , ν n ) ∈ R n + : ν i = 0 whenever x * i = 0} (5.7)
where row(A) = (ker A) ⊥ ⊆ R n denotes the row space of A [16, Ex. 5.2.6]. As a result, we see that x * is a solution of (VI) if and only if g(x * ) can be written in the form

g(x * ) - i∈A ν i e i ∈ row(A) (5.8)
for an ensemble of non-negative slackness coefficients ν i ≥ 0, i ∈ A, where

A ≡ A(x * ) = {i : x * i = 0} (5.9)
denotes the set of inequality constraints of (5.6) that are active at x * . With all this in mind, we will distinguish the following solution configurations (see also Fig. 2).

X NC(x * ) . x * -g(x * ) X NC(x * ) . x * -g(x * ) X NC(x * ) . x * -g(x * )
Figure 2: Different solution configurations: a non-extreme solution where g is sharp (left), an extreme solution where g is not sharp (center), and a sharp solution (right).

Definition 3 (Sharpness). Let x * ∈ X be a solution of (VI) with associated slackness coefficients ν i , i ∈ A, as per (5.8). The set of sharp (♯) and flat (♭) directions at x * are respectively defined as

A ♯ = {i ∈ A : ν i > 0} and A ♭ = {i ∈ A : ν i = 0}, (5.10)
and we say that g is sharp at

x * if A ♯ = A (or, equivalently, if A ♭ = ∅).
The sharpness of g at x * is then defined as ν * = min i∈A ♯ ν i (5.11) and, finally, if x * is an extreme point of X , we will say that x * is itself sharp.

The terminology "sharp" and "flat" alludes to the case where g is a gradient field, and is best illustrated by an example. To wit, let f (x 1 , x 2 ) = x 1 + 1 2 (x 2 -1) 2 for x 1 , x 2 ≥ 0, so f admits a (unique) global minimizer at x * = (0, 1). Applying Definition 3 to g = ∇f , we readily get A ♯ = {1} and A ♭ = {2}, reflecting the fact that f (x 1 , 1) exhibits a sharp minimum at 0 along x 1 whereas the landscape of f (0, x 2 ) is flat to first-order around 1 along x 2 . 5.3. Convergence rate analysis. We are now in a position to state and prove our refinement of Theorem 1 for linearly constrained problems. To that end, following [START_REF] Alvarez | Hessian Riemannian gradient flows in convex programming[END_REF], we will assume in the rest of this section that (BPM) is run with a Bregman regularizer h that is adapted to the polyhedral structure of X as per the definition below: Definition 4. Let X be a polyhedral domain of the form (5.6). Then, a Bregman regularizer h on X is said to be decomposable with kernel θ if

h(x) = n i=1 θ(x i ) for all x ∈ X .
(5.12)

for some continuous function θ : R + → R such that (i) θ ′′ (z) > 0 for all z > 0; and (ii) θ ′′ is locally Lipschitz on (0, ∞).

In addition to facilitating calculations, the notion of decomposability will further allow us to describe the convergence rate of the iterates of (BPM) near the boundary of X in finer detail. In fact, as it turns out, the speed of convergence along a given direction will actually be determined by the behavior of the derivative of the Bregman kernel θ near 0.

In this regard, there are two distinct regimes to consider. First, if lim x→0 + θ ′ (x) = -∞, it is straightforward to see that dom ∂h = ri X so, by Lemma A.2 in Appendix A, the iterates x t of (BPM) will remain in ri X for all t; in this case h is essentially smooth -or Legendre -in the sense of [START_REF] Rockafellar | Convex Analysis[END_REF]Chap. 26], and we will refer to it as steep. Otherwise, if θ ′ (0) exists and is finite, x t may reach the boundary of X in a finite number of iterations; we will refer to this case as non-steep. The key difference between these two regimes is that, in the non-steep case, the algorithm may achieve convergence in a finite number of steps (at least along certain directions). On the other hand, even though finite-time convergence is not possible in the steep regime, the algorithm's rate of convergence may still depend on the boundary behavior of θ. To illustrate this, we will consider the following concrete cases: Assumption 4. Let θ : R + → R be a kernel function as per Definition 4. Then θ ′ exhibits one of the following behaviors as x → 0 + : (a) Euclidean-like:

lim inf x→0 + θ ′ (x) > -∞. (b) Entropy-like: lim inf x→0 + [θ ′ (x) + log x] > -∞. (c) Power-like: lim inf x→0 + x p θ ′ (x) > -∞ for some p ∈ (0, 1).
Remark. Cases (a)-(c) above respectively mean that |θ ′ (x)| grows as O(1), O(|log x|) or O(1/x p ) as x → 0 + . Clearly, we have (a) =⇒ (b) =⇒ (c) so these cases are not exclusive; nonetheless, to avoid overloading the presentation, when we assume (b), we will tacitly imply that (a) does not also hold at the same time -and likewise for (c).

♦ With all this in hand, we proceed below to show that, in linearly constrained problems, (BPM) converges along sharp directions at x * at an accelerated rate relative to Theorem 1: sublinear rates may become linear, and linear rates transform to convergence in finite time.

Theorem 2. Suppose that (BPM) is run in a polyhedral domain with a decomposable regularizer as per Definition 4. Suppose further that Assumptions 1-4 hold, and that the method's step-size and initialization satisfy the requirements of Theorem 1. Then, for all i ∈ A ♯ , we have: (a) Under Assumption 4(a): there exists some T ≥ 1 such that x i,t = 0 for all t ≥ T (5.13a) (b) Under Assumption 4(b):

x i,t = O exp(-γν eff t/2) (5.13b)
(c) Under Assumption 4(c):

x i,t = O (γν eff t/2) -1/p (5.13c)
where

ν eff = ν * if g is sharp at x * (i.e., A ♭ = ∅), ν * /ϱ otherwise, (5.14)
and ϱ ≡ ϱ(A, b, x * ) ≥ 1 is a positive constant that depends only on X and x * .

In particular, if x * is sharp, we have the following immediate corollary of Theorem 2.

Corollary 2. If x * is sharp, then ∥x t -x * ∥ → 0 at a rate given by (5.13a), (5.13b), or (5.13c), depending respectively on whether Case (a), (b), or (c) of Assumption 4 holds.

Proof. First, note that x * is sharp if and only if span{e i : i ∈ A ♯ } + row(A) = R n : indeed, since X is a polyhedron, x * is extreme if and only if NC(x * ) has nonempty topological interior, and this, combined with (5.7) and the fact that A ♯ = A (since g is sharp at x * ),

Tsallis [q(1 -q)] -1 (x -x q ) O(1/t q/(2-q) ) O(1/t 1/(1-q) ) Hellinger - √ 1 -x 2 O(1/t 1/3 ) O(1/t 2 )
Table 2: proves our assertion. We thus conclude that, for all i = 1, . . . , n, there exist λ ij ∈ R, j ∈ A ♯ , such that e i -j∈A ♯ λ ij e j ∈ row(A), and hence, for all t = 1, 2, . . . , we have x i,t -x * i = j∈A ♯ λ ij x j,t . Our claim then follows from Theorem 2 and the fact that all norms are equivalent on R n .

Summary

■

To facilitate comparisons with the non-sharp regime, the guarantees of Theorem 2 are juxtaposed with those of Theorem 1 in Table 2. Beyond this comparison, Theorem 2 is the main result of this section so, before proving it, some further remarks are in order.

Remark 1 (Solution configurations). By construction (and the fact that X admits a Slater point), it is straightforward to verify that g is sharp at x * if and only if g(x * ) ∈ -ri(NC(x * )); likewise, x * is itself sharp if and only if g(x * ) ∈ -int(NC(x * )). As we noted in the proof of Corollary 2, the latter condition is equivalent to asking that span{e i : i ∈ A ♯ } + row(A) = R n , a condition which describes precisely the informal requirement that the sharp directions at x * suffice to characterize it. By contrast, if g is sharp at some non-extreme point x * , there exists some (nonzero) z ∈ TC(x * ) such that ⟨g(x * ), z⟩ = 0, indicating that the accelerated rates of Theorem 2 cannot be active along the residual direction z. We illustrate these distinct solution configurations in Fig. 2. ♦ Remark 2 (Tightness and the structure of X ). We underline that the dependence of ν eff on the structure of X in the second branch of (5.14) cannot be lifted. To see this, let

X = {x ∈ R 2 + : x 1 = εx 2 } (5.15)
i.e., A = [1 -ε] and ker A is spanned by the vector z = (ε, 1). Then, if we take g(x) = x-p with p 1 ≤ 0 and p 2 ≥ 0 (so that the origin is a solution), and we equip X with the Bregman regularizer induced by the entropic kernel θ(x) = x log x, a straightforward calculation shows that the iterates of (MD) satisfy the recursion

εθ ′ (x 1,t+1 ) + θ ′ (x 2,t+1 ) = εθ ′ (x 1,t ) + θ ′ (x 2,t ) + γ εx 1,t + x 2,t -γ⟨p, z⟩. (5.16)
Thus, letting χ t = x 2,t = x 1,t /ε, the above can be rewritten as

εθ ′ (εχ t+1 ) + θ ′ (χ t+1 ) = εθ ′ (εχ t ) + θ ′ (χ t ) + γ(ε 2 + 1)χ t -γ⟨p, z⟩ (5.17)
and hence, with θ ′ (x) = 1 + log x, we finally get

χ t+1 = χ t exp -γ ε 2 + 1 ε + 1 χ t + γ ⟨p, z⟩ ε + 1 . (5.18)
Now, if ⟨p, z⟩ < 0, we readily infer that χ t = x 2,t converges linearly to 0 at a rate of χ t ∼ exp(-|⟨p, z⟩|/(1 + ε) • γt) as predicted by Theorem 2. In particular, if g(0) = -p = (ν * , ν * ) with ν * > 0, the iterates of (MD) converge geometrically to zero with exponent γν * , which matches the estimate of Theorem 2 up to a factor of 1/2 in the exponent. On the contrary, if g(0) = -p = (ν * , 0), the term |⟨-p, z⟩| = εν * depends on the linear structure of X and can be arbitrarily bad as ε goes to zero. This illustrates why one cannot do away with the dependence on the linear structure of X when g is not sharp at x * . ♦ 5.4. Proof of Theorem 2. We now proceed to the proof of Theorem 2, beginning with two helper lemmas tailored to the polyhedral structure of X . The first is a book-keeping result regarding the subdifferentiability of h.

Lemma 5. Let h be a decomposable regularizer on X with kernel θ as per Definition 4. Then the domain of subdifferentiability of h is X h = {x ∈ X : x i ∈ dom ∂θ for all i = 1, . . . , n} and a continuous selection of ∂h is given by the expression

∇h(x) = n i=1 θ ′ (x i )e i for all x ∈ X h . (5.19)
Proof. See [START_REF] Rockafellar | Convex Analysis[END_REF]Thm 23.8], whose conditions are satisfied because X is polyhedral.

■

The second ingredient we will need is a separation result in the spirit of Farkas' lemma.

Lemma 6. Let X be a polyhedral domain of the form (5.6). Then, for all x * ∈ X , there exists P = P (A, b, x * ) ≥ 1 such that, for all I ⊆ A ≡ A(x * ), at least one of the following holds:

(a) I ̸ = ∅ and there exists i ∈ A \ I such that x i ≤ P max{x j : j ∈ I} for all x ∈ X .

(b) There exists z ∈ ker A such that ∥z∥ ≤ P , z i = 0 if i ∈ I and P ≥ z i ≥ 1 if i ∈ A \ I.

The proof of Lemma 6 is based on Farkas' lemma so we relegate it to Appendix A. Armed with all this, we can finally proceed to prove our result for linearly constrained problems.

Proof of Theorem 2. We will consider two main cases, namely lim x→0 + h ′ (x) = -∞ (the steep case) and lim x→0 + h ′ (x) > -∞ (the non-steep case). The steep regime will cover Cases (b) and (c) of Assumption 4, whereas the non-steep regime will account for Case (a).

Case 1: the steep regime. Note first that, without loss of generality, Cases (b) and (c) respectively imply that there exist C ∈ R and δ > 0 such that, for all x ∈ (0, δ), we have:

Under Assumption 4(b): θ ′ (x) ≥ log x -C (5.20a) Under Assumption 4(c): θ ′ (x) ≥ -Cx -p (5.20b)
With this in mind, let r > 0 be sufficiently small so that B r := {x ∈ X : ∥x -x * ∥ ≤ r} satisfies:

• B r ⊆ U ∩ B with B and U defined by (SOS) and (4.4) respectively.

• If x ∈ B r then x i < δ for all i ∈ A.

• If x ∈ B r , then ∥g(x) -g(x * )∥ * ≤ ν * /(2P ) with ν * by (5.11) and P by Lemma 6.

Now, recall that

Step 1 of the proof of Theorem 1 implies that the iterates of (BPM) will remain in B r for all t if x 1 is initialized sufficiently close to x * . Given this stability guarantee, we will construct below two sets I ♯ ⊆ A ♯ , I ♭ ⊆ A ♭ such that, for all i ∈ I := I ♯ ∪ I ♭ , we have

x i,t ≤ P |I| • exp(C + P ∥∇ h(x 1 )∥ * +P R-γν eff (t -1)/2) under Assumption 4(b), C[γν eff (t -1)/2 -P ∥∇ h(x 1 )∥ * -P R] -1 p + under Assumption 4(c), (5.21)
where R := sup x∈Br i / ∈A θ ′ (x i )e i * , P ≥ 1 is the constant given by Lemma 6 and ν eff is defined as in (5.14) with ϱ = P |A|.

Our construction proceeds inductively, starting with I ♯ = I ♭ = ∅, for which the stated property holds trivially. For the inductive step, if (5.21) holds for I ♯ ⊊ A ♯ and I ♭ ⊆ A ♭ , we will show that there exists some j ∈ A \ I such that (5.21) still holds for I ∪ {j}. Since the number of active constraints is finite, the progressive addition of these indices will allow us to reach I ♯ = A ♯ , thus proviing our initial claim.

To carry all this out, assume that

I ♯ ⊊ A ♯ , I ♭ ⊆ A ♭ ,
θ ′ (x i,t )z i + i∈I θ ′ (x i,t )z i = ⟨∇h(x 1 ), z⟩ -γ t-1 s=1 i∈A ν i z i + γ t-1 s=1 ⟨g(x * ) -g s+1/2 , z⟩ - i / ∈A θ ′ (x i,t )z i . (5.25)
Finally, by the properties we used to construct z, we further have .26) where the second inequality follows from how we chose B r at the beginning of the proof.

i∈A\I θ ′ (x i,t )z i ≤ P ∥∇h(x 1 )∥ * -γ(t -1) i∈A\I ν i z i + γ t-1 s=1 P ∥g(x * ) -g s+1/2 ∥ * + P i / ∈A θ ′ (x i,t )e i * ≤ P ∥∇h(x 1 )∥ * -γ(t -1) i∈A\I ν i z i + 1 2 γν * (t -1) + P R. ( 5
We conclude by distinguishing whether g is sharp at x * (i.e., if A ♭ = ∅ or not). Case 1.1: If A ♭ = ∅, we have A \ I = A ♯ \ I ♯ and ν i ≥ ν * for all i ∈ A \ I, so (5.26) gives i∈A\I θ ′ (x i,t )z i + γ(t -1)ν * z i ≤ P ∥∇h(x 1 )∥ * + 1 2 γν * (t -1) + P R.

(5.27)

Choosing the coordinate j ∈ A \ I corresponding to the smallest term in the sum on the left-hand side (LHS), we get

(θ ′ (x j,t ) + γ(t -1)ν * )(|A| \ I)z j ≤ P ∥∇h(x 1 )∥ * + 1 2 γν * (t -
′ (x j,t )z j ≤ P ∥∇h(x 1 )∥ * - γν * 2|A| (t -1) + P R . (5.31)
Dividing both sides by z j and using that it lies between 1 and P gives θ ′ (x j,t ) ≤ P ∥∇h(x 1 )∥ * -γν * 2P |A| (t -1) + P R (5.32)

Thus, in both cases, there exists j ∈ A \ I such that (5.32) holds, since P |A| ≥ 1. Therefore, combining this inequality with (5.20) we conclude that (5.21) holds for j, and since P ≥ 1, we conclude that we can augment I ♯ or I ♭ by j, depending on whether it belongs to A ♯ or A ♭ .

Case 2: the non-steep regime. The proof borrows the structure of the first case, though it is more direct. Take r > 0 small enough such that B r := {x ∈ X : ∥x -x * ∥ ≤ r} satisfies

• B r is included in U ∩ B, and if x ∈ B r , then ∥g(x) -g(x * )∥ * ≤ ν * 3P , • if x, x ′ ∈ B r then, ∥∇h(x) -∇h(x ′ )∥ * ≤ γν *
3P where P is given by Lemma 6. This is possible since ∇h is continuous at x * , • No other constraint x i = 0 with i / ∈ A becomes active in B r . As we have seen in the proof of Theorem 1, if (BPM) is initialized close enough to x * , then all the iterates x t and the half-iterates x t+1/2 for s = 1, 2, . . . are contained in B r .

As above, fix some t ≥ T . We will build sets I ♯ ⊆ A ♯ , I ♭ ⊆ A ♭ with the property that that for all i ∈ I ♯ ∪ I ♭ , x i,t = 0 .

(5.33)

Starting with I ♯ = I ♭ = ∅, Eq. (5.33) is trivially verified. Now, take I ♯ ⊊ A ♯ , I ♭ ⊆ A ♭ which satisfy the desired property and, as before, apply Lemma 6 with I := I ♯ ∪ I ♭ . If the first case of Lemma 6 holds, then I ̸ = ∅ and there exists j ∈ A such that,

x j,t ≤ P max(x i,t : i ∈ I) (5.34) which yields the result by adding i to I ♯ or I ♭ depending whether it belongs to A ♯ or A ♭ . Otherwise, if the second case holds, there is some z ∈ ker A such that ∥z∥ ≤ P , z i = 0 if i ∈ I and P ≥ z i ≥ 1 if i ∈ A \ I. For the sake of contradiction, assume that for all i ∈ A \ I, x i,t > 0. Showing that this results in a contradiction will give us an additional coordinate j ∈ A \ I for which x j,t = 0 that we will then add to I ♯ or I ♭ as in the first case. Now, let us determine the normal cone at x t . Since x t belongs to B X r (x * ), no other constraint other than the ones of I can become active, and these constraints are actually active by the definition of I and Eq. (5.33). Hence, the normal cone at x t (see Eq. (5.7)) becomes NC(x t ) = -i∈I ν i e i : (ν i ) i∈I ∈ (R + ) I + row(A) .

(5.35)

Taking a scalar product between the last inclusion of Lemma A.2 and z, we get that, ∇h(x t ) -∇h(x t-1 ) -γg t-1/2 , z = 0 .

( 5.36) This means, from the definition of B X r (x * ), that, γ⟨g(x * ), z⟩ = ⟨∇h(x t ) -∇h(x t-1 ), z⟩ + γ⟨g(x * ) -γg t-1/2 , z⟩ ≤ 2γν * 3 (5.37) However, by (5.8) and the properties of z, we also have

⟨g(x * ), z⟩ = i∈A\I ν i z i ≥ ν * (5.38)
which is in contradiction with (5.37). We may therefore iteratively add coordinates of A for which x i,t = 0; this completes the induction and our proof. ■

Concluding remarks

Our results indicate that Euclidean regularization leads to faster trajectory convergence rates near second-order sufficient (SOS) solutions. While this does not contradict the analysis of [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF] -which concerns the method's ergodic average and advocates the use of non-Euclidean regularizers in domains with a favorable geometry -it does run contrary to its spirit. We attribute the source of this discrepancy to the fact that Lipschitz continuity and second-order sufficiency are both norm-based conditions, so it is plausible to expect that norm-based regularizers would lead to better results. This raises the question of what the corresponding rate analysis would give in the case of Bregman-based variants of (LC) and (SOS), e.g., as in the recent works of [START_REF] Antonakopoulos | An adaptive mirror-prox algorithm for variational inequalities with singular operators[END_REF][START_REF] Bauschke | A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications[END_REF][START_REF] Lu | Relatively-smooth convex optimization by first-order methods and applications[END_REF]. We defer this analysis to future work.

Case (a). Assume that there exists (λ i ) i∈A\I ∈ (R + ) A\I not all zero and (µ i ) i∈I ∈ R I such that (A.18) holds. In this case, I must be non-empty since otherwise X would be reduced to {0} (see the first inclusion), violating the definition (5.6) of X . In addition, there is some i ∈ A \ I such that λ i > 0 and thus we have ∀x ∈ X , x i ≤ max (|λ j | : j ∈ I) λ i max(x j : j ∈ I) (A. [START_REF] Kakade | Regularization techniques for learning with matrices[END_REF] which corresponds to the first case of the lemma.

Case (b). For all (λ i ) i∈A\I ∈ (R + ) A\I not all zero and (µ i ) i∈I ∈ R I , (A.18) does not hold. To interpret this situation, we use the fact that X is of the general polyhedral form (5.6) so aff X = x * + ker A and x * always satisfies i∈A\I λ i x * i = i∈I µ i x * i = 0 so that Eq. (A.18) ⇐⇒ aff X ⊂

   x ∈ R n : i∈A\I λ i x i = i∈I µ i x i    ⇐⇒ ker A ⊂    x ∈ R n : i∈A\I λ i x i = i∈I µ i x i    ⇐⇒ i∈A\I λ i e i -

Figure 1 :

 1 Figure 1: The rate of convergence of (MD) in Examples 3.1-3.4. The Euclidean and shifted Hellinger regularizers lead to a geometric rate (see left figure); all other examples converge at a polynomial rate.

. 13 )

 13 Proof. Let c = 1 -a -b so c ≥ 0 by Assumption 3. Corollary 1 then yields

■

♦ Example 5 . 2 (

 52 Entropic regularization). Under the entropic regularizer of Example 3.2, and taking again g(x) = x + 1, the update rule (3.2) becomes

  of the accelerated rates of convergence observed along sharp directions depending on the Bregman kernel (cf. Definition 4). The Euclidean, entropic and Tsallis kernels are the prototypical examples of Cases (a)-(c) of Assumption 4.

i∈I µ i e-x 2 1 1 -x 2 ) 2 = G(x) 1 - 2 = 2 = 2 = 1 - 2 = 1 2 = - 1 +♦Example 5 . 4 (

 21212221212154 i ∈ row(A). (A.20) = x -γ(x + 1) √ 1 -x 2 + (x -γ(x + 1) √ x 2 + G(x) 2 , (B.1) with G(x) = x -γ(x + 1) √ 1 -x 2. Now, the behavior of G near x * = -1 can be approximated asG(x) = x -γ(x + 1) 3/2 (1 -x) 1/-1 + (x + 1) -γ(x + 1) 3/2 (2 -(x + 1)) 1/-1 + (x + 1) -√ 2γ(x + 1) 3/2 1 -1 4 (x + 1) + o(x + 1) = -1 + (x + 1) -√ 2γ(x + 1) 3/2 + √ 2γ 4 (x + 1) 5/2 + o (x + 1) 5/2 . (B.2)Another Taylor expansion then yieldsG(x) 2 = 1 -(x + 1) + √ 2γ(x + 1) 3/2 -√ 2γ 4 (x + 1) 5/2 + o (x + 1) 5/2 2(x + 1) + 2 √ 2γ(x + 1) 3/2 + (x + 1) 2 -√ 2γ 2 (x + 1) 5/2 + o (x + 1) 5/2 (B.3)so the denominator of Eq. (B.1) becomes1 -x 2 + G(x) 2 = (x + 1)(2 -(x + 1) + G(x) 2 + 1) 5/2 + o (x + 1) 5/2 . (B.4)Thus, plugging this expansion and Eq. (B.2) into Eq. (B.1) givesF (x) = -1 + (x + 1) -√ 2γ(x + 1) 3/2 + √ 2γ 4 (x + 1) 5/2 + o (x + 1) 5/2 1 + √ 2γ(x + 1) 3/2 -√ 2γ 4 (x + 1) 5/2 + o (x + 1) 5/2 = -1 + (x + 1) -√ 2γ(x + 1) 3/2 + √ 2γ 4 (x + 1) 5/2 + o (x + 1) (x + 1) -2 √ 2γ(x + 1) 5/2 + o (x + 1) 5/2 , (B.5)which gives our assertion when x * = -1. Three-dimensional simplex, continuing from p. 15

  and apply Lemma 6 to I = I ♯ ∪ I ♭ . If the first case of Lemma 6 holds, then I ̸ = ∅ and there exists i ∈ A such that x i,t ≤ P max j∈I x j,t(5.22)

			t-1 s=1 g s+1/2 ∈ row(A).
	Thus, taking the scalar product with z ∈ ker A yields
	n	t-1
	θ ′ (x i,t )z i = ⟨∇h(x 1 ), z⟩ -γ	⟨g s+1/2 , z⟩	(5.24)
	i=1	s=1
	so, after rearranging and invoking (5.8) to write ⟨g(x

so (5.21) still holds when i is appended to I ♯ or I ♭ . Otherwise, the second case of Lemma 6 holds and there exists z ∈ ker A with ∥z∥ ≤ P , z i = 0 if i ∈ I, and

P ≥ z i ≥ 1 if i ∈ A \ I. Since h is steep, this means that x t+1 = P xt (-γg t+1/2

) belongs to X h = ri X for all t = 1, 2, . . . , so NC(x t+1 ) is the affine hull of X , i.e., NC(x t+1 ) = row(A). Hence, Lemma A.2 guarantees ∇h(x t+1 ) -∇h(x t ) + γg t+1/2 ∈ row(A)

(5.23) so, telescoping from s = 1 to t -1, we get ∇h(x t ) -∇h(x 1 ) + γ * ), z⟩ = i∈A ν i z i , we get i∈A\I

  If A ♭ ̸ = ∅, then, since I ♯ ⊊ A ♯ , the intersection of A \ I and A ♯ is not empty so that, i∈A\I ν i z i ≥ i∈A\I ν i ≥ ν * and the inequality (5.26) above becomes,

	1) + P R	(5.28)
	and noting that |A \ I|z j ≥ 1 yields
	θ ′ (x j,t ) ≤ P ∥∇h(x 1 )∥ * -1 2 γν * (t -1) + P R .	(5.29)
	Case 1.2:

i∈A\I θ ′ (x i,t )z i ≤ P ∥∇h(x 1 )∥ * -1 2 γν * (t -1) + P R . (5.30)

Now, choosing j to be the coordinate in A \ I which minimizes the LHS and bounding |A \ I| by 1 and |A|, we get that θ

  ). We conclude our treatment of the simplex by showing that x 2,t ∼ x 2,t /x 3,t = Ω(1/t) if ν 2 = 0 but ν 1 > 0. To begin with, we have g 2 (x t ) = x 2,t = o(1) so, arguing as in the first part of the example, we readily get x 1,t+1 x 2,t+1 = x 1,t x 2,t exp(-γν 1 + o(1)), (B.6) so x 1,t /x 2,t converges to 0 at a geometric rate. Accordingly, the quantity x 3,t /x 2,t is /x 2,t vanishes as t → ∞. Hence, after telescoping, we conclude that ≤ 2γt + o(t), which in turn shows that x 2,t ∼ x 2,t /x 3,t = Ω(1/t), as claimed.

	bounded as
	x 3,t+1 x 2,t+1	=	x 3,t x 2,t	exp(γg 2 (x t ) -γg 3 (x t )) =	x 3,t x 2,t	exp(γx 2,t -γ(x 3,t -1))
		=	x 3,t x 2,t	exp(2γx 2,t + γx 1,t ) ≤	x 3,t x 2,t	exp 2γ	x 2,t x 3,t	+ γx 1,t	(B.7)
	Now, since both	x2,t x3,t and x 1,t go to zero,
		x 3,t+1 x 2,t+1	≤	x 3,t x 2,t	1 + 2γ	x 3,t x 2,t	+ γx 1,t + o 2	x 3,t x 2,t	+ x 1,t
				=	x 3,t x 2,t	+ 2γ + o(1) .	(B.8)
	since x 1,t x3,t x2,t ♦

Inequality constraints of the form Ax ≤ b can also be accommodated in(5.6) by introducing the associated slack variables s = b -Ax ≥ 0. Even though this leads to a more verbose presentation of X , the form (5.6) is much more convenient in terms of notational overhead, so we stick with the equality formulation throughout.

Acknowledgments

The authors gratefully acknowledge financial support by the French National Research Agency (ANR) in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02), the LabEx PERSYVAL (ANR-11-LABX-0025-01), MIAI@Grenoble Alpes (ANR-19-P3IA-0003), and the bilateral ANR-NRF grant ALIAS (ANR-19-CE48-0018-01). This work has also been partially supported by project MIS 5154714 of the National Recovery and Resilience Plan Greece 2.0 funded by the European Union under the NextGenerationEU Program. Part of this work was done while P. Mertikopoulos was visiting the Simons Institute for the Theory of Computing.

Appendix A. Auxiliary results

We provide here a series of basic properties, helper lemmas and auxiliary results that we use repeatedly in our paper.

A.1. Lemmas on numerical sequences. The first two results concern numerical sequences.

Lemma A.1. Consider two sequences of real numbers u t , α t ≥ 0, t = 1, 2, . . . , such that u t+1 ≤ u t -α t u 1+r t for some r > 0 and all t = 1, 2, . . .

Then, for all t = 1, 2, . . . , we have:

Proof. See [34, p. 46].

■

The second result that we prove here is a slight variant of the above lemma.

Lemma 1. Suppose that f : R + → R + admits the asymptotic expansion

for positive constants λ, r > 0. Then, for u 1 > 0 small enough, the sequence

Proof. By the assumption on f , there exists some ε > 0 such that

Note first that, if u 1 ≤ ε, Lemma A.1 readily implies that u t converges to 0 and that u t ≤ ε for all t. Moreover, if ε is small enough so that 1 -2λε r > 0 and u 1 is positive, this implies that all u t , for t = 1, 2, . . . , are positive. In particular, we consider the sequence u -r t , t = 1, 2, . . . , for which we get

Hence, u -r t ∼ rλt which gives the result.

■ A.2. Properties of Bregman divergences and the induced prox-mappings. We recall some basic properties of the Bregman divergence and the induced prox-mapping. Variants of these properties are fairly well known in the literature, so we omit their proofs and we refer to [START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF][START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF].

Lemma A.2. Let h be a Bregman regularizer on X and let ∇h be a continuous selection of ∂h. Then, for all x ∈ X h , x + ∈ X and y ∈ Y, we have:

where PC(x) = {y ∈ Y : ⟨y, x ′ -x⟩ ≤ 0 for all x ′ ∈ X } denotes the polar cone to X at x.

Lemma A.3 (3-point identity). For all p ∈ X and all x, x + ∈ X h , we have:

Lemma A.4 (Non-expansiveness). For all x ∈ X h and all y, y + ∈ Y we have:

The next two results that we provide consider the evolution of the Bregman divergence before and after a prox step (or two):

Then, for all p ∈ X , v ∈ PC(p), we have:

Then, after rearranging to isolate D(p, x + ), we get

where the inequality in the last line follows from Lemma A.2. Hence, given that ⟨v, x + -p⟩ ≤ 0 by the fact that v ∈ PC(p), we readily obtain

For the second inequality of the lemma, note that

where the penultimate inequality follows directly from Young's inequality and the last one from (4.1). Our assertion is then obtained by combining this last bound with (A.10). ■ Lemma 3. Let x + i = P x (y i ) for some x ∈ X h and y i ∈ Y, i = 1, 2. Then, for all p ∈ X and all v ∈ PC(p), we have:

To lower bound D(x + 2 , x), we use again Lemma 2 with p ← x + 2 and x + 1 = P x (y 1 ), which gives

and, after rearranging to isolate D(x + 2 , x) and substituting the resulting bound in (A.12),

Thus, by Young's inequality and the strong convexity of h, we finally obtain

and our proof is complete.

■ A.3. Legendre exponent for interior points. We discussed in Section 4.1 that β h (p) = 0 whenever p is an interior point. We formalize this intuition below.

Lemma A.5. Suppose that ∇h is locally Lipschitz continuous. Then β h (p) = 0 for all p ∈ X h ; in particular, β h (p) = 0 whenever p ∈ ri X .

Proof. Fix some p ∈ X h and suppose that ∇h is locally Lipschitz continuous. Then there exists a neighborhood U of p in X and some κ > 0 such that

Now, since ∇h(p) ∈ ∂h(p), we also have

for all x ∈ U ∩ X h . This shows that (4.3) holds with β = 0, i.e., β h (p) = 0.

■ A.4. A separation result. We now prove Lemma 6, which we restate below for convenience: Lemma 6. Let X be a polyhedral domain of the form (5.6). Then, for all x * ∈ X , there exists P = P (A, b, x * ) ≥ 1 such that, for all I ⊆ A ≡ A(x * ), at least one of the following holds:

(a) I ̸ = ∅ and there exists i ∈ A \ I such that x i ≤ P max{x j : j ∈ I} for all x ∈ X .

(b) There exists z ∈ ker A such that ∥z∥ ≤ P , z i = 0 if i ∈ I and

Proof. Our claim is trivial if I = A, so we will focus exclusively on the case I ⊊ A. The stated constant P = P (A, b, x * ) will then be obtained as the maximum of 1 and the constants we obtain for each possible I ⊊ A.

The proof consists in discussing whether there exists (λ i ) i∈A\I ∈ (R + ) A\I not all zero and (µ i ) i∈I ∈ R I such that the inclusion

holds. Case (a) considers the case when such coefficients exist, while Case (b) when this is not possible. Therefore, the fact that (A.18) does not hold for all (λ i ) i∈A\I ∈ (R + ) A\I not all zero and (µ i ) i∈I ∈ R I means that, the system,

with variables (λ i ) i∈A\I ∈ (R + ) A\I not all zero, (µ i ) i∈I ∈ R I , r ∈ R m , has no solution. Hence, by Motzkin's theorem on the alternative (see e.g., [8, Ex. 1.4.2]) 4 , this means that the system,

admits a solution z ∈ R n . Rescaling z and setting P to max(∥z∥, ∥z∥ ∞ ) then gives the second case.

■ Appendix B. Omitted calculations

In this appendix, we provide some computational details that were left out of the main text to streamline our presentation.

Example 3.4 (Hellinger distance, continuing from p. 6). We proceed to compute the Taylor expansion of F near x * = -1 for the shifted operator g(x) = x + 1. Indeed, in this case, the fixed point operator F is given by F (x) = P x (-γg(x)) = P x (-γ(x + 1)) 4 With the notations of [START_REF] Ben-Tal | Lectures on modern convex optimization: analysis, algorithms, and engineering applications[END_REF]Ex. 1.4.2], the lines of the matrix S are made of the e i for i ∈ A \ I and the lines of the matrix N are the e i for i ∈ I, -e i for i ∈ I, the lines of A and their opposite.