# Riemannian stochastic approximation algorithms

Mohammad Reza Karimi∗ ∗ Institute for Machine Learning, CAB G 65,
Universitaetstrasse 6, 8092 Zurich, Switzerland.
<EMAIL_ADDRESS>,  Ya-Ping Hsieh∗ ,
Panayotis Mertikopoulos⋄ ⋄ Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP,
LIG, 38000 Grenoble, France<EMAIL_ADDRESS>and  Andreas
Krause∗<EMAIL_ADDRESS>
###### Abstract.

We examine a wide class of stochastic approximation algorithms for solving
(stochastic) nonlinear problems on Riemannian manifolds. Such algorithms arise
naturally in the study of Riemannian optimization, game theory and optimal
transport, but their behavior is much less understood compared to the
Euclidean case because of the lack of a global linear structure on the
manifold. We overcome this difficulty by introducing a suitable Fermi
coordinate frame which allows us to map the asymptotic behavior of the
Riemannian Robbins–Monro (RRM) algorithms under study to that of an associated
deterministic dynamical system. 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, despite the significant complications
that arise due to the curvature and topology of the underlying manifold. We
showcase the flexibility of the proposed framework by applying it to a range
of retraction-based variants of the popular optimistic / extra-gradient
methods for solving minimization problems and games, and we provide a unified
treatment for their convergence.

###### Key words and phrases:

Stochastic approximation; Robbins-Monro algorithms; Riemannian methods;
optimization on manifolds.

###### 2020 Mathematics Subject Classification:

Primary 62L20, 37N40; secondary 90C15, 90C47, 90C48.

A one-page extended abstract of this paper was presented at COLT 2022.

##  1\. Introduction

Consider a nonlinear system of equations of the general form

Find $x^{\ast}\in\mathcal{M}$ such that $v(x^{\ast})=0$ (Root)

where $\mathcal{M}$ is a smooth $d$-dimensonal manifold and $v$ is a vector
field on $\mathcal{M}$. Root-finding problems of this type play a crucial role
in many areas of mathematical programming and data science, from Riemannian
optimization and game theory to reinforcement learning, signal processing, and
information theory. In particular, in addition to standard minimization
problems – that is, when $v=-\nabla f$ for some smooth function $f$ on
$\mathcal{M}$ – the general form of (Root) includes bilevel and saddle-point
problems, dynamic programming, games in normal form, and many other
equilibrium problems that arise in practice. For a range of applications and a
comprehensive introduction to the topic, see Absil et al. [1], Boumal [10],
and references therein.

The vast majority of methods for solving (Root) are iterative in nature, and
they rely on building a successively finer “model” function which is applied
to the last computed approximation of a root in order to get a new
approximation. Usually, this model function is based on the value of $v$ at a
candidate solution;111More sophisticated methods – like Newton’s algorithm –
rely on a first-order Taylor approximation of $v$ around an iterate; Halley’s
method employs a second-order model, and the hierarchy continues with the
general class of Householder methods that use $k$-th order derivatives to
build more precise polynomial models of $v$ around an iterate. However,
because these tensor methods involve the computation of higher-order
derivatives of $v$, they may become highly impractical for moderate-to-high
values of $d$. in many cases however, even this first-order model is too
costly or even impossible to compute – e.g., if
$v(x)=\operatorname{\mathbb{E}}[V(x;\omega)]$ for some random variable
$\omega$ with unknown distribution. In this case, a popular alternative is to
rely on stochastic approximation algorithms that only require oracle access to
a random – and possibly incomplete – approximation of $v$ at the queried
point.

In this general context, when $\mathcal{M}=\mathbb{R}^{d}$, the method of
choice for solving (Root) is the _Robbins–Monro_ (RM) algorithm

$X_{n+1}=X_{n}+\gamma_{n}V_{n}$ (RM)

where $\gamma_{n}>0$ is a variable step-size sequence and $V_{n}$ is a random
estimate of $v(X_{n})$. This method was introduced in the seminal papers of
Robbins and Monro [55] and Kiefer and Wolfowitz [33], and the first general
convergence results were obtained by Ljung [43, 44] for gradient problems.
This has subsequently led to substantial activity on the topic, with major
contributions by Benaïm and Hirsch [6], Benveniste et al. [8], Kushner and
Clark [37], Kushner and Yin [39], and many others. However, the linear
structure of $\mathbb{R}^{d}$ is deeply ingrained in all these works – and the
method itself – preventing its use for solving (Root) in a manifold setting –
e.g., the $d$-dimensional torus for a robotic arm with $d$ joints, Grassman or
Stiefel manifolds for robust principal component analysis, hyperbolic space
for text and graph embeddings, etc. Because of this, the applicability of
Robbins–Monro (RM) methods to general root-finding problems is significantly
more narrow than one might expect, even when $\mathcal{M}$ has a relatively
simple structure (like a matrix group or a Grassmannian).

#### Our contributions in the context of related work.

In view of the above, our main objective is to bridge the gap between
Euclidean and Riemannian stochastic approximation schemes by replacing the
“$+$” operation in (RM) with the Riemannian exponential map
$\exp_{X_{n}}(\cdot)$ on $\mathcal{M}$ – or, more generally (and often more
tractably), a _retraction_ on $\mathcal{M}$ based at $X_{n}$. In Riemannian
optimization, this approach was pioneered by Bonnabel [9], who examined the
case where $v$ is the Riemannian gradient of some objective function $f$.
Subsequent works [65, 61, 11, 15, 41, 63] expanded on the results of Bonnabel
[9] for Riemannian stochastic gradient descent, while similar results were
obtained in [22, 42, 7, 28] 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 (Root). A partial
extension to the non-gradient case was provided by a concurrent line of works
which examined the use of Riemannian extra-gradient methods under the
assumption of (geodesic) monotonicity [23, 60, 50, 21, 32, 14]. 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 Hilbert spaces [3].

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

  1. (1)

We introduce a _generalized_ _Riemannian Robbins–Monro_ template which
includes as special cases all methods mentioned above (Riemannian stochastic
gradient descent, extra-gradient, proximal point methods, etc.), as well as a
number of new stochastic approximation (SA) schemes for (Root).

  2. (2)

Under some mild technical conditions on $\mathcal{M}$, we show that the
sequence of generated points forms an “approximate solution” – an asymptotic
pseudotrajectory (APT) to be exact – of an associated deterministic dynamical
system (Theorem 2).

This stochastic approximation result extends the seminal theory of Benaïm and
Hirsch [6] for Euclidean Robbins–Monro schemes to a Riemannian setting, and
allows us to infer almost sure convergence of RRM schemes to the internally
chain-transitive (ICT) sets of the underlying Riemannian dynamics (cf.
Corollaries 1 and 2). In gradient and strictly monotone problems, these ICT
sets are precisely the roots of $v$, so we readily recover many of the
asymptotic convergence results mentioned above (often under much weaker
assumptions). In addition, as we show in Section 4, our framework applies to
several 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, techniques, and related approaches.

In the absence of a linear structure on $\mathcal{M}$, the major challenge we
have to overcome is the lack of a suitable coordinate frame within which to
analyze the trajectories of Riemannian SA algorithms. This reflects the
dichotomy that, unlike in the case of $\mathbb{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 [45], 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 [56], 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
$\mathcal{M}$ implies that the manifold is globally _flat_ , i.e., essentially
Euclidean [29]; this assumption is too restrictive for bona fide Riemannian
applications, so the analysis of [56] is not relevant for our purposes. An
additional issue is that the error bounds employed by Shah [56, p. 1131] rule
out vector fields with a rotational component – such as $v(x,y)=(-y,x)$ on
$\mathbb{R}^{2}$ – further limiting the applicability of their techniques to
the setting under consideration.

Finally, Durmus et al. [18, 19] also recently 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 [19] are necessarily ergodic in nature and
hence beyond our paper’s scope. The setting of [18] 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
function_ – the so-called “gradient-like” case [4] – so there is no overlap
with our analysis.

#### Basic notions and notation.

Throughout our paper, $\mathcal{M}$ will denote a $d$-dimensional,
geodesically complete, Riemannian manifold. We will write
$\langle\cdot,\cdot\rangle_{x}$, $x\in\mathcal{M}$, for its underlying metric,
$\lVert\cdot\rVert_{x}$ for the induced norm, and we will assume that the
sectional curvatures of $\mathcal{M}$ are bounded from above and below by
$K_{\textup{up}}$ and $K_{\textup{low}}$ respectively.

For any curve $\gamma$ on $\mathcal{M}$, the notation $\dot{\gamma}(t)$ will
always denote its velocity at time $t$. Given any pair of points
$x,x^{\prime}\in\mathcal{M}$ and a vector $v$ in the tangent space at $x$,
$\mathcal{T}_{x}\mathcal{M}$, we denote by $\Gamma_{x\to
x^{\prime}}\left(v\right)\in\mathcal{T}_{x^{\prime}}\mathcal{M}$ the vector
obtained by parallel transporting $v$ along the minimizing geodesic connecting
$x$ and $x^{\prime}$; if minimizing geodesics are not unique, $\Gamma_{x\to
x^{\prime}}\left(v\right)$ will be understood as parallel transport along any
of them. We also say that a vector field $v$ on $\mathcal{M}$ is
(geodesically) $L$-Lipschitz if, for all $x,x^{\prime}\in\mathcal{M}$, we have

$\lVert\Gamma_{x\to
x^{\prime}}\left(v(x)\right)-v(x^{\prime})\rVert_{x^{\prime}}=\lVert
v(x)-\Gamma_{x^{\prime}\to x}\left(v(x^{\prime})\right)\rVert_{x}\leq
L\operatorname{dist}(x,x^{\prime}),$ (1)

where $\operatorname{dist}(\cdot,\cdot)$ denotes the distance function on
$\mathcal{M}$ induced by $\langle\cdot,\cdot\rangle$. All vector fields in
this paper are assumed to be $L$-Lipschitz, complete and bounded, i.e.,
$G\coloneqq\sup_{x\in\mathcal{M}}\lVert v(x)\rVert_{x}<\infty$. Finally, the
Riemannian gradient of a smoooth function $f\colon\mathcal{M}\to\mathbb{R}$
will be denoted by $\nabla f$.

For a detailed account, we refer the reader to the masterful treatment of Lee
[40].

##  2\. Riemannian Robbins–Monro algorithms: Definitions and assumptions

###  2.1. The Riemannian Robbins–Monro template

We begin by discussing the basic template of _Riemannian Robbins–Monro_ (RRM)
algorithms. As we stated before, the main difference with their Euclidean
counterpart is that addition along “straight lines” is replaced with the
Riemannian exponential mapping. This leads to the abstract update rule

$X_{n+1}=\exp_{#1}(\gamma_{n}\left(v(X_{n})+W_{n}\right))$ (RRM)

where

  1. (1)

$X_{n}\in\mathcal{M}$ denotes the state of the algorithm at each iteration
counter $n=1,2,\dotsc$

  2. (2)

$W_{n}\in\mathcal{T}_{\scriptscriptstyle X_{n}}\mathcal{M}$ is an abstract
error term (described in detail below).

  3. (3)

$\gamma_{n}>0$ is the method’s step-size (also discussed below).

In the above, we tacitly assume that the error term $W_{n}$ is generated
_after_ $X_{n}$, so it is not adapted to the history
$\mathcal{F}_{n}\coloneqq\sigma(X_{1},\dotsc,X_{n})$ of $X_{n}$ – that is,
$W_{n}$ is random relative to $\mathcal{F}_{n}$. In addition, to differentiate
between “random” (zero-mean) and “systematic” (non-zero-mean) errors, it will
be convenient to further decompose $W_{n}$ as

$W_{n}=U_{n}+b_{n}$ (2)

where
$U_{n}=W_{n}-\operatorname{\mathbb{E}}[W_{n}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]$
captures the random, zero-mean part of $W_{n}$, while
$b_{n}=\operatorname{\mathbb{E}}[W_{n}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]$
represents the systematic component thereof.

To quantify all this, we will assume in the sequel that $U_{n}$ and $b_{n}$
are bounded as

$\operatorname{\mathbb{E}}[\lVert U_{n}\rVert_{\scriptscriptstyle
X_{n}}^{2}|\mathcal{F}_{n}]\leq\sigma_{n}^{2}\quad\textrm{and}\quad\operatorname{\mathbb{E}}[\lVert
b_{n}\rVert_{\scriptscriptstyle
X_{n}}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]\leq
B_{n}\quad\text{{(}a.s.{)}}$ (3)

where $\sigma_{n}$ and $B_{n}$ are to be construed as upper bounds on the
noise and bias of the error terms entering (RRM). Finally, for concision, we
will also write

$V_{n}=v(X_{n})+W_{n}$ (4)

so $V_{n}$ can be seen as a noisy – and potentially biased – estimator of
$v(X_{n})$ in (RRM). Obviously, $V_{n}=v(X_{n})$ whenever
$B_{n}=\sigma_{n}=0$; we will refer to this case as “deterministic”.

###  2.2. Stochastic approximation

Our main goal in the sequel will be to connect the asymptotic behavior of the
trajectories of (RRM) to a continuous-time dynamical system on $\mathcal{M}$.
To do so, in analogy with the ODE method for standard, Euclidean Robbins–Monro
schemes [39, 4], we will view (RRM) as an inexact (forward) Euler
discretization of the _Riemannian mean dynamics_

$\dot{\theta}(t)=v(\theta(t))$ (RMD)

and we will try to establish a measure of “closeness” between the discrete-
time sequence $X_{n}$ generated by (RRM) and the continuous-time orbits
$\theta(t)$ of (RMD).

In the Euclidean case, this is provided by taking an affine interpolation
$X(t)$ of $X_{n}$ that agrees with $X_{n}$ at all instances of the “effective
time” variable $\tau_{n}=\sum\nolimits_{k=1}^{n-1}\gamma_{k}$, viz.

$X(t)=X_{n}+\frac{t-\tau_{n}}{\tau_{n+1}-\tau_{n}}(X_{n+1}-X_{n})\quad\text{for
all $t\in[\tau_{n},\tau_{n+1}]$, $n=1,2,\dotsc$}$ (5)

Of course, this definition is not meaningful in a Riemannian setting because
of the lack of an affine structure on $\mathcal{M}$. Instead, noting that the
increments in (5) can be rewritten as
$(X_{n+1}-X_{n})/(\tau_{n+1}-\tau_{n})=V_{n}$, a natural Riemannian analogue
would be to follow the geodesic emanating from $X_{n}$ along $V_{n}$ until
reaching $X_{n+1}$. With this in mind, we define the _geodesic interpolation_
$X(t)$ of $X_{n}$ as

$X(t)=\exp_{#1}\left((t-\tau_{n})V_{n}\right)\quad\text{for all
$t\in[\tau_{n},\tau_{n+1}]$, $n\geq 1$},$ (GI)

so, by construction,  (itn) $X(\tau_{n})=X_{n}$for all $n$; and  (itn) each
segment of $X(t)$ is a geodesic.

Now, to compare $X(t)$ to the solution orbits of (RMD), let
$\Theta\colon\mathbb{R}_{+}\times\mathcal{M}\to\mathcal{M}$ denote the _flow_
of (RMD), i.e., $\Theta_{h}(x)$ is simply the position at time $h\geq 0$ of
the solution orbit of (RMD) that starts at $x\in\mathcal{M}$ (recall here that
$v$ is complete). We then have the following notion of “asymptotic closeness”:

######  Definition 1 (Benaïm and Hirsch, 1996).

We say that $X(t)$ is an _asymptotic pseudotrajectory_ (APT) of the mean
dynamics (RMD) if, for all $T>0$, we have

$\lim_{t\to\infty}\sup_{0\leq h\leq
T}\operatorname{dist}(X(t+h),\Theta_{h}(X(t)))=0.$ (6)

Intuitively, Definition 1 states that one cannot distinguish between the
geodesically interpolated process $X(t+h)$ and the orbit of (RMD) that starts
at $X(t)$ as $t\to\infty$. This is a highly non-trivial requirement, so much
of the analysis to follow hinges on establishing exactly this property; we
carry all this out in Section 3, where we also describe the precise connection
between the asymptotic behavior of (RRM) and that of (RMD).

###  2.3. Technical assumptions

We conclude this section with the basic assumptions that underlie the rest of
our paper. These are as follows:

######  Assumption 1 (Robbins–Monro step-sizes).

The step-size sequence $\gamma_{n}$, $n=1,2,\dotsc$, of (RRM) satisfies the
Robbins–Monro summability conditions $\sum_{n}\gamma_{n}=\infty$ and
$\sum_{n}\gamma_{n}^{2}<\infty$.

######  Assumption 2 (Error bounds).

The bounds on the noise and bias in (3) satisfy

$\sum_{n=1}^{\infty}\gamma_{n}^{2}\operatorname{\mathbb{E}}\left[\sigma_{n}^{2}\right]<\infty,\quad\sum_{n=1}^{\infty}\gamma_{n}\operatorname{\mathbb{E}}\left[B_{n}\right]<\infty,\quad\text{and}\quad
B_{n}\to 0\;\;\text{ {(}a.s.{)}.}$ (7)

Variants of the above assumptions are standard in the context of (Euclidean)
Robbins–Monro methods, cf. Kushner and Clark [37], Benaïm [4], and references
therein. Nonetheless, it is worth noting that Assumption 1 lies under the
explicit control of the algorithm designer, while Assumption 2 is _implicit_
and depends on the specific problem at hand – the mechanism providing access
to $v$, the specific form of (RRM), etc. In this regard, Assumption 2 is more
delicate than Assumption 1; we discuss this issue in detail in Section 4,
where we show that Assumption 2 is indeed satisfied for a wide range of
practical algorithms that adhere to the general template (RRM).

Moving forward, to exclude cases where (RRM) becomes unstable over time, a
standard practice in the literature is to assume that the sequence $X_{n}$ is
contained in a compact subset of $\mathcal{M}$, a property known as
_precompactness_. Formally, we have:

######  Assumption 3 (Precompactness).

The set of iterates $\\{X_{n}:n=1,2,\dotsc\\}$ has compact closure in
$\mathcal{M}$.

Albeit standard in the literature [37, 8, 4], Assumption 3 may be difficult to
verify if $\mathcal{M}$ is not itself compact (e.g., if
$\mathcal{M}=\mathbb{R}^{d}$). To account for this, we introduce in Section 3
a set of structural hypotheses on $\mathcal{M}$ and $v$ which guarantee that
Assumption 3 holds with probability $1$.

With all this in hand, our last blanket assumption is a technical requirement
that interfaces between the geometry of $\mathcal{M}$ and the dynamics of
(RMD). To state it, recall first that $x\in\mathcal{M}$ is said to be
_conjugate_ to $x^{\prime}=\exp_{#1}(v)$ along the geodesic $\exp_{#1}(tv)$ if
the exponential map $\exp_{x}(\cdot)$ is not a diffeomorphism in a
neighborhood of $v$ [40]. In addition, define the _Picard flow_
$\lambda\colon\mathbb{R}_{+}\to\mathcal{M}$ associated with $X(t)$ to be the
dynamical system

$\dot{\lambda}(h)=\Gamma_{X(t+h)\to\lambda(h)}\left(v(X(t+h))\right)$ (PFlow)

with initial condition ${\lambda}(0)=X(t)$.222The term “Picard flow” stems
from the fact that the integral $\int_{0}^{h}v(X(t+s))\>ds$ is the basic
iteration in Picard’s method of successive approximations for solving an ODE
in Euclidean spaces. In the case of (PFlow), the parallel transport is the
extra ingredient required to express the idea of “integrating $v$ along
$X(t)$”, so (PFlow) can be seen as a bona fide generalization of the Picard
iteration to Riemannian manifolds. Finally, given some $T,t>0$, consider the
sets

$\displaystyle\mathcal{C}_{\lambda}(t,T)$
$\displaystyle\coloneqq\left\\{h\in[0,T]:X(t+h)\textup{ is conjugate to
}\lambda(h)\right\\},$ (8a) $\displaystyle\mathcal{C}_{\theta}(t,T)$
$\displaystyle\coloneqq\left\\{h\in[0,T]:X(t+h)\textup{ is conjugate to
}\Theta_{h}(X(t))\right\\},$ (8b)

and let
$\mathcal{C}(t,T)=\mathcal{C}_{\lambda}(t,T)\cup\mathcal{C}_{\theta}(t,T)$. We
then make the following assumption:

######  Assumption 4 (Nowhere dense conjugates).

The set $\mathcal{C}(t,T)$ is nowhere dense in $[0,T]$ for sufficiently large
$t>0$ and for all $T>0$.

At first sight, Assumption 4 may appear quite opaque but it is otherwise quite
mild. Indeed, since the set of points conjugate to $X(t+h)$ is at most one-
dimensional [64], the only way that Assumption 4 can fail is if the curves
$X(t+h)$ and $\theta(h)$ or $\lambda(h)$ simultaneously traverse the same one-
dimensional submanifold of $\mathcal{M}$ – a fact which occurs with
probability $0$ if the distribution of $W_{n}$ is non-singular. It is also
worth noting that Assumption 4 is automatically satisfied on negatively curved
spaces by the Cartan–Hadamard theorem (cf. Proposition 1 in the next section)
and, finally, it is straightforward to verify Assumption 4 manually on many of
the manifolds that arise in practical applications – such as spheres,
Grassmannians, Stiefel manifolds and fixed-rank spectrahedra, cf. [40, 1, 46]
and references therein. We discuss all this in more detail in the next
section.

##  3\. Analysis and main results

###  3.1. Statement and discussion of main results

The connecting tissue between asymptotic pseudotrajectories and the mean
dynamics (RMD) is provided by the _internally chain-transitive_ (ICT) sets of
(RMD). These are defined as follows:

######  Definition 2 (Benaïm, 1999).

Let $\mathcal{S}$ be a nonempty compact subset of $\mathcal{M}$. Then:

  1. (1)

$\mathcal{S}$ is _invariant_ under (RMD) if
$\Theta_{t}(\mathcal{S})=\mathcal{S}$ for all $t\in\mathbb{R}$.

  2. (2)

$\mathcal{S}$ is an _attractor_ of (RMD) if it is invariant under (RMD) and
there exists a compact neighborhood $\mathcal{K}$ of $\mathcal{S}$ such that
$\lim_{t\to\infty}\operatorname{dist}(\Theta_{t}(x),\mathcal{S})=0$ uniformly
in $x\in\mathcal{K}$.

  3. (3)

$\mathcal{S}$ is _internally chain-transitive_ (ICT) if it is invariant and
$\Theta|_{\mathcal{S}}$ admits no attractors other than $\mathcal{S}$.

In words, ICT sets are the “terminal objects” of the dynamics (RMD): the
orbits that converge to an ICT set cannot be contained in a smaller subset
thereof, so ICT sets may be viewed as minimal connected periodic orbits up to
arbitrarily small numerical errors [4, Prop. 5.3]. The importance of this
property for the theory of stochastic approximation is owed to the following
theorem:

######  Theorem 1 (Benaïm and Hirsch, 1996).

Let $X(t)$, $t\geq 0$, be a precompact asymptotic pseudotrajectory of (RMD).
Then the limit set $\mathcal{L}(X)\coloneqq\operatorname*{\bigcap}_{t\geq
0}\operatorname{cl}\\{X(s):s\geq t\\}$ of $X$ is an ICT set of (RMD).

######  Remark.

Theorem 1 above was originally stated in the context of abstract metric
spaces; we have adapted it here to the Riemannian setting of (RMD) for
concision and concreteness. $\lozenge$

The limit set theorem above provides a fundamental link between APTs and the
long-run behavior of (RMD) as captured by its ICT sets. That being said, the
APT property itself may be difficult to verify from first principles, so the
application of Theorem 1 to Riemannian Robbins–Monro algorithms can be just as
difficult. In the Euclidean case ($\mathcal{M}=\mathbb{R}^{d}$), Benaïm and
Hirsch [6] address this issue via a series of criteria under which standard
(Euclidean) Robbins–Monro methods give rise to an APT of the associated mean
dynamics [4, Props. 4.1 and 4.2]. Unfortunately however, in a bona fide
Riemannian setting, these criteria cannot be used because they are
inextricably tied to the affine structure of $\mathbb{R}^{d}$; as a result, it
is not clear how to leverage Theorem 1 to obtain a theory of stochastic
approximation for Robbins–Monro methods on Riemannian manifolds. We tackle
this question below:

######  Theorem 2.

Suppose that Assumptions 1–4 hold. Then, with probability $1$, the geodesic
interpolation $X(t)$ of the sequence of iterates $X_{n}$ of (RRM) is an APT of
(RMD).

Theorem 2 plays a pivotal role in our paper and is the cornerstone for the
analysis to follow. In particular, by invoking Theorems 1 and 2 in tandem, we
readily obtain the following characterization of the limiting behavior of
(RRM):

######  Corollary 1.

Suppose that Assumptions 1–4 hold. Then $X_{n}$ converges to an ICT set of
(RMD) w.p.$1$.

An important consequence of Corollary 1 is that the analysis of the
stochastic, discrete-time system (RRM) boils down to that of the
_deterministic_ , _continuous-time_ system (RMD). In this way, Theorem 2
allows us to employ the same high-level strategies as the classic literature
on stochastic approximation: For example, if $v$ admits a potential or is
geodesically (strictly) monotone, it is easy to verify that the only ICT sets
of (RMD) are the roots of $v$ [57], so we readily recover the series of
asymptotic convergence results mentioned in the introduction. In Section 4, we
present a wide range of problems whose ICT sets coincide with the solutions of
(Root), and we further illustrate how Theorem 2 captures a series of
Riemannian stochastic approximation algorithms – old and new – in a unified
fashion.

The proof of Theorem 2 is fairly arduous, so we defer it to the end of this
section; instead, we proceed to discuss here in more detail the theorem’s
precompactness and conjugacy requirements (Assumptions 3 and 4 respectively).
First, as mentioned above, Assumption 3 is common in the stochastic
approximation literature but, in general, it is impractical to verify directly
from the primitives of the problem at hand – that is, the ambient manifold
$\mathcal{M}$ and the defining vector field $v$. Likewise, if the cut loci of
different points on $\mathcal{M}$ happen to have a complicated topological
structure, verifying Assumption 4 could also be impractical. To account for
all this, we provide below a set of structural conditions on $\mathcal{M}$ and
$v$ which guarantee that Assumptions 3 and 4 both hold with probability $1$:

######  Proposition 1.

Suppose that the following hypotheses hold:

  1. n(H1)

$\mathcal{M}$ is a Hadamard manifold with the Heine-Borel property.333A
Hadamard manifold is a simply connected Riemannian manifold with non-positive
sectional curvatures, and the Heine-Borel property simply posits that every
closed bounded subset thereof is compact.

  2. n(H2)

$v$ satisfies the weak coercivity condition

$\left\langle v(x),\nabla\operatorname{dist}^{2}(p,x)\right\rangle_{x}\leq 0$
(WC)

for some base point $p\in\mathcal{M}$ and for all
$x\in\mathcal{M}\setminus\mathbb{B}_{R}(p)$ outside a closed geodesic ball
$\mathbb{B}_{R}(p)\coloneqq\\{x\in\mathcal{M}:\operatorname{dist}(p,x)\leq
R\\}$ of radius $R>0$ and centered at $p$.

Then, with probability $1$, Assumptions 3 and 4 hold as stated – with
Assumption 4 only requiring 1.

######  Corollary 2.

Suppose that Assumptions 1 and 2 and Items 1 and 2 hold. Then $X_{n}$
converges to an ICT set of (RMD) w.p.$1$.

To streamline our discussion, we postpone the proof of Proposition 1 to
Appendix A. For now, we simply note that the Hadamard requirement 1 is fairly
common in applications to Riemannian optimization and includes Stiefel and
Grassmann manifolds (as homogeneous spaces), standard matrix manifolds,
hyperbolic spaces, etc. [13, 26].

In this regard, the most delicate requirement in Proposition 1 is Item 2. This
hypothesis may be viewed as a Riemannian relaxation of the coercivity
condition $\lim_{x\to\infty}\langle v(x),x\rangle/\lVert x\rVert_{2}=-\infty$,
which posits that the “inward-pointing” component of $v$ grows unbounded at
infinity, a property which is frequently used to ensure the stability of
Euclidean iterative algorithms [52, 20]. In our Riemannian setting, the role
of the radial field is played by the gradient of the squared distance function
$\nabla\operatorname{dist}^{2}(p,x)$, which is itself equal to $-2\log_{x}(p)$
if $p$ lies in the injectivity radius of $x$ (and hence for all
$x\in\mathcal{M}$ if $\mathcal{M}$ is Hadamard). In addition, it is important
to bear in mind that (WC) _does not_ impose any growth requirements on the
radial component of $v$; it only requires that $v$ does not have a consistent
outward-pointing component that could lead the process to diverge, so it is
significantly lighter in that respect (hence the adjective “weak”).

###  3.2. Proof of Theorem 2

Because the proof of Theorem 2 and the geometric scaffolding required are
quite delicate, we begin with a high-level outline outlining the main
difficulties and technical challenges involved. In brief, the basic obstacle
that we have to overcome is as follows:

  * •

On the one hand, we need a coordinate system to compare and compute distances
between different _points_ on $\mathcal{M}$. This can be done efficiently in
normal coordinates [40].

  * •

On the other hand, we also need to compare _vectors_ living on different
tangent spaces. In general, this comparison is very difficult to carry out in
a normal coordinate frame, but it is much easier in the _parallel frame
system_ (PFS) that we describe in detail in Section B.1.

Intuitively, the normal coordinate system is where _distances_ behave as if
the ambient space were Euclidean, and the parallel frame system is where
_vectors_ behave as in the Euclidean setting. Unfortunately however, the only
regime where these two systems can coexist is when $\mathcal{M}$ is _flat_ ,
i.e., the problem is “essentially Euclidean” to begin with. To circumvent this
obstacle, we take the following technical approach:

  1. (1)

Based on the notion of _Fermi coordinates_ [45] – which can be intuitively
understood as “normal coordinates along a curve”, cf. Fig. 1 – we begin by
constructing an _extended Fermi coordinate frame_ that allows us to focus on a
neighborhood of $X(t)$ containing all the information we need. [This is a
challenging, but otherwise mostly technical, step that does not affect the big
picture.]

  2. (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 incurs
an error term that is not present in the analysis of _Euclidean_ stochastic
approximation schemes, and which is difficult to control in regions of high
sectional curvature, cf. Eqs. 25 and 27.

  3. (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 invariably
introduces an additional error relative to the Euclidean analysis, cf. Eqs. 37
and 38.

  4. (4)

Serendipitously, these additional error terms can be managed without any
further assumptions, and a series of arguments in the spirit of Benaïm and
Hirsch [6] concludes our proof.

Figure 1. Illustration of an extended Fermi coordinate frame. The moving
observer $X(t+h)$ measures a curve $\lambda(h)$ using time-indexed geodesics
and an “inertial frame system”, i.e., frames obtained via parallel transport
from $X(t)$. For each $h$, $\tilde{\lambda}(h)\in\mathbb{R}^{d}$ is the normal
coordinate of $\lambda(h)\in\mathcal{M}$. The space-time map
$\tilde{\Theta}\colon\mathbb{R}_{+}\times\mathcal{M}\to\mathbb{R}_{+}\times\mathbb{R}^{d}$
is locally defined on a neighborhood containing $X(t+h)$ and $\lambda(h)$.

To formalize all this, we will need the following notions and results from
Riemannian geometry:

  * •

The concept of a _parallel frame system_.

  * •

A technical lemma by Fujita and Kotani [24] and Takahashi and Watanabe [59]
which bounds the distortion of velocities measured by a moving observer on a
manifold relative to flat spaces.

  * •

A comparison lemma to estimate the difference between parallel transport and
the pushforward of the Riemannian exponential map.

These elements are quite technical and involved, so we proceed directly to the
proof of Theorem 2 and we relegate all precise statements regarding the above
to Appendix B.

#### Step 1: Discrete-to-continuous time comparisons.

Following Benaïm [4], consider the “continuous-to-discrete” counter

$m(t)=\sup\\{n\geq 1:t\geq\tau_{n}\\}$ (9)

which measures the number of iterations required for the effective time
$\tau_{n}=\sum_{k=1}^{n-1}\gamma_{k}$ to reach a given timestamp $t\geq 0$.
Moroever, given an arbitrary sequence $A_{n}$, we will denote its piecewise-
constant interpolation as

$\bar{A}(t)=A_{n}\quad\text{for all $t\in[\tau_{n},\tau_{n+1})$, $n\geq 1$.}$
(10)

Using this notation, we may express the geodesic interpolation $X(t)$ of
$X_{n}$ in differential form as

$\displaystyle\dot{X}(t+h)$ $\displaystyle=\Gamma_{\bar{X}(t+h)\to
X(t+h)}\left(v\left(\bar{X}(t+h)\right)+\bar{W}(t+h)\right)$ (11)

and we further let

$\bar{\gamma}^{\ast}(t)\coloneqq\sup_{t\leq h\leq
t+T}\bar{\gamma}(h)\quad\text{and}\quad\bar{B}^{\ast}(t)\coloneqq\sup_{t\leq
h\leq t+T}\bar{B}(h)\quad\text{for all $t,h\geq 0$}.$ (12)

By the Stolz-Cesàro theorem, Assumptions 1 and 2 readily imply that
$\lim_{n\to\infty}B_{n}=0$, so we also have
$\lim_{t\to\infty}\bar{B}^{\ast}(t)=0$; we will use this fact freely in the
sequel.

As a last comparison step, we will also need a noise stability criterion. To
that end, let $\\{e_{k}(n)\\}_{k=1}^{d}$ be an arbitrary sequence of
orthonormal bases for $\mathcal{T}_{X_{n}}\mathcal{M}$, and let
$U_{n}^{\shortparallel}$ be the (Euclidean) noise vector composed of
components of the noise $U_{n}$ in the basis $\\{e_{k}(n)\\}_{k=1\dotsc d}$,
viz.

$U_{k,n}^{\shortparallel}\coloneqq\langle U_{n},e_{k}(n)\rangle_{X_{n}}.$ (13)

It is then easy to see that
$\operatorname{\mathbb{E}}[U_{n}^{\shortparallel}|\mathcal{F}_{n}]=0$, and,
moreover

$\operatorname{\mathbb{E}}\left[\left\lVert
U_{n}^{\shortparallel}\right\rVert_{2}^{2}\nonscript\,\middle|\nonscript\,\mathopen{}\mathcal{F}_{n}\right]=\operatorname{\mathbb{E}}\left[\lVert
U_{n}\rVert^{2}_{X_{n}}\nonscript\,\middle|\nonscript\,\mathopen{}\mathcal{F}_{n}\right]\leq\sigma^{2}_{n}$
(14)

by Assumption 2. Then, letting

$\Delta(t;T)\coloneqq\sup_{0\leq h\leq
T}\left\lVert\int_{t}^{t+h}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2},$
(15)

a classical argument by Benaïm [4, cf. Eq. (13) and onwards] readily gives

$\lim_{t\to\infty}\Delta(t;T)=0\quad\text{for all $T>0$ {(}a.s.{)}}.$ (16)

#### Step 2: Preliminary error bounds.

We first note that, by Assumption 3, $\sup_{t}r\left(X(t)\right)<\infty$ where
$r(\cdot)\coloneqq\operatorname{dist}(\cdot,p)$ is the radial function defined
in (WC). We claim that Assumption 3 also implies the boundedness of the Picard
flow. To see this, recall that the parallel transport is an isometry, so

$\displaystyle\left\lVert\dot{\lambda}(h)\right\rVert_{\lambda(h)}=\left\lVert\Gamma_{X(t+h)\to\lambda(h)}\left(v(X(t+h))\right)\right\rVert_{\lambda(h)}=\left\lVert
v(X(t+h))\right\rVert_{X(t+h)},$ (17)

and hence

$\sup_{0\leq h\leq T}\operatorname{dist}\left(\lambda(0),\lambda(h)\right)\leq
T\cdot\sup_{0\leq h\leq
T}\left\lVert\dot{\lambda}(h)\right\rVert_{\lambda(h)}<\infty,$ (18)

which implies $\sup_{0\leq h\leq T}r(\lambda(h))<\infty$. On the other hand,
the boundedness for the flow follows readily from the fact that $v$ and
$\operatorname{dist}$ are both $L$-Lipschitz,444Due to the smoothness of the
flow, the function $r\left(\theta(h)\right)$ is always differentiable in $h$
in the metric space sense [2], even though $\theta(h)$ might reach the cut
locus of $p$. so

$\displaystyle{\frac{\>d}{\>dh}r\left(\theta(h)\right)}$
$\displaystyle\leq\left\lVert\dot{\theta}(h)\right\rVert_{\theta(h)}=\left\lVert
v(\theta(h))\right\rVert_{\theta(h)}$ $\displaystyle\leq\left\lVert
v(\theta(h))-\Gamma_{p\to\theta(h)}\left(v(p)\right)\right\rVert_{\theta(h)}+\left\lVert\Gamma_{p\to\theta(h)}\left(v(p)\right)\right\rVert_{\theta(h)}$
$\displaystyle\leq
Lr\left(\theta(h)\right)+\left\lVert{v(p)}\right\rVert_{p}.$ (19)

An application of Grönwall’s inequality then gives $\sup_{0\leq h\leq
T}r\left(\theta(h)\right)<\infty$. Therefore, all computations in the sequel
can be restricted to a compact set, and we may assume without loss of
generality that $\bar{\gamma}^{\ast}(t)<1$ for all $t\geq 0$ and, in addition,
there exists some $R\equiv R({T,L,G})$ such that

$\displaystyle\max\left\\{\operatorname{dist}\left(X(t+h),\theta(h)\right),\operatorname{dist}\left(X(t+h),\lambda(h)\right)\right\\}\leq
R\quad\text{for all $h\in[0,T]$}.$ (20)

#### Step 3: Constructing the extended Fermi coordinates.

Let $\tilde{\mathcal{U}}$ be the neighborhood defined in Section B.2
restricted to $[t,t+T]$, i.e.,
$\tilde{\mathcal{U}}=\operatorname*{\bigcup}_{h=0}^{T}\mathcal{U}_{h}$ with
$\mathcal{U}_{h}$ given by (B.5). Clearly, $\tilde{\mathcal{U}}$ contains
$\\{X(t+h):h\in[0,T]\\}$; moreover, by construction, $\tilde{\Phi}$ carries a
system of orthonormal frames $\left\\{e_{k}(h)\right\\}_{k=1}^{d}$, one for
each $X(t+h)$. In what follows, all quantities will be expressed in these
frames unless explicitly mentioned otherwise.

Now, fix some $h_{0}\in[0,T]$ and let $\gamma_{\theta}$ and $\gamma_{\lambda}$
be two minimizing geodesics such that
$\gamma_{\theta}(0)=\gamma_{\lambda}(0)=X(t+h_{0})$,
$\gamma_{\theta}(1)=\theta(h_{0})$ and $\gamma_{\lambda}(1)=\lambda(h_{0})$.
Our first goal will be to extend $\tilde{\mathcal{U}}$ to an open set of
$\mathcal{M}$ that contains the geodesics $\gamma_{\theta}$ and
$\gamma_{\lambda}$, while retaining the exponential mapping as a local
diffeomorphism. This will serve a dual purpose:

  1. (1)

It enables us, for a fixed $h_{0}$, to consider the parallel frame systems at
$X(t+h_{0})$ and $\theta(h_{0})$ – or $X(t+h_{0})$ and $\lambda(h_{0})$,
depending on the context – so that we can easily compare the vector fields at
these points; see Section B.1.

  2. (2)

We want to apply Lemma B.1 to the curves $\theta(\cdot)$ and $\lambda$;
however, for $\dot{\tilde{\theta}}(h_{0})$ and $\dot{\tilde{\lambda}}(h_{0})$
to make sense, $\tilde{\mathcal{U}}$ must contain both curves for at least an
open time interval that includes $h_{0}$.

This is where Assumption 4 comes into play: Since the conjugate points arise
by definition when the exponential map ceases to be local diffeomorphisms
[40], it is reasonable to expect that, away from the time points where
$X(t+h)$ is conjugate to $\theta({h})$ or $\lambda({h})$, it is always
possible to extend $\tilde{\mathcal{U}}$ to include the geodesics connecting
$X(t+h)$ to $\theta({h})$ and $\lambda({h})$. In this regard, Assumption 4
simply posits that there cannot be “too many” such conjugate points.

To formalize this, fix some $h_{0}\notin\mathcal{C}(t,T)$ where
$\mathcal{C}(t,T)$ is defined as in (4), and assume also that
$t+h_{0}\neq\tau_{n}$ for all $n$ – i.e., $t+h_{0}$ is not a “kink point” of
(GI). Since the exponential mapping is a local diffeomorphism away from
conjugate points [40], it follows that
$\exp_{X(t+h_{0})}\colon\mathcal{T}_{X(t+h_{0})}\mathcal{M}\to\mathcal{M}$ is
a local diffeomorphism at $\tilde{\lambda}(h_{0})$ and
$\tilde{\theta}(h_{0})$, where $\tilde{\lambda}(h_{0})$ and
$\tilde{\theta}(h_{0})$ are the normal coordinates of $\lambda(h_{0})$ and
$\theta(h_{0})$ with center $X(t+h_{0})$ respectively. By the continuity of
the Picard flow and the frame system $\\{e_{k}(\cdot)\\}_{k=1}^{d}$, there
exists an open interval $(h_{\text{init}},h_{\text{fin}})$ containing $h_{0}$
such that, for all $h\in(h_{\text{init}},h_{\text{fin}})$,
$\exp_{X(t+h)}\colon\mathcal{T}_{X(t+h)}\mathcal{M}\to\mathcal{M}$ is a local
diffeomorphism at $\tilde{\lambda}(h)$ and $\tilde{\theta}(h)$, where
$\tilde{\lambda}(h)$ and $\tilde{\theta}(h)$ are the normal coordinates of
$\lambda({h})$ and $\theta({h})$ with center $X(t+h)$ respectively.

On that account, let $\gamma^{h}_{\theta}$ be a family of minimizing geodesics
such that $\gamma^{h}_{\theta}(0)=X(t+h)$ and
$\gamma^{h}_{\theta}(1)=\theta(h)$, and define $\gamma^{h}_{\lambda}$
similarly. Since both $\gamma^{h}_{\theta}$ and $\gamma^{h}_{\lambda}$ are
minimizing geodesics and $\theta(h)$ and $\lambda(h)$ are not conjugate to
$X(t+h)$, [40, Theorem 10.15] ensures that no point on $\gamma^{h}_{\theta}$
or $\gamma^{h}_{\lambda}$ is conjugate to $X(t+{h})$. Summing up, we have
shown that the exponential mapping is a local diffeomorphism at any point of
the set

$\\{\gamma^{h}_{\theta}\\}_{h\in(h_{\text{init}},h_{\text{fin}})}\cup\\{\gamma^{h}_{\lambda}\\}_{h\in(h_{\text{init}},h_{\text{fin}})}.$
(21)

The final step in our construction is to consider the union of all such
$(h_{\text{init}},h_{\text{fin}})$ for all $h_{0}\notin\mathcal{C}(t,T)$ and
$t+h_{0}\neq\tau_{n}$; we denote the set obtained in this way by
$\mathcal{H}$. More precisely, we claim that  itn) $\mathcal{H}$is a dense
open subset of $[0,T]$; and  itn) $\mathcal{H}$can be written as a countable
union of disjoint open intervals, i.e.,
$\mathcal{H}=\operatorname*{\bigcup}_{k}(h_{k},h_{k+1})$.  Indeed, the first
claim follows readily from Assumption 4 and the fact that the set
$\\{h_{0}:t+h_{0}=\tau_{n}\text{ for some }n\\}$ is countable. The second
claim is due to the compactness of $[0,T]$ and the fact that every nonempty
open interval in $\mathbb{R}$ contains a rational number.

In view of the above, it is possible to extend $\tilde{\mathcal{U}}$ and
$\tilde{\Phi}$ to an open set containing
$\\{\gamma^{h}_{\theta}\\}_{h\in(h_{\text{k}},h_{\text{k+1}})}\cup\\{\gamma^{h}_{\lambda}\\}_{h\in(h_{\text{k}},h_{\text{k+1}})},$
which, in turn, obviously contains
$\operatorname*{\bigcup}_{h\in(h_{\text{k}},h_{\text{k+1}})}\\{\theta(h),\lambda(h)\\}$.
We call such a pair $(\tilde{\Phi},\tilde{\mathcal{U}})$ the _extended Fermi
coordinate frame_ because it not only contains the central curve $h\to X(t+h)$
as in the classical case, but also $\theta(h)$ and $\lambda(h)$ for almost
every $h\in[0,T]$.

#### Step 4: Controlling the distance by decomposition.

From this point forward, we will assume that all computations take place in
the extended Fermi coordinate system $(\tilde{\Phi},\tilde{\mathcal{U}})$. By
the definition of $\theta(t)$, and given that $\tilde{\theta}(h)$ is the
normal coordinate of $\theta({h})$ with center $X(t+{h})$, we have, for all
$h\in\mathcal{H}$,

$\displaystyle\operatorname{dist}(X(t+h),\Theta_{h}(X(t)))$
$\displaystyle=\operatorname{dist}(X(t+h),\theta(h))=\lVert\tilde{\theta}(h)\rVert_{2}\leq\lVert\tilde{\theta}(h)-\tilde{\lambda}(h)\rVert_{2}+\lVert\tilde{\lambda}(h)\rVert_{2}.$
(22)

Since $\mathcal{H}$ is a dense open subset of $[0,T]$ and since it is a
countable union of open intervals, it follows that $\tilde{\theta}(h)$ and
$\tilde{\lambda}(h)$ are differentiable except on a set of measure zero. We
may thus write

$\displaystyle\lVert\tilde{\theta}(h)-\tilde{\lambda}(h)\rVert_{2}=\left\lVert\int_{0}^{h}\big{[}\dot{\tilde{\theta}}(u)-\dot{\tilde{\lambda}}(u)\big{]}\>du\right\rVert_{2}.$
(23)

Moreover, by Lemma B.1 and the definition of (RMD) and (PFlow), we have

$\displaystyle\dot{\tilde{\theta}}^{k}(u)$
$\displaystyle=\dot{\theta}_{k}(u)-\dot{X}_{k}(t+u)+\operatorname{\mathcal{O}}\big{(}\lVert\tilde{\theta}(u)\rVert_{2}^{2}\big{)}=\tilde{v}_{k}(u,\tilde{\theta}(u))-\dot{X}_{k}(t+u)+\operatorname{\mathcal{O}}\big{(}\lVert\tilde{\theta}(u)\rVert_{2}^{2}\big{)},$
(24a) $\displaystyle\dot{\tilde{\lambda}}^{k}(u)$
$\displaystyle=\dot{\lambda}_{k}(u)-\dot{X}_{k}(t+u)+\operatorname{\mathcal{O}}\big{(}\lVert\tilde{\lambda}(u)\rVert_{2}^{2}\big{)}=\tilde{\Lambda}_{k}\left(u,\tilde{\lambda}(u)\right)-\dot{X}_{k}(t+u)+\operatorname{\mathcal{O}}\big{(}\lVert\tilde{\lambda}(u)\rVert_{2}^{2}\big{)},$
(24b)

where $\dot{X}_{k}$ is defined in Lemma B.1 and
$\tilde{v}_{k}(u,\tilde{\theta}(u))$ and
$\tilde{\Lambda}_{k}\left(u,\tilde{\lambda}(u)\right)$ are, respectively, the
$k$-th components of the vectors $v\big{(}\theta(u)\big{)}$ and
$\Gamma_{X(t+u)\to\lambda(u)}\left(v(X(t+u))\right)$ in the frame induced by
the normal coordinate with center $X(t+u)$ and frame
$\left\\{e_{k}(u)\right\\}_{k=1}^{d}$. Now, denoting by
$\tilde{v}(u,\tilde{\theta}(u))$ the (Euclidean) vector with components
$\tilde{v}_{k}(u,\tilde{\theta}(u))$ and, by
$\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)$ the vector with components
$\tilde{\Lambda}_{k}\left(u,\tilde{\lambda}(u)\right)$, we may write

$\displaystyle\lVert\tilde{\theta}(h)-\tilde{\lambda}(h)\rVert_{2}\leq\Bigg{\|}\int_{0}^{h}\left(\tilde{v}(u,\tilde{\theta}(u))-\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)\right)\>du\Bigg{\|}_{2}+\int_{0}^{h}R_{1}(u)\>du$
(25)

where the remainder term $R_{1}(u)$ is of order
$\operatorname{\mathcal{O}}\left(\lVert\tilde{\theta}(u)\rVert_{2}^{2}+\lVert\tilde{\lambda}(u)\rVert_{2}^{2}\right)$.
Noting that, by (20),

$\displaystyle\lVert\tilde{\theta}(u)\rVert_{2}$
$\displaystyle=\operatorname{dist}\left(X(t+u),\Theta_{u}(X(t))\right)\leq R,$
(26a) $\displaystyle\lVert\tilde{\lambda}(u)\rVert_{2}$
$\displaystyle=\operatorname{dist}\left(X(t+u),\lambda(t+u)\right)\leq R,$
(26b)

we have
$\lVert\tilde{\theta}(u)\rVert_{2}^{2}+\lVert\tilde{\lambda}(u)\rVert_{2}^{2}\leq
R\left(\lVert\tilde{\theta}(u)\rVert_{2}+\lVert\tilde{\lambda}(u)\rVert_{2}\right)$,
we get
$R_{1}(u)=\operatorname{\mathcal{O}}_{R}\left(\lVert\tilde{\theta}(u)\rVert_{2}+\lVert\tilde{\lambda}(u)\rVert_{2}\right)$
where $\operatorname{\mathcal{O}}_{R}(\cdot)$ includes any constants depending
on $R$.

In the same vein, denoting by $\dot{\tilde{X}}(t+u)$ the Euclidean vector
whose $k$-th component is $\dot{X}_{k}(t+u)$, we have

$\displaystyle\lVert\tilde{\lambda}(h)\rVert_{2}\leq\Bigg{\|}\int_{0}^{h}\left(\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)-\dot{\tilde{X}}(t+u)\right)\>du\Bigg{\|}_{2}+\int_{0}^{h}R_{2}(u)\>du$
(27)

with
$R_{2}(u)=\operatorname{\mathcal{O}}\left(\lVert\tilde{\lambda}(u)\rVert_{2}^{2}\right)=\operatorname{\mathcal{O}}_{R}\left(\lVert\tilde{\lambda}(u)\rVert_{2}\right)$.

#### Step 5: From Fermi to parallel coordinates.

So far, we have reduced the proof to comparing the vectors in Eqs. 25 and 27.
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 Section
B.1) with a common base point $X(t+u)$, the vectors $v\big{(}\theta(u)\big{)}$
and $\Gamma_{X(t+u)\to\lambda(u)}\left(v(X(t+u))\right)$ possess some
favorable properties. To see this, parallel transport the frame
$\left\\{e_{k}(u)\right\\}_{k=1}^{d}$ along the geodesic from $X(t+u)$ to form
an orthonormal frame $\left\\{e^{\prime}_{k}(u)\right\\}_{k=1}^{d}$ of
$\mathcal{T}_{\lambda(u)}\mathcal{M}$, and consider the Euclidean vector
$\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)$ whose
components are defined as

$\displaystyle\tilde{\Lambda}^{\shortparallel}_{k}\left(u,\tilde{\lambda}(u)\right)$
$\displaystyle\coloneqq\langle\Gamma_{X(t+u)\to\lambda(u)}\left(v(X(t+u))\right),e^{\prime}_{k}(u)\rangle_{\lambda(u)}=\langle
v(X(t+u)),e_{k}(u)\rangle_{X(t+u)}.$ (28)

Similarly, parallel transport the frame $\left\\{e_{k}(u)\right\\}_{k=1}^{d}$
along the geodesic from $X(t+u)$ to form an orthonormal frame
$\left\\{e^{\prime\prime}_{k}(u)\right\\}_{k=1}^{d}$ of
$\mathcal{T}_{\theta(u)}\mathcal{M}$, and let
$\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\theta}(u)\right)\coloneqq\langle
v\big{(}\theta(u)\big{)},e^{\prime\prime}_{k}(u)\rangle_{\theta(u)}.$ Since
the parallel transport is an isometry, we get

$\displaystyle\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}^{\shortparallel}_{k}\left(u,\tilde{\lambda}(u)\right)$
$\displaystyle=\left\langle
v\big{(}\theta(u)\big{)},e^{\prime\prime}_{k}(u)\right\rangle_{\theta(u)}-\left\langle\Gamma_{X(t+u)\to\lambda(u)}\left(v(X(t+u))\right),e^{\prime}_{k}(u)\right\rangle_{\lambda(u)}$
$\displaystyle=\left\langle\Gamma_{\theta(u)\to
X(t+u)}\left(v\big{(}\theta(u)\big{)}\right),{e_{k}(u)}\right\rangle_{X(t+u)}-\left\langle{v(X(t+u))},e_{k}(u)\right\rangle_{X(t+u)}$
$\displaystyle=\left\langle\Gamma_{\theta(u)\to
X(t+u)}\left(v\big{(}\theta(u)\big{)}\right)-v(X(t+u)),{e_{k}(u)}\right\rangle_{X(t+u)},$
(29)

and hence

$\displaystyle\left\lVert\tilde{v}^{\shortparallel}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)\right\rVert_{2}$
$\displaystyle=\lVert\Gamma_{\theta(u)\to
X(t+u)}\left(v\big{(}\theta(u)\big{)}\right)-v(X(t+u))\rVert_{X(t+u)}$
$\displaystyle\leq
L\operatorname{dist}(\theta(u),X(t+u))=L\lVert\tilde{\theta}(u)\rVert_{2}$
(30)

where we have used the fact that $v$ is $L$-Lipschitz and that
$\tilde{\theta}(u)$ are normal coordinates with center $X(t+u)$.

We would therefore like to replace
$\tilde{v}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)$
with
$\tilde{v}^{\shortparallel}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)$
in (25). To this end, consider the difference in the $k$-th component of
$\tilde{v}\left(u,\tilde{\theta}(u)\right)$ and
$\tilde{v}^{\shortparallel}\left(u,\tilde{\theta}(u)\right)$

$\displaystyle{\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\theta}(u)\right)-\tilde{v}_{k}\left(u,\tilde{\theta}(u)\right)}=\left\langle
v\big{(}\theta(u)\big{)},e^{\prime\prime}_{k}(u)\right\rangle_{\theta(u)}-\left\langle
v\big{(}\theta(u)\big{)},\frac{\partial}{\partial
x_{k}}\Big{|}_{\theta(u)}\right\rangle_{\theta(u)}$ (31)

where $\frac{\partial}{\partial x_{i}}\Big{|}_{\theta(u)}$ is the $k$-th basis
in the frame induced by the normal coordinate with center $X(t+u)$ and frame
$\\{e_{k}(u)\\}_{k=1}^{d}$. More specifically, denote by $v$ the vector
$\sum_{k=1}^{d}\tilde{\theta}^{k}(u)e_{k}(u)\in\mathcal{T}_{X(t+u)}\mathcal{M}$
and consider the family of geodesics

$\gamma(s,s^{\prime})\coloneqq\exp_{#1}(s(v+s^{\prime}e_{k}(u)))$ (32)

so

$\frac{\partial}{\partial x_{i}}\Big{|}_{\theta(u)}=\frac{\partial}{\partial
s^{\prime}}\gamma(1,0)=\>d\exp_{X(t+u)}(v)(e_{k}(u)).$ (33)

Invoking Cauchy-Schwartz, Lemma B.2, and the fact that
$\\{e_{k}(u)\\}_{k=1}^{d}$ is orthonormal, we get

$\displaystyle\left\lvert\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\theta}(u)\right)-\tilde{v}_{k}\left(u,\tilde{\theta}(u)\right)\right\rvert$
$\displaystyle=\left\lvert\left\langle
v\big{(}\theta(u)\big{)},e^{\prime\prime}_{k}(u)-\frac{\partial}{\partial
x_{k}}\Big{|}_{\theta(u)}\right\rangle_{\theta(u)}\right\rvert$
$\displaystyle\leq\lVert
v\big{(}\theta(u)\big{)}\rVert_{\theta(u)}\cdot\left\lVert
e^{\prime\prime}_{k}(u)-\frac{\partial}{\partial
x_{k}}\Big{|}_{\theta(u)}\right\rVert_{\theta(u)}$ $\displaystyle\leq G\cdot
K_{\max}\cdot f_{-K_{\max}}(v)\cdot\lVert e_{k}(u)_{\perp}\rVert_{X(t+u)}$
$\displaystyle\leq G\cdot K_{\max}\cdot f_{-K_{\max}}(v)=G\cdot K_{\max}\cdot
f_{-K_{\max}}(\lVert\tilde{\theta}(u)\rVert_{2})$ (34)

where $K_{\max}=\max\\{\lvert K_{\textup{up}}\rvert,\lvert
K_{\textup{low}}\rvert\\}$. Since $(\sinh x)/x-1\leq\cosh x\leq e^{x}$ for all
$x\geq 0$, (B.9) and (B.10) yield

$f_{-K_{\max}}(\lVert\tilde{\theta}(u)\rVert_{2})\leq\frac{1}{K_{\max}}\exp\left(\sqrt{K_{\max}}\lVert\tilde{\theta}(u)\rVert_{2}\right)\leq\frac{e^{R\sqrt{K_{\max}}}}{R\cdot
K_{\max}}\cdot\lVert\tilde{\theta}(u)\rVert_{2}$ (35)

where the last inequality follows from
$\lVert\tilde{\theta}(u)\rVert_{2}=\operatorname{dist}(X(t+u),\Theta_{u}(X(t)))\leq
R$. Combining (34) and (35), we thus get

$\displaystyle\left\lvert\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\theta}(u)\right)-\tilde{v}_{k}\left(u,\tilde{\theta}(u)\right)\right\rvert$
$\displaystyle\leq\frac{Ge^{R\sqrt{K_{\max}}}}{R}\cdot\lVert\tilde{\theta}(u)\rVert_{2}.$
(36)

In short, we have shown that

$\tilde{v}\left(u,\tilde{\theta}(u)\right)=\tilde{v}^{\shortparallel}\left(u,\tilde{\theta}(u)\right)+R_{3}(u)$
(37)

where
$R_{3}(u)=\operatorname{\mathcal{O}}_{G,K_{\max},R}\left(\lVert\tilde{\theta}(u)\rVert_{2}\right)$
collects constants that depend on $G$, $K_{\max}$ and $R$. Exactly the same
computation shows that, for some
$R_{4}(u)=\operatorname{\mathcal{O}}_{G,K_{\max},R}\left(\lVert\tilde{\lambda}(u)\rVert_{2}\right)$,

$\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)=\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)+R_{4}(u).$
(38)

#### Step 6: Putting everything together.

With all these preliminaries in hand, we are finally in a position to complete
our proof.

###### Proof of Theorem 2.

We will proceed by bounding (25) and (27). Using (37), (38), and (30) in (25),
we obtain:

$\displaystyle\lVert\tilde{\theta}(h)-\tilde{\lambda}(h)\rVert_{2}$
$\displaystyle\leq\Bigg{\|}\int_{0}^{h}\left(\tilde{v}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)\right)\>du\Bigg{\|}_{2}+\int_{0}^{h}R_{1}(u)\>du$
$\displaystyle\leq\int_{0}^{h}\left\lVert\tilde{v}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)\right\rVert_{2}\>du+\int_{0}^{h}R_{1}(u)\>du$
$\displaystyle\leq\int_{0}^{h}\left\lVert\tilde{v}^{\shortparallel}\left(u,\tilde{\theta}(u)\right)-\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)\right\rVert_{2}\>du+\int_{0}^{h}(R_{1}+R_{3}+R_{4})(u)\>du$
$\displaystyle\leq
L\int_{0}^{h}\left\lVert\tilde{\theta}(u)\right\rVert_{2}\>du+\int_{0}^{h}(R_{1}+R_{3}+R_{4})(u)\>du.$
(39)

We next turn to (27). Our first task is to obtain an expression for
$\dot{\tilde{X}}(t+u)$, i.e., the Euclidean vector whose $k$-th component is
$\dot{{X}}_{k}(t+u)$. To this end, fix an iteration count $n$, and consider
all $u$ such that $t+u\in[\tau_{n},\tau_{n+1})$ (that is, consider only the
interpolation between $X_{n}$ and $X_{n+1}$). We claim that
$\dot{{X}}_{k}(t+u)$ is constant throughout all such $u$, and, in particular,
$\dot{{X}}_{k}(t+u)=\dot{{X}}_{k}(\tau_{n})$. This readily follows by noticing
that

  1. (1)

The curve $X(t)$ is a geodesic segment when restricted to
$[\tau_{n},\tau_{n+1})$; see (GI).

  2. (2)

The Fermi coordinates along $X(\cdot)$, when restricted to
$\\{X(s):s\in[\tau_{n},\tau_{n+1})\\}$, is simply a parallel frame system;
this follows from the fact that the frame $\\{e_{k}(u)\\}_{k=1}^{d}$ is
obtained from parallel transporting $\\{e_{k}(\tau_{n})\\}_{k=1}^{d}$ along
$X(\cdot)$ for all $u$ such that $t+u\in[\tau_{n},\tau_{n+1})$.

In this way, a simple calculation akin to (B.2) yields, for all such $u$,

$\displaystyle\dot{{X}}_{k}(t+u)$
$\displaystyle=\langle\dot{X}(t+u),e_{i}(u)\rangle_{X(t+u)}$
$\displaystyle=\left\langle\Gamma_{X(t+u)\to
X(\tau_{n})}\left(\left(\dot{X}(t+u)\right)\right),e_{i}(\tau_{n})\right\rangle_{X(\tau_{n})}$
$\displaystyle=\left\langle
v(X(\tau_{n}))+W_{n},e_{i}(\tau_{n})\right\rangle_{X(\tau_{n})}=\left\langle
v(X(\tau_{n}))+U_{n}+b_{n},e_{i}(\tau_{n})\right\rangle_{X(\tau_{n})}$ (40)

where we have used (11) and the definition of $W_{n}$ in the last equality.

Armed with the above, we can obtain a succinct expression for
$\dot{\tilde{X}}(t+u)$ as follows. First, let $\tilde{\bar{x}}(u)$ be the
normal coordinate of $\bar{X}(t+u)$ with center $X(t+u)$ (i.e.,
$(u,\tilde{\bar{x}}(u))=\tilde{\Phi}(u,\bar{X}(t+u))$), and define an
Euclidean vector $\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))$ by setting
its $k$-th component to

$\tilde{v}^{\shortparallel}_{k}\left(u,\tilde{\bar{x}}(u)\right)\coloneqq\left\langle
v\left(\bar{X}(t+u)\right),e_{i}(m(t+u))\right\rangle_{\bar{X}(t+u)}$ (41)

where the mapping $m(\cdot)$ is defined in (9). Define also the Euclidean
noise and bias vectors $U_{n}^{\shortparallel}$ and $b_{n}^{\shortparallel}$
by setting their components to

$\displaystyle U_{k,n}^{\shortparallel}$ $\displaystyle\coloneqq\left\langle
U_{n},e_{i}(m(t+u))\right\rangle_{\bar{X}(t+u)},$ (42a) $\displaystyle
b_{k,n}^{\shortparallel}$ $\displaystyle\coloneqq\left\langle
b_{n},e_{i}(m(t+u))\right\rangle_{\bar{X}(t+u)}.$ (42b)

Then (40) states precisely that

$\displaystyle\dot{\tilde{X}}(t+u)=\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))+\bar{U}^{\shortparallel}(t+u)+\bar{b}^{\shortparallel}(t+u).$
(43)

Substituting (43) into (27) and invoking (38), (12), and (15), we obtain

$\displaystyle\lVert\tilde{\lambda}(h)\rVert_{2}$
$\displaystyle\leq\left\lVert\int_{0}^{h}\left(\tilde{\Lambda}\left(u,\tilde{\lambda}(u)\right)-\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))-\bar{U}^{\shortparallel}(t+u)-\bar{b}^{\shortparallel}(t+u)\right)\>du\right\rVert_{2}+\int_{0}^{h}R_{2}(u)\>du$
$\displaystyle\leq\int_{0}^{h}\left\lVert\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)-\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))\right\rVert_{2}\>du$
$\displaystyle\qquad+\left\lVert\int_{0}^{h}\bar{U}^{\shortparallel}(t+u)\>du\right\rVert_{2}+\left\lVert\int_{0}^{h}\bar{b}^{\shortparallel}(t+u)\>du\right\rVert_{2}+\int_{0}^{h}[R_{2}(u)+R_{4}(u)]\>du$
$\displaystyle\leq\int_{0}^{h}\left\lVert\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)-\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))\right\rVert_{2}\>du$
$\displaystyle\qquad+\Delta(t,T)+\bar{B}^{\ast}(t)\cdot
h+\int_{0}^{h}[R_{2}(u)+R_{4}(u)]\>du.$ (44)

To bound the first term in (44), recall (28) and (41). An identical argument
leading to (B.4) shows that

$\displaystyle\left\lVert\tilde{\Lambda}^{\shortparallel}\left(u,\tilde{\lambda}(u)\right)-\tilde{v}^{\shortparallel}(u,\tilde{\bar{x}}(u))\right\rVert_{2}$
$\displaystyle=\left\lVert v(X(t+u))-\Gamma_{\bar{X}(t+u)\to
X(t+u)}\left(v(\bar{X}(t+u))\right)\right\rVert_{X(t+u)}$ $\displaystyle\leq
L\cdot\operatorname{dist}\left(\bar{X}(t+u),X(t+u)\right).$ (45)

Since $X(\cdot)$ is a (not necessarily minimizing) geodesic on $[m(t+u),t+u]$,
Eqs. 43 and 12 yield

$\displaystyle\operatorname{dist}(\bar{X}(t+u),X(t+u))$
$\displaystyle\leq\left\lVert\int_{m(t+u)}^{t+u}\dot{\tilde{X}}(s)\>ds\right\rVert_{2}$
$\displaystyle=\left\lVert\int_{m(t+u)}^{t+u}\tilde{v}^{\shortparallel}(s,\tilde{\bar{x}}(s))+\bar{U}^{\shortparallel}(s)+\bar{b}^{\shortparallel}(s)\>ds\right\rVert_{2}$
$\displaystyle\leq\left\lVert\int_{m(t+u)}^{t+u}\tilde{v}^{\shortparallel}(s,\tilde{\bar{x}}(s))\>ds\right\rVert_{2}$
$\displaystyle\qquad+\left\lVert\int_{m(t+u)}^{t+u}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}+\left\lVert\int_{m(t+u)}^{t+u}\bar{b}^{\shortparallel}(s)\>ds\right\rVert_{2}$
$\displaystyle\leq\left(G+\bar{B}^{\ast}(t)\right)\cdot\left(\int_{m(t+u)}^{t+u}\>ds\right)+\left\lVert\int_{m(t+u)}^{t+u}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}$
$\displaystyle\leq\left(G+\bar{B}^{\ast}(t)\right)\bar{\gamma}^{\ast}(t)+\left\lVert\int_{m(t+u)}^{t+u}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}.$
(46)

For $t$ large enough, we have $\bar{\gamma}^{\ast}(t)<1$, and hence

$\displaystyle\left\lVert\int_{m(t+u)}^{t+u}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}$
$\displaystyle\leq\left\lVert\int_{t-1}^{m(t+u)}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}+\left\lVert\int_{t-1}^{t+u}\bar{U}^{\shortparallel}(s)\>ds\right\rVert_{2}\leq
2\Delta(t-1,T+1).$ (47)

Combining Eqs. 44–47 then gives

$\displaystyle\lVert\tilde{\lambda}(h)\rVert_{2}$ $\displaystyle\leq
Lh[(G+\bar{B}^{\ast}(t))\bar{\gamma}^{\ast}(t)$
$\displaystyle\qquad+2\Delta(t-1,T+1)]+\Delta(t,T)+h\bar{B}^{\ast}(t)+\int_{0}^{h}[R_{2}(u)+R_{4}(u)]\>du$
$\displaystyle\leq
hC_{L,G}\left(\bar{B}^{\ast}(t)+\bar{\gamma}^{\ast}(t)+\Delta(t-1,T+1)\right)+\int_{0}^{h}[R_{2}(u)+R_{4}(u)]\>du$
(48)

for some constant $C_{L,G}$ that depends only on $L$ and $G$. Using (39) and
(48), we can then bound (22) as

$\displaystyle\operatorname{dist}(X(t+h),\Theta_{h}(X(t)))$
$\displaystyle=\lVert\tilde{\theta}(h)\rVert_{2}\leq\lVert\tilde{\theta}(h)\rVert_{2}+\lVert\tilde{\lambda}(h)\rVert_{2}$
(49)
$\displaystyle\leq\lVert\tilde{\theta}(h)-\tilde{\lambda}(h)\rVert_{2}+2\lVert\tilde{\lambda}(h)\rVert_{2}$
$\displaystyle\leq
L\int_{0}^{h}\left\lVert\tilde{\theta}(u)\right\rVert_{2}\>du+2hC_{L,G}\left(\bar{B}^{\ast}(t)+\bar{\gamma}^{\ast}(t)+\Delta(t-1,T+1)\right)$
$\displaystyle+3\int_{0}^{h}[R_{1}(u)+R_{2}(u)+R_{3}(u)+R_{4}(u)]\>du$ (50)

where
$[R_{1}(u)+R_{2}(u)+R_{3}(u)+R_{4}(u)]=\operatorname{\mathcal{O}}_{G,K_{\max},R}\left(\left\lVert\tilde{\theta}(h)\right\rVert_{2}+\left\lVert\tilde{\lambda}(h)\right\rVert_{2}\right)$.
Therefore, there exists a constant $C_{L,G,K_{\max},R}$ depending only on
$L,G,K_{\max}$, and $R$ such that we may bound (49) as

$\displaystyle\lVert\tilde{\theta}(h)\rVert_{2}+\lVert\tilde{\lambda}(h)\rVert_{2}$
$\displaystyle\leq
C_{L,G,K_{\max},R}\int_{0}^{h}\left(\left\lVert\tilde{\theta}(h)\right\rVert_{2}+\left\lVert\tilde{\lambda}(h)\right\rVert_{2}\right)\>du$
$\displaystyle+2hC_{L,G}\left(\bar{B}^{\ast}(t)+\bar{\gamma}^{\ast}(t)+\Delta(t-1,T+1)\right).$
(51)

Grönwall’s inequality then implies

$\displaystyle\lVert\tilde{\theta}(h)\rVert_{2}+\lVert\tilde{\lambda}(h)\rVert_{2}$
$\displaystyle\leq
2hC_{L,G}\left(\bar{B}^{\ast}(t)+\bar{\gamma}^{\ast}(t)+\Delta(t-1,T+1)\right)e^{h\cdot
C_{L,G,K_{\max},R}}.$ (52)

From (52), we may conclude that, with probability $1$, we have

$\displaystyle\lim_{t\to\infty}\sup_{t\leq h\leq T}$
$\displaystyle\operatorname{dist}(X(t+h),\Theta_{h}(X(t)))$
$\displaystyle\leq\lim_{t\to\infty}\sup_{t\leq h\leq
T}\left(\lVert\tilde{\theta}(h)\rVert_{2}+\lVert\tilde{\lambda}(h)\rVert_{2}\right)$
$\displaystyle\leq\lim_{t\to\infty}T\cdot
2C_{L,G}\left(\bar{B}^{\ast}(t)+\bar{\gamma}^{\ast}(t)+\Delta(t-1,T+1)\right)\cdot
e^{T\cdot C_{L,G,K_{\max},R}}=0$ (53)

since
$\lim_{t\to\infty}\bar{B}^{\ast}(t)=\lim_{t\to\infty}\bar{\gamma}^{\ast}(t)=0$
by assumption, and $\lim_{t\to\infty}\Delta\left(t-1,T+1\right)=0$ by (16),
both with probability $1$. ∎

##  4\. Applications and implications

Taking a step back, our goal in this section will be to show how a diverse
range of algorithms can be seen as special cases of (RRM), enabling in this
way the use of Theorems 2, 1 and 2 to deduce their convergence properties. To
simplify our presentation, we will make the following technical assumption:

######  Assumption 5.

The injectivity radius of $\mathcal{M}$ is bounded from below by $\delta>0$.

######  Remark.

Recall that the injectivity radius of $\mathcal{M}$ at $p$ is the radius of
the largest geodesic ball on which $\exp_{p}$ is a diffeomorphism, and the
injectivity radius of $\mathcal{M}$ is the infimum over all such radii [40].
In this regard, Assumption 5 serves to ensure that the exponential map is
invertible at consecutive iterates of (RRM) so no local topological
complications can arise; we will in fact prove in Proposition 2 that
$\log_{X_{n}}\coloneqq\exp_{X_{n}}^{-1}$ is well-defined for all sufficiently
large $n$ in the algorithms to follow. $\lozenge$

For concreteness, we will also assume that the algorithms considered in this
section have black-box access to $v$ via a _stochastic first-order oracle_
(SFO) [49]. Specifically, when called at $x\in\mathcal{M}$ with random seed
$\omega\in\Omega$, an stochastic first-order oracle (SFO) returns a random
vector $V(x;\omega)\in\mathcal{T}_{x}\mathcal{M}$ of the form

$V(x;\omega)=v(x)+\operatorname{Err}(x;\omega)$ (SFO)

where the error term
$\operatorname{Err}(x;\omega)\in\mathcal{T}_{x}\mathcal{M}$ is zero-mean and
bounded in $L^{2}$, i.e.,

$\operatorname{\mathbb{E}}[\operatorname{Err}(x;\omega)]=0\quad\text{and}\quad\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{2}_{x}]\leq\sigma^{2}\quad\text{for
all $x\in\mathcal{M}$.}$ (54)

Finally, we will also assume that all algorithms under study are run with a
step-size policy $\gamma_{n}$ such that

$\frac{A}{n}\leq\gamma_{n}\leq\frac{B}{n^{1/2}(\log
n)^{1/2+\varepsilon}}\quad\text{for some $A,B,\varepsilon>0$ and all
$n=1,2,\dotsc$}$ (55)

Clearly, Assumption 1 is satisfied automatically under (55).

######  Remark.

To facilitate comparisons with the Riemannian optimization literature, we will
sometimes refer to queries of the SFO $V(x;\omega)$ as (stochastic)
“gradients”; we stress however that, in general, _$v$ is not assumed to admit
a potential_. $\lozenge$

###  4.1. Basic examples

We now proceed to outline below a – highly incomplete – series of algorithms
that can be seen as special cases of the general template (RRM).

######  Algorithm 1 (RSGM).

The simplest _Riemannian stochastic gradient method_ (RSGM) [9] queries (SFO)
and proceeds as

$X_{n+1}=\exp_{X_{n}}\left(\gamma_{n}V(X_{n};\omega_{n})\right),$ (RSGM)

As such, (RSGM) is an RRM scheme with
$U_{n}=\operatorname{Err}(X_{n};\omega_{n})$ and $b_{n}=0$. $\lozenge$

######  Algorithm 2 (RPPM).

The (deterministic) _Riemannian proximal point method_ (RPPM) [22] is an
implicit update rule of the form

$\exp^{-1}_{X_{n+1}}(X_{n})=-\gamma_{n}v(X_{n+1}).$ (RPPM)

The RRM representation of (RPPM) is then obtained by taking
$b_{n}=\Gamma_{X_{n+1}\to X_{n}}\left(v(X_{n+1})\right)-v(X_{n})$ and
$U_{n}=0$ in the decomposition (2) of the error term $W_{n}$ of (RRM).
$\lozenge$

######  Algorithm 3 (RSEG).

Inspired by the original work of Korpelevich [35], the _Riemannian stochastic
extra-gradient_ (RSEG) method [60, 50] proceeds as

$\displaystyle X_{n+1/2}$
$\displaystyle=\exp_{#1}(\gamma_{n}V(X_{n};\omega_{n})),$ (RSEG)
$\displaystyle X_{n+1}$ $\displaystyle=\exp_{#1}(\Gamma_{X_{n+1/2}\to
X_{n}}\left(\gamma_{n}V(X_{n+1/2};\omega_{n+1/2})\right))$

where $\omega_{n}$ and $\omega_{n+1/2}$ are independent seeds for (SFO). Thus,
to recast (RSEG) in the RRM framework, simply take $U_{n}=\Gamma_{X_{n+1/2}\to
X_{n}}\left(\operatorname{Err}(X_{n+1/2};\omega_{n+1/2})\right)$ and
$b_{n}=\Gamma_{X_{n+1/2}\to X_{n}}\left(v(X_{n+1/2})\right)-v(X_{n})$.
$\lozenge$

######  Algorithm 4 (ROG).

Compared to (RSGM), the scheme (RSEG) involves two oracle queries per
iteration. Building on an original idea by Popov [53], the last oracle query
can be “recycled”, leading to the Riemannian optimistic gradient method

$\displaystyle X_{n+1/2}$
$\displaystyle=\exp_{#1}(\gamma_{n}V(X_{n-1/2};{\omega_{n-1}})),$ (ROG)
$\displaystyle X_{n+1}$ $\displaystyle=\exp_{#1}(\Gamma_{X_{n+1/2}\to
X_{n}}\left(\gamma_{n}V(X_{n+1/2};\omega_{n})\right)).$

With this in mind, (ROG) may be seen as a special case of (RRM) by taking
$U_{n}=\Gamma_{X_{n+1/2}\to
X_{n}}\left(\operatorname{Err}(X_{n+1/2};\omega_{n})\right)$ and
$b_{n}=\Gamma_{X_{n+1/2}\to X_{n}}\left(v(X_{n+1/2})\right)-v(X_{n})$.
$\lozenge$

In view of Theorems 2 and 1, the convergence analysis of Algorithms 1–4
essentially boils down to verifying Assumptions 1–4. Proposition 2 below does
much of the heavy lifting for this:

######  Proposition 2.

Suppose that Assumption 5 holds. If Algorithms 1–4 are run with oracle and
step-size parameters satisfying (54) and (55), then:

  1. (1)

Assumption 2 holds as stated.

  2. (2)

wp1, $X_{n+1}$ lies in the injectivity radius of $X_{n}$ for all sufficiently
large $n$, so Assumption 4 also holds as stated.

######  Corollary 3.

Suppose that Items 1 and 2 hold. If Algorithms 1–4 are run with oracle and
step-size parameters satisfying (54) and (55), then, with probability $1$, the
generated sequence $X_{n}$, $n=1,2,\dotsc$, converges to an ICT set of (RMD).

To streamline our discussion, we postpone the proof of Proposition 2 to
Appendix C.

###  4.2. Algorithmic variants and modifications

To increase the computational efficiency of Riemannian iterative schemes,
several important operations are routinely tacked on to the base algorithms
described in the previous section. In view of this, we proceed below to
illustrate a range of algorithmic variants and modifications that can be
readily incorporated in the general framework of (RRM).

#### Retraction-based methods.

When the exponential map is expensive to compute, a popular alternative is to
employ a so-called _retraction map_ [1, 10], defined here as a smooth mapping
$\mathcal{R}\colon\mathcal{T}\mathcal{M}\to\mathcal{M}$ that agrees with the
exponential map up to first order, i.e.,

${\mathcal{R}}_{x}(0)=x\quad\text{and}\quad\left.\frac{d}{dt}\right|_{t=0}{\mathcal{R}}_{x}(tv)=v\quad\text{for
all $(x,v)\in\mathcal{T}\mathcal{M}$}.$ (Rtr)

As it turns out, to replace the exponential map in Algorithms 1–4 with a
retraction, we only need to slightly strengthen our assumptions on the noise
in (SFO):

######  Proposition 3.

Suppose that the error term $\operatorname{Err}(x;\omega)$ of (SFO) is bounded
in $L^{4}$, i.e.,
$\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]\leq\kappa^{2}$
for some $\kappa>0$ and all $x\in\mathcal{M}$. Then Proposition 2 holds as
stated if the exponential map in Algorithms 1–4 is replaced by a retraction.

######  Corollary 4.

Suppose that Items 1 and 2 hold. If a retraction-based variant of Algorithms
1–4 is run with a step-size schedule satisfying (55) and an SFO with bounded
fourth moments, then, with probability $1$, the generated sequence $X_{n}$,
$n=1,2,\dotsc$, converges to an ICT set of (RMD).

######  Remark 1.

Albeit relatively light, the condition
$\sup_{x}\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]<\infty$
cannot be relaxed. This is true even in the Euclidean case with
${\mathcal{R}}_{x}(v)=x+\left(e^{v}-1\right)$ and
$\operatorname{Err}(x;\omega)\overset{\mathrm{law}}{=}\sqrt{X}-\operatorname{\mathbb{E}}[\sqrt{X}]$,
where $\operatorname{\mathbb{P}}\left(X>x\right)=x^{-3/2}$ if $x\geq 1$, and
$\operatorname{\mathbb{P}}\left(X>x\right)=1$ otherwise. $\lozenge$

######  Remark 2.

In many practical settings, the map ${\mathcal{R}}_{x}(v)$ satisfies the
stronger requirement of being a “second-order retraction”, i.e., it agrees
with $\exp$ up to _second_ order [1, 10]. In this case, the proof technique in
Section C.2 can be used to show that it suffices to have
$\sup_{x\in\mathcal{M}}\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{3}_{x}]<\infty$
in Proposition 3. $\lozenge$

#### Alternating variants.

The next set of variants concerns the case where $v$ is generated from a
$2$-player, min-max game (the $N$-player general-sum case is similar, but the
notation is more involved so we omit it).

To state it, consider a min-max game of the form
$\min_{y\in\mathcal{Y}}\max_{z\in\mathcal{Z}}\ell(y,z)$ played over the smooth
manifolds $\mathcal{Y}$ and $\mathcal{Z}$ (for the min and max player
respectively). Then, instead of performing _simultaneous_ updates for each
player in Algorithms 1–4, a common variant is to _alternate_ variables
according to the basic recursion

$\displaystyle Y_{n+1}$
$\displaystyle=\exp_{#1}(\gamma_{n}[v_{y}(Y_{n},Z_{n})+W_{y,n}])$ (RRM-alt)
$\displaystyle Z_{n+1}$
$\displaystyle=\exp_{#1}(\gamma_{n}[v_{z}(Y_{n+1},Z_{n})+W_{z,n}])$

where $(v_{y},v_{z})\coloneqq(-\nabla_{y}\ell,\nabla_{z}\ell)$, and $W_{y,n}$
and $W_{z,n}$ respectively denote the error terms entering (RRM) for the min
and max player respectively. Our next result shows that the alternating
variants of Algorithms 1–4 retain the convergence properties of their
simultaneous counterparts:

######  Proposition 4.

Suppose that Algorithms 1–4 are run with alternating updates as per (RRM-alt)
and oracle and step-size parameters satisfying (54) and (55). Then Proposition
2 holds as stated.

######  Corollary 5.

Suppose that Items 1 and 2 hold. If an alternating variant of Algorithms 1–4
is run with oracle and step-size parameters satisfying (54) and (55), then,
with probability $1$, the generated sequence $X_{n}$, $n=1,2,\dotsc$,
converges to an ICT set of (RMD).

######  Remark.

A simple but useful observation is that compositions of RRM schemes do not
change its asymptotic behavior: specifically, for any two RRM schemes
${\mathrm{RM}}_{1}$ and ${\mathrm{RM}}_{2}$, the update
$X_{n+1}={\mathrm{RM}}_{2}\circ{\mathrm{RM}}_{1}(X_{n})$ is equivalent to a
new RRM scheme $\tilde{X}_{n}$ where $\tilde{X}_{2n}=X_{n}$ and
$\tilde{X}_{2n+1}={\mathrm{RM}}_{1}(\tilde{X}_{2n})$. This allows us to “mix-
and-match” Propositions 2–4 to prove, for instance, the convergence of
alternating (RSEG) minimizer vs. retraction-based (RPPM) maximizer in
Riemannian two-player games. $\lozenge$

###  4.3. Implications for optimization and learning

We conclude this section with some concrete implications of our general theory
when $v$ is specialized to specific instances that arise in optimization and
learning theory. This allows us to extend and unify several existing results,
but also to obtain completely new ones altogether.

#### Optimization on manifolds.

Perhaps the most common task for learning on manifolds is the basic
minimization problem

$\min\nolimits_{x\in\mathcal{M}}f(x)$ (Opt)

where $f\colon\mathcal{M}\to\mathbb{R}$ is a $C^{d}$-smooth function (not
necessarily geodesically convex). In this case, applying our general theory to
$v=-\nabla f$ yields:

######  Proposition 5.

Suppose that Algorithms 1–4 are run against (Opt) with oracle and step-size
parameters satisfying (54) and (55). Assume further that the problem satisfies
n(itn) Assumptions 3and 5; or, alternatively,  n(itn) Items 1and 2.  Then,
with probability $1$, the induced sequence of iterates $X_{n}$,
$n=1,2,\dotsc$, converges to a component of critical points of $f$ on which
$f$ is constant.

If, in addition,
$\sup_{x}\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]<\infty$,
the above conclusions apply to all retraction-based variants of Algorithms
1–4.

###### Proof.

By Corollary 1, Propositions 2 and 3 (or Propositions 3 and 4 for the
retraction-based case), it follows that $X_{n}$ converges to an ICT set of
(RMD) with probability $1$. Moreover, by Sard’s theorem [57], the set of
critical values of $f$ has empty interior, so Proposition 6.4 of Benaïm [4]
implies that every ICT of (RMD) is contained in a set of critical points of
$f$ on which $f$ is constant. Our assertion then follows by combining the
above. ∎

Proposition 5 contains as a special case the analysis of Bonnabel [9] who
established the almost sure convergence of Algorithm 1 (and its retraction-
based variants) under the assumptions that  n(itn) $X_{n}$is precompact
(Assumption 3);  n(itn) the injectivity radius of $\mathcal{M}$ is uniformly
bounded from below (Assumption 5); and  n(itn) the oracle $V(x;\omega)$ is
uniformly bounded in both $x\in\mathcal{M}$ and $\omega\in\Omega$, i.e.,
$\operatorname{ess}\sup_{x,\omega}\lVert V(x;\omega)\rVert_{x}<\infty$.  In
this regard, Proposition 5 not only relaxes the bounded oracle requirement of
Bonnabel [9], but it provides a straightforward way to establish the
convergence of a wide array of Riemannian algorithms which cannot otherwise be
covered by the tailor-made analysis of Bonnabel [9]. We are not aware of a
comparable convergence result for Algorithms 2–4 (or their retraction-based
variants).

In addition, we note that Proposition 5 includes the celebrated _natural
policy gradient_ (NPG) method of Kakade [31] for reinforcement learning –
which, as noted by Bonnabel [9], is a particular case of retraction-based
Algorithm 1. We stress that our result applies not only to “vanilla” natural
policy gradient (NPG) methods, but also to its optimistic/extra-gradient
variants. To our knowledge, there are no comparable results in the NPG
literature.

#### Games on manifolds.

We next move on to the game setting, i.e., when $\mathcal{M}$ decomposes as a
product of the form
$\mathcal{M}=\mathcal{M}_{1}\times\dotsm\times\mathcal{M}_{N}$ for some
$N\in\mathbb{N}$ and, likewise, the $i$-th component of
$v=(v_{1},\dotsc,v_{N})$ is of the form
$v_{i}=\operatorname{\nabla}_{x_{i}}u_{i}(x_{1},\dotsc,x_{N})$, where
$u_{i}\colon\mathcal{M}\to\mathbb{R}$ is the payoff function of the $i$-th
player – for applications and a detailed discussion, cf. Ratliff et al. [54]
and references therein.

There are many solution concepts of games on manifolds that can be seen as a
natural generalization of Nash equilibria. In this work, we will consider the
so-called _Nash–Stampacchia equilibria_ (NSE) [36, 47]:

######  Definition 3.

We say that
$x^{\ast}\coloneqq(x^{\ast}_{1},\dotsc,x^{\ast}_{N})\in\mathcal{M}$ is a
_Nash–Stampacchia_ equilibrium of $v$ if

$\langle
v(x^{\ast}),\exp^{-1}_{x^{\ast}}(x)\rangle_{x^{\ast}}\coloneqq\sum\nolimits_{i=1}^{N}\langle
v_{i}(x^{\ast}),\exp^{-1}_{x^{\ast}_{i}}(x_{i})\rangle_{x^{\ast}_{i}}\geq
0,\quad\text{for all $x\in\mathcal{M}$}.$

Our first result below concerns the convergence of Algorithms 2–4 in a general
class of (Riemannian) monotone games known as $\alpha$-accretive games [62]:

######  Proposition 6 (Riemannian $\alpha$-accretive games).

Let $v=[-\nabla_{x_{i}}\ell_{i}]$ be an $\alpha$-accretive game field, i.e.,
all $r\geq 0$, we have

$(1+\alpha
r)\operatorname{dist}(x,x^{\prime})\leq\operatorname{dist}\left(\exp_{x}(rv(x)),\exp_{x^{\prime}}(rv(x^{\prime}))\right),\quad\alpha>0.$

Then Algorithms 1–4, as well as their alternating / retraction-based versions,
converge to the game’s set of Nash–Stampacchia equilibria.

###### Proof.

This immediately follows by combining Propositions 2–4 with the main result of
[36]. ∎

To the best of our knowledge, most of the algorithms we consider are new in
the setting of Riemannian monotone games except for the _deterministic_
gradient and extra-gradient methods [60, 50, 21, 32, 14].

While being quite general, accretivity is a strong, convexity-like assumption
about the games. In our next result, we prove general convergence for a class
of _non-convex_ games.

######  Proposition 7 (Riemannian potential games).

Let $v=[-\nabla_{x_{i}}\ell_{i}]$ be a game field associated with a Riemannian
potential game [34, 48]. Then Algorithms 1–4, as well as their alternating /
retraction-based versions, converge to the critical points of the game
potential.

###### Proof.

Simply combine Propositions 2–4. ∎

For Riemannian potential games, the convergence of the continuous-time
dynamics (RMD) is well known, but we are not otherwise aware of a similar
result for stochastic, discrete-time RRM methods. Our theory bridges this gap
by showing that the same guarantees are in fact achieved by a wide array of
RRM schemes – not only by their continuous-time ancestor.

#### Limit cycles in Riemannian manifolds.

We conclude this section by showing that, in complement to the pointwise
convergence results above, our theory can also be used to derive convergence
to limit cycles that arise in more general Riemannian settings.

######  Example 4.1 (Cycling of RRM schemes).

The following example is taken from [17]: Consider the vector field on
$\mathbb{S}^{2}\coloneqq\\{x\in\mathbb{R}^{3}:x^{2}_{1}+x_{2}^{2}+x_{3}^{3}=1\\}$,
defined by

$v(x)=\begin{bmatrix}-x_{2}\\\ x_{1}\\\
0\end{bmatrix}+\left(x_{3}^{2}-\frac{1}{4}\right)\begin{bmatrix}x_{1}x_{3}\\\
x_{2}x_{3}\\\ -x_{1}^{2}-x_{2}^{2}\end{bmatrix}.$ (56)

Then it is known that the associated (RMD) _cycles_ in the sense that the ICT
sets of $v$ in (56) contain attracting periodic orbits; see Fig. 2. Our
Propositions 2–4 then imply that any RRM scheme driven by $v$ also cycles. To
the best of our knowledge, this is the first rigorous example of a cycling
problem for Riemannian stochastic approximation in the literature.

Figure 2. The flow of (RMD) induced by the game field (56).

##  5\. Concluding remarks

Our theory provides a unified analysis for the convergence of Riemannian
Robbins–Monro schemes that might seem vastly different from each other at
first sight; instead, by verifying certain simple criteria on the error terms
$W_{n}$ in Assumption 2, our analysis shows that we can study the
deterministic dynamics (RMD) to directly infer the algorithm’s long-run
behavior. At the same time, our results offer but a glimpse of the flexibility
of (RRM), and several important research directions remain open:

  1. (1)

In many applications (especially to game theory and sequential online
learning), oracle access to $v$ may be out of reach, and one would need to
employ _zeroth-order_ – or _bandit_ – optimization methods. It this case, a
key question that arises is whether a Riemannian _Kiefer–Wolfowitz_ algorithm
can be analyzed as a special case of (RRM) – and, if so, whether there are any
fundamental differences relative to the Euclidean setting.

  2. (2)

The diminishing step-size assumption is indispensable for our analysis, which
covers many practically relevant settings. However, another common strategy
involves _constant_ step-sizes, and this is not covered by our theory. It
would thus be interesting to see if the techniques for analyzing constant
step-size SA schemes [38, 39] can be generalized to a manifold setting.

  3. (3)

Finally, several Riemannian algorithms are known to escape undesirable
solutions [15, 58]. We conjecture that the general avoidance theory of [51, 5,
12, 27] can be likewise extended to Riemannian manifolds; if true, this would
imply that many iterative Riemannian methods (including retraction-based ones)
converge with probability $1$ _only_ to local minimizers.

We defer the study of the above questions to future work.

## 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 Proposition 1. To that end, let
$p\in\mathcal{M}$ be as in (WC), define the radial distance function
$r(x)\coloneqq\operatorname{dist}(x,p)$, and let $k(x)\coloneqq r^{2}(x)/2$.
The following theorem makes it clear under which assumptions $k$ is smooth and
provides a control on its Hessian.

######  Theorem 3 (Jost [30, Theorem 6.6.1]).

With notation as above, suppose that the exponential map $\exp_{p}$ is a
diffeomorphism on $\\{v\in\mathcal{T}_{x}\mathcal{M}:\lVert
v\rVert_{x}\leq\rho\\}$. Moreover, suppose that the sectional curvature of
$\mathcal{M}$ is nonpositive and bounded below by $-\kappa^{2}$ on a geodesic
ball of radius $\rho>0$ centered at $p$. Then $k$ is smooth on the punctured
ball $\mathbb{B}_{\rho}^{\ast}(p)\coloneqq\\{x\neq
p:\operatorname{dist}(x,p)\leq\rho\\}$ and, in particular,

$\operatorname{\nabla}k(x)=-\log_{p}(x)$ (A.1)

In addition, $\lVert\operatorname{\nabla}k(x)\rVert=r(x)$, and

$\operatorname{Hess}k(x)[v,v]\leq\kappa r(x)\operatorname{ctgh}(\kappa
r(x))\lVert v\rVert_{x}^{2},\quad\text{for all
$x\in\mathbb{B}_{\rho}^{\ast}(p)$ and $v\in\mathcal{T}_{x}\mathcal{M}$}.$
(A.2)

######  Remark.

Since $\mathcal{M}$ is simply connected and complete – by Item 1 – we may take
$\rho=\infty$ in the theorem above [30, Corollary 6.9.1].

Now, to proceed, recall that $v$ satisfies the weak coercivity condition (WC)
if, for some $R>0$, we have

$\langle v(x),\operatorname{\nabla}k(x)\rangle_{x}\leq 0$ (A.3)

whenever $x\in\mathcal{M}$ lies outside a geodesic $\mathbb{B}_{R}(p)$
contained in the interior of $\mathcal{M}$. 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

$\Phi(x)=f(r(x))$ (A.4)

where $f$ is a $\mathcal{C}^{\infty}$ non-negative function with $f(x)=0$ for
all $x\leq R$ and satisfies

$0\leq f^{\prime}(x)\leq C_{1},\quad f^{\prime\prime}(x)\leq C_{2}$ (A.5)

for all $x\geq R$. Moreover, we require $f(x)=\Omega(x)$ as $x\to\infty$ so
that controlling $f$ implies control of $x$ (for a concrete example of such a
function, cf. Lemma A.2). Our first result is that $\Phi=f\circ r$ has a
bounded Hessian and is smooth.

######  Lemma A.1.

Let $\Phi$ be defined as above. Then $\Phi$ is negatively correlated with $v$
in the sense that

$\langle\operatorname{\nabla}\Phi(x),v(x)\rangle_{x}\leq 0\quad\text{for all
$x\in\mathcal{M}$}.$ (A.6)

Moreover, there exists a constant $C>0$ such that
$\operatorname{Hess}\Phi(x)[v,v]\leq C\lVert v\rVert_{x}^{2}$, and hence

$\Phi(x^{\prime})\leq\Phi(x)+\langle\operatorname{\nabla}\Phi(x),\log_{x}(x^{\prime})\rangle+\frac{C}{2}\operatorname{dist}^{2}(x,x^{\prime})\quad\text{for
all $x,x^{\prime}\in\mathcal{M}$}.$ (A.7)

###### Proof.

We begin by noting that the gradient of $\Phi$ is given by

$\operatorname{\nabla}\Phi(x)=\begin{cases}0&\quad\text{if $r(x)\leq R$},\\\
\frac{f^{\prime}(r(x))}{r(x)}\operatorname{\nabla}k(x)&\quad\text{if
$r(x)>R$}.\end{cases}$ (A.8)

By assumption, $f^{\prime}(r(x))/r(x)\geq 0$ so
$\langle\operatorname{\nabla}\Phi(x),v(x)\rangle_{x}$ has the same sign as
$\langle\operatorname{\nabla}k(x),v(x)\rangle_{x}$ if $r(x)>R$ and is
otherwise zero if $r(x)\leq R$. We thus conclude that $\Phi$ and $v$ are
negatively correlated, as claimed.

Now, to compute the Hessian of $\Phi$, notice that
$\operatorname{Hess}\Phi(x)[v,v]=\langle\operatorname{\nabla}_{v}\operatorname{\nabla}\Phi(x),v\rangle_{x}$.
Hence,

$\displaystyle\operatorname{Hess}\Phi(x)[v,v]$
$\displaystyle=\operatorname{\nabla}_{v}\frac{f^{\prime}(r(x))}{r(x)}\cdot\langle\operatorname{\nabla}k(x),v\rangle_{x}+\frac{f^{\prime}(r(x))}{r(x)}\langle\operatorname{\nabla}_{v}\operatorname{\nabla}k(x),v\rangle_{x}$
$\displaystyle=\underbrace{\left\langle\operatorname{\nabla}\frac{f^{\prime}(r(x))}{r(x)},v\right\rangle\cdot\langle\operatorname{\nabla}k(x),v\rangle_{x}}_{\mathrm{I}}+\underbrace{\frac{f^{\prime}(r(x))}{r(x)}\,\operatorname{Hess}k(x)[v,v]}_{\mathrm{II}}.$
(A.9)

Here we use the same notation for directional derivative of a scalar function
and the covariant derivative. With this in mind, the first step in computing
$\mathrm{I}$ is the observation that

$\operatorname{\nabla}\frac{f^{\prime}(r(x))}{r(x)}=\left(f^{\prime\prime}(r(x))-\frac{f^{\prime}(r(x))}{r(x)}\right)\frac{1}{r^{2}(x)}\operatorname{\nabla}k(x),$
(A.10)

and hence

$\mathrm{I}=\left(f^{\prime\prime}(r(x))-\frac{f^{\prime}(r(x))}{r(x)}\right)\frac{1}{r^{2}(x)}\langle\operatorname{\nabla}k(x),v\rangle_{x}^{2}\leq\frac{C_{2}}{r^{2}(x)}\lVert\operatorname{\nabla}k(x)\rVert_{x}^{2}\lVert
v\rVert_{x}^{2}=C_{2}\,\lVert v\rVert_{x}^{2}.$ (A.11)

For $\mathrm{II}$, as $x\operatorname{ctgh}x\leq 1+x$ for $x\geq 0$, we obtain

$\mathrm{II}\leq\frac{f^{\prime}(r(x))}{r(x)}(1+\kappa r(x))\lVert
v\rVert_{x}^{2}\leq C_{1}\left(1/R+\kappa\right)\lVert v\rVert_{x}^{2}.$
(A.12)

Summing up everything, we obtain

$\operatorname{Hess}\Phi(x)[v,v]\leq(C_{2}+C_{1}/R+C_{1}\kappa)\lVert
v\rVert_{x}^{2}\eqqcolon C\lVert v\rVert_{x}^{2},$ (A.13)

that is, $\Phi$ has bounded Hessian. Moreover, $\Phi$ is smooth as a
composition of smooth functions. Let $x,x^{\prime}\in\mathcal{M}$ be
arbitrary, and let $\gamma:[0,1]\to\mathcal{M}$ be a geodesic connecting the
two. By Taylor’s remainder theorem, there exists some $t\in(0,1)$ such that

$\Phi(x^{\prime})=\Phi(x)+\langle\operatorname{\nabla}\Phi(x),\dot{\gamma}(0)\rangle_{x}+\frac{1}{2}\operatorname{Hess}\Phi(\gamma(t))[\dot{\gamma},\dot{\gamma}].$
(A.14)

Thus, invoking (A.13) and noting that
$\lVert\dot{\gamma}\rVert=\operatorname{dist}(x,x^{\prime})$ and
$\dot{\gamma}(0)=\log_{x}(x^{\prime})$, we obtain (A.7) and our proof is
complete. ∎

We now proceed to the main argument, where we show how to use $\Phi$ to
control the iterates $X_{n}$. Letting $\Phi_{n}=\Phi(X_{n})$ and using Lemma
A.1 we get

$\displaystyle\Phi_{n+1}$
$\displaystyle=\Phi\left(\exp_{X_{n}}(\gamma_{n}V_{n})\right)$
$\displaystyle\leq\Phi_{n}+\gamma_{n}\langle\operatorname{\nabla}\Phi(X_{n}),V_{n}\rangle_{X_{n}}+\frac{C\gamma_{n}^{2}}{2}\lVert
V_{n}\rVert_{X_{n}}^{2}$
$\displaystyle\leq\Phi_{n}+\gamma_{n}\langle\operatorname{\nabla}\Phi(X_{n}),U_{n}+b_{n}\rangle_{X_{n}}+\frac{3C\gamma_{n}^{2}}{2}\big{[}\lVert
v(X_{n})\rVert_{X_{n}}^{2}+\lVert U_{n}\rVert_{X_{n}}^{2}+\lVert
b_{n}\rVert_{X_{n}}^{2}\big{]}$ (A.15)

where the second line follows from the negative correlation of $\Phi$ and $v$,
the definition (4) of $V_{n}$, and the Cauchy-Schwarz inequality. Conditioning
on $\mathcal{F}_{n}$ and taking expectations, and invoking Cauchy-Schwarz and
the fact that
$\lVert\operatorname{\nabla}\Phi(X_{n})\rVert\leq\frac{C_{1}}{r(X_{n})}\lVert\operatorname{\nabla}k(X_{n})\rVert=C_{1}$,
we obtain:

$\operatorname{\mathbb{E}}[\Phi_{n+1}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]\leq\Phi_{n}+\gamma_{n}C_{1}B_{n}+\tfrac{3}{2}C\gamma_{n}^{2}\big{[}G^{2}+B_{n}^{2}+\sigma_{n}^{2}\big{]},$
(A.16)

where we have bounded the second moments by their respective upper bounds.

To proceed, let
$\varepsilon_{n}=\gamma_{n}C_{1}B_{n}+(3/2)C\gamma_{n}^{2}\big{[}G^{2}+B_{n}^{2}+\sigma_{n}^{2}\big{]}$
denote the “residual” term in (A.16). Notice that

$\sum_{n=1}^{\infty}\varepsilon_{n}\leq
C_{1}\sum_{n=1}^{\infty}\gamma_{n}B_{n}+\frac{3C}{2}\sum_{n=1}^{\infty}\gamma_{n}^{2}(G^{2}+B_{n}^{2}+\sigma_{n}^{2}),$
(A.17)

and hence, by Assumption 2 and the dominated convergence theorem, we infer
that $\operatorname{\mathbb{E}}[\sum_{n}\varepsilon_{n}]<\infty$.

Next, consider the auxiliary process
$\hat{\Phi}_{n}=\Phi_{n}+\operatorname{\mathbb{E}}[\sum_{k=n}^{\infty}\varepsilon_{k}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]$,
adapted to the same filtration. By (A.16), we have
$\operatorname{\mathbb{E}}[\hat{\Phi}_{n+1}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]\leq\Phi_{n}+\operatorname{\mathbb{E}}[\sum_{k=n}^{\infty}\varepsilon_{k}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]=\hat{\Phi}_{n-1}$,
i.e., $\hat{\Phi}_{n}$ is a supermartingale with respect to $\mathcal{F}_{n}$.
This shows that
$\operatorname{\mathbb{E}}[\hat{\Phi}_{n}]\leq\operatorname{\mathbb{E}}[\hat{\Phi}_{1}]<\infty$,
i.e., $\hat{\Phi}_{n}$ is uniformly bounded in $L^{1}$. Hence, by Doob’s
supermartingale convergence theorem [25, Theorem 2.5], it follows that
$\hat{\Phi}_{n}$ converges with probability $1$ to some finite random limit
$\hat{\Phi}_{\infty}$. In turn, since $\sum_{n}\varepsilon_{n}<\infty$ (a.s.),
this implies that
$\Phi_{n}=\hat{\Phi}_{n}-\operatorname{\mathbb{E}}[\sum_{k=n}^{\infty}\varepsilon_{n}\nonscript\,|\nonscript\,\mathopen{}\mathcal{F}_{n}]$
also converges to some (random) finite limit (a.s.). From this and the fact
that $\Phi_{n}=\Omega(r(X_{n}))$, we deduce $\limsup_{n}r(X_{n})<\infty$ as
claimed.

######  Lemma A.2.

Let $h\colon\mathbb{R}\to\mathbb{R}$ be the function

$h(x)=\begin{cases}0&\quad\text{if $x\leq 0$}\\\
\frac{e^{-1/x}}{e^{-1/x}+e^{-1/(1-x)}}&\quad\text{if $x\in(0,1)$}\\\
1&\quad\text{otherwise}\end{cases}$ (A.18)

and, for $R>0$, let $f(x)=\int_{0}^{x}h(s-R)\>ds.$ Then $f$ is $C^{\infty}$
and it satisfies the conditions (A.5) with $C_{1}=1$ and $C_{2}=2$. In
addition, one has $f(x)\geq x-(R+1)$, and hence $f(x)=\Omega(x)$.

###### Proof.

As $h(x)\in[0,1]$, we obtain that $f^{\prime}(x)\in[0,1]$. By a
straightforward computation, one observes that the first derivative of $h$ is
bounded as $0\leq h^{\prime}(x)\leq 2$, so
$f^{\prime\prime}(x)=h^{\prime}(x-R)\leq 2.$ Then, to complete our proof,
simply notice that, for $x\geq R+1$, we have
$f(x)=\int_{0}^{x}h(s-R)\>ds\geq\int_{R+1}^{x}1\>ds=x-(R+1)$, as claimed. ∎

##  Appendix B Geometric preliminaries for the proof of Theorem 2

In this appendix, we collect the necessary technical prerequisites for the
proof of Theorem 2.

###  B.1. The parallel frame system

We begin with the definition of the frame system that we use to compare
vectors on different tangent spaces. To that end, fix any two points
$x,x^{\prime}\in\mathcal{M}$, and consider two arbitrary vectors
$v\in\mathcal{T}_{x}\mathcal{M}$ and
$w\in\mathcal{T}_{x^{\prime}}\mathcal{M}$. There is a convenient frame system
(i.e., a set of bases for $\mathcal{T}_{x}\mathcal{M}$ and
$\mathcal{T}_{x^{\prime}}\mathcal{M}$) for comparing $v$ and $w$, defined as
follows: Pick an arbitrary orthonormal frame
$\left\\{e_{k}\right\\}_{k=1}^{d}$ for $\mathcal{T}_{x}\mathcal{M}$. Since the
parallel transport map is an isometry, the vectors
$\\{e^{\prime}_{k}\\}_{k=1}^{d}\coloneqq\left\\{\Gamma_{x\to
x^{\prime}}\left(e_{k}\right)\right\\}_{k=1}^{d}$ form an orthonormal basis
for $\mathcal{T}_{x^{\prime}}\mathcal{M}$. Now consider the components of
$v,w$ in these two frames, namely

$\displaystyle v_{k}^{\shortparallel}\coloneqq\langle
v,e_{k}\rangle_{x}\;\text{and}\;w_{k}^{\shortparallel}\coloneqq\langle
w,e^{\prime}_{k}\rangle_{x^{\prime}}\quad\text{for all $k=1,\dotsc,d$}.$ (B.1)

We shall call (B.1) the _parallel frame system_ for vectors $v$ and $w$ (the
dependence on the initial frame $\left\\{e_{k}\right\\}_{k=1}^{d}$ is
suppressed).

By virtue of parallel transport, in the parallel frame system we have

$\displaystyle\Gamma_{x\to
x^{\prime}}\left(v\right)_{k}^{\shortparallel}=\langle\Gamma_{x\to
x^{\prime}}\left(v\right),e^{\prime}_{k}\rangle_{x^{\prime}}=\langle
v,\Gamma_{x^{\prime}\to
x}\left(e^{\prime}_{k}\right)\rangle_{x}=v_{k}^{\shortparallel}.$ (B.2)

In the same vein, we have $\Gamma_{x^{\prime}\to
x}\left(w\right)_{k}^{\shortparallel}=w_{k}^{\shortparallel}$. Thus, since
$\left\\{e_{k}\right\\}_{k=1}^{d}$ is orthonormal, we may write:

$\lVert w-\Gamma_{x\to x^{\prime}}\left(v\right)\rVert_{x}=\lVert
v-\Gamma_{x^{\prime}\to x}\left(w\right)\rVert_{x^{\prime}}=\lVert
v^{\shortparallel}-w^{\shortparallel}\rVert_{2}.$ (B.3)

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

$\lVert v^{\shortparallel}(x^{\prime})-v^{\shortparallel}(x)\rVert_{2}\leq
L\operatorname{dist}(x,x^{\prime})\quad\text{for all
$x,x^{\prime}\in\mathcal{M}$}.$ (B.4)

###  B.2. The Fermi coordinate system

For any $h$, let
$\mathcal{U}_{h}\subset\mathcal{T}_{X(t+h)}\mathcal{M}\simeq\mathbb{R}^{d}$ be
a neighborhood of 0 on which the mapping

$\exp_{X(t+h)}\colon\mathcal{U}_{h}\to\mathcal{M}$ (B.5)

is a diffeomorphism between $\mathcal{U}_{h}$ and
$\exp_{X(t+h)}\left(\mathcal{U}_{h}\right)$. It is well-known that such a
neighborhood exists, and that the exponential map $\exp_{X(t+h)}$ along with
an arbitrary orthonormal frame at $\mathcal{T}_{X(t+h)}\mathcal{M}$ induces a
local coordinate system on $\exp_{X(t+h)}\left(\mathcal{U}_{h}\right)$, called
the _normal coordinate frame with center $X(t+h)$_ [40]. Normal coordinates
are best-suited for comparing _distances_ of points on manifolds: for
instance, if $\tilde{x}^{\prime}$ is the normal coordinate of $x^{\prime}$
with center $x$, then
$\operatorname{dist}(x,x^{\prime})=\lVert\tilde{x}^{\prime}\rVert_{2}$.

The _Fermi coordinate system_ (FCS) [45], roughly speaking, is a system of
normal coordinates “along a curve”. To define it, fix an arbitrary orthonormal
frame $\left\\{e_{k}(0)\right\\}_{k=1}^{d}$ for
$\mathcal{T}_{X(t)}\mathcal{M}$. We can obtain a system of orthonormal frames
$\left\\{e_{k}(h)\right\\}_{k=1}^{d}$ by parallel transporting
$\left\\{e_{k}(0)\right\\}_{k=1}^{d}$ from $\mathcal{T}_{X(t)}\mathcal{M}$ to
$\mathcal{T}_{X(t+h)}\mathcal{M}$ along the curve $h\mapsto X(t+h)$. Let
$\mathcal{U}_{h}\subset\mathcal{T}_{X(t+h)}\mathcal{M}$ be a neighborhood of 0
defined as in (B.5), and set
$\tilde{\mathcal{U}}\coloneqq\operatorname*{\bigcup}_{h}\\{\exp_{X(t+h)}\left(\mathcal{U}_{h}\right)\\}\subset\mathcal{M}$.
Finally, consider the mapping

$\tilde{\Phi}\colon\mathbb{R}_{+}\times\tilde{\mathcal{U}}\to\mathbb{R}_{+}\times\mathbb{R}^{d}$
(B.6)

that sends a point $(h,x)\in\mathbb{R}_{+}\times\tilde{\mathcal{U}}$ to
$(h,\tilde{x})\in\mathbb{R}_{+}\times\mathbb{R}^{d}$, where $\tilde{x}$ is the
normal coordinate of $x$ with center $X(t+h)$ and frame
$\left\\{e_{k}(h)\right\\}_{k=1}^{d}$. By virtue of the normal coordinates, we
know that $\tilde{\Phi}$ is a diffeomorphism between
$\mathbb{R}_{+}\times\tilde{\mathcal{U}}$ and a neighborhood of
$\mathbb{R}_{+}\times\\{0\\}$. The mapping $\tilde{\Phi}$ and its inverse is
called the _Fermi coordinate system_ along the curve $h\mapsto X(t+h)$. In the
sequel, we will abuse the notation and simply call it the Fermi coordinates
along $X$.

The following property of the Fermi coordinate system plays a key role in our
analysis and we will use it freely:

######  Lemma B.1 (Takahashi and Watanabe, 1981; Fujita and Kotani, 1982).

Let $\gamma$ be a differentiable curve on $\mathcal{M}$ such that
$\gamma(h)\in\mathcal{U}_{h}$ for all $h\in\mathbb{R}_{+}$, and let
$\tilde{\gamma}$ be the curve of $\gamma$ in the Fermi coordinate system (FCS)
along $X$ (i.e., $(h,\tilde{\gamma}(h))=\tilde{\Phi}(h,\gamma(h))$. Then

$\dot{\tilde{\gamma}}_{k}(h)=\dot{\gamma}_{k}(h)-\dot{X}_{k}(t+h)+\operatorname{\mathcal{O}}\big{(}\lVert\tilde{\gamma}(h)\rVert_{2}^{2}\big{)}$
(B.7)

where $\dot{X}_{k}(t+h)\coloneqq\langle\dot{X}(t+h),e_{k}(h)\rangle_{X(t+h)}$
is the $k$-th component of $\dot{X}(t+h)$ in the frame
$\left\\{e_{k}(h)\right\\}_{k=1}^{d}$, and $\dot{\gamma}_{k}(h)$ is the $k$-th
component of $\dot{\gamma}(h)$ in the (possibly non-orthonormal) frame induced
by the normal coordinate system with center $X(t+h)$ and frame
$\left\\{e_{k}(h)\right\\}_{k=1}^{d}$.

###  B.3. Comparing the differential of $\exp$ and parallel transport

As will become clear in the proof, the parallel frame system (PFS) is
convenient for comparing _vectors_ at different points, whereas the FCS 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.6) is simultaneously a PFS for all points nearby
the curve $X(\cdot)$, then the underlying manifold $\mathcal{M}$ must be flat;
i.e., the Riemannian curvature tensor vanishes everywhere [29].

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 [41, Theorem 3.12] or [16, Proposition A.1]:

######  Lemma B.2 (Comparing ${\mathrm{d}}\exp$ and parallel transport).

Let $\mathcal{M}$ be a Riemannian manifold whose sectional curvatures are in
the interval $[K_{\textup{low}},K_{\textup{up}}]$, and let $K=\max\left(\lvert
K_{\textup{low}}\rvert,\lvert K_{\textup{up}}\rvert\right)$. For
$v\in\mathcal{T}_{x}\mathcal{M}$, consider the geodesic
$\gamma(t)=\exp_{x}(tv)$. If $\gamma$ is defined and has no interior conjugate
point on the interval $[0,1]$, then

$\forall w\in\mathcal{T}_{x}\mathcal{M},\quad\quad\lVert
T_{v}(w)-\Gamma_{v}(w)\rVert_{\gamma(1)}\leq K\cdot
f_{K_{\textup{low}}}(\lVert v\rVert_{x})\cdot\lVert w_{\perp}\rVert_{x}$ (B.8)

where $w_{\perp}\coloneqq w-\frac{\langle v,w\rangle_{x}}{\langle
v,v\rangle_{x}}v$ is the component of $w$ orthogonal to $v$,
$T_{v}={\mathrm{d}}\exp_{x}(v)$ is the differential of the exponential map,
and $\Gamma_{v}$ denotes the parallel transport along $\gamma$ from
$\gamma(0)$ to $\gamma(1)$. The function $f_{K_{\textup{low}}}$ in (B.8) is
defined as

$f_{K_{\textup{low}}}(a)=\begin{cases}\frac{r^{2}}{6}&\quad\text{if
$K_{\textup{low}}=0$},\\\
r^{2}\left(1-\frac{\sin({a}/{r})}{{a}/{r}}\right)&\quad\text{if
$K_{\textup{low}}=\frac{1}{r^{2}}>0$},\\\
r^{2}\left(\frac{\sinh({a}/{r})}{{a}/{r}}-1\right)&\quad\text{if
$K_{\textup{low}}=-\frac{1}{r^{2}}<0$}.\end{cases}$ (B.9)

Moreover, the function $f_{K_{\textup{low}}}$ is dominated by the case
$K_{\textup{low}}<0$: For all $K\geq\lvert K_{\textup{low}}\rvert$ and
$a\in\mathbb{R}_{+}$,

$f_{K_{\textup{low}}}(a)\leq f_{-K}(a).$ (B.10)

##  Appendix C Proofs for Section 4

###  C.1. Proof of Proposition 2

We prove Proposition 2 in this appendix. To this end, we first provide a
convenient lemma which shows that, almost surely, the effect of the noise is
asymptotically annihilated by the step-size:

######  Lemma C.1.

Under the assumptions in Proposition 2, we have, with probability 1,

$\lim_{n\to\infty}\left\lVert\gamma_{n}V(X_{n};\omega_{n})\right\rVert_{X_{n}}=0.$
(C.1)

###### Proof.

By definition,

$\displaystyle\left\lVert\gamma_{n}V(X_{n};\omega_{n})\right\rVert_{X_{n}}$
$\displaystyle=\left\lVert\gamma_{n}v(X_{n})+\operatorname{Err}(X_{n};\omega_{n})\right\rVert_{X_{n}}$
$\displaystyle\leq\gamma_{n}G+\gamma_{n}\left\lVert\operatorname{Err}(X_{n};\omega_{n})\right\rVert_{X_{n}}.$
(C.2)

The first term goes to 0 by choice of step-sizes. To control the second term,
note that by Chebyshev’s inequality and (54), we have

$\operatorname{\mathbb{P}}\left({\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert_{X_{n}}\geq\sqrt{n\log^{1+\frac{\epsilon}{2}}n}}\right)\leq\frac{\sigma^{2}}{n\log^{1+\frac{\epsilon}{2}}n}$
(C.3)

where $\varepsilon$ is the same as in our choice of step-size in Proposition
2. In turn, this implies that

$\sum_{n=2}^{\infty}\operatorname{\mathbb{P}}\left({\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert_{X_{n}}\geq\sqrt{n\log^{1+\frac{\epsilon}{2}}n}}\right)<\infty$

so, by the Borel-Cantelli lemma, we have
$\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert_{X_{n}}=\operatorname{\mathcal{O}}\left(\sqrt{n\log^{1+\frac{\epsilon}{2}}n}\right)$
with probability $1$. Hence, by our assumptions for the method’s step-size, we
get

$\gamma_{n}\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert_{X_{n}}=\operatorname{\mathcal{O}}\left(\frac{\sqrt{n\log^{1+\frac{\epsilon}{2}}n}}{\sqrt{n\log^{1+\varepsilon}n}}\right)=\operatorname{\mathcal{O}}\left(\frac{1}{\log^{\frac{\varepsilon}{4}}n}\right)$

which, combined with (C.2), implies (C.1). ∎

We are now ready for the full proof.

The second claim of Proposition 2 is a direct consequence of Lemma C.1 and
Assumption 5. As for the first, note that $\sum_{n}\gamma_{n}^{2}<\infty$ by
our choice of step-sizes so Assumption 1 holds by construction and we are left
to establish the summability conditions (7) of Assumption 2. To this end, by
Assumption 1, (54), and Lemma C.1, it suffices to show that $\lVert
b_{n}\rVert_{X_{n}}=\operatorname{\mathcal{O}}(\gamma_{n}\cdot\text{noise})$
where “noise” is a query from (SFO) at the appropriate state. We proceed
method-by-method:

####  Algorithm 1: RSGM.

For (RSGM), we have $W_{n}=U_{n}=\operatorname{Err}(X_{n};\omega_{n})$ and
$b_{n}=0$, so (7) follows from the stated assumptions for (SFO).

####  Algorithm 2: RPPM.

For (RPPM), we have $U_{n}=0$ and

$\displaystyle\lVert b_{n}\rVert_{X_{n}}$
$\displaystyle=\lVert\Gamma_{X_{n+1}\to
X_{n}}\left(v(X_{n+1})\right)-v(X_{n})\rVert_{X_{n}}$ $\displaystyle\leq
L\operatorname{dist}\left(X_{n},X_{n+1}\right)=\gamma_{n}L\lVert
v(X_{n})\rVert_{X_{n}}\leq\gamma_{n}LM=\operatorname{\mathcal{O}}(\gamma_{n})$
(C.4)

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: RSEG.

For the Riemannian stochastic extra-gradient (RSEG) algorithm, the recurrence
(RSEG) gives $U_{n}=\Gamma_{X_{n+1/2}\to
X_{n}}\left(\operatorname{Err}(X_{n+1/2};\omega_{n+1/2})\right)$ so

$\operatorname{\mathbb{E}}\big{[}\left\lVert
U_{n}\right\rVert^{2}_{X_{n}}\big{]}=\operatorname{\mathbb{E}}\big{[}\left\lVert\operatorname{Err}(X_{n+1/2};\omega_{n+1/2})\right\rVert^{2}_{X_{n+1/2}}\big{]}\leq\sigma^{2}$
(C.5)

by (54) and the fact that the parallel transport map is a linear isometry.
Finally, for the second part of (7), the definition of (RSEG) yields

$\lVert b_{n}\rVert=\lVert\Gamma_{X_{n+1/2}\to
X_{n}}\left(v(X_{n+1/2})\right)-v(X_{n})\rVert_{X_{n}}\leq
L\operatorname{dist}\left(X_{n+1/2},X_{n}\right)=\gamma_{n}L\lVert
V(X_{n};\omega_{n})\rVert_{X_{n}}$ (C.6)

so our claim follows from the assumptions on $V$.

####  Algorithm 4: ROG.

For the ROG algorithm, the recurrence (ROG) gives $U_{n}=\Gamma_{X_{n+1/2}\to
X_{n}}\left(\operatorname{Err}(X_{n};\omega_{n+1/2})\right)$ and
$b_{n}=\Gamma_{X_{n+1/2}\to X_{n}}\left(v(X_{n+1/2})\right)-v(X_{n})$, so
Assumption 2 can be checked exactly as in the case of Algorithm 3 above.

###  C.2. Proof of Proposition 3

By definition, ${\mathcal{R}}_{x}(v)$ is a smooth map and hence satisfies
$\lim_{v\to 0}{\mathcal{R}}_{x}(v)=x$. Then Lemma C.1 readily implies that
$X_{n+1}$ lies in the injectivity radius of $X_{n}$ with probability 1 for $n$
large enough.

We first consider the retraction-based Algorithm 1:

$\displaystyle X_{n+1}$
$\displaystyle={\mathcal{R}}_{X_{n}}(\gamma_{n}V(X_{n};\omega_{n})).$ (C.7)

Let $\tilde{v}_{n}\in\mathcal{T}_{X_{n}}\mathcal{M}$ be the vector such that
$\exp_{X_{n}}\left(\gamma_{n}\tilde{v}_{n}\right)=X_{n+1}$, i.e.,

$\gamma_{n}\tilde{v}_{n}=\exp^{-1}_{X_{n}}\Big{(}{\mathcal{R}}_{X_{n}}(\gamma_{n}V(X_{n};\omega_{n}))\Big{)}.$
(C.8)

Then (C.7) is an RRM scheme with $W_{n}=\tilde{v}_{n}-v(X_{n})$ where
$\tilde{v}_{n}$ is defined in (C.8). We will show that, under the assumption
$\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]<\infty$,
the following holds with probability 1:

$\displaystyle
b_{n}=\operatorname{\mathbb{E}}\left[W_{n}\nonscript\,\middle|\nonscript\,\mathopen{}\mathcal{F}_{n}\right]\to
0,\quad\sup_{n}\operatorname{\mathbb{E}}\left[\lVert
W_{n}\rVert^{2}_{X_{n}}\right]<\infty$ (C.9)

which obviously implies (7).

Consider the curve $c(t)\coloneqq{\mathcal{R}}_{X_{n}}(tV(X_{n};\omega_{n}))$.
By Lemma C.1, for $n$ large enough, $c(t)$ lies in the injectivity radius of
$X_{n}$ almost surely for all $t\in[0,\gamma_{n}]$. Let $\hat{c}(t)$ be the
smooth curve of $c(t)$ in the normal coordinate with base $X_{n}$ and an
arbitrary orthonormal frame, and let $\hat{X}_{n+1}$ be the normal coordinate
of $X_{n+1}$. Also, let $\tilde{v}_{n}^{\mathrm{N}}$ be the (Euclidean) vector
of $\tilde{v}_{n}$ expanded in the chosen orthonormal basis, and define
$V^{\mathrm{N}}(X_{n};\omega_{n})$ and
$\operatorname{Err}^{\mathrm{N}}(X_{n};\omega_{n})$ similarly. By definition,
$\hat{X}_{n+1}$ is nothing but $\gamma_{n}\tilde{v}_{n}^{\mathrm{N}}$.

Since $X_{n}=c(0)$ and $X_{n+1}=c(\gamma_{n})$, by properties of a retraction
map we must have

$\displaystyle\gamma_{n}\tilde{v}_{n}^{\mathrm{N}}$
$\displaystyle=\hat{c}(\gamma_{n})$
$\displaystyle=\hat{c}(0)+\gamma_{n}\dot{\hat{c}}(0)+\operatorname{\mathcal{O}}\left(\gamma_{n}^{2}\lVert\dot{\hat{c}}(0)\rVert_{2}^{2}\right)$
$\displaystyle=\gamma_{n}V^{\mathrm{N}}(X_{n};\omega_{n})+\operatorname{\mathcal{O}}\left(\gamma_{n}^{2}\lVert
V(X_{n};\omega_{n})\rVert_{X_{n}}^{2}\right)$
$\displaystyle\eqqcolon\gamma_{n}V^{\mathrm{N}}(X_{n};\omega_{n})+\gamma_{n}\tilde{b}_{n}$
(C.10)

where $\tilde{b}_{n}=\operatorname{\mathcal{O}}\left(\gamma_{n}\lVert
V(X_{n};\omega_{n})\rVert_{X_{n}}^{2}\right)$. Therefore, since
$\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]<\infty$
for all $x\in\mathcal{M}$, we have

$\operatorname{\mathbb{E}}\left[\left\lVert
W_{n}\right\rVert_{X_{n}}^{2}\right]=\operatorname{\mathbb{E}}\left[\left\lVert\operatorname{Err}^{\mathrm{N}}(X_{n};\omega_{n})+\tilde{b}_{n}\right\rVert_{2}^{2}\right]<\infty.$
(C.11)

On the other hand,

$\left\lVert
b_{n}\right\rVert_{X_{n}}=\left\lVert\operatorname{\mathbb{E}}\left[W_{n}\nonscript\,\middle|\nonscript\,\mathopen{}\mathcal{F}_{n}\right]\right\rVert_{X_{n}}=\left\lVert\operatorname{\mathbb{E}}\left[\tilde{b}_{n}\nonscript\,\middle|\nonscript\,\mathopen{}\mathcal{F}_{n}\right]\right\rVert_{2}=\operatorname{\mathcal{O}}\left(\gamma_{n}\lVert
V(X_{n};\omega_{n})\rVert_{X_{n}}^{2}\right).$ (C.12)

By Chebyshev’s inequality and the fact that
$\operatorname{\mathbb{E}}[\lVert\operatorname{Err}(x;\omega)\rVert^{4}_{x}]\leq\kappa^{2}<\infty$,
we have

$\operatorname{\mathbb{P}}\left({\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert^{2}_{X_{n}}\geq\sqrt{n\log^{1+\frac{\epsilon}{2}}n}}\right)\leq\frac{\kappa^{2}}{n\log^{1+\frac{\epsilon}{2}}n}$

where $\varepsilon$ is the same as in our choice of step-size in Proposition
2. Using an calculation identical to Lemma C.1, we conclude that

$\gamma_{n}\lVert\operatorname{Err}(X_{n};\omega_{n})\rVert^{2}_{X_{n}}=\operatorname{\mathcal{O}}\left(\frac{\sqrt{n\log^{1+\frac{\epsilon}{2}}n}}{\sqrt{n\log^{1+\varepsilon}n}}\right)=\operatorname{\mathcal{O}}\left(\frac{1}{\log^{\frac{\varepsilon}{4}}n}\right)$

which concludes the proof of (C.9).

For the retraction-based variants of Algorithms 2–4, by the above analysis, we
may replace ${\mathcal{R}}_{x}(\gamma_{n}V(\cdot;\omega_{n}))$ with
$\exp_{x}\left(\gamma_{n}\left(V(\cdot;\omega_{n})+\tilde{b}_{n}\right)\right)$
where $\tilde{b}_{n}\to 0$ almost surely. The rest is the same as in the proof
of Proposition 2. $\blacksquare$

###  C.3. Proof of Proposition 4

Similar to Proposition 2, Lemma C.1 guarantees that all geodesics are
minimizing and invertible. Hence, by the $L$-Lipschitzness of $v$ and (RRM-
alt), we have, with probability 1,

$\displaystyle\left\lVert
v_{z}(Y_{n+1},Z_{n})-v_{z}(Y_{n},Z_{n})\right\rVert_{z_{n}}$
$\displaystyle\leq L\operatorname{dist}\left(\begin{bmatrix}Y_{n+1}\\\
Z_{n}\end{bmatrix},\begin{bmatrix}Y_{n}\\\ Z_{n}\end{bmatrix}\right)$
$\displaystyle=L\gamma_{n}\left\lVert
v_{y}({Y_{n}},Z_{n})+U_{y,n}+b_{y,n}\right\rVert_{Y_{n}}$ $\displaystyle\to 0$
(C.13)

by Lemma C.1 and Proposition 2. Therefore, we may rewrite (RRM-alt) as

$\displaystyle
X_{n+1}=\exp_{X_{n}}\left(\gamma_{n}\begin{bmatrix}v_{y}(Y_{n},Z_{n})+U_{y,n}+b_{y,n}\\\
v_{z}({Y_{n}},Z_{n})+U_{z,n}+b^{\prime}_{z,n}\end{bmatrix}\right)$ (C.14)

where $b^{\prime}_{z,n}=b_{z,n}+v_{z}(Y_{n+1},Z_{n})-v_{z}(Y_{n},Z_{n})$. By
(C.13) and Proposition 2, $b^{\prime}_{z,n}$ satisfies (7), so our proof is
complete. $\blacksquare$

## References

  * Absil et al. [2008] Absil PA, Mahony R, Sepulchre R (2008) Optimization Algorithms on Matrix Manifolds. Princeton University Press
  * Ambrosio et al. [2005] Ambrosio L, Gigli N, Savaré G (2005) Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media
  * Bauschke and Combettes [2017] Bauschke HH, Combettes PL (2017) Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd edn. Springer, New York, NY, USA
  * Benaïm [1999] Benaïm M (1999) Dynamics of stochastic approximation algorithms. In: Azéma J, Émery M, Ledoux M, Yor M (eds) Séminaire de Probabilités XXXIII, Lecture Notes in Mathematics, vol 1709, Springer Berlin Heidelberg, pp 1–68
  * Benaïm and Hirsch [1995] Benaïm M, Hirsch MW (1995) Dynamics of Morse-Smale urn processes. Ergodic Theory and Dynamical Systems 15(6):1005–1030
  * Benaïm and Hirsch [1996] Benaïm M, Hirsch MW (1996) Asymptotic pseudotrajectories and chain recurrent flows, with applications. Journal of Dynamics and Differential Equations 8(1):141–176
  * Bento et al. [2017] Bento GC, Ferreira OP, Melo JG (2017) Iteration-complexity of gradient, subgradient and proximal point methods on Riemannnian manifolds. Journal of Optimization Theory and Applications 173(2):548–562
  * Benveniste et al. [1990] Benveniste A, Métivier M, Priouret P (1990) Adaptive Algorithms and Stochastic Approximations. Springer
  * Bonnabel [2013] Bonnabel S (2013) Stochastic gradient descent on Riemannian manifolds. IEEE Trans Autom Control 58(9):2217–2229
  * Boumal [2022] Boumal N (2022) An introduction to optimization on smooth manifolds. https://www.nicolasboumal.net/book/
  * Boumal et al. [2019] Boumal N, Absil PA, Cartis C (2019) Global rates of convergence for nonconvex optimization on manifolds. IMA Journal of Numerical Analysis 39(1):1–33
  * Brandière and Duflo [1996] Brandière O, Duflo M (1996) Les algorithmes stochastiques contournent-ils les pièges ? Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 32(3):395–427
  * Bridson and Haefliger [1999] Bridson MR, Haefliger A (1999) Metric Spaces of Non-Positive Curvature. Springer
  * Chen et al. [2021] Chen J, Liu S, Chang X (2021) Modified tseng’s extragradient methods for variational inequality on Hadamard manifolds. Applicable Analysis 100(12):2627–2640
  * Criscitiello and Boumal [2019] Criscitiello C, Boumal N (2019) Efficiently escaping saddle points on manifolds. Advances in Neural Information Processing Systems 32:5987–5997
  * Criscitiello and Boumal [2020] Criscitiello C, Boumal N (2020) An accelerated first-order method for non-convex optimization on manifolds. arXiv preprint arXiv:200802252
  * Drivas and Mailybaev [2021] Drivas TD, Mailybaev AA (2021) ‘life after death’in ordinary differential equations with a non-lipschitz singularity. Nonlinearity 34(4):2296
  * Durmus et al. [2020] Durmus A, Jiménez P, Moulines É, Said S, Wai HT (2020) Convergence analysis of Riemannnian stochastic approximation schemes. arXiv preprint arXiv:200513284
  * Durmus et al. [2021] Durmus A, Jiménez P, Moulines É, Said S (2021) On Riemannnian stochastic approximation schemes with fixed step-size. In: International Conference on Artificial Intelligence and Statistics, PMLR, pp 1018–1026
  * Facchinei and Pang [2003] Facchinei F, Pang JS (2003) Finite-Dimensional Variational Inequalities and Complementarity Problems. Springer Series in Operations Research, Springer
  * Fan et al. [2020] Fan J, Qin X, Tan B (2020) Tseng’s extragradient algorithm for pseudomonotone variational inequalities on Hadamard manifolds. Applicable Analysis pp 1–14
  * Ferreira and Oliveira [2002] Ferreira O, Oliveira P (2002) Proximal point algorithm on Riemannnian manifolds. Optimization 51(2):257–270
  * Ferreira et al. [2005] Ferreira OP, Pérez LL, Németh SZ (2005) Singularities of monotone vector fields and an extragradient-type algorithm. Journal of Global Optimization 31(1):133–151
  * Fujita and Kotani [1982] Fujita T, Kotani Si (1982) The onsager-machlup function for diffusion processes. Journal of Mathematics of Kyoto University 22(1):115–130
  * Hall and Heyde [1980] Hall P, Heyde CC (1980) Martingale Limit Theory and Its Application. Probability and Mathematical Statistics, Academic Press, New York
  * Hosseini and Sra [2020] Hosseini R, Sra S (2020) Recent advances in stochastic Riemannnian optimization. In: Handbook of Variational Methods for Nonlinear Geometric Data, Springer
  * Hsieh et al. [2021] Hsieh YP, Mertikopoulos P, Cevher V (2021) The limits of min-max optimization algorithms: Convergence to spurious non-critical sets. In: ICML ’21: Proceedings of the 38th International Conference on Machine Learning
  * Huang and Wei [2021] Huang W, Wei K (2021) Riemannnian proximal gradient methods. Mathematical Programming pp 1–43
  * Iliev [2006] Iliev BZ (2006) Handbook of normal frames and coordinates, vol 42. Springer Science & Business Media
  * Jost [2017] Jost J (2017) Riemannian geometry and geometric analysis, 7th edn. Springer, DOI https://doi.org/10.1007/978-3-319-61860-9
  * Kakade [2001] Kakade SM (2001) A natural policy gradient. Advances in neural information processing systems 14
  * Khammahawong et al. [2020] Khammahawong K, Kumam P, Chaipunya P, Yao JC, Wen CF, Jirakitpuwapat W (2020) An extragradient algorithm for strongly pseudomonotone equilibrium problems on Hadamard manifolds. Thai Journal of Mathematics 18(1):350–371
  * Kiefer and Wolfowitz [1952] Kiefer J, Wolfowitz J (1952) Stochastic estimation of the maximum of a regression function. The Annals of Mathematical Statistics 23(3):462–466
  * Klavins and Koditschek [2002] Klavins E, Koditschek DE (2002) Phase regulation of decentralized cyclic robotic systems. The International Journal of Robotics Research 21(3):257–275
  * Korpelevich [1976] Korpelevich GM (1976) The extragradient method for finding saddle points and other problems. Èkonom i Mat Metody 12:747–756
  * Kristály [2014] Kristály A (2014) Nash-type equilibria on riemannian manifolds: a variational approach. Journal de Mathématiques Pures et Appliquées 101(5):660–688
  * Kushner and Clark [1978] Kushner HJ, Clark DS (1978) Stochastic Approximation Methods for Constrained and Unconstrained Systems. Springer
  * Kushner and Huang [1981] Kushner HJ, Huang H (1981) Asymptotic properties of stochastic approximations with constant coefficients. SIAM Journal on Control and Optimization 19(1):87–105
  * Kushner and Yin [1997] Kushner HJ, Yin GG (1997) Stochastic approximation algorithms and applications. Springer-Verlag, New York, NY
  * Lee [1997] Lee JM (1997) Riemannian Manifolds: an Introduction to Curvature. No. 176 in Graduate Texts in Mathematics, Springer
  * Lezcano-Casado [2020] Lezcano-Casado M (2020) Curvature-dependant global convergence rates for optimization on manifolds of bounded geometry. arXiv preprint arXiv:200802517
  * Li et al. [2009] Li C, López G, Martín-Márquez V (2009) Monotone vector fields and the proximal point algorithm on Hadamard manifolds. Journal of the London Mathematical Society 79(3):663–683
  * Ljung [1977] Ljung L (1977) Analysis of recursive stochastic algorithms. IEEE Trans Autom Control 22(4):551–575
  * Ljung [1978] Ljung L (1978) Strong convergence of a stochastic approximation algorithm. Annals of Statistics 6(3):680–696
  * Manasse and Misner [1963] Manasse F, Misner CW (1963) Fermi normal coordinates and some basic concepts in differential geometry. Journal of mathematical physics 4(6):735–745
  * Massart et al. [2019] Massart E, Hendrickx JM, Absil PA (2019) Curvature of the manifold of fixed-rank positive-semidefinite matrices endowed with the bures–wasserstein metric. In: International Conference on Geometric Science of Information, Springer, pp 739–748
  * Mertikopoulos et al. [2022] Mertikopoulos P, Hsieh YP, Cevher V (2022) Learning in games from a stochastic approximation viewpoint. https://arxiv.org/abs/2206.03922
  * Muhammad and Egerstedt [2005] Muhammad A, Egerstedt M (2005) Decentralized coordination with local interactions: Some new directions. In: Cooperative Control, Springer, pp 153–170
  * Nesterov [2004] Nesterov Y (2004) Introductory Lectures on Convex Optimization: A Basic Course. No. 87 in Applied Optimization, Kluwer Academic Publishers
  * Neto et al. [2016] Neto JC, Santos P, Soares P (2016) An extragradient method for equilibrium problems on Hadamard manifolds. Optimization Letters 10(6):1327–1336
  * Pemantle [1990] Pemantle R (1990) Nonconvergence to unstable points in urn models and stochastic aproximations. Annals of Probability 18(2):698–712
  * Phelps [1993] Phelps RR (1993) Convex Functions, Monotone Operators and Differentiability, 2nd edn. Lecture Notes in Mathematics, Springer-Verlag
  * Popov [1980] Popov LD (1980) A modification of the Arrow–Hurwicz method for search of saddle points. Mathematical Notes of the Academy of Sciences of the USSR 28(5):845–848
  * Ratliff et al. [2016] Ratliff LJ, Burden SA, Sastry SS (2016) On the characterization of local Nash equilibria in continuous games. IEEE Trans Autom Control 61(8):2301–2307
  * Robbins and Monro [1951] Robbins H, Monro S (1951) A stochastic approximation method. Annals of Mathematical Statistics 22:400–407
  * Shah [2021] Shah SM (2021) Stochastic approximation on riemannian manifolds. Applied Mathematics & Optimization 83(2):1123–1151
  * Sternberg [1999] Sternberg S (1999) Lectures on differential geometry, vol 316. American Mathematical Soc.
  * Sun et al. [2019] Sun Y, Flammarion N, Fazel M (2019) Escaping from saddle points on riemannian manifolds. Advances in Neural Information Processing Systems 32
  * Takahashi and Watanabe [1981] Takahashi Y, Watanabe S (1981) The probability functionals (onsager-machlup functions) of diffusion processes. In: Stochastic Integrals, Springer, pp 433–463
  * Tang and Huang [2012] Tang Gj, Huang Nj (2012) Korpelevich’s method for variational inequality problems on Hadamard manifolds. Journal of Global Optimization 54(3):493–509
  * Tripuraneni et al. [2018] Tripuraneni N, Flammarion N, Bach F, Jordan MI (2018) Averaging stochastic gradient descent on Riemannnian manifolds. In: Conference On Learning Theory, PMLR, pp 650–687
  * Wang et al. [2010] Wang J, López G, Martín-Márquez V, Li C (2010) Monotone and accretive vector fields on riemannian manifolds. Journal of optimization theory and applications 146(3):691–708
  * Wang et al. [2021] Wang X, Tu Z, Hong Y, Wu Y, Shi G (2021) No-regret online learning over Riemannnian manifolds. In: Thirty-Fifth Conference on Neural Information Processing Systems
  * Warner [1965] Warner FW (1965) The conjugate locus of a Riemannian manifold. American Journal of Mathematics 87(3):575–604
  * Zhang and Sra [2016] Zhang H, Sra S (2016) First-order methods for geodesically convex optimization. In: Conference on Learning Theory, PMLR, pp 1617–1638

  *[RRM]: Riemannian Robbins–Monro
  *[RM]: Robbins–Monro
  *[SA]: stochastic approximation
  *[APT]: asymptotic pseudotrajectory
  *[ICT]: internally chain-transitive
  *[ODE]: ordinary differential equation
  *[APTs]: asymptotic pseudotrajectory
  *[w.p.$1$]: with probability 1
  *[PFS]: parallel frame system
  *[SFO]: stochastic first-order oracle
  *[RSGM]: Riemannian stochastic gradient method
  *[RPPM]: Riemannian proximal point method
  *[RSEG]: Riemannian stochastic extra-gradient
  *[NPG]: natural policy gradient
  *[NSE]: Nash–Stampacchia equilibrium
  *[FCS]: Fermi coordinate system
  *[ROG]: Riemannian optimistic gradient