Title: Online and Stochastic Optimization beyond Lipschitz Continuity: A Riemannian Approach

Abstract: Motivated by applications to machine learning and imaging science, we study a class of online and stochastic optimization problems with loss functions that are not Lipschitz continuous; in particular, the loss functions encountered by the optimizer could exhibit gradient singularities or be singular themselves. Drawing on tools and techniques from Riemannian geometry, we examine a Riemann-Lipschitz (RL) continuity condition which is tailored to the singularity landscape of the problem's loss functions. In this way, we are able to tackle cases beyond the Lipschitz framework provided by a global norm, and we derive optimal regret bounds and last iterate convergence results through the use of regularized learning methods (such as online mirror descent). These results are subsequently validated in a class of stochastic Poisson inverse problems that arise in imaging science.

1. Introduction:

The surge of recent breakthroughs in machine learning and artificial intelligence has reaffirmed the prominence of first-order methods in solving large-scale optimization problems. One of the main reasons for this is that the computation of higher-order derivatives of functions with thousands -if not millions -of variables quickly becomes prohibitive; another is that gradient calculations are typically easier to distribute and parallelize, especially in large-scale problems. In view of this, first-order methods have met with prolific success in many diverse fields, from machine learning and signal processing to wireless communications, nuclear medicine, and many others [10, 34, 37] . This success is especially pronounced in the field of online optimization, i.e., when the optimizer faces a sequence of time-varying loss functions f t , t = 1, 2, . . . , one at a time -for instance, when drawing different sample points from a large training set [11, 35] . In this general framework, first-order methods have proven extremely flexible and robust, and the attained performance guarantees are well known to be optimal [1, 11, 35] . Specifically, if the optimizer faces a sequence of G-Lipschitz convex losses, the incurred min-max regret after T rounds is Ω(GT 1/2 ), and this bound can be achieved by inexpensive first-order methods -such as online mirror descent and its variants [11, 35, 36, 41] .

Nevertheless, in many machine learning problems (support vector machines, Poisson inverse problems, quantum tomography, etc.), the loss landscape is not Lipschitz continuous, so the results mentioned above do not apply. Thus, a natural question that emerges is the following: Is it possible to apply online optimization tools and techniques beyond the standard Lipschitz framework? And, if so, how?

Our approach and contributions. Our point of departure is the observation that Lipschitz continuity is a property of metric spaces -not normed spaces. Indeed, in convex optimization, Lipschitz continuity is typically stated in terms of a global norm (e.g., the Euclidean norm), but such a norm is de facto independent of the point in space at which it is calculated. Because of this, the standard Lipschitz framework is oblivious to the finer aspects of the problem's loss landscape -and, in particular, any singularities that may arise at the boundary of the problem's feasible region. On the other hand, in general metric spaces, this is no longer the case: the distance between two points is no longer given by a global norm, so it is much more sensitive to the geometry of the feasible region. For this reason, if the (Riemannian) distance dist(x, x ) between two points x and x becomes larger and larger as the points approach the boundary of the feasible region, a condition of the form | f (x)f (x )| = O(dist(x , x)) may still hold even if f becomes singular at the boundary.

We leverage this observation by introducing the notion of Riemann-Lipschitz (RL) continuity, an extension of "vanilla" Lipschitz continuity to general spaces endowed with a Riemannian metric. We show that this metric can be chosen in a principled manner based solely on the singularity landscape of the problem's loss functions -i.e., their growth rate at infinity and/or the boundary of the feasible region. Subsequently, using a similar mechanism to choose a Riemannian regularizer, we provide an optimal O(T 1/2 ) regret guarantee through the use of regularized learning methods -namely, "follow the regularized leader" (FTRL) and online mirror descent (OMD).

Our second contribution concerns an extension of this framework to stochastic programming. First, in the context of stochastic convex optimization, we show that an online-to-batch conversion yields an O(T -1/2 ) value convergence rate. Second, motivated by applications to nonconvex stochastic programming (where averaging is not a priori beneficial), we also establish the convergence of the method's last iterate in a class of nonconvex problems satisfying a weak secant inequality. Finally, we supplement our theoretical analysis with numerical experiments in Poisson inverse problems.

Related work. To the best of our knowledge, the first treatment of a similar question was undertaken by Bauschke et al. [3] who focused on deterministic, offline convex programs ( f t = f for all t) without a Lipschitz smoothness assumption (i.e., Lipschitz continuity of the gradient, as opposed to Lipschitz continuity of the objective). To tackle this issue, Bauschke et al. [3] introduced a second-order "Lipschitz-like" condition of the form ∇ 2 f β∇ 2 h for some suitable Bregman function h, and they showed that Bregman proximal methods achieve an O(1/T ) value convergence rate in offline convex problems with perfect gradient feedback.

Always in the context of deterministic optimization, Bolte et al. [8] extended the results of Bauschke et al. [3] to unconstrained non-convex problems and established trajectory convergence to critical points for functions satisfying the Kurdyka-Łojasiewicz (KL) inequality. In a slightly different vein, Lu et al. [25] considered functions that are also strongly convex relative to the Bregman function defining the Lipschitz-like condition for the gradients, and they showed that mirror descent achieves a geometric convergence rate in this context. Finally, in a very recent preprint, Hanzely et al. [17] examined the rate of convergence of an accelerated variant of mirror descent under the same Lipschitz-like smoothness assumption. Importantly, all these works concern offline, deterministic optimization problems with perfect gradient feedback and regularity assumptions that cannot be exploited in an online optimization setting (such as the KL inequality). Beyond offline, deterministic optimization problems, Lu [24] established the ergodic convergence of mirror descent in stochastic non-adversarial convex problems under a "relative continuity" condition of the form ∇ f (x) ≤ G inf x √ 2D(x , x)/ xx (with D denoting the divergence of an underlying "reference" Bregman function h). More recently, Hanzely and Richtárik [16] examined the performance of stochastic mirror descent under a combination of relative strong convexity and relative smoothness / Lipschitz-like conditions, and established a series of convergence rate guarantees that mirror the corresponding rates for ordinary (Euclidean) stochastic gradient descent. Except for trivial cases, these conditions are not related to Riemann-Lipschitz continuity, so there is no overlap in our results our methodology.

Finally, in a very recent paper, Bécigneul and Ganea [5] established the convergence of a class of adaptive Riemannian methods in geodesically convex problems (extending in this way classical results for AdaGrad to a manifold setting). Importantly, the Riemannian methodology of [5] involves the exponential mapping of the underlying metric and focuses on geodesic convexity, so it concerns an orthogonal class of problems. The only overlap would be in the case of flat Riemannian manifolds: however, even though the manifolds we consider here are topologically simple, they are not flat. 1 In view of this, there is no overlap with the analysis and results of [5] .

2. Problem setup:

We begin by presenting the core online optimization framework that we will consider throughout the rest of our paper. This can be described by the following sequence of events:

1. At each round t = 1, 2, . . . , the optimizer chooses an action X t from a convex -but not necessarily closed or compact -subset X of an ambient normed space V d . 2. The optimizer incurs a loss f t (X t ) based on some (unknown) convex loss function f t : X → . 3. The optimizer updates their action and the process repeats.

Remark 1. For posterity, we note that if X is not closed, f t (or its derivatives) could become singular at a residual point x ∈ bd(X ) \ X ; in particular, we do not assume here that f t admits a smooth extension to the closure cl(X ) of X (or even that it is bounded over bounded subsets of X ).

In this broad framework, the most widely used figure of merit is the minimization of the agent's regret. Formally, the regret of a policy X t ∈ X , t = 1, 2, . . . , is defined as

EQUATION

for all x ∈ X . We then say that the policy X t leads to no regret if Reg x (T ) = o(T ) for all x ∈ X .

In addition to convexity, the standard assumption in the literature for the problem's loss functions is Lipschitz continuity, i.e.,

| f t (x ) -f t (x)| ≤ G t x -x (LC)

for some G t ≥ 0, t = 1, 2, . . . , and for all x, x ∈ X . Under (LC), if the agent observes at each stage t an element v t of ∂ f t (X t ), straightforward online policies based on gradient descent enjoy a bound of the form Reg x (T ) = O( ḠT T 1/2 ), with Ḡ2 11, 35, 41] . In particular, if G ≡ lim sup T →∞ ḠT < ∞ (e.g., if each f t is G-Lipschitz continuous over X ), we have the bound

T = T -1 T t=1 G 2 t [

Reg x (T ) = O(GT 1/2 ) (2)

which is well known to be min-max optimal in this setting [1] .

A note on notation. Throughout our paper, we make a clear distinction between V and its dual, and we use Dirac's notation v|x for the duality pairing between v ∈ V * and x ∈ V (not to be confused with the notation •, • for a scalar product on V). Also, unless mentioned otherwise, all notions of boundary and interior should be interpreted in the relative (as opposed to topological) sense. We also make the blanket assumption that the subdifferential ∂ f t of f t admits a continuous selection

∇ f t (x) ∈ ∂ f t (x) for all x ∈ dom ∂ f t ≡ {x ∈ X : ∂ f t (x) ∅}.

3 Riemann-Lipschitz continuity Despite its generality, (LC) may fail to hold in a wide range of problems and applications, ranging from support vector machines to Poisson inverse problems, quantum tomography, etc. [3, 8, 25] . The loss functions of these problems exhibit singularities at the boundary of the feasible region, so the standard regret analysis cited above no longer applies. Accordingly, our first step will be to introduce a family of local norms • x , x ∈ X , such that a variant of (LC) holds even if the derivatives ∂ f /∂x i of f blow up near the boundary of X .

To achieve this, we will employ the notion of a Riemannian metric. This is simply a positiondependent scalar product on V, i.e., a continuous assignment of bilinear pairings •, • x , x ∈ X , satisfying the following conditions for all z, z ∈ V and all x ∈ X :

1. Symmetry: z, z x = z , z x .

2. Positive-definiteness: z, z x ≥ 0 with equality if and only if z = 0.

More concretely, in the standard basis {e i } d i=1 of d , we define the metric tensor of •, • x as the matrix g(x) ∈ d×d with components g i j (x) = e i , e j x i, j = 1, . . . , d.

The norm of z ∈ V at x ∈ X is then defined as

EQUATION

In this way, a Riemannian metric allows us to measure lengths and angles between displacement vectors at each x ∈ X ; for illustration, we provide some notable examples below:

Example 1 (Euclidean geometry). The ordinary Euclidean metric on X = d is g(x) = I. This yields the standard expressions z 2 x = d i=1 z 2 i and z, z x = d i=1 z i z i , both independent of x. Example 2 (Hyperbolic geometry). The Poincaré metric on the positive orthant X = d ++ is

EQUATION

leading to the local norm z 2 x = d i=1 z 2 i /x 2 i . Under (5), d ++ can be seen as a variant of Poincaré's half-space model for hyperbolic geometry [22] ; this will become important later.

Given a Riemannian metric on X , the length of a curve γ : [0, 1] → X is defined as L g [γ] = 1 0 γ(s) γ(s) ds, and the Riemannian distance between x 1 , x 2 ∈ X is given by dist

EQUATION

Under this definition, it is natural to introduce the following Riemannian variant of (LC):

Definition 1. We say that f : X → is Riemann-Lipschitz continuous relative to g if

| f (x ) -f (x)| ≤ G dist g (x,

x ) for some G ≥ 0 and all x, x ∈ X .

Albeit simple to state, (RLC) may be difficult to verify because it requires the computation of the distance function dist g of g -which, in turn, relies on geodesic calculations to identify the shortest possible curve between two points. Nevertheless, if f is differentiable, Proposition 1 below provides an alternative characterization of Riemann-Lipschitz continuity which is easier to work with:

Proposition 1. Suppose that f : X → is differentiable. Then, (RLC) holds if and only if

grad f (x) x ≤ G for all x ∈ X . (RLC)

Remark 2. In the above, the Riemannian gradient grad f (x) of f at x is defined as follows: First, let Z = span{xx : x, x ∈ X } denote the tangent hull of X . Then, grad f (x) ∈ Z is defined by the characteristic property f (x; z) = grad f (x), z x for all z ∈ Z.

Existence and uniqueness of grad f (x) is due to the fact that g(x) is positive-definite -and, hence, invertible [22] . In particular, if X is full-dimensional (so Z = V), we have:

EQUATION

)

The proof of Proposition 1 requires the introduction of further tools from Riemannian geometry; seeing as these notions are not used anywhere else in our paper, we relegate it to the appendix. Instead, we close this section with a simple example of a singular function which is nonetheless Riemann-Lipschitz continuous:

EQUATION

Thus, although f is not Lipschitz continuous in the standard sense, it is Riemann-Lipschitz continuous relative to g; we will revisit this example in our treatment of Poisson inverse problems in Section 6.

More generally, Example 3 suggests the following rule of thumb: if f exhibits a gradient singularity of the form

|∂ i f (x)| = O(φ(x)) at some residual point x ∈ cl(X ) \ X of X , taking g i j (x) = φ(x) 2 δ i j gives grad f (x) 2 x = φ(x) -2 d i=1 [∂ i f (x)] 2 = O(1)

. On that account, f is Riemann-Lipschitz continuous, even though its derivative is singular; we find this heuristic particularly appealing because it provides a principled choice of Riemannian metric under which f satisfies (RLC).

4. Algorithms:

In this section, we present the algorithms that we will study in the sequel: "follow the regularized leader" (FTRL) and online mirror descent (OMD). Both methods have been widely studied in the literature in the context of "vanilla" Lipschitz continuity; however, beyond this basic setting, treating FTRL/OMD in the Riemannian framework of the previous section is an intricate affair that requires several conceptual modifications. For this reason, we take an in-depth look into both methods below. 2

4.1. Regularization:

We begin with the idea of regularization through a suitable penalty function. In our Riemannian setting, we adapt this notion as follows: Definition 2. Let g be a Riemannian metric on X and let h : V → be a proper lower semicontinuous (l.s.c.) convex function with dom h = X . 3 We say that h is a Riemannian regularizer on X if:

1. The subdifferential of h admits a continuous selection, i.e., a continuous function ∇h such that ∇h(x) ∈ ∂h(x) for all x ∈ X • ≡ dom ∂h.

2. h is strongly convex relative to g, i.e.,

EQUATION

)

for some K > 0 and all x ∈ X • , x ∈ X .

The Bregman divergence induced by h is then defined for all p ∈ X , x ∈ X • as

EQUATION

There are two points worth noting in the above definition. First, the domain of h is all of X , but this need not be the case for the subdifferential ∂h of h: by convex analysis arguments [33, Chap. 26], we have ri X ⊆ X • ≡ dom ∂h ⊆ X . To connect the two, we will say that h is a Riemann-Legendre regularizer when X • = ri X and D(p, x n ) → 0 whenever x n → p.

Second, strong convexity in (11) is defined relative to the underlying Riemannian metric. If the norm in (11) does not depend on x, we recover the standard definition; however, the dependence of the second-order term in (11) on g can change the landscape significantly. Lemma 1 and Example 5 below provide an illustration of this interplay between g and h: Lemma 1. A Riemannian regularizer h is K-strongly convex relative to g if and only if

EQUATION

The proof of Lemma 1 follows from a rearrangement of ( 11) so we omit it. Instead, we present below some examples of Riemannian regularizers:

Example 4. Let X = [0, 1] d , and consider the so-called Burg entropy h(x) = -d i=1 log x i . It is easy to see that h(x) is strongly convex relative to the standard Euclidean norm • 2 . Moreover, we have

EQUATION

and, by Taylor's theorem with Lagrange remainder, we readily get

EQUATION

)

where z x = d i=1 z 2 i /x i denotes the so-called Shahshahani norm on X (i.e., h is also strongly convex relative to • x ). This regularizer has been used extensively in the setting of Poisson inverse problems and plays a central role in the analysis of Bauschke et al. [3] , Hanzely and Richtárik [16] , He et al. [18] , and Lu et al. [25] . For completeness, we revisit it in Section 6.

Example 5. Let X and g be as in Example 3, and let h

(x) = (1 + r 2 )/ d i=1 x i with r 2 = d i=1 x 2 i .

Then, a tedious (but otherwise straightforward) algebraic calculation gives

EQUATION

i.e., h is strongly convex relative to g. By contrast, due to the singularity of g at 0, it is easy to check that the Euclidean regularizer h(x) = (1/2) d i=1 x 2 i is not strongly convex relative to g.

4.2. Algorithms and feedback structure:

With these preliminaries in hand, we begin with the FTRL algorithm, which we state here as follows:

X t+1 = arg min x∈X t s=1 f s (x) + γ -1 h(x) . (FTRL)

In the above, γ > 0 is a step-size parameter whose role is discussed below; as for the existence of the arg min, this is justified by the lower semicontinuity and strong convexity of h together with the fact that dom h = X (so the minimum cannot be attained in the residual set cl(X ) \ X of X ).

In terms of feedback, FTRL assumes that the optimizer has access to all the loss functions encountered up to a given round (except, of course, for the current one). In many cases of practical interest, this assumption is too restrictive and, instead, the optimizer only has access to a first-order oracle for each f t . To model this feedback structure, we assume that once X t has been chosen, the optimizer receives an estimate v t of ∇ f t (X t ) satisfying the following hypotheses:

EQUATION

In the above, • * denotes the dual norm of • X t , i.e., the Riemannian norm at the point X t where the oracle was called (we suppress here the index X t and write • * instead of • X t , * to lighten the notation). In particular, the oracle feedback v t may fail to be bounded in L 2 relative to any global norm on V * ; as such, (17) is considerably weaker than the standard L 2 -boundedness assumption for global norms. Finally, in terms of measurability, the expectation in ( 17) is conditioned on the history F t of X t up to stage t; since v t is generated randomly from X t , it is not F t -measurable.

To proceed, the main idea of mirror descent is to replace f s (x) in (FTRL) with the first-order surrogate

f s (x) ← f s (X s ) + ∇ f s (X s )|x -X s .

In this way, substituting v s for the estimate of ∇ f s (X s ) received at stage s, we obtain the linearized FTRL scheme

EQUATION

)

To rewrite this process in recursive form, introduce the auxiliary (dual) variable

EQUATION

so Y t+1 = -γ t s=1 v s , and hence

EQUATION

Therefore, letting

EQUATION

denote the so-called "mirror map" of the method, we obtain the following incarnation of the online mirror descent (OMD) algorithm:

Y t+1 = Y t -γv t X t+1 = Q(Y t+1 ). (OMD)

This version of OMD is also known as "dual averaging" [26, 29, 30, 38] or "lazy mirror descent" [35] ; for a "greedy" variant, see [6, 27, 28] and references therein.

5 Analysis and results

5.1. Regret analysis:

We begin by stating our main results for the regret minimization properties of FTRL and OMD. Throughout this section, we make the following blanket assumptions:

1. Both algorithms are initialized at the "prox-center" x c = arg min h of X and are run with (constant) step-size α/T 1/2 for some α > 0 chosen by the optimizer.

2. The t-th stage loss function f t : X → is convex and satisfies (RLC) with constant G t .

3. The optimizer's aggregate loss T t=1 f t attains its minimum value at some x * ∈ X .

The purpose of the last assumption is to avoid cases where the infimum of a loss function is not attained within the problem's feasible region (such as e -x over + ). We then have:

Theorem 1. Let Reg(T ) ≡ Reg x * (T ), Ḡ2 T = T -1 T t=1 G 2 t , and M2 T = T -1 T t=1 M 2 t . Then:

a) The FTRL algorithm enjoys the regret bound

EQUATION

b) The OMD algorithm with noisy feedback of the form (17) enjoys the mean regret bound

EQUATION

In particular, if

sup t G t < ∞, sup t M t < ∞, both algorithms guarantee O( √ T ) regret.

Remark 3. We emphasize here that the O( 4 Of course, if H < ∞ and G (or M) is known to the optimizer, ( 22) can be optimized further by tuning α. When these constants are unknown, achieving an optimized constant by means of an adaptive step-size policy is an important question, but one which lies beyond the scope of this paper.

√ T ) regret bound above is achieved even if X is unbounded or if the "Bregman depth" H ≡ sup x∈X D(x, x c ) = sup h -inf h of X is infinite.

The main idea behind the proof of Theorem 1 is to relate the Riemannian structure of X to the Bregman regularization framework underlying (FTRL) and (OMD). A first such link is provided by the Bregman divergence (12) ; however, because of the primal-dual interplay between X t ∈ X and Y t ∈ V * , the Bregman divergence is not sufficiently adapted. To overcome this difficulty, we employ the Fenchel coupling between a target point p ∈ X and y ∈ V * , defined here as

EQUATION

with h * (y) = max x∈X { y|xh(x)} denoting the convex conjugate of h. As we show in the appendix, the Fenchel coupling (which is non-negative by virtue of Young's inequality) enjoys the key property

EQUATION

It is precisely this primal-dual inequality which allows us to go beyond the standard Lipschitz framework: compared to (primal-primal) inequalities of a similar form for global norms [2, 21, 27, 30, 40] , the distinguishing feature of ( 24) is the advent of the Riemannian norm v x, * . Thanks to the intricate connection between g and h, the second-order term in (24) can be controlled even when the received gradient is unbounded relative to any global norm, i.e., even if the objective is singular.

The main obstacle to achieve this is that the underlying Riemannian metric g, the Fenchel coupling Φ and the Bregman divergence D (all state-dependent notions of distance) need not be compatible with one another. That this is indeed the case is owed to Lemma 1: tethering the Riemannian norm in (13) to the second argument of the Bregman divergence instead of the first (or any other point in-between) plays a crucial role in deriving (24) . Any other relation between g and h along these lines is not amenable to analyzing (FTRL) or (OMD) in this framework.

5.2. Applications to stochastic optimization:

The second part of our analysis concerns stochastic optimization problems of the form

minimize f (x) = ¾[F(x; ω)] subject to x ∈ X (Opt)

with the expectation taken over some model sample space Ω. Our first result here is as follows:

Theorem 2. Assume that f is convex and Riemann-Lipschitz continuous in mean square, i.e.,

sup x ¾[ ∇F(x; ω) 2 x, * ] ≤ M 2 for some M > 0. If (OMD)

is run for T iterations with a constant step-size of the form α/

√:

T and stochastic gradients v t = ∇F(X t ; ω t ) generated by an i.i.d. sequence ω t ∈ Ω, we have

EQUATION

)

where XT = (1/T ) T t=1 X t is the "ergodic average" of X t and D c = inf x * ∈arg min f D(x * , x c ) < ∞ denotes the Bregman distance of the prox-center x c of X to arg min f .

The key novelty in Theorem 2 is that the optimal O(T -1/2 ) convergence rate of OMD is maintained even if the stochastic gradients of F become singular at residual points x ∈ cl(X ) \ X . As with the regret guarantee of Theorem 1, this is achieved by the intricate three-way relation between the landscape of f , the underlying Riemannian metric g (which is tailored to the singularity profile of the latter), and the Riemannian regularizer h. The proof of Theorem 2 likewise relies on an online-to-batch conversion of the regret guarantees of (OMD) for the sequence of stochastic gradients ∇F(•; ω t ) of f ; the details can be found in the appendix.

To go beyond the ergodic guarantees of Theorem 2, we also analyze below the convergence of the "last iterate" of OMD, i.e., the actual sequence of generated points X t . This is of particular interest for non-convex problems where ergodic convergence results are of limited value (because Jensen's inequality no longer applies). To obtain global convergence results in this setting, we focus on a class of functions which satisfy a weak secant inequality of the form

inf{ ∇ f (x)|x -x * : x * ∈ arg min f, x ∈ K} > 0 (SI)

for every closed subset K of X that is separated by neighborhoods from arg min f . Variants of this condition have been widely studied in the literature and include non-convex functions with complicated ridge structures [9, 13, 19, 20, 23, 31, 39, 40] . In this very general setting, we have: Theorem 3. Assume f satisfies (SI) and is Riemann-Lipschitz continuous in L 2 . Suppose further that arg min f is bounded and (OMD) is run with a sequence of stochastic gradients v t = ∇F(X t ; ω t ), a Riemann-Legendre regularizer h, and a variable step-size γ t such that ∞ t=1 γ t = ∞, ∞ t=1 γ 2 t < ∞. Then, with probability 1, X t converges to some (possibly random) x * ∈ arg min f . The proof for Theorem 3 hinges on combining (quasi-)supermartingale convergence results with the basic inequality (24); we detail the proof in the paper's appendix. Seeing as (SI) holds trivially for (pseudo-)convex functions, we only note here that Theorem 3 complements Theorem 2 in an important way: the convergence of X t implies that of Xt , so the convergence of ¾[ f ( Xt )] to min f follows immediately from Theorem 3; however, the rate of convergence (25) doesn't. In practice, the ergodic average converges to interior minimizers faster than the last iterate but lags behind when tracking boundary points and/or in non-convex landscapes; we explore this issue in Section 6 below.

������� ���� ��� �������������� ��� ��������������                                   • • • • • • ••• • •• • •• • • • • • • • • • • • • • • • • • • • • • • • ◼ ◼ ◼ ◼ ◼ ◼◼◼◼ ◼◼◼ ◼ ◼◼ ◼◼ ◼◼◼ ◼ ◼ ◼◼ ◼ ◼◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◼ ◻ ◻ ◻ ◻ ◻ ◻◻◻ ◻ ◻◻◻ ◻ ◻◻ ◻◻ ◻◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◻ ◆ ◆ ◆ ◆ ◆◆◆ ◆ ◆◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆◆ ◆ ◇ ◇ ◇ ◇ ◇ ◇◇ ◇◇ ◇ ◇ ◇◇ ◇ ◇◇ ◇ ◇◇ ◇◇ ◇◇ ◇ ◇◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇ ◇

6. Numerical experiments in Poisson inverse problems:

For the purposes of validation, we proceed with an application of our algorithmic results to a broad class of Poisson inverse problems that arise in tomography problems. Referring the reader to the appendix for the details, the objective of interest here is the Poisson likelihood loss (generalized Kullback-Leibler divergence):

EQUATION

where u ∈ m + is a vector of Poisson data observations (e.g., pixel intensities) and H ∈ m×d is an ill-conditioned matrix representing the data-gathering protocol. Since the generalized KL objective of (26) exhibits an O(1/x) singularity at the boundary of the orthant, we consider the Poincaré metric g(x) = diag(1/x 1 , . . . , 1/x d ) under which the KL divergence is Riemann-Lipschitz continuous (cf. Example 2). Going back to Example 3, a suitable Riemannian regularizer for this metric is h(x) = m i=1 1/x 2 i , which is 1-strongly convex relative to g. We then run the induced mirror descent algorithm with an online-to-batch conversion mechanism as described in Section 5.2. For reference purposes, we call the resulting process Riemannian mirror descent (RMD).

Subsequently, we ran RMD on a Poisson denoising problem for a 384 × 384 test image contaminated with Poisson noise (so d ≈ 10 5 in this case). For benchmarking, we also ran a fast variant of the widely used Lucy-Richardson (LR) algorithm [7] , and the recent composite mirror prox (CMP) method of [18] ; all methods were run with stochastic gradients and the same minibatch size. Because of the "dark area" gradient singularities when [Hx] j → 0, Euclidean stochastic gradient methods oscillate without converging, so they are not reported. As we see in Fig. 1 , the RMD process provides the sharpest reconstruction of the original. In particular, after an initial warm-up phase, the last iterate of Riemannian mirror descent consistently outperforms the LR algorithm by 7 orders of magnitude, and CMP by 3. We also note that the Poisson likelihood loss decreases faster under the last iterate of RMD relative to the different algorithmic variants that we tested, exactly because of the hysteresis effect that is inherent to ergodic averaging.

Overall, we note that the introduction of an additional degree of freedom (the choice of Bregman function and that of the local Riemannian norm), makes RMD a particularly flexible and powerful paradigm for loss models with singularities. We find these results particularly encouraging for further investigations on the interplay between Riemannian geometry and Bregman-proximal methods.

7. Concluding remarks:

Owing to its connections with machine learning (support vector machines, Poisson inverse problems, quantum tomography, etc.), venturing beyond Lipschitz continuity is a fruitful research direction that has recently generated considerable interest in the literature. Depending on the type of continuity or smoothness encountered (Lipschitz continuity of the objective or Lipschitz continuity of the objective's gradients), the results can be significantly different, and it is not a priori clear which surrogate smoothness/continuity condition would be the most appropriate for any given problem. The present paper provides a complementary, Riemannian-geometric viewpoint which we feel can be fairly promising for the design of efficient optimization algorithms in this context. The precise characterization of the interplay between the different continuity conditions considered in the literature is an important open issue which we leave for future work.

A Riemann-Lipschitz continuity:

In this appendix, our main goal is to prove Proposition 1, i.e., the equivalence between ( 7) and (RLC) when f is differentiable.

To do so, we first need to introduce the notion of a geodesic, i.e., a length-minimizing curve that attains the infimum inf γ L[γ] over all piecewise smooth curves joining two points x 1 , x 2 ∈ U. That such a curve exists and is unique in our setting is a basic fact of Riemannian geometry [22] . Moreover, given a tangent vector z ∈ Z, this leads to the definition of the exponential mapping exp :

U × Z → U so that (x, z) → exp x (z) = γ z (1), (A.1)

where γ z denotes the unique geodesic emanating from x with initial velocity vector γz = z. We then have exp x (tz) = γ z (t) for all t and, moreover, for sufficiently small r > 0, the restriction of exp x to a ball of radius r in Z is a diffeomorphism onto its image in U. The largest positive number i x such that the above holds for all r < i x is then known as the injectivity radius of U at x [22] .

Our proof of the equivalence between ( 7) and (RLC) follows a simplified version of the approach of Ferreira [14] who, to our knowledge, was the first to discuss the concept of proximal subgradients in Riemannian manifolds. To that end, fix some x ∈ U, let z = grad f (x) and consider the ray emanating from x

γ(t) = exp x (tz/ z x ). (A.2)

Since γ is a geodesic, we readily obtain dist g (x, γ(t)) = t for all sufficiently small t. Also, by construction, we have exp -1 x γ(t) = tz/ z x . Hence, with f convex, it follows that, for some constant α > 0 and for sufficiently small positive δ < i x , we have

f (γ(t)) -f (x) ≥ z, exp -1 x (γ(t)) (A.3) -α dist g (x, γ(t)) 2 , (A.4)

where we used the local topological equivalence of the Riemannian topology and the standard topology of d , and the fact that exp z (t) is a diffeomorphism for sufficiently small δ > 0 -and hence, for all δ < i x . Thus, if f is also Riemann-Lipschitz continuous in the sense of (7), we will also have

Gt ≥ f (γ(t)) -f (x) ≥ z, tz/ z x x -αt 2 . (A.5)

Thus, by isolating the leftmost and rightmost hand sides, dividing by t, and taking the limit t → 0, we get grad f (x

EQUATION

as was to be shown.

To establish the converse, assume that (RLC) holds, fix x, x ∈ X pick K > G and a sufficiently small δ > 0, and consider the Riemannian distance majorant w(x) = dist g (x, x ) if dist g (x, x ) < δ, and

w(x) = K dist g (x, x ) + ε 2 /(δ -ε) when δ < dist g (x, x ) < 2δ, with ε = dist g (x, x ) -δ.

A simple calculation then shows that grad w(x) x ≥ K > G. It is also straightforward to show that the minimum of f + w is attained at x , so

f (x ) = f (x ) + w(x ) ≤ f (x) + w(x) ≤ f (x) + K dist g (x, x ). (A.7)

Then, by interchanging x and x above, we obtain

EQUATION

B Properties of mirror mappings and the Fenchel coupling:

We begin by recalling and clarifying some of the notational conventions used in the paper. First, let V d be a finite-dimensional real space; then, its dual space will be denoted by Y ≡ V * , and we write y|x for the duality pairing between y ∈ Y and x ∈ V. Also, if • is a norm on V, the dual norm on Y is defined as y * ≡ sup{ y|x : x ≤ 1}.

Given an extended-real-valued convex function f : X → ∪ {∞}, we will write dom f ≡ {x ∈ V : f (x) < ∞} for its effective domain. The subdifferential of f at x ∈ dom f is then defined as ∂ f (x) ≡ {y ∈ Y : f (x )f (x) + y|xx ≥ 0 for all x ∈ V} and the domain of subdifferentiability of f is dom ∂ f ≡ {x ∈ dom f : ∂ f ∅}. Finally, assuming it exists, the directional derivative of f at x along z ∈ V is defined as f (x; z) ≡ d/dt| t=0 f (x + tz). We will then say that f is differentiable at x ∈ dom f if there exists ∇ f (x) ∈ Y such that ∇ f (x)|z = f (x; z) for all vectors of the form z = xx, x ∈ dom f . With these notational conventions at hand, we proceed to prove some auxiliary results and estimates that are used throughout the analysis of Section 5. To recall the basic setup, we assume throughout what follows that h is a Riemannian regularizer in the sense of Definition 2. The convex conjugate h * : Y → of h is then defined as

h * (y) = sup x∈X { y|x -h(x)}. (B.1)

Since h is K-strongly convex relative to g, it is also strongly convex relative to the Euclidean norm (recall here that g(x) µI). As a result, the supremum in (B.1) is always attained, and h * (y) is finite for all y ∈ Y [4] . Moreover, by standard results in convex analysis [33, Chap. 26] , h * is differentiable on Y and its gradient satisfies the identity

∇h * (y) = arg max x∈X { y|x -h(x)}. (B.2)

Thus, recalling the definition of the mirror map Q : Y → X (cf.. Section 4):

EQUATION

Together with the prox-mapping induced by h, all these notions are related as follows: Lemma B.1. Let h be a Riemannian regularizer on X . Then, for all x ∈ dom ∂h and all y, v ∈ Y, we have:

a) x = Q(y) ⇐⇒ y ∈ ∂h(x). (B.5a) b) x + = Q(∇h(x) + v) ⇐⇒ ∇h(x) + v ∈ ∂h(x + ) (B.5b)

Finally, if x = Q(y) and p ∈ X , we have

∇h(x)|x -p ≤ y|x -p . (B.6)

Remark. Note that (B.5b) directly implies that ∂h(x + ) ∅, i.e., x + ∈ dom ∂h for all v ∈ Y. An immediate consequence of this is that the update rule x + = Q(∇h(x) + v) is well-posed, i.e., it can be iterated in perpetuity.

Proof of Lemma B.1. To prove (B.5a), note that x solves (B.2) if and only if y -∂h(x) 0, i.e., if and only if y ∈ ∂h(x). Eq. (B.5b) is then obtained in the same manner.

For the inequality (B.6), it suffices to show it holds for all p ∈ X • ≡ dom ∂h (by continuity). To do so, let

φ(t) = h(x + t(p -x)) -[h(x) + y|x + t(p -x) ].

(B.7) Since h is strongly convex relative to g and y ∈ ∂h(x) by (B.5a), it follows that φ(t) ≥ 0 with equality if and only if t = 0. Moreover, note that ψ(t) = ∇h(x + t(px)) -y|px is a continuous selection of subgradients of φ. Given that φ and ψ are both continuous on [0, 1], it follows that φ is continuously differentiable and φ = ψ on [0, 1]. Thus, with φ convex and φ(t) ≥ 0 = φ(0) for all t ∈ [0, 1], we conclude that φ (0) = ∇h(x) -y|px ≥ 0, from which our claim follows.

As we mentioned earlier, much of our analysis revolves around a "primal-dual" divergence between a target point p ∈ X and a dual vector y ∈ Y, called the Fenchel coupling. Following [26] , this is defined as follows for all p ∈ X , y ∈ Y:

Φ(p, y) = h(p) + h * (y) -y|p . (B.8)

The following lemma illustrates some basic properties of the Fenchel coupling:

Lemma B.2. Let h be a Riemannian regularizer on X with convexity modulus K. Then, for all p ∈ X and all y ∈ Y, we have:

1. Φ(p, y) = D(p, Q(y)) if Q(y) ∈ X • (but not necessarily otherwise). 2. If x = Q(y), then Φ(p, y) ≥ K 2 x -p 2 x

Proof. For our first claim, let x = Q(y). Then, by definition we have:

Φ(p, y) = h(p) -y|Q(y) -h(Q(y)) -y|p = h(p) -h(x) -y|p -x . (B.9)

Since y ∈ ∂h(x), we have h (x; px) = y|px whenever x ∈ X • , thus proving our first claim. For our second claim, working in the previous spirit we get that:

Φ(p, y) = h(p) -h(x) -y|p -x (B.10)

Thus, we obtain the result by recalling the strong convexity assumption for h with respect to the Riemannian norm • x .

We continue with some basic relations connecting the Fenchel coupling relative to a target point before and after a gradient step. The basic ingredient for this is a primal-dual analogue of the so-called "three-point identity" for Bregman functions [12] : Lemma B.3. Let h be a regularizer on X . Fix some p ∈ X and let y, y + ∈ Y. Then, letting x = Q(y), we have

Φ(p, y + ) = Φ(p, y) + Φ(x, y + ) + y + -y|x -p . (B.11)

Proof. By definition, we get: With all this at hand, we have the following key estimate: Proposition B.1. Let h be a Riemannian regularizer on X with convexity modulus K, fix some p ∈ X , let x = Q(y) for some y ∈ Y. Then, for all v ∈ Y, we have: and hence, after rearranging:

Φ(p, y + ) = h(p) + h * (y + ) -y + |p Φ(p, y) = h(p) + h * (y) -y|p . (B.

EQUATION

Φ(p, y + v) = Φ(p, y) -Φ(Q(y + v), y) + v|Q(y + v) -p = Φ(p, y) -Φ(Q(y + v), y) + v|x -p + v|Q(y + v) -x (B.16)

By Young's inequality [33] , we also have

EQUATION

x, * (B.17

)

Our claim then follows by the fact that Φ(Q(y + v), y)

≥ K 2 Q(y + v) -x 2

x (cf. Lemmas 1 and B.2).

C Analysis of FTRL:

Our goal here is to prove the regret bound (22a) of Theorem 1. The starting point of our analysis is the following basic bound: Lemma C.1 (35, Lemma 2.3). The sequence of actions generated by (FTRL) satisfies

Reg x (T ) ≤ h(x) -h(X 1 ) + T t=1 [ f t (X t ) -f t (X t+1 )]. (C.1)

Importantly, the above bound does not require any Lipschitz continuity or strong convexity assumptions, so it applies to our setting "as is". The importance of Riemann-Lipschitz continuity lies in the following: Lemma C.2. If f is convex and Riemann-Lipschitz continuous with constant G, then:

f (x) -f (x ) ≤ G x -x x for all x, x ∈ X • . (C.2)

Proof. By the convexity of f , we have:

f (x) -f (x ) ≤ ∇ f (x)|x -x = grad f (x), x -x x ≤ grad f (x) x x -x x ≤ G x -x x (C.3)

where the first line follows from the definition of the Riemannian gradient of f , the second one from the Cauchy-Schwartz inequality, and the last from Proposition 1.

With these preliminary results at hand, we obtain the following basic bound for FTRL: Proposition C.1. Suppose that (FTRL) is run against a sequence of loss convex loss functions with assumptions as in Section 5. Then, for all x ∈ X , we have:

Reg x (T ) ≤ D(x, x c ) γ + 2γ K T t=1 G 2 t (C.4)

Proof. Our proof is patterned after Shalev-Shwartz [35] , but with an important difference regarding the use of Riemannian norms and Riemann-Lipschitz continuity. To begin, let

F t (x) = t-1 s=1 f s (x) + h(x) γ (C.5)

denote the "cumulative" loss faced by the optimizer up to roun t -1, including the regularization penalty. By the definition of the FTRL policy, X t ∈ arg min F t (x), so ∇F t (X t )|x -X t ≥ 0 and, likewise, ∇F t+1 (X t+1 )|x -X t+1 ≥ 0 for all x ∈ X . Furthermore, since h is K-strongly convex relative to g, F t and F t+1 will be (K/γ)-strongly convex relative to g.

Putting all this together, we obtain:

F t (X t+1 ) ≥ F t (X t ) + K 2γ X t+1 -X t 2 X t+1 (C.6a) F t+1 (X t ) ≥ F t+1 (X t+1 ) + K 2γ X t -X t+1 2 X t (C.6b)

and hence, after summing the above inequalities:

f t (X t ) -f t (X t+1 ) ≥ K 2γ X t+1 -X t 2 X t+1 + K 2γ X t -X t+1 2 X t ≥ K 2γ X t+1 -X t 2 X t . (C.7)

On the other hand, Lemma C.2 gives

f t (X t ) -f t (X t+1 ) ≤ G t X t -X t+1 X t (C.8)

so, combining the last two inequalities, we get:

K 2γ X t+1 -X t X t ≤ G t . (C.9)

Therefore, plugging this back into (C.8) yields

f t (X t ) -f t (X t+1 ) ≤ 2γG 2 t K , (C.10)

and our result obtains from Lemma C.1.

The proof of (22a) then follows by applying Proposition C.1 with a step-size of the prescribed form.

D Ergodic analysis of OMD:

This appendix is devoted to the proof of our main regret bound for (OMD). We begin by recalling ther recursive definition of the (lazy) OMD method:

Y t+1 = Y t -γv t X t+1 = Q(Y t+1 ) (OMD)

with Q defined as in Appendix B and oracle feedback subject to the hypotheses (17) . We may then write the oracle feedback received by the optimizer at time t as v t = ∇ f t (X t ) + U t+1 ; hence, by the unbiasedness assumption (17a), it follows that ¾[U t+1 | F t ] = 0, i.e., U t is a martingale difference sequence (MDS) relative to F t . Now, applying Proposition B.1 to (OMD), we get:

Φ(x * , Y t+1 ) ≤ Φ(x * , Y t ) -γ v t |X t -x * + γ 2 2K v t 2 X t , * = Φ(x * , Y t ) + γ ∇ f t (X t )|x * -X t -γ U t+1 |X t -x * + γ 2 2K v t 2 X t , * . (D.1)

Hence, after rearranging and telescoping, we obtain

Reg(T ) ≤ T t=1 ∇ f t (X t )|X t -x * ≤ D(x * , x c ) γ + T t=1 ξ t+1 + γ 2K T t=1 v t 2 X t , * (D.2)

where, in the last line, we used the definition of the Riemannian dual norm • * ≡ • x * , * , and we set ξ t+1 = U t+1 |x * -X t . Our result then follows by taking expectations on both sides.

E Last-iterate analysis of OMD:

In this last section, we will present the convergence analysis for the last iterate of (OMD) to arg min f . We begin by recalling two important results from probability theory. The first is a version of the law of large numbers for martingale difference sequences that are bounded in L 2 [15] : Theorem E. 

on the set ∞ t=1 β -2 t ¾[ζ 2 t | F t-1 ] < ∞.

The second is a convergence result for quasi-supermartingales due to Robbins and Sigmund [32] : Lemma E.1. Let (F t ) t∈ be a non-decreasing sequence of σalgebras. Let (α t ) t∈ , (θ t ) t∈ nonnegative F tmeasurable random variables, (η t ) t∈ is an F tmeasurable non-negative summable random variable and the following inequality holds:

¾[α t+1 | F t ] ≤ α t -θ t + η t almost surely (E.2)

Then, (α t ) t∈ converges almost surely towards a [0, ∞)-valued random variable.

By the unbiasedness assumption for U t , we have ¾

[ξ t+1 | F t ] = ¾[U t+1 | F t ]|X t -x * = 0. Moreover, for all x * ∈ arg min f , we have ∞ t=1 γ 2 t ¾[ξ t+1 | F t ] ≤ ∞ t=1 γ 2 t X t -x * 2 X t ¾[U t+1 | F t ] ≤ ∞ t=1 γ 2 t Φ(x * , Y t ) ¾[U t+1 | F t ] < ∞ (E.9)

where the last (strict) inequality is obtained due to the finite mean square property, the boundness of Φ(x * , Y t ) and the fact that ∞ t=1 γ 2 t < ∞. Thus, we can apply the law of large numbers for L 2martingales stated above and conclude that β -1 t t s=1 γ s ξ s+1 converges to 0 almost surely. On the other hand, for the term S t+1 = t s=1 γ 2 s v s 2 X t , * , since v s+1 is F s -measurable for all s = 1, 2 . . . , t -1 we have:

¾[S t+1 | F t ] = ¾         t-1 i=1 γ 2 t v i 2 x i , * + γ 2 t v t 2 X t , * F t         = S t + γ 2 t ¾ v t 2 X t , * F t ≥ S t (E.10)

so S t is a submartingale with respect to F t . Furthermore, by the law of total expectation, we also get:

¾[S t+1 ] = ¾[¾[S t+1 | F t ]] ≤ σ 2 t i=1 γ 2 i ≤ σ 2 ∞ t=1 γ 2 t < ∞, (E.11)

implying that S t is bounded in L 1 . Thus, due to Doob's submartingale convergence theorem [15] , we coclude that S t converges to some (almost surely finite) random variable S ∞ so lim t→∞ S t+1

β t = 0 with probability 1.

Hence, obtaining a maximum likelihood estimate for x leads to the archetypal Poisson inverse problem minimize f (x) ≡ D KL (u, Hx),

EQUATION

where D KL (p, q) = m j=1 [p j log(p j /q j ) + q jp j ] denotes the generalized KL divergence on m + . For an extensive review of Poisson inverse problems, we refer the reader to Bertero et al. [7] .

In many cases of practical interest, measurements arrive in distinct batches over time -e.g., as sequential optical sections in microscopy and tomography. Moreover, due to the large numbers of pixels/voxels involved (a typical range of values for m is between 10 6 and 10 7 ), gradients of f are very costly to compute; as such, optimization methods that rely on accurate gradient data are difficult to apply in this setting. Accordingly, a natural workaround to this obstacle is to exploit the online nature of the measurement process, model (PIP) as an online optimization problem, and then to use an online-to-batch conversion to get a candidate solution [35] .

On the downside, this online optimization analysis crucially requires the loss functions faced by the optimizer to be Lipschitz continuous. However, this assumption does not hold for (PIP): if f j (x) = -u j log(u j /(Hx) j ) denotes the singular part of the KL divergence for the j-th sample, we readily get so ∇ f j x is bounded under this modified norm. This is the principal motivation for considering the Riemannian mirror descent method defined with respect to this metric and the regularizer presented in Example 5.

EQUATION

F.2 Details on the experiments:

In the rest of this appendix, we discuss in more detail the algorithms tested in Section 6. The algorithms we considered are 1. The accelerated Lucy-Richardson algorithm, as presented in [7] and corresponding to OMD with the entropic regularizer h(x) = i x i log x i . 2. The composite mirror prox (CMP) of He et al. [18] , corresponding to OMD with an extra gradient step and the log-barrier (Burg) regularizer h(x) =i log x i of Example 4. 3. The Riemannian mirror descent (RMD) algorithm detailed in Section 6, corresponding to the Poincaré-like regularizer of Example 5.

All algorithms were run with stochastic gradients drawn with the same minibatch size (n = 256) and a step-size of the form γ t ∝ 1/ √ t (corresponding to the stochastic variant of each algorithm). For comparison purposes, we harvested each algorithm's last generated sample ("last iterate") as well as the corresponding ergodic average (defined here as XT = T t=1 γ t X t T t=1 γ t ). Overall, the algorithms' last generated sample provided consistently better results than the ergodic average. The ground truth and the evolution of the Poisson likelihood loss are all reported in Fig. 1 . Remark. We should note here that the method of He et al. [18] can be seen as an "extra-gradient" version of the NoLips algorithm of Bauschke et al. [3] and the "relative stochastic gradient descent" scheme of Hanzely and Richtárik [16] (the difference between the last two being the step-size policy). In our experiments, Burg mirror descent with and without an extra-gradient step behaved similarly, with the extra-gradient version (CMP) performing slightly better. To minimize clutter, and because we are already comparing RMD to CMP above, we do not report this extra set of numerical experiments.

:

For example, the open unit simplex endowed with the Shahshahani metric is isometric to the positive orthant of a sphere with the round metric (cf. Section 3). This space has constant positive curvature and the geodesics are portions of great circles, so the two analyses are very different in that case.

For convenience, we tacitly assume in what follows that g(x) µI for some µ > 0 and all x ∈ X . Since g 0, this can always be achieved without loss of generality by replacing g by g + µI.

Following standard convex analysis terminology, "proper" means here that h +∞ while lower semicontinuous refers to the property that lim inf x→x 0 h(x) ≥ h(x 0 ) for all x 0 ∈ V.

To see this, simply note that D(x, x c ) = h(x)h(x c ) -∇h(x c )|xx c < ∞ for all x ∈ X = dom h (recall also that, since x c = arg min h, we have 0 ∈ ∂h(x c ) so x c ∈ dom ∂h).

Figures and Tables:

Figure 1: Reconstruction of the Lena test image from a sample contaminated with Poisson noise. Left to right: (a) the contaminated sample; (b) CMP reconstruction; (c) RMD reconstruction; and (d) Poisson likelihood loss at each iteration. The RMD process provides a sharper definition of image features relative to the CMP algorithm (which is the second-best).

Then, by subtracting the above we get:Φ(p, y + ) -Φ(p, y) = h(p) + h * (y + ) -y + |ph(p)h * (y) + y|p = h * (y + )h * (y) -y + -y|p = h * (y + ) -y|Q(y) + h(Q(y)) -y + -y|p = h * (y + ) -y|x + h(x) -y + -y|p = h * (y + ) + y + -y|x -y + |x + h(x) -y + -y|p = Φ(x, y + ) + y + -y|xp (B.13)and our proof is complete.

Proof. By the three-point identity (B.11), we getΦ(p, y) = Φ(p, y + v) + Φ(Q(y + v), y) + y -(y + v)|Q(y + v)p (B.15)

Hx) j .(F.3) This shows that the gradient of f j exhibits an O(1/x) singularity at the boundary of d + , so f cannot be Lipschitz under any global norm on d .As suggested by Example 3, this singularity can be lifted by considering the local normz 2 x = (x 1 + • • • + x d ) 2 d i=1 zfor all v ∈ d . (F.4) In this case, we have ∇ f j (x

1. Let Y t = t i=1 ζ i be a martingale and β t a non-decreasing positive sequence such that lim t→∞ β t = ∞. Then, lim

Table content: <table><tr><td>t→∞</td><td>Y t /β t = 0 almost surely</td><td>(E.1)</td></tr></table>

Peer review (iclr2020):
The paper establishes optimal regret bounds of the order O(sqrt{T}) for Follow The Regularised Leader (FTRL) and Online Mirror Descent (OMD) for convex loss functions and potentials (a.k.a. Riemannian regularizers) that are, respectively, Lipschitz continuous and strongly convex with respect to a given Riemannian metric. These conditions naturally generalize the classical conditions typically considered in the literature, which are defined with respect to a global norm and, as such, are not well-suited to problems where the loss functions and its gradient present singularities at the boundary of the feasibility region. The authors suggest a principled way to choose both the Riemannian metric and the potential function based on the singularity landscape of the gradient of the loss function. Via standard online-to-batch conversion, the authors also address the offline setting and give O(1/sqrt{T}) error bounds for ergodic averages in convex problems and for last iterates in non-convex problems satisfying a weak secant inequality. The authors include numerical experiments involving a Poisson inverse problem. The paper is well-written, with a very clean narrative highlighting the main ideas and results. To the best of my knowledge, the literature review is complete and rightly highlights the fact that most results in the literature on Riemannian mirror descent methods have so far primarily addressed offline deterministic problems with exact oracle gradients. The contribution of this work lies not only in the focus on online and noisy setting but also on establishing natural results upon natural generalizations of well-known conditions in standard (non-Riemannian) settings. The techniques used are extensions of the classical theory and follow quite naturally, with the exception of the non-trivial primal-dual inequality (20). QUESTIONS/SUGGESTIONS: 1) Can the authors be more explicit about how the (OMD) equations are derived from (16). While this is standard, I feel currently there is a bit of a jump in the narrative—which otherwise is very good. 2) The workhorse behind the established results is the primal-dual inequality (20) which relies on the introduction of the Fenchel coupling. Can the authors be more explicit about the use of this inequality, and what makes the Riemann generalization difficult in general? In particular, can the authors comment on the applicability of this inequality (or similar) to the smooth setting? 3) Typo: sometimes the notation seems to be used instead of the notation . See, for instance, equation (7) in Definition 1 and the definition of the set in Remark 2. AFTER REBUTTAL: I thank the authors for addressing my questions.