Online and Stochastic Optimization beyond Lipschitz Continuity: A Riemannian Approach
Year: 2020

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.

text: 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 ft, 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 ( ft = 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 4 β∇2h 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 infx′ √ 2D(x′, x)/‖x′ − x‖ (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 Xt from a convex – but not necessarily closed or compact – subset X of an ambient normed space V d.
2. The optimizer incurs a loss ft(Xt) based on some (unknown) convex loss function ft : X → . 3. The optimizer updates their action and the process repeats.
Remark 1. For posterity, we note that if X is not closed, ft (or its derivatives) could become singular at a residual point x ∈ bd(X ) \ X ; in particular, we do not assume here that ft 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 Xt ∈ X , t = 1, 2, . . . , is defined as
Regx(T ) = T∑
t=1
[ ft(Xt) − ft(x)], (1)
for all x ∈ X . We then say that the policy Xt leads to no regret if Regx(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., | ft(x′) − ft(x)| ≤ Gt‖x′ − x‖ (LC) for some Gt ≥ 0, t = 1, 2, . . . , and for all x, x′ ∈ X . Under (LC), if the agent observes at each stage t an element vt of ∂ ft(Xt), straightforward online policies based on gradient descent enjoy a bound of the form Regx(T ) = O(ḠT T 1/2), with Ḡ2T = T−1 ∑T t=1 G 2 t [11, 35, 41]. In particular, if G ≡ lim supT→∞ ḠT < ∞ (e.g., if each ft is G-Lipschitz continuous over X ), we have the bound
Regx(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 ∂ ft of ft admits a continuous selection ∇ ft(x) ∈ ∂ ft(x) for all x ∈ dom ∂ ft ≡ {x ∈ X : ∂ ft(x) , ∅}.

heading: 3 Riemann–Lipschitz continuity
text: 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 /∂xi 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 :
1For 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.
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 {ei}di=1 of d, we define the metric tensor of 〈·, ·〉x as the matrix g(x) ∈ d×d with components
gi j(x) = 〈ei, e j〉x i, j = 1, . . . , d. (3)
The norm of z ∈ V at x ∈ X is then defined as
‖z‖2x ≡ 〈z, z〉2x = d∑
i, j=1
gi j(x)ziz j = z>g(x)z. (4)
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‖2x = ∑d i=1 z 2 i and 〈z, z′〉x = ∑d i=1 ziz ′ i , both independent of x.
Example 2 (Hyperbolic geometry). The Poincaré metric on the positive orthant X = d++ is
g(x) = diag(1/x21, . . . , 1/x 2 d), (5)
leading to the local norm ‖z‖2x = ∑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 Lg[γ] =∫ 1 0 ‖γ̇(s)‖γ(s) ds, and the Riemannian distance between x1, x2 ∈ X is given by
distg(x1, x2) = infγ Lg[γ]. (6)
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 distg(x, x′) for some G ≥ 0 and all x, x′ ∈ X . (7)
Albeit simple to state, (RLC) may be difficult to verify because it requires the computation of the distance function distg 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{x′ − x : 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 . (8) 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:
[grad f (x)]i = ∑d
j=1 g(x)−1i j ∂ j f (x). (9)
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:
Example 3. Let X = [0, 1]d\{0} (so X is convex but neither open nor closed) and let f (x) = − log(a>x) for some positive vector a ∈ d++. If we take gi j(x) = δi j/(x1 + · · · + xd)2, we get
‖grad f (x)‖2x = ∑d i=1 a 2 i · (∑d i=1 xi )2 (a>x)2 ≤ ∑d i=1 a 2 i (min j a j)2 . (10)
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 gi j(x) = φ(x)2δi j gives ‖grad f (x)‖2x = φ(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.,
h(x′) ≥ h(x) + 〈∇h(x)|x′ − x〉 + 12 K‖x ′ − x‖2x (11)
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 D(p, x) = h(p) − h(x) − 〈∇h(x)|p − x〉. (12)
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 riX ⊆ X ◦ ≡ dom ∂h ⊆ X . To connect the two, we will say that h is a Riemann–Legendre regularizer when X ◦ = riX and D(p, xn)→ 0 whenever xn → 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
D(p, x) ≥ 12 K‖p − x‖ 2 x. (13)
2For 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.
3Following standard convex analysis terminology, “proper” means here that h , +∞ while lower semicontinuous refers to the property that lim infx→x0 h(x) ≥ h(x0) for all x0 ∈ V .
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) = −∑di=1 log xi. It is easy to see that h(x) is strongly convex relative to the standard Euclidean norm ‖·‖2. Moreover, we have
D(p, x) = d∑
i=1
[ pi xi − log pi xi − 1 ] (14)
and, by Taylor’s theorem with Lagrange remainder, we readily get
D(p, x) ≥ 1 2 d∑ i=1 (pi − xi)2 xi = ‖p − x‖2x (15)
where ‖z‖x = ∑d i=1 z 2 i /xi 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 + r2)/ ∑d
i=1 xi with r 2 = ∑d i=1 x 2 i . Then,
a tedious (but otherwise straightforward) algebraic calculation gives
D(p, x) ≥ d∑
i=1
(xi − pi)2 ( ∑d j=1 x j)2 = ‖p − x‖2x (16)
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:
Xt+1 = arg min x∈X
{∑t s=1 fs(x) + γ−1h(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 ft. To model this feedback structure, we assume that once Xt has been chosen, the optimizer receives an estimate vt of ∇ ft(Xt) satisfying the following hypotheses:
a) Unbiasedness: [vt | Ft] = ∇ ft(Xt). (17a) b) Finite mean square: [‖vt‖2∗ | Ft] ≤ M2t . (17b)
In the above, ‖·‖∗ denotes the dual norm of ‖·‖Xt , i.e., the Riemannian norm at the point Xt where the oracle was called (we suppress here the index Xt and write ‖·‖∗ instead of ‖·‖Xt ,∗ to lighten the notation). In particular, the oracle feedback vt may fail to be bounded in L2 relative to any global norm on V∗; as such, (17) is considerably weaker than the standard L2-boundedness assumption for global norms. Finally, in terms of measurability, the expectation in (17) is conditioned on the history Ft of Xt up to stage t; since vt is generated randomly from Xt, it is not Ft-measurable. To proceed, the main idea of mirror descent is to replace fs(x) in (FTRL) with the first-order surrogate fs(x)← fs(Xs) + 〈∇ fs(Xs)|x − Xs〉. In this way, substituting vs for the estimate of ∇ fs(Xs) received at stage s, we obtain the linearized FTRL scheme
Xt+1 = arg minx∈X { γ ∑t s=1〈vs|x〉 + h(x) } . (18)
To rewrite this process in recursive form, introduce the auxiliary (dual) variable
Yt+1 = Yt − γvt (19)
so Yt+1 = −γ ∑t s=1 vs, and hence
Xt+1 = arg min x∈X {h(x) − 〈Yt+1|x〉} = arg max x∈X {〈Yt+1|x〉 − h(x)}. (20)
Therefore, letting Q(y) = arg maxx∈X {〈y|x〉 − h(x)} (21)
denote the so-called “mirror map” of the method, we obtain the following incarnation of the online mirror descent (OMD) algorithm:
Yt+1 = Yt − γvt Xt+1 = Q(Yt+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” xc = 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 ft : X → is convex and satisfies (RLC) with constant Gt. 3. The optimizer’s aggregate loss ∑T t=1 ft 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 ) ≡ Regx∗ (T ), Ḡ2T = T−1 ∑T t=1 G 2 t , and M̄ 2 T = T −1 ∑T t=1 M 2 t . Then:
a) The FTRL algorithm enjoys the regret bound Reg(T ) ≤ D(x∗, xc)
α + 2αḠ2T K √T . (22a) b) The OMD algorithm with noisy feedback of the form (17) enjoys the mean regret bound
[Reg(T )] ≤ D(x∗, xc)
α + αM̄2T 2K √T . (22b) In particular, if supt Gt < ∞, supt Mt < ∞, both algorithms guarantee O( √ T ) regret. Remark 3. We emphasize here that the O( √
T ) regret bound above is achieved even if X is unbounded or if the “Bregman depth” H ≡ supx∈X D(x, xc) = sup h − inf h of X is infinite.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.
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 Xt ∈ X and Yt ∈ 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
Φ(p, y) = h(p) + h∗(y) − 〈y|p〉 for all p ∈ X , y ∈ V∗, (23) 4To see this, simply note that D(x, xc) = h(x) − h(xc) − 〈∇h(xc)|x − xc〉 < ∞ for all x ∈ X = dom h (recall
also that, since xc = arg min h, we have 0 ∈ ∂h(xc) so xc ∈ dom ∂h).
with h∗(y) = maxx∈X {〈y|x〉 − h(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
Φ(p, y − γv) ≤ Φ(p, y) − 〈γv|Q(y) − p〉 + γ 2
2K ‖v‖2Q(y),∗. (24)
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., supx [‖∇F(x;ω)‖2x,∗] ≤ M2 for some M > 0. If (OMD) is run for T iterations with a constant step-size of the form α/ √ T and stochastic gradients vt = ∇F(Xt;ωt) generated by an i.i.d. sequence ωt ∈ Ω, we have
[ f (X̄T )] ≤ min f + [
Dc α + αM2 2K
] 1 √
T (25)
where X̄T = (1/T ) ∑T t=1 Xt is the “ergodic average” of Xt and Dc = infx∗∈arg min f D(x ∗, xc) < ∞ denotes the Bregman distance of the prox-center xc 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 Xt. 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 L2. Suppose further that arg min f is bounded and (OMD) is run with a sequence of stochastic gradients vt = ∇F(Xt;ω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, Xt 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 Xt implies that of X̄t, so the convergence of [ f (X̄t)] 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):
f (x) = ∑m
j=1
[ u j log
u j (Hx) j
+ (Hx) j − u j ]
(26)
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/x1, . . . , 1/xd) 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 ≈ 105 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.
Acknowledgments
The authors gratefully acknowledge financial support from the French National Research Agency (ANR) under grants ORACLESS (ANR–16–CE33–0004–01) and ELIOT (ANR-18-CE40-0030), as well as the FAPESP 2018/12579-7 project. This work also benefited from financial support by MIAI Grenoble Alpes (Multidisciplinary Institute in Artificial Intelligence).
References [1] Abernethy, Jacob, Peter L. Bartlett, Alexander Rakhlin, Ambuj Tewari. 2008. Optimal strategies and
minimax lower bounds for online convex games. COLT ’08: Proceedings of the 21st Annual Conference on Learning Theory.
[2] Balandat, Maximilian, Walid Krichene, Claire Tomlin, Alexandre Bayen. 2016. Minimizing regret on reflexive Banach spaces and Nash equilibria in continuous zero-sum games. NIPS ’16: Proceedings of the 30th International Conference on Neural Information Processing Systems.
[3] Bauschke, Heinz H., Jérôme Bolte, Marc Teboulle. 2017. A descent lemma beyond Lipschitz gradient continuity: First-order methods revisited and applications. Mathematics of Operations Research 42(2) 330–348.
[4] Bauschke, Heinz H., Patrick L. Combettes. 2017. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. 2nd ed. Springer, New York, NY, USA.
[5] Bécigneul, Gary, Octavian-Eugen Ganea. 2019. Riemannian adaptive optimization methods. ICLR ’19: Proceedings of the 2019 International Conference on Learning Representations.
[6] Beck, Amir, Marc Teboulle. 2003. Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters 31(3) 167–175.
[7] Bertero, Mario, Patrizia Boccacci, Gabriele Desiderà, Giuseppe Vicidomini. 2009. Image deblurring with Poisson data: from cells to galaxies. Inverse Problems 25(12) 123006.
[8] Bolte, Jérôme, Shoham Sabach, Marc Teboulle, Yakov Vaisbourd. 2018. First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems. SIAM Journal on Optimization 28(3) 2131–2151.
[9] Bottou, Léon. 1998. Online learning and stochastic approximations. On-line learning in neural networks 17(9) 142.
[10] Bubeck, Sébastien. 2015. Convex optimization: Algorithms and complexity. Foundations and Trends in Machine Learning 8(3-4) 231–358.
[11] Bubeck, Sébastien, Nicolò Cesa-Bianchi. 2012. Regret analysis of stochastic and nonstochastic multi-armed bandit problems. Foundations and Trends in Machine Learning 5(1) 1–122.
[12] Chen, Gong, Marc Teboulle. 1993. Convergence analysis of a proximal-like minimization algorithm using Bregman functions. SIAM Journal on Optimization 3(3) 538–543.
[13] Facchinei, Francisco, Jong-Shi Pang. 2003. Finite-Dimensional Variational Inequalities and Complementarity Problems. Springer Series in Operations Research, Springer.
[14] Ferreira, Orizon P. 2006. Proximal subgradient and a characterization of Lipschitz function on Riemannian manifolds. Journal of Mathematical Analysis and Applications 313 587–597.
[15] Hall, P., C. C. Heyde. 1980. Martingale Limit Theory and Its Application. Probability and Mathematical Statistics, Academic Press, New York.
[16] Hanzely, Filip, Peter Richtárik. 2018. Fastest rates for stochastic mirror descent methods. https: //arxiv.org/abs/1803.07374.
[17] Hanzely, Filip, Peter Richtárik, Lin Xiao. 2018. Accelerated Bregman proximal gradient methods for relatively smooth convex optimization. https://arxiv.org/abs/1808.03045.
[18] He, Niao, Zaid Harchaoui, Yichen Wang, Le Song. 2016. Fast and simple optimization for Poisson likelihood models. https://arxiv.org/abs/1608.01264.
[19] Jiang, Houyuan, Huifu Xu. 2008. Stochastic approximation approaches to the stochastic variational inequality problem. IEEE Trans. Autom. Control 53(6) 1462–1475.
[20] Karimi, Hamed, Julie Nutini, Mark Schmidt. 2016. Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition. https://arxiv.org/abs/1608.04636.
[21] Krichene, Walid. 2016. Continuous and discrete dynamics for online learning and convex optimization. Ph.D. thesis, Department of Electrical Engineering and Computer Sciences, University of California, Berkeley.
[22] Lee, John M. 1997. Riemannian Manifolds: an Introduction to Curvature. No. 176 in Graduate Texts in Mathematics, Springer.
[23] Ljung, Lennart. 1978. Strong convergence of a stochastic approximation algorithm. Annals of Statistics 6(3) 680–696.
[24] Lu, Haihao. 2017. "Relative-continuity" for non-Lipschitz non-smooth convex optimization using stochastic (or deterministic) mirror descent. https://arxiv.org/abs/1710.04718.
[25] Lu, Haihao, Robert M. Freund, Yurii Nesterov. 2018. Relatively-smooth convex optimization by first-order methods and applications. SIAM Journal on Optimization 28(1) 333–354.
[26] Mertikopoulos, Panayotis, Zhengyuan Zhou. 2019. Learning in games with continuous action sets and unknown payoff functions. Mathematical Programming 173(1-2) 465–507.
[27] Nemirovski, Arkadi Semen, Anatoli Juditsky, Guanghui Lan, Alexander Shapiro. 2009. Robust stochastic approximation approach to stochastic programming. SIAM Journal on Optimization 19(4) 1574–1609.
[28] Nemirovski, Arkadi Semen, David Berkovich Yudin. 1983. Problem Complexity and Method Efficiency in Optimization. Wiley, New York, NY.
[29] Nesterov, Yurii. 2007. Dual extrapolation and its applications to solving variational inequalities and related problems. Mathematical Programming 109(2) 319–344.
[30] Nesterov, Yurii. 2009. Primal-dual subgradient methods for convex problems. Mathematical Programming 120(1) 221–259.
[31] Nevel’son, M. B., Rafail Z. Khasminskii. 1976. Stochastic Approximation and Recursive Estimation. American Mathematical Society, Providence, RI.
[32] Robbins, Herbert, David Sigmund. 1971. A convergence theorem for nonnegative almost supermartingales and some applications. J. S. Rustagi, ed., Optimizing Methods in Statistics. Academic Press, New York, NY, 233–257.
[33] Rockafellar, Ralph Tyrrell. 1970. Convex Analysis. Princeton University Press, Princeton, NJ. [34] Scutari, Gesualdo, Francisco Facchinei, Daniel Pérez Palomar, Jong-Shi Pang. 2010. Convex optimization,
game theory, and variational inequality theory in multiuser communication systems. IEEE Signal Process. Mag. 27(3) 35–49.
[35] Shalev-Shwartz, Shai. 2011. Online learning and online convex optimization. Foundations and Trends in Machine Learning 4(2) 107–194.
[36] Shalev-Shwartz, Shai, Yoram Singer. 2007. Convex repeated games and Fenchel duality. Advances in Neural Information Processing Systems 19. MIT Press, 1265–1272.
[37] Sra, Suvrit, Sebastian Nowozin, Stephen J. Wright. 2012. Optimization for Machine Learning. MIT Press, Cambridge, MA, USA.
[38] Xiao, Lin. 2010. Dual averaging methods for regularized stochastic learning and online optimization. Journal of Machine Learning Research 11 2543–2596.
[39] Zhang, Hui, Wotao Yin. 2013. Gradient methods for convex minimization: Better rates under weaker conditions. https://arxiv.org/abs/1303.4645.
[40] Zhou, Zhengyuan, Panayotis Mertikopoulos, Nicholas Bambos, Stephen Boyd, Peter W. Glynn. 2017. Stochastic mirror descent for variationally coherent optimization problems. NIPS ’17: Proceedings of the 31st International Conference on Neural Information Processing Systems.
[41] Zinkevich, Martin. 2003. Online convex programming and generalized infinitesimal gradient ascent. ICML ’03: Proceedings of the 20th International Conference on Machine Learning. 928–936.

heading: A Riemann–Lipschitz continuity
text: 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 x1, x2 ∈ 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) 7→ expx(z) = γz(1), (A.1) where γz denotes the unique geodesic emanating from x with initial velocity vector γ̇z = z. We then have expx(tz) = γz(t) for all t and, moreover, for sufficiently small r > 0, the restriction of expx to a ball of radius r in Z is a diffeomorphism onto its image in U . The largest positive number ix such that the above holds for all r < ix 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) = expx(tz/‖z‖x). (A.2)
Since γ is a geodesic, we readily obtain distg(x, γ(t)) = t for all sufficiently small t. Also, by construction, we have exp−1x γ(t) = tz/‖z‖x. Hence, with f convex, it follows that, for some constant α > 0 and for sufficiently small positive δ < ix, we have
f (γ(t)) − f (x) ≥ 〈z, exp−1x (γ(t))〉 (A.3) − α distg(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 expz(t) is a diffeomorphism for sufficiently small δ > 0 – and hence, for all δ < ix. 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 − αt2. (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)‖x = ‖z‖x ≤ G, (A.6) 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) = distg(x, x′) if distg(x, x′) < δ, and w(x) = K distg(x, x′) + ε2/(δ − ε) when δ < distg(x, x′) < 2δ, with ε = distg(x, x′) − δ. A simple calculation then shows that ‖gradw(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 distg(x, x′). (A.7)
Then, by interchanging x and x′ above, we obtain | f (x) − f (x′)| ≤ K distg(x, x′). Since K > G has been chosen arbitrarily, (7) follows.
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|x′ − x〉 ≥ 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 = x′ − x, 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): Q(y) = arg max
x∈X {〈y|x〉 − h(x)}, (B.3)
we readily get Q(y) = ∇h∗(y). (B.4)
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) 3 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(p− x))− y|p− x〉 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|p − x〉 ≥ 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) ≥ K2 ‖x − p‖2x
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; p − x) = 〈y|p − x〉 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:
Φ(p, y+) = h(p) + h∗(y+) − 〈y+|p〉 Φ(p, y) = h(p) + h∗(y) − 〈y|p〉. (B.12)
Then, by subtracting the above we get:
Φ(p, y+) − Φ(p, y) = h(p) + h∗(y+) − 〈y+|p〉 − h(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|x − p〉 (B.13)
and our proof is complete.
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:
Φ(p, y + v) ≤ Φ(p, y) + 〈v|x − p〉 + 1 2K ‖v‖2x,∗ (B.14)
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) and hence, after rearranging:
Φ(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
〈v|Q(y + v) − x〉 ≤ K 2 ‖Q(y + v) − x‖2x + 1 2K ‖v‖2x,∗ (B.17)
Our claim then follows by the fact that Φ(Q(y+ v), y) ≥ K2 ‖Q(y+ v)− x‖2x (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
Regx(T ) ≤ h(x) − h(X1) + T∑
t=1
[ ft(Xt) − ft(Xt+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:
Regx(T ) ≤ D(x, xc)
γ + 2γ K T∑ t=1 G2t (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
Ft(x) = t−1∑ s=1 fs(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, Xt ∈ arg min Ft(x), so 〈∇Ft(Xt)|x − Xt〉 ≥ 0 and, likewise, 〈∇Ft+1(Xt+1)|x− Xt+1〉 ≥ 0 for all x ∈ X . Furthermore, since h is K-strongly convex relative to g, Ft and Ft+1 will be (K/γ)-strongly convex relative to g.
Putting all this together, we obtain:
Ft(Xt+1) ≥ Ft(Xt) + K 2γ ‖Xt+1 − Xt‖2Xt+1 (C.6a)
Ft+1(Xt) ≥ Ft+1(Xt+1) + K 2γ ‖Xt − Xt+1‖2Xt (C.6b)
and hence, after summing the above inequalities:
ft(Xt) − ft(Xt+1) ≥ K 2γ ‖Xt+1 − Xt‖2Xt+1 + K 2γ ‖Xt − Xt+1‖2Xt
≥ K 2γ ‖Xt+1 − Xt‖2Xt . (C.7)
On the other hand, Lemma C.2 gives
ft(Xt) − ft(Xt+1) ≤ Gt‖Xt − Xt+1‖Xt (C.8)
so, combining the last two inequalities, we get:
K 2γ ‖Xt+1 − Xt‖Xt ≤ Gt. (C.9)
Therefore, plugging this back into (C.8) yields
ft(Xt) − ft(Xt+1) ≤ 2γG2t
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:
Yt+1 = Yt − γvt Xt+1 = Q(Yt+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 vt = ∇ ft(Xt) + Ut+1; hence, by the unbiasedness assumption (17a), it follows that [Ut+1 | Ft] = 0, i.e., Ut is a martingale difference sequence (MDS) relative to Ft. Now, applying Proposition B.1 to (OMD), we get:
Φ(x∗,Yt+1) ≤ Φ(x∗,Yt) − γ〈vt |Xt − x∗〉 + γ2
2K ‖vt‖2Xt ,∗
= Φ(x∗,Yt) + γ〈∇ ft(Xt)|x∗ − Xt〉 − γ〈Ut+1|Xt − x∗〉 + γ2
2K ‖vt‖2Xt ,∗. (D.1)
Hence, after rearranging and telescoping, we obtain
Reg(T ) ≤ T∑
t=1
〈∇ ft(Xt)|Xt − x∗〉 ≤ D(x∗, xc)
γ + T∑ t=1 ξt+1 + γ 2K T∑ t=1 ‖vt‖2Xt ,∗ (D.2)
where, in the last line, we used the definition of the Riemannian dual norm ‖·‖∗ ≡ ‖·‖x∗,∗, and we set ξt+1 = 〈Ut+1|x∗ − Xt〉. 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 L2 [15]: Theorem E.1. Let Yt = ∑t i=1 ζi be a martingale and βt a non-decreasing positive sequence such that limt→∞ βt = ∞. Then, lim t→∞ Yt/βt = 0 almost surely (E.1)
on the set ∑∞
t=1 β −2 t [ζ 2 t | Ft−1] < ∞.
The second is a convergence result for quasi-supermartingales due to Robbins and Sigmund [32]: Lemma E.1. Let (Ft)t∈ be a non-decreasing sequence of σ− algebras. Let (αt)t∈ , (θt)t∈ nonnegative Ft− measurable random variables, (ηt)t∈ is an Ft− measurable non-negative summable random variable and the following inequality holds:
[αt+1 | Ft] ≤ αt − θt + ηt almost surely (E.2) Then, (αt)t∈ converges almost surely towards a [0,∞)-valued random variable.
An application of this lemma leads us to the following result which is of independent interest:
Proposition E.1. Let Xt be the sequence of iterates generated by (OMD) run with a step-size sequence γt such that ∑∞ t=1 γ 2 t < ∞ and a stochastic oracle as in the statement of Theorems 2 and 3. Then, for all x∗ ∈ arg min f , Φ(x∗,Yt) converges with probability 1.
Proof. Let x∗ ∈ arg min f . Recalling our main estimation:
Φ(x∗,Yt+1) ≤ Φ(x∗,Yt) − γt〈vt |Xt − x∗〉x + γ2t 2K ‖vt‖2Xt ,∗ (E.3)
and taking conditional expectations on both sides, we get due to Ft− measurability arguments:
[Φ(x∗,Yt+1)|Ft] ≤ Φ(x∗,Yt) − γt〈vt |Xt − x∗〉x + γ2t 2K [‖vt‖2Xt ,∗|Ft]. (E.4)
Since, (2K)−1 ∑∞
t=1 γ 2 t [‖vt‖2Xt ,∗|Ft] ≤ M(2K) −1 ∑∞ t=1 γ 2 t < ∞. Thus, by applying the above we get
the result.
Having this at hand, we can establish the following proposition:
Proposition E.2. Let Xt be the sequence of iterates generated by (OMD) with assumptions as in Theorem 3. Then, for all x∗ ∈ arg min f , the sequence ‖Xt − x∗‖Xt is bounded with probability 1.
Proof. Recalling our main estimation and taking condition expectations on both sides, we get:
[Φ(x∗,Yt+1) | Ft] ≤ Φ(x∗,Yt) − γt〈vt |Xt − x∗〉x + γ2t 2K [‖vt‖2Xt ,∗|Ft] (E.5)
Hence, by the above corollary, we have that the sequence Φ(x∗,Yt) converges with probability 1 for all x∗ ∈ arg min f . Thus, it is also bounded with probability 1 for all x∗. We then get
‖Xt − x∗‖2Xt ≤ 2 K Φ(x∗,Yt) (E.6)
which concludes our proof.
We continue by showing that Xt possesses a subsequence that converges to arg min f :
Proposition E.3. Let Xt be the sequence of iterates generated by (OMD) with assumptions as in Theorem 3. Then, with probability 1, there exists a (possibly random) subsequence of Xt which converges to arg min f .
Proof. Assume to the contrary that, with positive probability, the sequence Xt generated by (OMD) admits no limit points in arg min f . Conditioning on this event, there exists a (nonempty) closed set C ⊂ X which is separated by neighborhoods from arg min f and is such that Xt ∈ C for all suffiently large t. Then, by relabeling Xt if necessary, we can assume without loss of generality that Xt ∈ C for all t ∈ . Thus, by Proposition B.1, we get:
Φ(x∗,Yt+1) ≤ Φ(x∗,Yt) − γt〈vt |Xt − x∗〉 + γ2t 2K ‖vt‖2Xt ,∗
= Φ(x∗,Yt) − γt〈∇ f (Xt)|Xt − x∗〉 − γt〈Ut+1|Xt − x∗〉 + γ2t 2K ‖vt‖2Xt ,∗ ≤ Φ(x∗,Yt) − γtδ(C) + γtξt+1 + γ2t 2K ‖vt‖2Xt ,∗ (E.7)
where in the last line we set δ(C) = inf{〈∇ f (x)|x − x∗〉 : x∗ ∈ arg min f , x ∈ C} > 0 (by (SI)), Ut+1 = vt − ∇ f (Xt), ξt+1 = −〈Ut+1|Xt − x∗〉 and βt = ∑t i=1 γi. Thus, by telescoping and factorizing we get:
Φ(x∗,Yt+1) ≤ Φ(x∗,Y1) − βt δ(C) − ∑ts=1 γsξs+1βt − ∑t s=1 γ 2 s‖vs‖2Xs,∗ 2Kβt  (E.8)
By the unbiasedness assumption for Ut, we have [ξt+1 | Ft] = 〈 [Ut+1 | Ft]|Xt − x∗〉 = 0. Moreover, for all x∗ ∈ arg min f , we have
∞∑ t=1 γ2t [ξt+1 | Ft] ≤ ∞∑ t=1 γ2t ‖Xt − x∗‖2Xt [Ut+1 | Ft] ≤ ∞∑ t=1 γ2t Φ(x ∗,Yt) [Ut+1 | Ft] < ∞ (E.9)
where the last (strict) inequality is obtained due to the finite mean square property, the boundness of Φ(x∗,Yt) and the fact that ∑∞ t=1 γ 2 t < ∞. Thus, we can apply the law of large numbers for L2−
martingales stated above and conclude that β−1t ∑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‖vs‖2Xt ,∗, since vs+1 is Fs-measurable for all s = 1, 2 . . . , t − 1 we have:
[S t+1 | Ft] =  t−1∑
i=1
γ2t ‖vi‖2xi,∗ + γ 2 t ‖vt‖2Xt ,∗ ∣∣∣∣∣∣∣ Ft  = S t + γ2t [‖vt‖2Xt ,∗ ∣∣∣ Ft] ≥ S t (E.10)
so S t is a submartingale with respect to Ft. Furthermore, by the law of total expectation, we also get:
[S t+1] = [ [S t+1 | Ft]] ≤ σ2 t∑
i=1
γ2i ≤ σ2 ∞∑
t=1
γ2t < ∞, (E.11)
implying that S t is bounded in L1. Thus, due to Doob’s submartingale convergence theorem [15], we coclude that S t converges to some (almost surely finite) random variable S∞ so limt→∞ S t+1βt = 0 with probability 1.
Now, by letting t → ∞ in (E.8), we get Φ(x∗,Yt)→ −∞, a contradiction. Going back to our original assumption, this shows that there exists a subsequence of Xt which converges to arg min f with probability 1, as claimed.
With all this at hand, we proceed to the proof of our last-iterate convergence result:
Proof of Theorem 3. By the boundedness (and hence compactness) of arg min f , Proposition E.3 implies that, with probability 1, there exists some x∗ ∈ arg min f such that Xtk → x∗ for some (possibly random) subsequence Xtk of Xt. By the Riemann–Legendre property of h, it follows that Φ(x
∗,Ytk ) = D(x∗, Xtk ) → 0 as k → ∞, implying in turn that limt→∞ D(x∗, Xt) = 0 (by Proposition E.1). Since D(x∗, Xt) ≥ K‖Xt − x∗‖2Xt ≥ µ‖Xt − x ∗‖2, we conclude that Xt → x∗, and our proof is complete.
F Applications to Poisson inverse problems and numerical experiments
F.1 Detailed statement of the problem
The class of Poisson inverse problems that we consider stem from linear systems of the form
u = Hx + z (F.1)
where
• x ∈ d+ is the object under study (a signal, image, . . . ). • u ∈ m++ is the observed data (usually m d). • The kernel matrix H ∈ m×d+ is a representation of the data-gathering protocol and is
typically highly ill-conditioned (e.g., a Toeplitz matrix in the case of image deconvolution problems). • z ∈ m is the noise affecting the measurements.
When data points are obtained by means of a counting process, measurements can be modeled as Poisson random variables of the form u j ∼ Pois(Hx) j. Then, up to an additive constant, the log-likelihood of x ∈ d given an observation u ∈ m++ will be
L(x; u) = − m∑
j=1
[ u j log
u j (Hx) j
+ (Hx) j − u j ] . (F.2)
Hence, obtaining a maximum likelihood estimate for x leads to the archetypal Poisson inverse problem minimize f (x) ≡ DKL(u,Hx), subject to x ∈ d+, (PIP)
where DKL(p, q) = ∑m
j=1[p j log(p j/q j) + q j − p 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 106 and 107), 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
∂ f j ∂xi = u jH ji (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 norm ‖z‖2x = (x1 + · · · + xd)2 ∑d i=1 z 2 i for all v ∈ d. (F.4)
In this case, we have
‖∇ f j(x)‖2x = d∑
i=1
u2j H2jix 2 i
(Hx)2j =
u2j ∑d i=1 H 2 jix
2 i[∑d i=1 H jixi ]2 = O(u2j ), (F.5)
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.
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 xi log xi. 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 xi 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 X̄T = ∑T t=1 γtXt /∑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.

Review summaries:
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.

This paper investigates online and stochastic convex optimization problems in which the objective function is not Lipschitz continuous. The originality of this study lies in the use of Riemannian geometry. Specifically, the standard condition of Lipschitz continuity is replaced with a more general condition involving Riemannian distances and called Riemann-Lipschitz Continuity (RLC). Based on an appropriate definition of Riemannian regularizer and a generalization of Fenchel coupling to Riemannian geometry, the authors provide $O(\sqrt T)$ regret (resp. risk) bounds for the online (resp. stochastic) mirror descent algorithm, under the Riemann-Lipchitz condition. The performance of the algorithm is validated on Poisson inverse problems.