Mohammad Reza Karimi ⋆,

Ya-Ping Hsieh
email: yaping.hsieh@inf.ethz.ch

AND Andreas Krause
email: krausea@ethz.ch.

THE DYNAMICS OF RIEMANNIAN ROBBINS-MONRO ALGORITHMS

Many important learning algorithms, such as stochastic gradient methods, are often deployed to solve nonlinear problems on Riemannian manifolds. Motivated by these applications, we propose a family of Riemannian algorithms generalizing and extending the seminal stochastic approximation framework of Robbins and Monro [60]. Compared to their Euclidean counterparts, Riemannian iterative algorithms are much less understood due to the lack of a global linear structure on the manifold. We overcome this difficulty by introducing an extended Fermi coordinate frame which allows us to map the asymptotic behavior of the proposed Riemannian Robbins-Monro (RRM) class of algorithms to that of an associated deterministic dynamical system under very mild assumptions on the underlying manifold. In so doing, we provide a general template of almost sure convergence results that mirrors and extends the existing theory for Euclidean Robbins-Monro schemes, albeit with a significantly more involved analysis that requires a number of new geometric ingredients. We showcase the flexibility of the proposed RRM framework by using it to establish the convergence of a retraction-based analogue of the popular optimistic / extra-gradient methods for solving minimization problems and games, and we provide a unified treatment for their convergence.

Introduction

Background and motivation. Consider a nonlinear system of equations of the general form Find z * ∈ M such that V (z * ) = 0 (NLS)

where M is a smooth manifold and V is a vector field on M. Root-finding problems of this type play a critical role in many areas of mathematical programming and learning theory, from Riemannian optimization and game theory to reinforcement learning and optimal control -e.g., when designing the optimal path of a robotic arm or employing natural gradient methods [START_REF] Kakade | A natural policy gradient[END_REF] over smooth statistical manifolds.

In this paper, we are interested in the case where V is stochastic, i.e., V (z) = E[V(z; ω)] for some random variable ω with unknown distribution. In this case, when M = R d , the method of choice for solving (NLS) is the Robbins-Monro (RM) algorithm

Z n+1 = Z n + γ n V(Z n ; ω n ) (RM)
where γ n > 0 is a variable step-size sequence and ω n is an i.i.d. sequence of samples (equivalently, this could be thought of as accessing V via a stochastic black-box oracle). This method was introduced in the seminal papers of Robbins and Monro [START_REF] Robbins | A stochastic approximation method[END_REF] and Kiefer and Wolfowitz [START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF], and the first general convergence results were obtained by Ljung [START_REF] Ljung | Analysis of recursive stochastic algorithms[END_REF][START_REF] Ljung | Strong convergence of a stochastic approximation algorithm[END_REF] for gradient problems, i.e., V = -∇f for some potential function f on M. This led to substantial activity on the topic, with major contributions by Kushner and co-authors [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF][START_REF] Kushner | Stochastic Approximation Methods for Constrained and Unconstrained Systems[END_REF][START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF], and many others. However, the linear structure of R d is deeply ingrained in all these works -as well as the method's very definition -preventing its use to solve many important stochastic approximation (SA) problems on manifolds: the d-dimensional torus for a robotic arm with d joints, the Grassman or Stiefel manifolds for robust principal component analysis, the hyperbolic spaces for text and graph embedding, etc.

Our contributions in the context of related work. Our paper seeks to lift this limitation by replacing the "+" operation in (RM) with the Riemannian exponential map exp Zn (•) on M -or, more generally (and often more tractably), a retraction on M based at Z n . In Riemannian optimization, this approach was pioneered by Bonnabel [START_REF] Bonnabel | Stochastic gradient descent on Riemannnian manifolds[END_REF] who examined the case where V is the Riemannian gradient of some objective function f . Subsequent works [START_REF] Boumal | Global rates of convergence for nonconvex optimization on manifolds[END_REF][START_REF] Criscitiello | Efficiently escaping saddle points on manifolds[END_REF][START_REF] Lezcano-Casado | Curvature-dependant global convergence rates for optimization on manifolds of bounded geometry[END_REF][START_REF] Tripuraneni | Averaging stochastic gradient descent on Riemannnian manifolds[END_REF][START_REF] Wang | No-regret online learning over Riemannnian manifolds[END_REF][START_REF] Zhang | First-order methods for geodesically convex optimization[END_REF] expanded on the results of [START_REF] Bonnabel | Stochastic gradient descent on Riemannnian manifolds[END_REF] for Riemannian stochastic gradient descent, while similar results were obtained in [START_REF] Bento | Iteration-complexity of gradient, subgradient and proximal point methods on Riemannnian manifolds[END_REF][START_REF] Ferreira | Proximal point algorithm on Riemannnian manifolds[END_REF][START_REF] Huang | Riemannnian proximal gradient methods[END_REF][START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF] for Riemannian proximal point methods.

All these works focus exclusively on the case where V is a gradient field, so they do not apply to general, non-gradient instances of (NLS) which are crucial for min-max problems, games and multi-agent learning problems. A partial extension to the non-gradient case was provided by a line of works [START_REF] Chen | Modified tseng's extragradient methods for variational inequality on Hadamard manifolds[END_REF][START_REF] Fan | Tseng's extragradient algorithm for pseudomonotone variational inequalities on Hadamard manifolds[END_REF][START_REF] Ferreira | Singularities of monotone vector fields and an extragradienttype algorithm[END_REF][START_REF] Khammahawong | An extragradient algorithm for strongly pseudomonotone equilibrium problems on Hadamard manifolds[END_REF][START_REF] Neto | An extragradient method for equilibrium problems on Hadamard manifolds[END_REF][START_REF] Gj | Korpelevich's method for variational inequality problems on Hadamard manifolds[END_REF], which examined the use of Riemannian extra-gradient methods under the assumption of (geodesic) monotonicity. This is a strong, convexity-type assumption which posits that V globally points towards its (necessarily connected) root system in a suitable, geodesic sense; convergence is then obtained following a similar line of reasoning as in the case of monotone operator theory in Euclidean spaces [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF].

Our paper does not make any such assumptions and directly examines the dynamics of Riemannian Robbins-Monro methods for general vector fields V , gradient and non-gradient, monotone and non-monotone alike. In this regard, our main contributions can be summarized as follows:

(1) We introduce a generalized Riemannian Robbins-Monro template which includes as special cases all methods mentioned above (Riemmanian stochastic gradient descent, extra-gradient, proximal point methods, etc.), as well as a number of new SA schemes for (NLS).

(2) We show that, under mild technical conditions on M, the sequence of generated points forms an "approximate solution" -an asymptotic pseudotrajectory (APT) to be exact -of an associated deterministic dynamical system (Theorem 1), and converges with probability 1 to the so-called internally chain-transitive (ICT) sets thereof (Theorem 2).

In gradient and strictly monotone problems, these ICT sets are precisely the roots of V [63], so we immediately recover many of the asymptotic convergence results mentioned above (often under much weaker assumptions). In addition, our framework applies to several interesting settings beyond gradient or monotone systems -such as ordinal potential games, supermodular games, and cooperative dynamics -and covers a significantly wider class of SA schemes.

Tools and techniques. In the absence of a linear structure on M, the major challenge we have to overcome is the lack of a suitable coordinate frame with which to analyze the trajectories of Riemannian SA algorithms. This reflects the dichotomy that, unlike the case of R d , points and vectors on manifolds obey fundamentally different rules and have to be compared using different moving frames. To circumvent this obstacle, we introduce an extended Fermi coordinate frame inspired by Manasse and Misner [START_REF] Manasse | Fermi normal coordinates and some basic concepts in differential geometry[END_REF], and we use it to prove that Riemannian SA schemes enjoy similar error bounds as in Euclidean spaces, up to some high-order terms that vanish in the long run. The aggregation and propagation of these errors can then be controlled using arguments from martingale limit theory which ultimately yield the convergence properties mentioned above.

A concurrent approach to establish the APT property in Riemannian SA schemes is due to Shah [START_REF] Shah | Stochastic approximation on riemannian manifolds[END_REF], who assumes the existence of a local diffeomorphism mapping geodesic interpolations to linear interpolations in a Euclidean space. However, the existence of such a diffeomorphism on every point of M implies that the manifold is globally flat, i.e., essentially Euclidean [START_REF] Iliev | Handbook of normal frames and coordinates[END_REF]; this assumption is far too restrictive for bona fide Riemannian applications, so the analysis of [START_REF] Shah | Stochastic approximation on riemannian manifolds[END_REF] is not relevant for our purposes. An additional issue is that the error bounds employed by Shah [61, p. 1131] rule out vector fields with a rotational componentsuch as V (x, y) = (-y, x) on R 2 -further limiting the applicability of their techniques to our setting.

Finally, the recent papers by Durmus et al. [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF][START_REF] Durmus | On Riemannnian stochastic approximation schemes with fixed step-size[END_REF] also consider a generic version of Robbins-Monro schemes, with both vanishing and constant step-sizes. The analysis of the latter type of schemes cannot lead to convergence with probability 1, so the results of [START_REF] Durmus | On Riemannnian stochastic approximation schemes with fixed step-size[END_REF] are necessarily ergodic in nature and hence out of our paper's scope. The setting of [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF] is closer in spirit to our own, and it also accounts for the effects of bias in the queries to V ; however, the results obtained therein concern dynamics that admit a Lyapunov functionthe so-called "gradient-like" case [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] -so there is no overlap with our analysis.

Notation and preliminaries

Throughout our paper, M will denote a d-dimensional, geodesically complete, smooth manifold equipped with a Riemannian metric ⟨•, •⟩ z with its induced norm ∥•∥ z ; we refer the reader to [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF] for standard definitions and notations (such as tangent spaces T z M). The Riemannian gradient will be simply denoted by ∇, and the Euclidean 2-norm is denoted by ∥•∥ 2 .

For any curve γ, the notation γ will always denote the derivative with respect to time. Given any points z, z ′ ∈ M and a vector v ∈ T z M, we denote by Γ z→z ′ (v) ∈ T z ′ M the vector obtained by parallel transporting v along the minimizing geodesic connecting z and z ′ ; if the minimizing geodesics are not unique, then Γ z→z ′ (v) is understood as the parallel transport along any of them.

We also say that a vector field

V on M is (geodesically) L-Lipschitz if, for all z, z ′ ∈ M, ∥Γ z→z ′ (V (z)) -V (z ′ )∥ z ′ = ∥V (z) -Γ z ′ →z (V (z ′ ))∥ z ≤ L dist(z, z ′ ),
where dist denotes the Riemannian distance induced by ⟨•, •⟩ z . All vector fields in this paper are assumed to be L-Lipschitz and bounded, i.e., G := sup z∈M ∥V (z)∥ z < ∞.

Riemannian Robbins-Monro schemes: Definition, dynamics, and convergence

We begin our analysis by introducing the core algorithmic template of Riemannian Robbins-Monro schemes. The material related to asymptotic pseudotrajectories is introduced in Section 3.2; Section 3.3 states all the assumptions and discusses their generality, and our main results are presented in Sections 3.4-3.5.

3.1.

The Riemannian Robbins-Monro template. As we stated before, the main idea behind the Riemannian Robbins-Monro (RRM) template is to replace addition along "straight lines" in Euclidean SA schemes with the Riemannian exponential mapping. This leads to the recursive update

Z n+1 = exp Zn (γ n (V (Z n ) + W n ))
(RRM) where (1) Z n ∈ M denotes the state of the algorithm at each iteration counter n = 1, 2, . . .

(2) W n ∈ T Zn M is an abstract error term (described in detail below).

(3) γ n > 0 is the method's step-size.

In the above, the error term W n is generated after Z n so W n is not adapted to the natural filtration F n := σ(Z 1 , . . . , Z n ) of Z n (i.e., W n is random relative to F n ). We will also write

Vn = V (Z n ) + W n (1)
so Vn can be seen as a noisy estimator of V (Z n ). To differentiate between "random" (zeromean) and "systematic" (non-zero-mean) errors, it will be convenient to further decompose W n as

W n = U n + b n ( 2
)
where

b n = E[W n | F n ]
represents the systematic component and U n = W n -b n captures the random, zero-mean part. To quantify all this, we will consider the following descriptors for W n : a) Bias:

B n = E[∥b n ∥ Zn | F n ] (3a) b) Variance: σ 2 n = E[∥U n ∥ 2 Zn |F n ] (3b)
Below, when B n = σ n = 0, we refer to the model as deterministic.

3.2. Mean dynamics and asymptotic pseudotrajectories. In analogy to the ODE method of stochastic approximation for ordinary, Euclidean Robbins-Monro schemes [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF], we will view (RRM) as a noisy Euler discretization of the Riemannian mean dynamics

θ(t) = V (θ(t)). (RMD)
To make this analogy precise, we will require a measure of "closeness" between the iterates of (RRM) and the integral curves of (RMD). To this end, let τ n = n-1 k=1 γ k denote the "effective time" that has elapsed till the n-th iteration of (RRM), and define the (continuous-time) geodesic interpolation Z(t) of Z n as

Z(t) = exp Zn t -τ n τ n+1 -τ n (γ n [V (Z n ) + W n ]) (Int) for all t ∈ [τ n , τ n+1 ], n ≥ 1 (so Z(τ n ) = Z n for all n).
To compare Z(t) to the solution orbits of (RMD), we will further consider the flow Θ : R + × M → M of (RMD), where Θ h (z) is simply the orbit of (RMD) at time h ∈ R + with initial condition z(0) = z ∈ M. We then have the following notion of "asymptotic closeness" due to Benaïm and Hirsch [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF]:

Definition 1 (Asymptotic pseudotrajectories). We say that Z(t) is ((a.s.)) an asymptotic pseudotrajectory (APT) of (RMD) if, for all T > 0, we have

lim t→∞ sup 0≤h≤T dist(Z(t + h), Θ h (Z(t))) = 0 almost surely. (4)
Intuitively, (4) states that, as t → ∞, one cannot distinguish between the interpolation Z(t + h) and the Riemannian mean dynamics with initial condition Z(t). The rigorous connection between (RMD) and (RRM) provided by this notion will become clear in Section 3.4.

3.3.

Technical assumptions and requirements. We now introduce the basic assumptions of our model and discuss their generality.

Step-size, noise and bias assumptions. Following the literature on stochastic approximation [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF], we will make the following standard assumptions for the step-size of (RRM): Assumption 1. The step-size sequence γ n , n = 1, 2, . . . of (RRM) satisfies the Robbins-Monro summability conditions n γ n = ∞ and n γ 2 n < ∞. Our second assumption concerns the error terms that appear in (RRM):

Assumption 2. The bounds on the noise and bias in (3a)-(3b) satisfy

B n → 0 and sup n σ 2 n < ∞ ((a.s.) ) (5)
Of the above, Assumption 1 is explicit and is controlled by the algorithm designer. By contrast, Assumption 2 depends on the primitives of the problem (the law providing measurements of V , the specific setup of (RRM), etc.), so it is more delicate; we discuss it in detail in Section 4, where we show that it is satisfied for a range of important instantiations of (RRM).

Boundedness and stability assumptions. To exclude problems where the iterates of stochastic approximation (SA) methods may escape to infinity, a standard practice in the literature is to make this into an explicit assumption. In our manifold setting, this can be stated as follows:

Assumption 3. For some (and hence, all) p ∈ M, the sequence Z n generated by (RRM) satisfies

sup n dist(p, Z n ) < ∞ with probability 1. (6)
Instead of simply imposing Assumption 3 as is commonly done in the literature [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF][START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF], we propose below a weaker condition that ensures [START_REF] Benaïm | Dynamics of Morse-Smale urn processes[END_REF], and which is general enough to cover most of the envisioned applications. To state it, we observe first that M typically falls under one of the following two categories:

(1) The manifold M is compact: This is the case for problems with simple manifold constraints such as spheres, balls, simplices, hypercubes, etc. Another important example is the Grassmannian manifold endowed with certain natural metrics; see [START_REF] Wong | Differential geometry of grassmann manifolds[END_REF][START_REF] Wong | Sectional curvatures of grassmann manifolds[END_REF].

(2) The sectional curvatures of M are non-positive: The most important such instances are matrix manifolds [START_REF] Lin | Riemannnian geometry of symmetric positive definite matrices via cholesky decomposition[END_REF][START_REF] Sra | Conic geometric optimization on the manifold of positive definite matrices[END_REF] and hyperbolic spaces [START_REF] Nickel | Poincaré embeddings for learning hierarchical representations[END_REF]; additional examples can be found in [START_REF] Bridson | Metric spaces of non-positive curvature[END_REF].

In the first case, Assumption 3 is trivially satisfied. For the second case, we propose the following explicit structural replacement:

Assumption 3 ′ . The sectional curvatures of M are non-positive and uniformly bounded from below by -K min > -∞, and the vector field V is weakly asymptotically coercive, i.e.,

V (z), ∇ dist 2 (p, z) ≤ 0 (WAC)
for some (and hence all) p ∈ M and for all z with sufficiently large dist(p, z). In addition, we assume that the noise and bias bounds in (3a)-(3b) satisfy

∞ n=1 E[γ n B n ] < ∞ and ∞ n=1 E γ 2 n B 2 n + σ 2 n < ∞. (7)
Remark 1. We shall see later in the proof of Proposition 2 that, for most algorithms, [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF] is weaker than (5) and hence can be discarded; for completeness, we also note that the differentiability of the radial function r 2 (z) = dist 2 (p, z) follows from Cartan's theorem [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF].

The inspiration for (WAC) comes from the Euclidean case where a standard criterion to ensure stability is the coercivity requirement lim z→∞ ⟨V (z), z⟩/∥z∥ 2 = -∞ [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Phelps | Convex Functions, Monotone Operators and Differentiability[END_REF]. This imposes a stringent growth condition on V which, mutatis mutandis, (WAC) relaxes significantly: under (WAC), it suffices if V does not have a persistent "outward-pointing" component.

In view of all this, our first result shows that, in simply connected manifolds, Assumption 3 can be replaced by Assumption 3 ′ (which, as we discussed, is much weaker in general):

Proposition 1. If M is simply connected, Assumption 3 ′ implies Assumption 3.
The most important step in the proof of Proposition 1 is to construct a suitable energy function which remains bounded under (WAC). However, due to the nonlinear structure of Riemannian manifolds, the resulting construction is fairly involved, so we defer it to Appendix A. Conjugacy assumptions. A key notion in the analysis to follow is the so-called Picard flow λ : R + → M associated with the interpolation Z(t), defined here as the solution of the following system of ODEs:

λ(h) = Γ Z (t+h)→λ(h) (V (Z(t + h))) (PFlow)
with initial condition λ(0) = Z(t). The term "Picard flow" stems from the fact that, in Euclidean spaces, the integral h 0 V (Z(t + s)) ds is the basic iteration in Picard's method of successive approximations for solving ODEs. In the case of (PFlow), the parallel transport is the extra ingredient required to express the idea of "integrating V along Z(t)", so (PFlow) can be seen as a natural generalization of the Picard iteration to Riemannian manifolds. Now, recall that two points z ∈ M and z ′ := exp z (v) are said to be conjugate along a minimizing geodesic if exp z (•) is minimizing but fails to be a locally diffeomorphism in a neighborhood of v [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF]. Define

C λ (t, T ) := {h ∈ [0, T ] | Z(t + h) is conjugate to λ(h)}, C θ (t, T ) := {h ∈ [0, T ] | Z(t + h) is conjugate to θ(h)},
and let C(t, T ) = C λ (t, T ) C θ (t, T ). We then make the following technical assumption: Assumption 4. For all t and T , C(t, T ) is a nowhere dense subset of [0, T ] with probability 1.

At first sight, Assumption 4 may appear somewhat opaque but, in practice, it is very mild and is satisfied by most applications. Indeed, since the set of points conjugate to Z(t + h) is at most one-dimensional [START_REF] Warner | The conjugate locus of a riemannian manifold[END_REF], the only way to violate Assumption 4 is if the curves Z(t + •) and θ(•) (or λ(•)) simultaneously traverse the same one-dimensional submanifold, which happens with probability 0 if the distribution of the noise U n is non-singular. It is also worth noting that Assumption 4 is trivially true on negatively curved spaces by the Cartan-Hadamard theorem [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF] and, finally, it is straightforward to verify Assumption 4 manually on the manifolds that arise most often in practical applications (such as spheres, balls, Grassmannians, or fixed-rank spectrahedra), cf. [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF], [START_REF] Massart | Curvature of the manifold of fixed-rank positive-semidefinite matrices endowed with the bures-wasserstein metric[END_REF] and references therein.

3.4.

From discrete to continuous: RRM schemes and APTs. We are finally in a position to state and discuss our main results concerning the convergence properties of (RRM). Our first result describes the way in which (RRM) can be approximated by (RMD): Theorem 1. Suppose that Assumptions 1-4 hold. Then, with probability 1, the geodesic interpolation Z(t) of the sequence of iterates of (RRM) is an APT of (RMD).

The proof of Theorem 1 is fairly involved and the geometric scaffolding required is quite delicate, so we first provide here a high-level outline. In brief, the main obstacle that we The moving observer Z(t + h) measures a curve λ(h) using time-indexed geodesics and an "inertial frame system", i.e., frames obtained via parallel transport from Z(t). For each

h, λ(h) ∈ R d is the normal coordinate of λ(h) ∈ M. The space-time map Φ : R + × M → R + × R d is locally defined on a neighborhood containing Z(t + h) and λ(h).
have to overcome in order to extend Euclidean APT theory to the manifold setting is the following dichotomy:

• In view of definition (4), we need a coordinate system to compare the distance between two points. This is typically done in normal coordinates [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF]; on the other hand, comparing vectors that belong to different tangent spaces in normal coordinates is a formidable task.

• Since the flow Θ t (z) is driven by a vector field, we will also need to compare vectors living on different tangent spaces. To this end, the most convenient framework is the parallel frame system, described in detail in Appendix B.1. Unfortunately, the parallel frame system cannot coexist with a normal coordinate system unless M is flat, i.e., the manifold's Riemannian curvature tensor vanishes everywhere [START_REF] Iliev | Handbook of normal frames and coordinates[END_REF].

In words, the normal coordinate system is where distances behave as if the space were Euclidean, and the parallel frame system is where vectors behave as in the Euclidean case; unfortunately, the only setting where these two descriptions come together is in flat manifolds, which are "almost" Euclidean to begin with. To circumvent this hurdle, we take the following technical path:

(1) Based on the idea of Fermi coordinates [START_REF] Manasse | Fermi normal coordinates and some basic concepts in differential geometry[END_REF], which can be intuitively understood as "normal coordinates along a curve" (see Fig. 1), we begin by constructing an extended Fermi coordinate frame that allows us to focus on a neighborhood of Z(t) containing all the information we need. [This is a challenging but otherwise purely technical step and can be safely omitted without losing the high-level picture.]

(2) Using the extended Fermi coordinates constructed above, we can reduce the task of comparing the distance between two Riemannian curves to comparing several Euclidean, albeit individually intractable, vector fields. This step will incur an error term as compared to the Euclidean analysis; see [START_REF] Criscitiello | An accelerated first-order method for non-convex optimization on manifolds[END_REF] and [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF].

(3) To obtain expressions of vector fields that are more amenable to computation, we will switch from the extended Fermi coordinates to the parallel frame system and bound the difference between the two. This step will pick up another error term as compared to the Euclidean analysis, cf. Eqs. ( 27) and ( 28).

(4) Serendipitously -and, perhaps, surprisingly -these additional error terms can be managed without any further assumptions, and a more involved analysis in the spirit of [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] concludes our proof.

Proof of Theorem 1. We will need the following ingredients:

• The concept of parallel frame system, reviewed in Appendix B.1.

• A key lemma by Fujita and Kotani [START_REF] Fujita | The onsager-machlup function for diffusion processes[END_REF] and Takahashi and Watanabe [START_REF] Takahashi | The probability functionals (onsager-machlup functions) of diffusion processes[END_REF], which bounds the distortion of velocities measured by a moving observer on a manifold relative to flat spaces; see Lemma B.1.

• A comparison lemma between the differential of exponential map and parallel transport, recited in Lemma B.2.

We now commence the proof. For ease of reading, we break down our proof into several steps:

Step 1: Discrete-to-continuous transformations; noise stability. Recall the "effective time" τ n = n-1 k=1 γ k as the time that has elapsed at the n-th iteration of the discrete-time process Z n ; recall also the definition (Int) of the continuous-time interpolation Z(t) of Z n as

Z(t) = exp Zn t -τ n τ n+1 -τ n (γ n [V (Z n ) + W n ]) . (Int)
We will further require the "continuous-to-discrete" correspondence

m(t) = sup{n ≥ 1 : t ≥ τ n } (8)
which measures the number of iterations required for the effective time τ n of the process to reach the timestamp t.

Given an arbitrary sequence A n (which can be numbers or points/vectors, either Euclidean or Riemannian), we will denote its piecewise-constant interpolation as

A(t) = A n for all t ∈ [τ n , τ n+1 ), n ≥ 1. (9)
Using this notation, we can express (Int) in a differential form:1

Ż(t + h) = Γ Z(t+h)→Z(t+h) V Z(t + h) + W (t + h) . (10)
Recall also the definition (RMD) and (PFlow):

θ(h) = V (θ(h)), (Flow)
λ(h) = Γ Z (t+h)→λ(h) (V (Z(t + h))) , (PFlow)
both with initial condition Z(t). We also set

γ ⋆ (t) := sup t≤h≤t+T γ(h), B ⋆ (t) := sup t≤h≤t+T B(h). (11)
We will also need a noise stability criterion. Let {e k (n)} d k=1 be an arbitrary sequence of orthonormal basis for T Zn M, and let U n, be the (Eucldiean) noise vector composed of components of the noise U n in the basis {e

k (n)} d k=1 : U k n, := ⟨U n , e k (n)⟩ Zn .
Then it is easy to see that E[U n, |F n ] = 0, and

E ∥U n, ∥ 2 2 |F n = E ∥U n ∥ 2 Zn |F n ≤ sup n σ 2 n < ∞ by Assumption 2. Therefore, {U n, } ∞
n=1 is a sequence of martingale Euclidean noise vectors that is with zero mean and finite variance. The convergence of such noise sequences is a well-studied subject in stochastic approximation: Define

∆(t; T ) := sup 0≤h≤T t+h t U (u) du 2 . ( 12
)
It is classical (see e.g., [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF]) that For all T, lim t→∞ ∆(t; T ) = 0 (a.s.)

Step 2: Everything is bounded. We first note that, by Proposition 1, sup t r(Z(t)) < ∞ where r(•) is the radial function defined in Assumption 3 ′ . We claim that Proposition 1 also implies the boundedness of the Picard flow. To see this, simply note that, since the parallel transport is an isometry,

λ(h) λ(h) = ∥Γ Z (t+h)→λ(h) (V (Z(t + h)))∥ λ(h) = ∥V (Z(t + h))∥ Z(t+h) , so that sup 0≤h≤T dist(λ(0), λ(h)) ≤ T • sup 0≤h≤T λ(h) λ(h) < ∞,
which implies sup 0≤h≤T r(λ(h)) < ∞. On the other hand, the boundedness for the flow follows readily from the L-Lipschitzness of V and the 1-Lipschitzness of dist:

2 d dh r(θ(h)) ≤ θ(h) θ(h) = ∥V (θ(h))∥ θ(h) ≤ ∥V (θ(h)) -Γ p→θ(h) (V (p))∥ θ(h) + ∥Γ p→θ(h) (V (p))∥ θ(h) ≤ Lr(θ(h)) + ∥V (p)∥ p .
An application of the Gronwall's inequality then concludes sup 0≤h≤T r(θ(h)) < ∞. Therefore, all computations in the sequel can be restricted to a compact set. We may thus further assume that the sectional curvatures K are bounded. Since γ n → 0, for t large enough we may assume γ ⋆ (t) < 1.

To summarize, for some

-∞ < K low ≤ K up < ∞, we have γ ⋆ (t) < 1, K low ≤ K ≤ K up . (14)
In addition, there exists an R = R(T, L, G) depending on T , L, and G such that

sup 0≤h≤T dist(Z(t + h), θ(h)), dist(Z(t + h), λ(h)) ≤ R < ∞. ( 15
)
2 Due to the smoothness of the flow, the function r(θ(h)) is always differentiable in h in the metric space sense [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF], even though θ(h) might reach the cut locus of p.

Step 3: Constructing the extended Fermi coordinates. Let V be the neighborhood defined in Appendix B.2 restricted to [t, t + T ] (i.e., V = T h=0 U h where U h is defined in (B.4)); evidently, V contains {Z(t + h) : h ∈ [0, T ]}. Recall the Fermi coordinate system Φ in (B.5) restricted to [t, t + T ]. Recall also that Φ carries a system of orthonormal frames {e k (h)} d k=1 , one for each Z(t + h). Below, we will always work in these frames with no explicit mention.

Fix ȟ ∈ [0, T ]. Let γ θ and γ λ be two minimizing geodesics such that γ θ (0) = γ λ (0) = Z(t + ȟ), γ θ (1) = θ( ȟ) and γ λ (1) = λ( ȟ). Our first goal is to extend V to an open set of M that contains the geodesics γ θ and γ λ , while retaining the exponential mapping as local diffeomorphisms. This will serve a dual purpose:

(1) It enables us, for a fixed ȟ, to consider the parallel frame systems at Z(t + ȟ) and θ( ȟ) (or Z(t + ȟ) and λ( ȟ)), so that we can easily compare the vector fields at these points; see Appendix B.1. (2) We want to apply Lemma B.1 to the curves θ(•) and λ(•); in order to make sense of θ( ȟ) or λ( ȟ), V needs to contain both curves for at least an open time interval that includes ȟ. This is where Assumption 4 comes into play: Since the conjugate points are precisely where the exponential map ceases to be local diffeomorphisms [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF], it is reasonable to expect that, away from the time points where Z(t + h) is conjugate to θ(h) or λ(h), it is always possible to extend V to include the geodesics connecting Z(t + h) to θ(h) and λ(h). Assumption 4 then simply posits that there cannot be "too many" such conjugate time points.

Consider any ȟ / ∈ C(t, T ) where C(t, T ) is defined in (4), and assume also that t + ȟ ̸ = τ n for all n (i.e., t + ȟ is not a "corner" of (Int)). Since the exponential mapping is a local diffeomorphism away from conjugate points [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF], it follows that exp Z(t+ ȟ) (•) : T Z (t+ ȟ) M → M is a local diffeomorphism at λ( ȟ) and θ( ȟ), where λ( ȟ) and θ( ȟ) are the normal coordinates of λ( ȟ) and θ( ȟ) with center Z(t + ȟ). By continuity of the flow/Picard flow and frame system {e k (•)}

H = k (h k , h k+1 ).
The first claim follows readily from Assumption 4 and the fact that the set { ȟ : t + ȟ = τ n for some n} is countable. To prove the second claim, simply note that, since all open intervals of R contain at least one rational number, it is impossible for an open set to be the union of uncountably many disjoint open intervals.

To summarize, we can extend V and Φ to an open set including

{γ h θ } h∈(h k ,h k+1 ) {γ h λ } h∈(h k ,h k+1 ) ,
which obviously contains h∈(h k ,h k+1 ) {θ(h), λ(h)}. We call such a pair (V, Φ) the extended Fermi coordinate, since it not only contains the central curve h → Z(t + h) as in the classical case, but also θ(h) and λ(h) for almost every h ∈ [0, T ].

Step 4: Controlling the distance by decomposition. From now on, all computation happens within the scope of the extended Fermi coordinate (V, Φ). By definition of θ(•), and since θ(h) is the normal coordinate of θ(h) with center Z(t + h), we have, for all h ∈ H,

dist(Z(t + h), Θ h (Z(t))) = dist(Z(t + h), θ(h)) = ∥ θ(h)∥ 2 ≤ ∥ θ(h) -λ(h)∥ 2 + ∥ λ(h)∥ 2 . ( 16
)
Since H is a dense open subset of [0, T ] and since it is a countable union of open intervals, it follows that θ(h) and λ(h) are differentiable except for a measure zero set. We can thus write

∥ θ(h) -λ(h)∥ 2 = h 0 θ(u) -λ(u) du 2 .
By Lemma B.1, (RMD), and (PFlow), we have

θk (u) = θk (u) -Żk (t + u) + O ∥ θ(u)∥ 2 2 = Ṽ i (u, θ(u)) -Żk (t + u) + O ∥ θ(u)∥ 2 2 , λk (u) = λk (u) -Żk (t + u) + O ∥ λ(u)∥ 2 2 = Λk u, λ(u) -Żk (t + u) + O ∥ λ(u)∥ 2 2 ,
where Ṽ k (u, θ(u)) and Λk u, λ(u) are, respectively, the k-th components of the vectors V θ(u) and Γ Z(t+u)→λ(u) (V (Z(t + u))) in the frame induced by the normal coordinate with center Z(t + u) and frame {e k (u)} d k=1 , and Żk (•) is defined in Lemma B.1. Denoting by Ṽ (u, θ(u)) the (Euclidean) vector with components Ṽ i (u, θ(u)) and, by Λ u, λ(u) the vector with components Λk u, λ(u) , we may thus write

∥ θ(h) -λ(h)∥ 2 ≤ h 0 Ṽ (u, θ(u)) -Λ u, λ(u) du 2 + h 0 R 1 (u)du (17)
where the remainder term

R 1 (u) is of order O ∥ θ(u)∥ 2 2 + ∥ λ(u)∥ 2 2 .
Noting that, by [START_REF] Chen | Modified tseng's extragradient methods for variational inequality on Hadamard manifolds[END_REF],

∥ θ(u)∥ 2 = dist (Z(t + u), Θ u (Z(t))) ≤ R, ∥ λ(u)∥ 2 = dist (Z(t + u), λ(t + u)) ≤ R, we have ∥ θ(u)∥ 2 2 +∥ λ(u)∥ 2 2 ≤ R ∥ θ(u)∥ 2 + ∥ λ(u)∥ 2 , and hence R 1 (u) = O R ∥ θ(u)∥ 2 + ∥ λ(u)∥ 2 where O R (•) hides constants depending on R.
In the same vein, denoting by Ż(t + u) the Euclidean vector whose k-th component is Żk (t + u), we have

∥ λ(h)∥ 2 ≤ h 0 Λ u, λ(u) -Ż(t + u) du 2 + h 0 R 2 (u)du (18)
where

R 2 (u) = O ∥ λ(u)∥ 2 2 = O R ∥ λ(u)∥ 2 .
Step 5: Switching coordinate systems: From Fermi to parallel. So far, we have reduced the proof to comparing between the vectors in ( 17) and [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF]. However, these vectors are not amenable to further computation as they are expressed in the frames induced by the normal coordinates, and these frames may not even be orthonormal.

On the other hand, when expressed in the parallel frame system (see Appendix B.1) with a common base point Z(t + u), the vectors V θ(u) and Γ Z (t+u)→λ(u) (V (Z(t + u))) possess some favorable properties. To see this, parallel transport the frame {e k (u)} d k=1 along the geodesic from Z(t + u) to form an orthonormal frame {e ′ k (u)} d k=1 of T λ(u) M, and consider the Euclidean vector Λ u, λ(u) whose components are defined as

Λk u, λ(u) := ⟨Γ Z (t+u)→λ(u) (V (Z(t + u))) , e ′ k (u)⟩ λ(u) = ⟨V (Z(t + u)), e k (u)⟩ Z(t+u) . (19)
Similarly, parallel transport the frame {e k (u)}

Since the parallel transport is an isometry, we get

Ṽ k u, θ(u) -Λk u, λ(u) = V θ(u) , e ′′ k (u) θ(u) -⟨Γ Z (t+u)→λ(u) (V (Z(t + u))) , e ′ k (u)⟩ λ(u) = Γ θ(u)→Z (t+u) V θ(u) , e k (u) Z(t+u) -⟨V (Z(t + u)), e k (u)⟩ Z(t+u) = Γ θ(u)→Z (t+u) V θ(u) -V (Z(t + u)), e k (u) Z(t+u) , whence Ṽ u, θ(u) -Λ u, λ(u) 2 = ∥Γ θ(u)→Z (t+u) V θ(u) -V (Z(t + u))∥ Z(t+u) ≤ L dist(θ(u), Z(t + u)) = L∥ θ(u)∥ 2 (21)
where we have used the L-Lipschitzness of V and the fact that θ(u) is the normal coordinate with center Z(t + u).

We would like to therefore replace Ṽ u, θ(u

) -Λ u, λ(u) with Ṽ u, θ(u) -Λ u, λ(u)
in [START_REF] Criscitiello | An accelerated first-order method for non-convex optimization on manifolds[END_REF]. To this end, consider the difference in the k-th component of Ṽ u, θ(u) and Ṽ u, θ(u) :

Ṽ k u, θ(u) -Ṽ k u, θ(u) = V θ(u) , e ′′ k (u) θ(u) -V θ(u) , ∂ ∂z k θ(u) θ(u) (22)
where

∂ ∂zi θ(u)
is the k-th basis in the frame induced by the normal coordinate with center

Z(t + u)
Then

∂ ∂z i θ(u) = ∂ ∂s ′ γ(1, 0) = d exp Z(t+u) (v)(e k (u)). (24)
Invoking Cauchy-Schwartz, ( 14), Lemma B.2, and using the fact that {e k (u)} d k=1 is orthonormal, we get

Ṽ k u, θ(u) -Ṽ k u, θ(u) = V θ(u) , e ′′ k (u) - ∂ ∂z k θ(u) θ(u) ≤ ∥V θ(u) ∥ θ(u) • e ′′ k (u) - ∂ ∂z k θ(u) θ(u) ≤ G • K max • f -Kmax (v) • ∥e k (u) ⊥ ∥ Z(t+u) ≤ G • K max • f -Kmax (v) = G • K max • f -Kmax (∥ θ(u)∥ 2 ) (25)
where

K max = max (|K up |, |K low |).
Since a standard Taylor expansion argument shows that sinh(x) x -1 ≤ cosh(x) ≤ e x for all x ≥ 0, by (B.7), (B.8), and ( 14), we get

f -Kmax (∥ θ(u)∥ 2 ) ≤ 1 K max exp K max ∥ θ(u)∥ 2 ≤ e R √ Kmax R • K max • ∥ θ(u)∥ 2 (26)
where the last inequality follows from ∥ θ(u)∥ 2 = dist(Z(t + u), Θ u (Z(t))) ≤ R. Combining ( 25) and ( 26), we thus get

Ṽ k u, θ(u) -Ṽ k u, θ(u) ≤ Ge R √ Kmax R • ∥ θ(u)∥ 2 .
In short, we have shown that

Ṽ u, θ(u) = Ṽ u, θ(u) + R 3 (u). (27)
Here,

R 3 (u) = O G,Kmax,R ∥ θ(u)∥ 2 ,
where O G,Kmax,R (•) hides constants that depend on G, K max and R. Exactly the same computation shows that, for some

R 4 (u) = O G,Kmax,R ∥ λ(u)∥ 2 , Λ u, λ(u) = Λ u, λ(u) + R 4 (u). ( 28
)
Step 6: Concluding the proof. We will proceed by bounding ( 17) and [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF]. Using ( 27), [START_REF] Iliev | Handbook of normal frames and coordinates[END_REF], and ( 21) in [START_REF] Criscitiello | An accelerated first-order method for non-convex optimization on manifolds[END_REF], we obtain:

∥ θ(h) -λ(h)∥ 2 ≤ h 0 Ṽ u, θ(u) -Λ u, λ(u) du 2 + h 0 R 1 (u)du ≤ h 0 Ṽ u, θ(u) -Λ u, λ(u) 2 du + h 0 R 1 (u) du ≤ h 0 Ṽ u, θ(u) -Λ u, λ(u) 2 du + h 0 (R 1 + R 3 + R 4 )(u) du ≤ L h 0 θ(u) 2 du + h 0 (R 1 + R 3 + R 4 )(u) du. ( 29
)
We next turn to [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF]. Our first task is to obtain an expression for Ż(t + u); recall that this is the Euclidean vector whose k-th component is Żk (t + u). To this end, fix an iteration count n, and consider all u such that t + u ∈ [τ n , τ n+1 ) (that is, consider only the interpolation between Z n and Z n+1 ). We claim that Żk (t + u) is constant throughout all such u, and, in particular, Żk (t + u) = Żk (τ n ). This readily follows by noticing that 1) The curve Z(•) is a geodesic segment when restricted to [τ n , τ n+1 ); see (Int).

2) The Fermi coordinate along Z(•), when restricted to {Z(s)

: s ∈ [τ n , τ n+1 )}, is simply a parallel frame system: The frame {e k (u)} d k=1 is obtained from parallel transporting {e k (τ n )} d k=1 along Z(•) for all u such that t + u ∈ [τ n , τ n+1
). Thus, a simple calculation akin to (B.2) yields, for all such u,

Żk (t + u) = ⟨ Ż(t + u), e i (u)⟩ Z(t+u) = Γ Z (t+u)→Z (τn ) Ż(t + u) , e i (τ n ) Z(τn) = ⟨V (Z(τ n )) + W n , e i (τ n )⟩ Z(τn) = ⟨V (Z(τ n )) + U n + b n , e i (τ n )⟩ Z(τn) (30)
where we have used [START_REF] Bonnabel | Stochastic gradient descent on Riemannnian manifolds[END_REF] and the definition of W n in the last equality.

Armed with the above, we can obtain a succinct expression for Ż(t+u) as follows. Let z(u) be the normal coordinate of Z(t + u) with center Z(t + u) (i.e., (u, z(u)) = Φ(u, Z(t + u))), and define a Euclidean vector Ṽ (u, z(u)) by setting its k-th component to

Ṽ k (u, z(u)) := V Z(t + u) , e i (m(t + u)) Z(t+u) (31)
where the mapping m(•) is defined in [START_REF] Bento | Iteration-complexity of gradient, subgradient and proximal point methods on Riemannnian manifolds[END_REF]. Define also the Euclidean noise and bias vectors U n, and b n, by setting their components to

U k n, := ⟨U n , e i (m(t + u))⟩ Z(t+u) , b k n, := ⟨b n , e i (m(t + u))⟩ Z(t+u) .
Then [START_REF] Kakade | A natural policy gradient[END_REF] states precisely that

Ż(t + u) = Ṽ (u, z(u)) + U (t + u) + b (t + u). (32)
Substituting ( 32) into [START_REF] Durmus | Convergence analysis of Riemannnian stochastic approximation schemes[END_REF] and invoking ( 28), [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF], and ( 12), we obtain

∥ λ(h)∥ 2 ≤ h 0 Λ u, λ(u) -Ṽ (u, z(u)) -U (t + u) -b (t + u) du 2 + h 0 R 2 (u)du ≤ h 0 Λ u, λ(u) -Ṽ (u, z(u)) 2 du + h 0 U (t + u) du 2 + h 0 b (t + u) du 2 + h 0 (R 2 + R 4 )(u)du ≤ h 0 Λ u, λ(u) -Ṽ (u, z(u)) 2 du + ∆(t, T ) + B ⋆ (t) • h + h 0 (R 2 + R 4 )(u)du. (33)
To bound the first term in [START_REF] Klavins | Phase regulation of decentralized cyclic robotic systems[END_REF], recall ( 19) and ( 31). An identical argument leading to (B.3) shows that

Λ u, λ(u) -Ṽ (u, z(u)) 2 = V (Z(t + u)) -Γ Z(t+u)→Z(t+u) V (Z(t + u)) Z(t+u) ≤ L • dist Z(t + u), Z(t + u) . (34)
Since Z(•) is a (not necessarily minimizing) geodesic on [m(t + u), t + u], we have, by ( 32), [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF], and (

ds 2 ≤ G + B ⋆ (t) • t+u m(t+u) ds + t+u m(t+u) U (s) ds 2 ≤ G + B ⋆ (t) γ ⋆ (t) + t+u m(t+u) U (s) ds 2 . ( 35
)
For t large enough, we have γ ⋆ (t) < 1, and hence

t+u m(t+u) U (s) ds 2 ≤ m(t+u) t-1 U (s) ds 2 + t+u t-1 U (s) ds 2 ≤ 2∆(t -1, T + 1). (36)
Combing ( 34), [START_REF] Kristály | Nash-type equilibria on riemannian manifolds: a variational approach[END_REF], and ( 36), (33) then becomes

∥ λ(h)∥ 2 ≤ L • h • G + B ⋆ (t) γ ⋆ (t) + 2∆(t -1, T + 1) + ∆(t, T ) + B ⋆ (t) • h + h 0 (R 2 + R 4 )(u)du ≤ h • C L,G B ⋆ (t) + γ ⋆ (t) + ∆(t -1, T + 1) + h 0 (R 2 + R 4 )(u)du (37)
for some constant C L,G that depends only on L and G. Using ( 29) and (37), we can then bound [START_REF] Criscitiello | Efficiently escaping saddle points on manifolds[END_REF] as

dist(Z(t + h), Θ h (Z(t))) = ∥ θ(h)∥ 2 ≤ ∥ θ(h)∥ 2 + ∥ λ(h)∥ 2 (38)
≤ ∥ θ(h) -λ(h)∥ 2 + 2∥ λ(h)∥ 2 ≤ L h 0 θ(u) 2 du + h • 2C L,G B ⋆ (t) + γ ⋆ (t) + ∆(t -1, T + 1) + 3 h 0 (R 1 + R 2 + R 3 + R 4 )(u) du (39)
where

(R 1 + R 2 + R 3 + R 4 )(u) = O G,Kmax,R θ(h) 2 + λ(h) 2
. Therefore, there exists a constant C L,G,Kmax,R depending only on L, G, K max , and R such that we may bound [START_REF] Kushner | Asymptotic properties of stochastic approximations with constant coefficients[END_REF] as

∥ θ(h)∥ 2 + ∥ λ(h)∥ 2 ≤ C L,G,Kmax,R h 0 θ(h) 2 + λ(h) 2 du + h • 2C L,G B ⋆ (t) + γ ⋆ (t) + ∆(t -1, T + 1) .
Grönwall's inequality then implies

∥ θ(h)∥ 2 + ∥ λ(h)∥ 2 ≤ h • 2C L,G B ⋆ (t) + γ ⋆ (t) + ∆(t -1, T + 1) • e h•C L,G,Kmax,R . (40)
From ( 40), we may conclude the proof as follows:

lim t→∞ sup t≤h≤T dist(Z(t + h), Θ h (Z(t))) ≤ lim t→∞ sup t≤h≤T ∥ θ(h)∥ 2 + ∥ λ(h)∥ 2 ≤ lim t→∞ T • 2C L,G B ⋆ (t) + γ ⋆ (t) + ∆(t -1, T + 1) • e T •C L,G,Kmax,R = 0 a.s.
since lim t→∞ B ⋆ (t) = lim t→∞ γ ⋆ (t) = 0 a.s. by assumption, and lim t→∞ ∆ (t -1, T + 1) = 0 a.s. by [START_REF] Boumal | Global rates of convergence for nonconvex optimization on manifolds[END_REF]. ■ 3.5. The Limits of Riemannian Robbins-Monro schemes. Now, having established that the iterates of (RRM) constitute an APT of (RMD), we proceed to explore the formal link between asymptotic pseudotrajectories and the solution orbits of (RMD). To this end, we will need the important notion of an internally chain-transitive (ICT) set.

Definition 2. Let S be a nonempty compact subset of M and let Θ be a flow on M. Then:

(1) S is invariant if Θ t (S) = S for all t ∈ R, i.e., Θ t (z) ∈ S for all z ∈ S.

(2) S is an attractor of Θ if it is invariant and there exists a compact neighborhood K of S such that lim t→∞ dist(Θ t (z), S) = 0 uniformly in z ∈ K. (3) S is internally chain-transitive if it is invariant and Θ| S admits no other attractors than S.

In words, ICT sets can be seen as a "terminal object" for the dynamics (RMD): the ensemble of orbits that converges to an ICT set will not be ultimately contained in a smaller subset thereof. Our next result shows that, with probability 1, all limit points of (RRM) lie in some such set: Theorem 2. If Assumptions 1-4 hold, then Z n converges almost surely to an ICT set of (RMD).

The proof of Theorem 2 leverages Theorem 1 in an essential way and employs a foundational result of Benaïm and Hirsch [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF] to establish the link with the ICT sets of (RMD):

Proof of Theorem 2. By Theorem 1, Z n generates APTs of the mean dynamics (RMD). Now, let L = t≥0 cl(Z(t, ∞)) be the limit set of Z(t), i.e., the set of limit points of convergent sequences Z(t n ) with lim n t n = ∞. Our claim then follows by the limit set theorem of Benaïm and Hirsch [7, Theorem 8.2], which holds for an arbitrary metric (even non-Riemannian) space.

■

An important message of Theorem 2 is that we can reduce the convergence analysis of (RRM) to the study of a deterministic, continuous-time dynamical system which is significantly simpler than the original stochastic, discrete-time system. In this sense, Theorem 2 delivers the same high-level message as the classic SA literature for Euclidean Robbins-Monro (RM) schemes. In particular, if V admits a potential or is strictly monotone, it is easy to verify that the only ICT sets of (RMD) are the roots of V [63], so we readily recover the series of asymptotic convergence results mentioned in the introduction. In the general case, the roots of V are always contained in its ICT sets, but the inclusion may be strict; however, this again becomes an equivalence under the standard assumptions which have been used in the Euclidean RM literature to obtain global convergence -cooperative dynamics, global asymptotic stability, etc. [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF][START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF]. In the next section, we shall see that Theorems 1-2 further capture a series of Riemannian algorithms -old and new -in a unified fashion.

Applications and implications

We now proceed to show how a wide array of Riemannian algorithms can be seen as special cases of (RRM) -enabling in this way the tandem use of Theorems 1 and 2 to deduce their convergence properties. To simplify the presentation, we will make the following technical assumption: Assumption 5. The injectivity radius of M is uniformly bounded from below by δ > 0.

The injectivity radius at a point p on M is the largest radius for which the exponential map at p is a diffeomorphism, and the injectivity radius of M is the infimum over all such radii [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF]. In this regard, Assumption 5 serves to ensure that exp is invertible at consecutive iterates of (RRM) so as to avoid local topological complications; we will later prove in Proposition 2 that exp -1 is well-defined in all the algorithms below. 4.1. Specific algorithms. Throughout this section, we adopt a black-box setup [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF] with stochastic first-order oracle (SFO) feedback. Specifically, when called at z ∈ M with random seed ω ∈ Ω, an SFO returns a random vector V (z; θ) ∈ T z M of the form

V (z; θ) = V (z) + Err(z; θ) (SFO)
where the error term Err(z; θ) ∈ T z M is zero-mean and finite-variance with respect to the Riemannian metric:

∀z ∈ M, E[Err(z; θ)] = 0 and E[∥Err(z; θ)∥ 2 z ] ≤ σ 2 . ( 41
)
To facilitate comparisons with the Riemannian optimization literature, we will abusively refer to queries of V as "gradients"; we stress, however, that V is not assumed to admit a potential.

Algorithm 1 (Riemannian stochastic gradient methods). The simplest Riemannian stochastic gradient method (RSGM) queries an SFO and proceeds as:

Z n+1 = exp Zn (γ n V (Z n ; θ n )) , (R-SGM)
As such, (R-SGM) admits a straightforward RRM representation by taking W n = Err(Z n ; θ n ). ¶ Algorithm 2 (Riemannian proximal point methods). The (deterministic) Riemannian proximal point method (RPPM) is an implicit update rule of the form:

exp -1 Z n+1 (Z n ) = -γ n V (Z n+1 ). (R-PPM)
The RRM representation of (R-PPM) is obtained by taking

W n = U n + b n with b n = Γ Z n+1 →Zn (V (Z n+1 )) -V (Z n )
and U n = 0 (recall here that Γ z→z ′ (v) denotes parallel transport along minimizing geodesics). ¶

Algorithm 3 (Riemannian stochastic extra-gradient). Inspired by the original work of Korpelevich [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF], the Riemannian stochastic extra-gradient (RSEG) method iterates as

Z n+1/2 = exp Zn (γ n V (Z n ; θ n )) , Z n+1 = exp Zn Γ Z n+1/2 →Zn γ n V (Z n+1/2 ; θ n+1/2 ) (R-SEG)
where θ n and θ n+1/2 are independent random seeds for (SFO). To recast (R-SEG) in the RRM framework, simply take

U n = Γ Z n+1/2 →Zn Err(Z n+1/2 ; θ n+1/2 ) and b n = Γ Z n+1/2 →Zn V (Z n+1/2 ) - V (Z n ). ¶
Algorithm 4 (Riemannian optimistic gradient). Compared to (R-SGM), the scheme (R-SEG) involves two oracle queries per iteration. Inspired by Popov [START_REF] Popov | A modification of the Arrow-Hurwicz method for search of saddle points[END_REF], we can instead "recycle" the last oracle query, leading to the Riemannian optimistic gradient method

Z n+1/2 = exp Zn γ n V (Z + n-1 ; θ n-1 ) , Z n+1 = exp Zn Γ Z n+1/2 →Zn γ n V (Z n+1/2 ; θ n ) .
(R-OG) (R-OG) can be seen as a special case of (RRM) by taking

U n = Γ Z n+1/2 →Zn Err(Z n+1/2 ; θ n ) and b n = Γ Z n+1/2 →Zn V (Z n+1/2 ) -V (Z n ). ¶
In light of Theorems 1-2, to analyze the convergence of Algorithms 1-4 under a specific choice of step-size, we can simply verify Assumptions 1-4. More precisely, it suffices to show that ( 5) and ( 7) hold with probability 1. This is the purpose of the next proposition: Proposition 2. Under Assumption 5 and a step-size strategy such that A/n ≤ γ n ≤ B/ n(log n) 1+ε for some A, B, ε > 0, the following hold for Algorithms 1-4:

(1) With probability 1, Z n+1 lies in the injectivity radius of Z n for n large enough.

(2) Equations ( 5) and (7) hold with probability 1. As a result, under Assumptions 3 ′ -4, Algorithms 1-4 generate APTs of (RMD) and converge to its ICT sets.

Before proving Proposition 2, we first provide a convenient lemma which shows that, almost surely, the effect of the noise is asymptotically annihilated by the step-size: Lemma 1. Under the assumptions in Proposition 2, we have, with probability 1,

lim n→∞ ∥γ n V (Z n ; θ n )∥ Zn = 0. ( 42
)
Proof of Lemma 1. By definition,

∥γ n V (Z n ; θ n )∥ Zn = ∥γ n V (Z n ) + Err(Z n ; θ n )∥ Zn ≤ γ n G + γ n ∥Err(Z n ; θ n )∥ Zn . ( 43
)
The first term goes to 0 by choice of step-sizes. To control the second term, note that by Chebyshev's inequality and ( 41), we have

P ∥Err(Z n ; θ n )∥ Zn ≥ n log 1+ ϵ 2 n ≤ σ 2 n log 1+ ϵ 2 n ( 44
)
where ε is the same as in our choice of step-size in Proposition 2. In turn, this implies that

∞ n=2 P ∥Err(Z n ; θ n )∥ Zn ≥ n log 1+ ϵ 2 n < ∞ so, by the Borel-Cantelli lemma, we have ∥Err(Z n ; θ n )∥ Zn = O n log 1+ ϵ 2 n with proba- bility 1.
Hence, by our assumptions for the method's step-size, we get

γ n ∥Err(Z n ; θ n )∥ Zn = O   n log 1+ ϵ 2 n n log 1+ε n   = O 1 log ε 4 n
which, combined with [START_REF] Li | Monotone vector fields and the proximal point algorithm on Hadamard manifolds[END_REF], implies [START_REF] Lezcano-Casado | Curvature-dependant global convergence rates for optimization on manifolds of bounded geometry[END_REF]. ■

We are now ready to prove Proposition 2.

Proof of Proposition 2. The first claim of Proposition 2 is a direct consequence of Lemma 1 and Assumption 5. To prove the second claim, note that n γ 2 n < ∞ by our choice of step-sizes. Thus, in order to prove [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF], it suffices to show

E[B n ] = E[∥b n ∥] = O(γ n ) and E[σ 2
n ] ≤ σ 2 for some constant σ.

We next proceed method-by-method: Algorithm 1: Riemannian stochastic gradient descent(/ascent). For (R-SGM), we have W n = U n = Err(Z n ; θ n ) and b n = 0, so both ( 5) and ( 7) follow from the stated assumptions for (SFO). Algorithm 2: Riemannian proximal point method. For (R-PPM), we have U n = 0 and

∥b n ∥ Zn = ∥Γ Z n+1 →Zn (V (Z n+1 )) -V (Z n )∥ Zn ≤ L dist(Z n , Z n+1 ) = γ n L∥V (Z n )∥ Zn ≤ γ n LM = O(γ n )
where we have used the L-Lipschitzness and G-boundedness of V , and the distance-minimizing property of exp within the injectivity radius. Algorithm 3: Riemannian stochastic extra-gradient. For (R-SEG), we have U n = Γ Z n+1/2 →Zn Err(Z n+1/2 ; θ n+1/2 ) so that

E σ 2 n = E ∥U n ∥ 2 Zn = E Err(Z n+1/2 ; θ n+1/2 ) 2 Z n+1/2 ≤ σ 2
by ( 41) and the fact that the parallel transport map is a linear isometry.

To verify [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF], by the definition of (R-SEG), we have

∥b n ∥ = ∥Γ Z n+1/2 →Zn V (Z n+1/2 ) -V (Z n )∥ Zn ≤ L dist Z n+1/2 , Z n = γ n L∥V (Z n ; θ n )∥ Zn → 0
almost surely by Lemma 1.

Algorithm 4: Riemannian optimistic gradient. For (R-OG), we have [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF]. The bias term can then be bounded exactly as in the case of Algorithm 3. ■ 4.2. Operations that preserve Riemannian Robbins-Monro schemes. To increase the efficiency of Riemannian iterative schemes, several important operations have been routinely performed on top of the base algorithms in Section 4.1. In this section, we show that the RRM template incorporates these operations in a unified fashion.

U n = Γ Z n+1/2 →Zn Err(Z n ; θ n+1/2 ) and b n = Γ Z n+1/2 →Zn V (Z n+1/2 ) -V (Z n ), so E[σ 2 n ] ≤ σ 2 again holds by
Retraction-Based Methods. When the exponential map is expensive to compute, a popular alternative in practice is to use the retraction map [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Boumal | An introduction to optimization on smooth manifolds[END_REF], which is defined as a smooth mapping R • (•) : T M → M that mimics the exponential map up to first order: For all

(z, v) ∈ T M, R z (0) = z, d dt R z (tv) t=0 = v. (Rtr)
It turns out that, to replace the exponential map in Algorithms 1-4 with retraction, we only need to pay a small price on the noise assumption:

Proposition 3.
If, in addition to the assumptions in Proposition 2, the SFO has finite fourth moment:

E[∥Err(z; θ)∥ 4 z ] ≤ κ 2 < ∞, then replacing exp Zn (•) with R Zn (•)
in Algorithms 1-4 results in another RRM scheme satisfying (5) and (7) with probability 1. As a result, under Assumptions 3 ′ -4, they generate APTs of (RMD) and converge to its ICT sets. Remark 3. In many practical settings, the map R z (v) satisfies the stronger notion of "secondorder retraction" [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF][START_REF] Boumal | An introduction to optimization on smooth manifolds[END_REF]. In this case, the proof technique below shows that it suffices to have E[∥Err(z; θ)∥ 3 z ] < ∞ in Proposition 3. Proof of Proposition 3. By definition, R z (v) is a smooth map and hence satisfies lim v→0 R z (v) = z. Then Lemma 1 readily implies that Z n+1 lies in the injectivity radius of Z n with probability 1 for n large enough.

We first consider the retraction-based Algorithm 1:

Z n+1 = R Zn (γ n V (Z n ; θ n )). (45)
Let Ṽn ∈ T Zn M be the vector such that exp Zn γ n Ṽn = Z n+1 , i.e.,

γ n Ṽn = exp -1 Zn R Zn (γ n V (Z n ; θ n )) . (46)
Then ( 45) is an RRM scheme with

W n = Ṽn -V (Z n )
where Ṽn is defined in [START_REF] Ljung | Strong convergence of a stochastic approximation algorithm[END_REF]. We will show that, under the assumption E[∥Err(z; θ)∥ 4 z ] < ∞, the following holds with probability 1:

b n = E[W n | F n ] → 0, sup n E ∥W n ∥ 2 Zn < ∞ (47)
which obviously implies ( 5) and [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF]. Consider the curve c(t) := R Zn (tV (Z n ; θ n )). By Lemma 1, for n large enough, c(t) lies in the injectivity radius of Z n almost surely for all t ∈ [0, γ n ]. Let ĉ(t) be the smooth curve of c(t) in the normal coordinate with base Z n and an arbitrary orthonormal frame, and let Ẑn+1 be the normal coordinate of Z n+1 . Also, let Ṽ N n be the (Euclidean) vector of Ṽn expanded in the chosen orthonormal basis, and define V N (Z n ; θ n ) and Err N (Z n ; θ n ) similarly. By definition, Ẑn+1 is nothing but γ n Ṽ N n . Since Z n = c(0) and Z n+1 = c(γ n ), by properties of a retraction map we must have

γ n Ṽ N n = ĉ(γ n ) = ĉ(0) + γ n ċ(0) + O γ 2 n ∥ ċ(0)∥ 2 2 = γ n V N (Z n ; θ n ) + O γ 2 n ∥V (Z n ; θ n )∥ 2 Zn =: γ n V N (Z n ; θ n ) + γ n bn (48)
where bn

= O γ n ∥V (Z n ; θ n )∥ 2 Zn . Therefore, since E[∥Err(z; θ)∥ 4 z ] < ∞ for all z ∈ M, we have E ∥W n ∥ 2 Zn = E Err N (Z n ; θ n ) + bn 2 2 < ∞.
On the other hand,

∥b n ∥ Zn = ∥E[W n | F n ]∥ Zn = E bn F n 2 = O γ n ∥V (Z n ; θ n )∥ 2
Zn . By Chebyshev's inequality and the fact that

E[∥Err(z; θ)∥ 4 z ] ≤ κ 2 < ∞, we have P ∥Err(Z n ; θ n )∥ 2 Zn ≥ n log 1+ ϵ 2 n ≤ κ 2 n log 1+ ϵ
2 n where ε is the same as in our choice of step-size in Proposition 2. Using an calculation identical to Lemma 1, we conclude that

γ n ∥Err(Z n ; θ n )∥ 2 Zn = O   n log 1+ ϵ 2 n n log 1+ε n   = O 1 log ε 4 n
which concludes the proof of [START_REF] Macua | Learning parametric closed-loop policies for markov potential games[END_REF].

For retraction-based Algorithms 2-4, by the above analysis, we may replace

R z (γ n V (•; θ n ))
with exp z γ n V (•; θ n ) + bn where bn → 0 almost surely. The rest is the same as in the proof of Proposition 2. ■ Alternation. Consider a Riemannian (zero-sum) two-player game:

min max X∈M1,Y ∈M2 ℓ(X, Y )
where M 1 and M 2 are smooth manifolds. Instead of performing simultaneous updates as in Algorithms 1-4, a common variant is to alternately update the min/max variables:

X n+1 = exp Xn (γ n [V X (X n , Y n ) + U X,n + b X,n ]) Y n+1 = exp Yn (γ n [V Y (X n+1 , Y n ) + U Y,n + b Y,n ]) (alt-R-RM)
where (V X , V Y ) := (-∇ X ℓ, ∇ Y ℓ), U X,n and b X,n are the noise and bias of the RRM scheme employed by the X variable, and similarly for U Y,n and b Y,n . It is easy to see that alternating Algorithms 1-4 results in another convergent RRM scheme:

Proposition 4. Let Z n := (X n , Y n ) where (X n , Y n
) is generated by using Algorithms 1-4 in (alt-R-RM). Under the same assumptions as in Proposition 2, Z n is itself an RRM scheme satisfying (5) and (7) with probability 1. As a result, under Assumptions 3 ′ -4, Z n generates an APT of (RMD) and converges to its ICT sets.

Remark 4. It is easy to generalize Proposition 4 to any number of players and any type of (not necessarily zero-sum) games.

Remark 5. A simple but useful observation is that compositions of RRM schemes merely "slow down" the algorithm but do not change its asymptotic behavior: For any two RRM schemes RM 1 and RM 2 , the update Z n+1 = RM 2 • RM 1 (Z n ) is equivalent to a new RRM scheme Zn where Z2n = Z n and Z2n+1 = RM 1 ( Z2n ). This allows us to "mix-and-match" Propositions 2-4 to prove, for instance, the convergence of alternating (R-SEG) minimizer vs. retraction-based (R-PPM) maximizer in Riemannian two-player games.

Proof of Proposition 4. Similar to Proposition 2, Lemma 1 guarantees that all geodesics are minimizing and invertible. Hence, by the L-Lipschitzness of V and (alt-R-RM), we have, with probability 1,

∥V Y (X n+1 , Y n ) -V Y (X n , Y n )∥ Yn ≤ L dist X n+1 Y n , X n Y n = Lγ n ∥V X (X n , Y n ) + U X,n + b X,n ∥ Xn → 0 (49)
by Lemma 1 and Proposition 2. Therefore, we may rewrite (alt-R-RM) as Riemannian non-convex optimization. When V (z) = -∇f (z) for some sufficiently smooth function f : M → R, Sard's theorem implies that the ICT sets consist of connected components of critical points [63]. Therefore, all the above-mentioned algorithms converge to critical points of f . This is a vast generalization of the main result of Bonnabel [START_REF] Bonnabel | Stochastic gradient descent on Riemannnian manifolds[END_REF] which, in the context of this paper, can only cover (retraction-based) Algorithm 1 and bounded (SFO).

Z n+1 = exp Zn γ n V X (X n , Y n ) + U X,n + b X,n V Y (X n , Y n ) + U Y,n + b ′ Y,n where b ′ Y,n = b Y,n + V Y (X n+1 , Y n ) -V Y (X n , Y n ).
Riemannian monotone games. Consider the game setting: V = [-∇ zi ℓ i ] where ℓ i : M = i M i → R denotes the loss of the i-th player. If V is monotone [START_REF] Wang | Monotone and accretive vector fields on riemannian manifolds[END_REF] and sufficiently smooth, then combining our results with [START_REF] Kristály | Nash-type equilibria on riemannian manifolds: a variational approach[END_REF] shows that all the above mentioned algorithms converge to the Nash-Stampacchia equilibria, a natural generalization of Nash equilibria to Riemannian manifolds. To the best of our knowledge, most of the algorithms we consider are new in this setting except for the deterministic gradient method and extra-gradient mentioned in the Introduction.

Riemannian non-convex potential games. Solving general non-convex games remains an open problem even for Euclidean spaces. However, a special class of games, known as the Riemannian potential games, has found interesting applications [START_REF] Klavins | Phase regulation of decentralized cyclic robotic systems[END_REF][START_REF] Muhammad | Decentralized coordination with local interactions: Some new directions[END_REF]. For such games, the continuous-time (RMD) has been shown to enjoy many desirable properties such as convergence to critical points or local stability [START_REF] Ratliff | Characterization and computation of local nash equilibria in continuous games[END_REF][START_REF] Ratliff | Genericity and structural stability of non-degenerate differential nash equilibria[END_REF][START_REF] Ratliff | On the characterization of local nash equilibria in continuous games[END_REF], whereas we are not aware of any similar result on stochastic RRM schemes. Our theory bridges this gap by showing that the same guarantees in [START_REF] Ratliff | Characterization and computation of local nash equilibria in continuous games[END_REF][START_REF] Ratliff | Genericity and structural stability of non-degenerate differential nash equilibria[END_REF][START_REF] Ratliff | On the characterization of local nash equilibria in continuous games[END_REF] are in fact achieved by a wide array of RRM schemes.

Markov potential games via Natural Gradient Descent. An interesting class of games that models collaborative behavior in multi-agent reinforcement learning is the Markov potential games (MPGs) [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF][START_REF] Macua | Learning parametric closed-loop policies for markov potential games[END_REF][START_REF] Marden | State based potential games[END_REF]. These games can be solved via a variety of learning schemes, notably the natural policy gradient (NPG) method of Kakade [START_REF] Kakade | A natural policy gradient[END_REF], which has found wide empirical success.

As noticed by Bonnabel [START_REF] Bonnabel | Stochastic gradient descent on Riemannnian manifolds[END_REF], NPG, as a particular case of the Natural Gradient Descent method [START_REF] Amari | Natural gradient works efficiently in learning[END_REF], is an important instance of retraction-based algorithms. In our framework, NPG is simply a retraction-based Algorithm 1. Since MPGs can be reduced to jointly optimizing a global objective as in Riemannian potential games [START_REF] Macua | Learning parametric closed-loop policies for markov potential games[END_REF], under mild assumptions on the noise, our theory together with [START_REF] Ratliff | Characterization and computation of local nash equilibria in continuous games[END_REF][START_REF] Ratliff | Genericity and structural stability of non-degenerate differential nash equilibria[END_REF][START_REF] Ratliff | On the characterization of local nash equilibria in continuous games[END_REF] then not only implies the convergence of plain NPG, but also its optimistic/extra-gradient variants. To our knowledge, these are new schemes that might be of independent interest.

Concluding remarks

Motivated by applications to stochastic approximation and optimization problems on manifolds, our paper proposed a general class of Riemannian iterative algorithms à la Robbins-Monro, and we showed that it captures a wide array of existing and novel algorithms. Our theory provides a unified analysis for the convergence of RRM schemes that might seem vastly different from each other at first sight: By verifying certain simple criteria on the error terms W n in Assumption 2, we can instead study (RMD) to infer the algorithm's long-run behavior.

At the same time, we should stress that our results offer but a glimpse of the flexibility of (RRM). Specifically, the following important questions are left open (and deferred to future work):

(1) We have only analyzed stochastic first-order methods. In practice, however, evaluating Riemannian (stochastic) gradients can be out of reach in online / sequential learning problems. To remedy this issue, various zeroth-order or bandit algorithms, which only require function evaluations, have been proposed. It is, however, unclear whether a Riemannian analogue of the Kiefer-Wolfowitz algorithm can analyzed viea the Robbins-Monro template as in the Euclidean case [START_REF] Hsieh | The limits of min-max optimization algorithms: Convergence to spurious non-critical sets[END_REF]. [START_REF] Amari | Natural gradient works efficiently in learning[END_REF] The assumption on the diminishing step-sizes is indispensable to our analysis, which covers many practically relevant settings. However, another common strategy in practice is the constant step-size rule, which is not covered by our theory. It is thus interesting to see if the techniques for analyzing constant step-size SA schemes can be generalized to Riemannian manifolds [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF][START_REF] Kushner | Asymptotic properties of stochastic approximations with constant coefficients[END_REF][START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF]. (3) Finally, several Riemannian algorithms have been designed to avoid undesirable solutions such as non-minimizing critical points [START_REF] Criscitiello | Efficiently escaping saddle points on manifolds[END_REF][START_REF] Sun | Escaping from saddle points on riemannian manifolds[END_REF]. We conjecture that the avoidance theory of Benaïm and Hirsch [START_REF] Benaïm | Dynamics of Morse-Smale urn processes[END_REF] can be extended to Riemannian manifolds; this, if true, would imply that many iterative Riemannian methods, including the retraction-based ones, converge only to local minimizers ((a.s.)).

Theorem 1 (Jost [29, Theorem 6.6.1]). Let p ∈ M be arbitrary, and suppose that the exponential map exp p is a diffeomorphism on {v ∈ T p M : ∥v∥ ≤ ρ}. Moreover, suppose the sectional curvature is nonpositive and bounded below by -κ 2 on the ball of radius ρ around p. Let r(•) and k(•) be defined as above. Then k(•) is smooth on B(p, ρ) := {z ̸ = p : dist(z, p) ≤ ρ} and ∇ k(z) = -exp -1 p z. In addition, ∥∇ k(z)∥ = r(z), and one has the following control on Hessian of k(•):

Hess k(z)[v, v] ≤ κr(z) ctgh(κr(z))∥v∥ 2 z , for all z ∈ B(p, ρ) and v ∈ T z M.
Since M is simply-connected and complete, we may take ρ = ∞ in the theorem above [29, Corollary 6.9.1].

Recall that the vector field V satisfies the weak-coercivity condition if for some R > 0 and some p ∈ M (and hence all p ∈ M) we have for all z ∈ M outside of B(p, R),

⟨V (z), ∇ k(z)⟩ z ≤ 0.
Our proof relies on constructing a suitable "energy function" that serves as an easy-to-control proxy for the distance of the iterates of (RRM) from the origin. This function will be of the form Φ(z) = f (r(z)) where f is a C ∞ non-negative function with f (x) = 0 for all x ≤ R and satisfies

0 ≤ f ′ (x) ≤ C 1 , f ′′ (x) ≤ C 2 (A.1)
for all x ≥ R. Moreover, we require f (x) = Ω(x) as x → ∞ so that controlling f implies control of x. One example of such functions is given in Lemma A.2. We first show that Φ has a bounded Hessian and is smooth.

Lemma A.1. Let Φ be defined as above. Then, Φ is negatively correlated with V :

⟨∇ Φ(z), V (z)⟩ z ≤ 0, ∀z ∈ M,
and there is some constant C such that Hess Φ

(z)[v, v] ≤ C∥v∥ 2 z . Moreover, Φ is C-smooth, that is, Φ(z 2 ) ≤ Φ(z 1 ) + ⟨∇ Φ(z 1 ), exp -1 z1 z 2 ⟩ + C 2 dist 2 (z 1 , z 2 ).
Proof. We first compute the gradient of Φ:

∇ Φ(z) = 0 r(z) ≤ R, f ′ (r(z)) r(z) ∇ k(z) r(z) > R.
By assumption, f ′ (r(z))/r(z) ≥ 0, and hence, ⟨∇ Φ(z), V (z)⟩ z has the same sign as ⟨∇ k(z), V (z)⟩ z if r(z) > R, and is zero if r(z) ≤ R. Thus, Φ and V are negatively correlated.

To compute the Hessian of Φ, notice that Hess Φ

(z)[v, v] = ⟨∇ v ∇ Φ(z), v⟩ z . Hence, Hess Φ(z)[v, v] = ∇ v f ′ (r(z)) r(z) • ⟨∇ k(z), v⟩ z + f ′ (r(z)) r(z) ⟨∇ v ∇ k(z), v⟩ z = ∇ f ′ (r(z)) r(z) , v • ⟨∇ k(z), v⟩ z a ○ + f ′ (r(z)) r(z) Hess k(z)[v, v] b ○
.

Here we use the same notation for directional derivative of a scalar function and the covariant derivative. We first compute a ○. We have

∇ f ′ (r(z)) r(z) = f ′′ (r(z)) - f ′ (r(z)) r(z) 1 r 2 (z) ∇ k(z), hence, a ○ = f ′′ (r(z)) - f ′ (r(z)) r(z) 1 r 2 (z) ⟨∇ k(z), v⟩ 2 z ≤ C 2 r 2 (z) ∥∇ k(z)∥ 2 z ∥v∥ 2 z = C 2 ∥v∥ 2 z . For b ○, as x ctgh x ≤ 1 + x for x ≥ 0, we obtain b ○ ≤ f ′ (r(z)) r(z) (1 + κr(z))∥v∥ 2 z ≤ C 1 1 R + κ ∥v∥ 2 z .
Summing up everything, we obtain

Hess Φ(z)[v, v] ≤ (C 2 + C 1 /R + C 1 κ)∥v∥ 2 z =: C∥v∥ 2 z
, that is, Φ has bounded Hessian. Moreover, Φ is smooth as a composition of smooth functions. Let z 1 , z 2 ∈ M be arbitrary, and let γ : [0, 1] → M be a geodesic connecting the two. By Taylor's remainder theorem, there exists some t ∈ (0, 1) such that

Φ(z 2 ) = Φ(z 1 ) + ⟨∇ Φ(z 1 ), γ(0)⟩ z1 + 1 2 Hess Φ(γ(t))[ γ, γ].
Thus, by the Hessian upper bound, and noticing that ∥ γ∥ = dist(z 1 , z 2 ) and γ(0) = exp -1 z1 z 2 , we get

Φ(z 2 ) ≤ Φ(z 1 ) + ⟨∇ Φ(z 1 ), exp -1 z1 z 2 ⟩ + C 2 dist 2 (z 1 , z 2 ),
as desired. ■

We now proceed to the main argument, where we show how to use Φ to control the iterates Z n . Letting Φ n = Φ(Z n ) and using Lemma A.1 we get

Φ n+1 = Φ exp Zn (γ n Vn ) ≤ Φ n + γ n ⟨∇ Φ(Z n ), Vn ⟩ Zn + Cγ 2 n 2 ∥ Vn ∥ 2 Zn ≤ Φ n + γ n ⟨∇ Φ(Z n ), U n + b n ⟩ Zn + 3Cγ 2 n 2 ∥V (Z n )∥ 2 Zn + ∥U n ∥ 2 Zn + ∥b n ∥ 2
Zn where the second line follows from negative correlation of Φ and V , the definition (1) of Vn , and the Cauchy-Schwartz inequality. Conditioning on F n and taking expectations, and invoking Cauchy-Schwartz and the fact that ∥∇ Φ

(Z n )∥ ≤ C1 r(Zn) ∥∇ k(Z n )∥ = C 1 , we obtain: E[Φ n+1 | F n ] ≤ Φ n + γ n C 1 B n + 3 2 Cγ 2 n G 2 + B 2 n + σ 2 n , (A.2)
where we have bounded the second moments by their respective upper bounds.

To proceed, let

ε n = γ n C 1 B n + (3/2)Cγ 2 n G 2 + B 2 n + σ 2 n denote the "residual" term in (A.2). Notice that ∞ n=1 ε n ≤ C 1 ∞ n=1 γ n B n + 3C 2 ∞ n=1 γ 2 n (G 2 + B 2 n + σ 2 n ),
and hence, by [START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF] and dominated convergence theorem, we have that E

[ n ε n ] < ∞. Next, consider the auxiliary process Φn = Φ n + E[ ∞ k=n ε k | F n ], adapted to the same filtration. By (A.2), we have E[ Φn+1 | F n ] ≤ Φ n + E[ ∞ k=n ε k | F n ] = Φn-1 , i.e.
, Φn is a supermartingale with respect to F n . This shows that E[ Φn ] ≤ E[ Φ1 ] < ∞, i.e., Φn is uniformly bounded in L 1 . Hence, by Doob's supermartingale convergence theorem [25, Theorem 2.5], it follows that Φn converges with probability 1 to some finite random limit Φ∞ . In turn, since n ε n < ∞ (a.s.), this implies that

Φ n = Φn -E[ ∞ k=n ε n | F n ]
also converges to some (random) finite limit (a.s.). From this and the fact that Φ n = Ω(r(Z n )), we deduce lim sup n r(Z n ) < ∞ as claimed.

Lemma A.2. Let h : R → R be the function

h(x) =      0 x ≤ 0 e -1/x e -1/x +e -1/(1-x) x ∈ (0, 1) 1 x ≥ 1 . It is easy to see that h is C ∞ . For R > 0 define f (x) = x 0 h(s -R) ds.
Then f is C ∞ , satisfies the conditions (A.1) with C 1 = 1 and C 2 = 2. In addition, one has f (x) ≥ x -(R + 1), and hence f (x) = Ω(x).

Proof. As h(x) ∈ [0, 1], we obtain that f ′ (x) ∈ [0, 1]. By straight-forward computation, one observes that h has bounded first derivative 0 ≤ h ′ (x) ≤ 2. Hence,

f ′′ (x) = h ′ (x -R) ≤ 2.
Notice that for x ≥ R + 1, f (x) = We shall call (B.1) the parallel frame system for vectors v and v ′ (the dependence on the initial frame {e k } d k=1 is suppressed). By virtue of parallel transport, in the parallel frame system we have

Γ z→z ′ (v) k = ⟨Γ z→z ′ (v) , e ′ k ⟩ z ′ = ⟨v, Γ z ′ →z (e ′ k )⟩ z = v k . (B.2)
In the same vein, we have Γ z ′ →z (v ′ ) k = v ′k . Thus, since {e k } d k=1 is orthonormal, we may write:

∥v ′ -Γ z→z ′ (v)∥ z = ∥v -Γ z ′ →z (v ′ )∥ z ′ = ∥v -v ′ ∥ 2 .
In other words, in the parallel frame system, the comparison of vectors living on different tangent spaces is reduced to simply comparing the Euclidean norms of their components. For instance, the L-Lipschitzness of V can be rephrased as: ∀z, z ′ ∈ M, ∥V (z ′ ) -V (z)∥ 2 ≤ L dist(z, z ′ ).

(B.3) B.2. The Fermi coordinate system. For any h, let U h ⊂ T Z (t+h) M ≃ R d be a neighborhood of 0 on which the mapping exp Z(t+h) (•) :

U h → M (B.4)
is a diffeomorphism between U h and exp Z(t+h) (U h ). It is well-known that such a neighborhood exists, and that the exponential map exp Z(t+h) along with an arbitrary orthonormal frame at T Z (t+h) M induces a local coordinate system on exp Z(t+h) (U h ), called the normal coordinate with center Z(t + h) [START_REF] Lee | Riemannnian manifolds: an introduction to curvature[END_REF]. Normal coordinates are best suitable for comparing distances of points on manifolds. For instance, if z′ is the normal coordinate of z ′ with center z, then dist(z, z ′ ) = ∥z ′ ∥ 2 . The Fermi coordinate [START_REF] Manasse | Fermi normal coordinates and some basic concepts in differential geometry[END_REF], roughly speaking, is a system of normal coordinates "along a curve". To define it, fix an arbitrary orthonormal frame {e k (0)} d k=1 for T Z(t) M. We can obtain a system of orthonormal frames {e k (h)} . By virtue of the normal coordinate, we know that Φ is a diffeomorphism between R + × V and a neighborhood of R + × {0}. The mapping Φ and its inverse is called the Fermi coordinate system along the curve h → Z(t + h). In the sequel, we will abuse the notation and simply call it the Fermi coordinate along Z(•).

The following property of the Fermi coordinate system plays a key role in our analysis: [START_REF] Fujita | The onsager-machlup function for diffusion processes[END_REF][START_REF] Takahashi | The probability functionals (onsager-machlup functions) of diffusion processes[END_REF]). Let γ be a differentiable curve on M such that γ(h) ∈ U h for all h ∈ R + , and let γ be the curve of γ in the Fermi coordinate system along Z(•) (i.e., (h, γ(h)) = Φ(h, γ(h)). Then Comparing the differential of exp and parallel transport. As will become clear in the proof, the parallel frame system is convenient for comparing vectors at different points, whereas the Fermi coordinate system is best suitable for comparing the distance between curves, both features being essential to our proof. There is, however, a dichotomy: It is known that if the Fermi coordinate in (B.5) is simultaneously a parallel frame system for all points nearby the curve Z(•), then the underlying manifold M must be flat; i.e., the Riemannian curvature tensor vanishes everywhere [START_REF] Iliev | Handbook of normal frames and coordinates[END_REF].

Lemma B.1 ([
Therefore, we would like to work with parallel frame and Fermi coordinate systems separately, and compare the difference between the two whenever needed. To this end, we will need the following technical lemma, whose proof can be found in [START_REF] Lezcano-Casado | Curvature-dependant global convergence rates for optimization on manifolds of bounded geometry[END_REF]Theorem 3.12] For v ∈ T z M, consider the geodesic γ(t) = exp z (tv). If γ is defined and has no interior conjugate point on the interval [0, 1], then

∀v ′ ∈ T z M, ∥T v (v ′ ) -Γ v (v ′ )∥ γ(1) ≤ K • f K low (∥v∥ z ) • ∥v ′ ⊥ ∥ z (B.6)
where v ′ ⊥ := v ′ -⟨v,v ′ ⟩z ⟨v,v⟩z v is the component of v ′ orthogonal to v, T v = d exp z (v) is the differential of the exponential map, and Γ v denotes the parallel transport along γ from γ(0) to γ(1). The function f K low in (B.6) is defined as

f K low (a) =       

Figure 1 :

 1 Figure1: Illustration of an extended Fermi coordinate frame. The moving observer Z(t + h) measures a curve λ(h) using time-indexed geodesics and an "inertial frame system", i.e., frames obtained via parallel transport from Z(t). For each h, λ(h) ∈ R d is the normal coordinate of λ(h) ∈ M. The space-time map Φ : R + × M → R + × R d is locally defined on a neighborhood containing Z(t + h) and λ(h).

d

  k=1 along the geodesic from Z(t + u) to form an orthonormal frame {e ′′ k (u)} d k=1 of T θ(u) M, and define Ṽ k u, θ(u) := ⟨V θ(u) , e ′′ k (u)⟩ θ(u) .

  and frame {e k (u)} d k=1 . More specifically, denote by v the vector d k=1 θk (u)e k (u) ∈ T Z (t+u) M and consider a family of geodesics γ(s, s ′ ) := exp Z(t+u) s(v + s ′ e k (u)) .

  , dist Z(t + u), Z(t + u)

Remark 2 .

 2 The condition E[∥Err(z; θ)∥ 4 z ] < ∞ cannot be relaxed as shown by the Euclidean example of R z (v) = z + (e v -1) and Err(z; θ) law = √ Y -E √ Y , where P(Y > y) = y -3/2 if y ≥ 1, and P(Y > y) = 1 otherwise.

4 . 3 .

 43 [START_REF] Marden | State based potential games[END_REF] and Proposition 2, b ′ Y,n satisfies (5) and[START_REF] Benaïm | Asymptotic pseudotrajectories and chain recurrent flows, with applications[END_REF], concluding the proof.■ Implications. We collect some immediate consequences of the results in Sections 4.1-4.2.

x 0 hB. 1 .

 01 (s -R) ds ≥ x R+1 1 ds = x -(R + 1). ■ Appendix B. Preliminaries for the Proof of Theorem 1In this appendix, we collect the necessary tools for the proof of Theorem 1. The parallel frame system. Fix any two points z, z ′ ∈ M, and consider two arbitrary vectors v ∈ T z M and v ′ ∈ T z ′ M. There is a convenient frame system (i.e., a set of bases for T z M and T z ′ M) for comparing v and v ′ , defined as follows: Pick an arbitrary orthonormal frame {e k } d k=1 for T z M. Since the parallel transport map is an isometry, the vectors {e ′ k } d k=1 := {Γ z→z ′ (e k )} d k=1 form an orthonormal basis for T z ′ M. Consider the component vectors of v, v ′ in these two frames: v ∈ R d , v k := ⟨v, e k ⟩ z ; v ′ ∈ R d , v ′k := ⟨v ′ , e ′ k ⟩ z ′ . (B.1)

d

  k=1 by parallel transporting {e k (0)}d k=1 from T Z (t) M to T Z (t+h) M along the curve h → Z(t + h). Let U h ⊂ T Z (t+h)M be a neighborhood of 0 defined as in (B.4), and set V := h {exp Z(t+h) (U h )} ⊂ M. Consider the mappingΦ : R + × V → R + × R d (B.5) by sending a point (h, z) ∈ R + × V to (h, z) ∈ R + × R d ,where z is the normal coordinate of z with center Z(t + h) and frame {e k (h)} d k=1

  γk (h) = γk (h) -Żk (t + h) + O ∥γ(h)∥ 2 2 . Here, Żk (t + h) := ⟨ Ż(t + h), e k (h)⟩ Z(t+h) is the k-th component of Ż(t + h) in the frame {e k (h)}d k=1 , and γk (h) is the k-th component of γ(h) in the (possibly non-orthonormal) frame induced by the normal coordinate with center Z(t + h) and frame {e k (h)} d k=1 . B.3.

  or [17, Proposition A.1]: Lemma B.2 (Comparing d exp and parallel transport). Let M be a Riemannian manifold whose sectional curvatures are in the interval [K low , K up ], and let K = max (|K low |, |K up |).

r 2 6 ,

 6 K low = 0 r 2 1 -sin(a/r) a/r , K low = 1 r 2 > 0 r 2 sinh(a/r) a/r -1 , K low = -1 r 2 < 0 . (B.7)Moreover, the function f K low is dominated by the caseK low < 0: For all K ≥ |K low | and a ∈ R + , f K low (a) ≤ f -K (a). (B.8)

  In short, we have shown that the exponential mapping is a local diffeomorphism at any point in the set:{γ h θ } h∈(hinit,h fin ) {γ h λ } h∈(hinit,h fin ). The final step in our construction is to consider the union of all such (h init , h fin ) for all ȟ / ∈ C(t, T ) and t + ȟ ̸ = τ n ; we denote the so obtained set by H. We claim that • H is a dense open subset of [0, T ], and • H can be written as a countable union of disjoint open intervals:

d k=1 , there exists an open interval (h init , h fin ) containing ȟ such that, for all h ∈ (h init , h fin ), exp Z(t+h) (•) : T Z (t+h) M → M is a local diffeomorphism at λ(h) and θ(h), where λ(h) and θ(h) are the normal coordinates of λ(h) and θ(h) with center Z(t + h). Let γ h θ be a family of minimizing geodesics such that γ h θ (0) = Z(t + h) and γ h θ (1) = θ(h), and define γ h λ similarly. Since both γ h θ and γ h λ are minimizing geodesics and since θ(h) and λ(h) are not conjugate to Z(t + h), [40, Theorem 10.15] ensures that no point on γ h θ or γ h λ is conjugate to Z(t + h).

In the literal sense; it has nothing to do with differential p-forms.

Acknowledgments

This research was supported by the SNSF grant 407540_167212 through the NRP 75 Big Data program, and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme grant aggreement No 815943. YPH acknowledges funding through an ETH Foundations of Data Science (ETH-FDS) postdoctoral fellowship. Part of this work was done while PM was visiting the Simons Institute for the Theory of Computing. PM is also grateful for financial support by the French National Research Agency (ANR) in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02), the LabEx PERSYVAL (ANR-11-LABX-0025-01), MIAI@Grenoble Alpes (ANR-19-P3IA-0003), and the grant ALIAS (ANR-19-CE48-0018-01).

Appendix A. Stability Analysis

The purpose of this appendix is to prove:

Fix an arbitrary (base) point p ∈ M playing the role of the origin, and define the radial distance r(z) := dist(z, p) and let k(z) := 1 2 r 2 (z). The following theorem makes it clear under which assumptions is k(•) smooth and provides a control on its Hessian.