{'fraction_non_alphanumeric': 0.04929725746170448, 'fraction_numerical': 0.014186450492182976, 'mean_word_length': 4.2055076037813395, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tracking and regret bounds for online zeroth-order Euclidean and Riemannian optimisation

We study numerical optimisation algorithms that use zeroth-order information to minimise time-varying geodesically-convex cost functions on Riemannian manifolds. In the Euclidean setting, zeroth-order algorithms have received a lot of attention in both the time-varying and time-invariant cases. However, the extension to Riemannian manifolds is much less developed. We focus on Hadamard manifolds, which are a special class of Riemannian manifolds with global nonpositive curvature that offer convenient grounds for the generalisation of convexity notions. Specifically, we derive bounds on the expected instantaneous tracking error, and we provide algorithm parameter values that minimise the algorithm's performance. Our results illustrate how the manifold geometry in terms of the sectional curvature affects these bounds. Additionally, we provide dynamic regret bounds for this online optimisation setting. To the best of our knowledge, these are the first regret bounds even for the Euclidean version of the problem. Lastly, via numerical simulations, we demonstrate the applicability of our algorithm on an online Karcher mean problem.


Introduction
Time-varying optimisation problems are popular in the machine learning community under the framework known as online convex optimisation (OCO) [32]. OCO is a promising methodology for modelling sequential tasks, and the main goal is developing algorithms that can track trajectories of the optimisers of the time-varying optimisation problem (up to asymptotic error bounds). OCO can be regarded as an iterative game between a player and an adversary. At each iteration k ∈ N 0 , the player (algorithm) selects a decision x k from a convex set K, and the adversary reveals a convex function f k : K → R. The player subsequently suffers an instantaneous loss f k (x k ). Particularly, in this time-varying context, it is often of interest to solve the following sequence of optimisation problems, min x∈K⊂R n f k (x).
(1) set of M, and each f k : M → R, k ∈ N 0 , is geodesically L-smooth and geodesically strongly convex (see Definitions 2 and 5 further below). Note that (2) is the extension of the Euclidean setting (1) to the Riemannian case. It is not uncommon that, in the aforementioned applications, explicit expressions for the cost functions f k may not be available, or their gradients are too costly to compute. For instance, in online system identification [53], f k is a complicated function with respect to x that depends on the output of a dynamical system, and it is evaluated through a real-world experiment. Moreover, in robust space tracking [23], when the problem is large-scale, evaluating the full gradient at every iteration can be excessively costly since the cost function in such problems is often expressible as a sum of several component functions, each depending only on a subset of measurements. Therefore, we consider a gradient-free setting where only function evaluations can be obtained via an oracle. Each cost function is thus seen as a black-box with time-varying input-output map. The zerothorder setting has been vastly studied in the Euclidean literature since it can be found in many applications where derivatives are either unavailable, or too expensive to compute, see e.g. [41,42] for machine learning, [35] for online controller tuning, [17] for deep neural networks, and [20] for mobile fog computing.
Formally, we aim to generate solutions to (2) using random gradient-free iterates of the form where P X denotes the projection that maps a point x ∈ M to P X (x) ∈ X ⊂ M such that dist(x, P X (x)) < dist(x, y), for all y ∈ X \{P X (x)}, and P X (x) = x for x ∈ X . The positive constant α k denotes the step size, g η,k + is the oracle or gradient estimator, and Exp x k (·) denotes the exponential mapping which we formally define in Section 2 below. We consider an extension of the recently proposed zeroth-order oracle for optimisation over Riemannian manifolds in [39] to make it suitable for our time-varying setting (2). Particularly, we define the oracle as where k + := k + 1/2, k ∈ N 0 and η > 0 corresponds to the oracle's precision. Some important discussions are in place behind the subscript k + and also the way u is constructed. We note that the oracle is essentially a two-point estimate of the directional derivative, and the subscript k + captures the time-varying nature of the underlying process. That is, our framework permits the cost function to change between the two evaluations. For example, consider an estimation problem where there is a continuously varying dynamical system, so from the moment we sample the cost function at x to the moment we evaluate it at Exp x (ηu), the cost function would have already changed by the underlying process, and this change is represented by k + . For the forthcoming analysis, let us define F := {f k : M → R|k ∈ N 0 ∪ {j + 1/2|j ∈ N 0 }}. Now, with respect to how u is constructed, we follow the same approach as [39]. That is, we let u = P u 0 ∈ T x M, where u 0 is normally distributed as per u 0 ∼ N (0, I n ) ∈ R n , and P ∈ R n×n is the orthogonal projection matrix onto the tangent space T x M of M at x.
Remark 1 Note that we assume that the cost functions f k are defined on submanifolds embedded in the Euclidean space, which makes it efficient to sample u from the associated tangent space. As also discussed in [39], we emphasise that the above gradient estimation methodology could be generally applied to other manifolds; however, efficiently sample Gaussian random objects on the tangent space of general manifolds is not an easy task [34]. Future work will focus on relaxing the way we sample the random vector u by using results such as the ones in [51].

More related work
In the general field of optimisation over Riemannian manifolds, and for the time-invariant-or offline-setting, several results have been proposed in the literature, and we briefly review them below. For instance, [11] provided convergence rates for deterministic Riemannian gradient descent and smooth cost functions. Stochastic algorithms were also considered for smooth Riemannian optimisation in [9,36,64,66,69]. Particularly, [9] extended the classical stochastic gradient descent algorithms to the Riemannian case, and provided convergence results. The authors in [66] introduced the variance-reduced RSVRG method and considered Riemannian optimisation of finite sums of geodesically smooth functions. The work [36] proposed a Riemannian stochastic recursive gradient algorithm (RSRG) which provides notable computational advantages in comparison to RSVRG. The work [69] introduced the Riemannian SPIDER method for non-convex Riemannian optimisation as a simple and efficient extension of the Euclidean SPIDER counterpart. Lastly, [64] studied stochastic projection-free methods for constrained optimisation of smooth functions on Riemannian manifolds. The stochastic Riemannian Frank-Wolfe methods for nonconvex and geodesically convex problems are introduced. For non-smooth cost functions, Riemannian subgradient methods have been proposed in [40], manifold ADMM methods in [37], manifold proximal gradient (ManPG) methods in [19], manifold proximal point algorithms (ManPPA) in [18], and stochastic ManPG in [63].
None of the aforementioned works have considered the zeroth-order setting, in which the oracle makes available only cost function values as opposed to first-order or second-order information.
To the best of our knowledge, zeroth-order Riemannian optimisation for time-invariant cost functions have been considered in [16,25,28,39]. Particularly, [28] presented the extended Riemannian stochastic derivative-free optimisation (RSDFO) algorithm, and proved it converges in finitely many steps in compact connected Riemannian manifolds. The authors in [16] extended the derivative-free optimisation method by [49] to Riemannian manifolds, but did not provide any complexity or convergence results. Just recently, [39] provided the first complexity results for both deterministic and stochastic zeroth-order Riemannian optimisation. Their zeroth-order methods rely on an estimator of the Riemannian gradient based on a modification of the Gaussian smoothing technique from the seminal work by Nesterov [45]. In [39], the authors illustrated that the proposed zeroth-order method has a comparable performance to its first-order counterparts in many applications of interest such as matrix approximation, k-PCA, sparse PCA, and the Karcher mean problem. Lastly, Dreisigmeyer in [25] studies direct search methods over general reductive homogenous spaces using maps from the tangent space to the manifold. We emphasise that these direct search methods use multiple mesh points per iteration, which is different to this paper, since we only evaluate the cost function twice at each step and take one action only.
We note that the existing works on zeroth-order Riemannian optimisation listed above do not consider the online setting (2) in which the cost-function is allowed to be time-varying. This problem, however, has been widely studied for M ≡ R n , see e.g. [6,14,21,23,65], in which it is assumed that cost function evaluations are carried out simultaneously when computing two-point estimates of the gradient. Later on, [56] relaxed this assumption and allowed the cost function to change between function evaluations, which added an extra modelling layer that better respects the time-varying nature of the problem.

Contributions
Our contributions are threefold: • We extend the OCO framework presented in [56] from the Euclidean setting to the case where the cost function is defined on a Hadamard manifold. Our proposed algorithm uses an extension of the zeroth-order oracle recently presented in [39] that allows for the function to be time-varying.
• We provide asymptotic bounds on the expected instantaneous tracking error, which to the best of our knowledge, are the first error bounds for online zeroth-order optimisation on Riemannian manifolds available in the literature. Our results illustrate how the manifold geometry-in terms of the sectional curvature-influences the performance of the algorithm. In addition, we provide explicit choices for the algorithm parameters-step size and oracle's precision-that minimise the performance of the algorithm.
• Lastly, we provide dynamic regret bounds for our time-varying setting. These are the first available dynamic regret bounds for (2), which were not previously available even for the Euclidean case.

Preliminaries
In this section, we present a brief introduction on the basics of manifold optimization. For a more in-depth revision we refer the reader to [1]. A smooth manifold is a topological manifold M with a globally defined differentiable structure. At any point x on a smooth manifold, tangent vectors are defined as the tangents of parametrised curves passing through x. The tangent space T x M of a manifold M at x is defined as the set of all tangent vectors at the point x. Tangent vectors on manifolds generalise the notion of a directional derivative. Formally, we can define the tangent space T x M as follows The tangent bundle of a differentiable manifold M is the manifold T M that assembles all the tangent vectors in M. As a set, it is the disjoint union of all the tangent planes, i.e. T M := ⊔ x∈M T x M. The dimension of a manifold M, denoted as d, is the dimension of the Euclidean space that the manifold is locally homeomorphic to. In particular, the dimension of the tangent space is always equal to the dimension of the manifold. A Riemannian manifold is a couple (M, g), where M is a smooth manifold equipped with a smoothly varying inner product (Riemannian metric) on the tangent space at every point, i.e. g(·, ·) := ·, · x : T x M × T x M → R. Without loss of generality, when the Riemannian metric is clear from the context, we simply talk about "the Riemannian manifold M". Throughout, we assume the Levi-Civita connection is associated with (M, g).
Definition 1 (Riemannian gradient) Suppose f is a smooth function on the Riemannian manifold M. The Riemannian gradient gradf (x) is defined as the unique element of A Riemannian submanifold M of a Riemanninan manifold N is a submanifold of N equipped with the Riemannian metric inherited from N . Since in this paper we assume that M is a Riemannian submanifold embedded in R n , M is equipped with the Riemannian metric inherited from R n . We thus write the inner product on the tangent space ·, · x at every point x ∈ M as ·, · x = ·, · , with the right-hand side being the Euclidean inner-product. Consequently, the Riemannian gradient in Definition 1 becomes the projection of its Euclidean gradient onto the tangent space, that is gradf (x) = P TxM (∇f (x)), where ∇f (x) denotes the Euclidean gradient of f at x.
We now introduce a family of local parametrisations often called retractions. In a nutshell, retractions allow us to move on a manifold (i.e. move in the direction of a tangent vector) while staying on the manifold. Formally, we say that a retraction mapping R x is a smooth mapping from T x M to M such that R x (0) = x, and the differential at 0 is an identity mapping, i.e. d dt R x (tξ) t=0 = ξ, for all ξ ∈ T x M. We refer to the latter as the local rigidity condition. In other words, for every tangent vector ξ ∈ T x M, the curve γ ξ : t → R x (tξ) satisfiesγ ξ (0) = ξ. Moving along the curve γ ξ is thought of as moving in the direction ξ while constrained to the manifold M. In particular, the exponential mapping Exp x is a retraction that generates geodesics. A geodesic is a curve representing in some sense the shortest path between two points in a Riemannian manifold.
Throughout this paper, we assume that M is a Hadamard manifold, as introduced in Section 1, i.e. complete, simply connected, and has nonpositive sectional curvature everywhere. Complete refers to the domain of the exponential mapping being the whole tangent bundle T M, and simply connected means there are no circular paths that cannot be shrunk to a point. Hadamard manifolds have strong properties. For instance, there exists a unique geodesic between any two points on M, and the exponential map is globally invertible at any point, where · is the norm associated with the Riemannian metric, which in our setting corresponds to the Euclidean norm as discussed above. On a Hadamard manifold, the notion of parallel transport provides a way to transport a vector along a geodesic. It is defined as the operator Γ y Definition 2 (Geodesically L-smoothness) A differentiable function f ∈ F is said to be geodesically L-smooth if there exists L ≥ 0 such that the following inequality holds for all x, y ∈ M, where we recall that dist(x, y) denotes the geodesic distance between x and y, and Γ x y is the parallel transport from y to x.
It can be shown that if f is geodesically L-smooth, then for any x, y ∈ M we have [67,69] We now provide convexity definitions on Hadamard manifolds. We emphasize that these definitions would be more subtle for general Riemannian manifolds, which can have several geodesics between two points, see e.g. [60].

Definition 3 (Geodesically convex set)
A set X ⊂ M is said to be geodesically convex if for any x, y ∈ X , the unique shortest length geodesic connecting x to y lies entirely in X .
where γ is the unique geodesic connecting x to y.
Definition 5 (Geodesic strong convexity) A function f : X → R is said to be geodesically σ-strongly convex if there exists σ ∈ (0, L] such that for any x, y ∈ X , The random object u in (4) is computed as in [39], i.e. u = P u 0 , where u 0 ∼ N (0, I n ) and P is the orthogonal projection matrix onto the tangent space. Therefore, as in [39], we define where ν is the normalising constant, see Section 2.2 and Lemma 2.1 in [39] for more details. Sometimes we will omit the subscript u 0 for simplicity.

Online optimisation using zeroth-order Riemannian oracles
Before presenting our results we state the underlying assumptions.

Assumption 1
(a) Every f ∈ F in (2) is geodesically L-smooth and geodesically strongly convex as per Definitions 2 and 5, respectively.
The above assumptions are standard in the context of online convex optimisation, and have been vastly used in the literature when M ≡ R n , see e.g. [23,44,56,68]. Particularly, (a) generalises assumptions such as Lipschitz gradient smoothness and strong convexity (to ensure a unique minimiser) often adopted in convex optimisation over R n . Item (b) bounds the change in minimiser, also called path variation in the literature [65], and (c) bounds the variation between consecutive cost functions. These assumptions on bounded variations essentially avoid two consecutive problems being arbitrarily different. For most tracking applications these assumptions are naturally satisfied, see e.g. [22,23]. For instance, in the context of robust subspace tracking, they ensure that the underlying target subspace does not change abruptly [62]. Note that δ → 0 corresponds to the special case where the cost-function does not vary between evaluations and thus g η,k + (x, u) = g η,k (x, u) in (4). Lastly, (d) assumes a lower bound in the curvature of M, and (e) an upper bound in the diameter of X , which are common assumptions in Riemannian optimisation see e.g. [9,67], and [39].
In online convex optimisation, there are usually two measures of interest to assess performance of the algorithms: the tracking error [23,56] and regret [32]. Henceforth, we study both measures and thus provide performance guarantees for the proposed algorithm (3) by means of upper bounds on both the tracking error and regret. The tracking error essentially captures how well the algorithm (3) can follow (or track ) the optimums of the time-varying optimisation problem (2) as k grows. On the other hand, regret (or in this case dynamic regret [44]), captures the accumulation of losses over the time horizon.

Tracking error bounds
The main objective of this section is to provide performance guarantees for algorithm (3) by means of upper bounds on the tracking error. That is, characterise how well (3) can track the optimisers of the time-varying problem (2). To that end, we define the tracking error as and the estimation error as The objects e k andē k are also known as pre-update optimality gap and post-update optimality gap, respectively, see e.g. [23].
To state our main tracking result, we require some intermediate steps which we present in the following. First, a common technique for analysing optimisation algorithms is to write the estimation error in terms of the tracking error by using the law of cosines in the Euclidean space. Unfortunately, this equality does not exist for general non-linear spaces, and in fact, there are no corresponding analytical expressions. However, in [67], a trigonometric distance bound for Alexandrov spaces with curvature bounded below was proposed. Alexandrov spaces are length spaces with curvature bound and form a generalisation of Riemannian manifolds with sectional curvature bounded below. The result uses the properties of geodesic triangles, and it can be used as an analogue to the law of cosines given its fundamental nature. For our setting, this inequality is formalised in the lemma below.
Proof See Appendix B.
Note that for the Euclidean setting, we set the curvature κ = 0 in Lemma 1, which leads to ζ = 1.
The next intermediate step is to derive a relation between the conditional expectations of e k+1 andē k . By Assumption 1(b) and the triangle inequality for the geodesic distance dist(·, ·), we can write The last required intermediate step is the proposition below, where we present bounds related to the zeroth-order oracle that are essential to show our optimiser tracking result.
Proof See Appendix C. Proposition 1 is the extension of Proposition 2.1 by [39] to the time-varying case. As also noted in [39] for the offline case, we see that the oracle in (4) is a biased estimator of the gradient in the online case, and the difference between them can be bounded as in Proposition 1. Moreover, note that our bounds have extra terms (w.r.t. [39]) that depend on δ, which is the upper bound on the cost function variation by means of Assumption 1(c). We emphasise that Proposition 1 recovers the oracle bounds by [39] for δ → 0 (time-invariant case).
We are now in a position to state a general result that illustrates how the conditional expectation of the tracking error evolves in time for any given step size α k . This is the main tool required to show that algorithm (3) can track the optimisers of (2).
We can see that the conditional expectation of the tracking error depends on ψ(e k ) which in turn depends on the parameters of the problem such as the Lipschitz constant L, manifold curvature κ, step size α k , oracle's precision η, and manifold dimension d. For instance, we can see that ψ increases with L, however, its dependence on other parameters such as the step size and oracle's precision is not trivial. However, it turns out that if we pick a constant step size, we can obtain a simpler expression for the expected tracking error. Particularly, if we choose a constant step size in the interval (0, σ/(2L 2 (d + 4)ζ(κ, R)), we can show that the expected tracking error remains bounded for k → ∞, which is the main objective of this section. This is formalised in the below corollary.
Corollary 1 shows that the expected value of the tracking error converges to a ball of radius ∆, as long as we pick a constant step size α k = α in algorithm (3) in the interval α ∈ 0, σ/(2L 2 (d + 4)ζ(κ, R)) . That is, the algorithm will track the time-varying optimisers of (2) with an error of ∆ as k → ∞.
The upper bound ∆ is used as a performance metric of the algorithm. Therefore, by minimising ∆, we aim to improve the tracking performance of the algorithm.
Remark 2 It is worth noticing that ∆ depends on the manifold geometry through κ, as opposed to the Euclidean counterpart by [56]. Moreover, ∆ depends on the intrinsic dimension d of the manifold, and not on the Euclidean ambient space dimension n which could be considerably larger. This is due to the fact that we work directly on the manifold and perform appropriate extensions of notions such as geodesic convexity, exponential maps, etc. Therefore, it would be more costly to work in the larger ambient Euclidean space, and our Riemannian method should be used preferably, unless a specific structure in the larger space can be considerably exploited to simplify calculations. The reported dependence of ∆ on the manifold curvature is consistent with recent literature. For instance, in Riemannian SVRG algorithms, the convergence rate of the algorithm depends on the manifold curvature as observed by [66], see also [67] for similar conclusions on subgradient methods.

Remark 3 (Convergence of the algorithm)
We can provide a complexity bound for the algorithm in terms of the number of iterations that takes the expected value of the tracking error to converge to a ball of radius ∆ + ǫ, where ǫ is an arbitrary positive scalar. Particularly, from the proof of Corollary 1, we have that (see (20)), With this inequality, we can show that there exists K ∈ N such that for all k ≥ K, E [e k ] ≤ ∆ + ǫ. In fact, such K satisfies It is assumed in (11) that (1 − ρ)E [e 0 ] − 2V > 0, otherwise K = 1.
In brief, we have shown that algorithm (3) can track optimisers of the time-varying optimisation problem (2) up to an error of ∆. We emphasise that this holds for any choice of constant step size in the interval α ∈ 0, σ/(2L 2 (d + 4)ζ(κ, R)) , and any choice of oracle's precision η. Indeed, ∆ may be large depending on the underlying parameters. Therefore, it is crucial to find parameters that minimise ∆, which is what we do in the following. Particularly, we find specific expressions for the step size α and oracle's precision η such that the performance metric ∆ is minimised.
Proof See Appendix F. As a summary of the results of this section, we essentially stated that, if we use algorithm (3) with choices of step size α and oracle's precision η as per Theorem 2, then the algorithm will track the optimisers of the time-varying optimisation problem (2) with performance ∆, which is in fact optimal for this choice of algorithm parameters.

Regret bounds
To cope with changing environments, another popular measure of performance for OCO algorithms is the so-called dynamic regret [31,65], which compares the cumulative loss of the learner/player to a sequence of optimal solutions. Consequently, we consider the following regret definitions, as counterparts to our tracking and estimation errors respectively, Reg Track.
Reg Est.
The main goal is to achieve sublinear regret. We say that an algorithm performs well if its regret is sublinear as a function of T , since this implies that, on the average, the algorithm performs as well as a clairvoyant who selects the minimiser at each step. Before providing the dynamic regret bounds, we impose an extra assumption, which is standard in the literature for constrained OCO, see e.g. [29,33,44,68].
Assumption 2 For all f ∈ F and x ∈ X , ∃G > 0 such that gradf (x) ≤ G.
Note that this assumption is widely used in convex optimisation, and it holds here since X is bounded and the cost functions are geodesically smooth.
For the above definitions of regret, we present the following bounds.
Theorem 3 If the step-size α k and oracle's precision η k are chosen as whereĀ and T k = 2 m for k ∈ [2 m − 1, 2 m+1 − 2], m ∈ N 0 . Then, the tracking and estimation regrets satisfy, for any T ≥ 1, Reg Track.
Reg Est.
We have showed in (14) that, disregarding V T , the proposed zeroth-order algorithm (3) achieves sublinear regret, which highlights its performance. This error term V T relates to the change in minimisers (or path variation) given the time-varying nature of the problem at hand. Obviously, if V T is sub-linear in T , we can see that the presented regret bounds would be sub-linear in T as well. It is therefore required that path variations diminish with T to achieve sublinear regret. This holds if the target being tracked slows down over time or eventually stops [31].

Numerical example
To validate our results, we apply our zeroth-order algorithm to the problem of computing the Karcher mean of a collection of symmetric positive definite (SPD) matrices, also known as Riemannian centre of mass or Fréchet mean [7]. This problem appears in a number of applications such as medical imaging [27], image segmentation [50], signal estimation [38], and particle filtering [10]. Note that the Karcher mean is guaranteed to exist and be unique on a Hadamard manifold, see e.g. [5].
We consider a time-varying version of the Karcher mean problem which arises in online scenarios. For instance, we may want to find a central representative for a collection of online noisy measurements of a moving object. Formally, we consider that the measurements are N SPD matrices of dimension m × m that become available at each time k, which we denote by {A k,1 , . . . , A k,N }.
The manifold of SPD matrices is defined as M := {X ∈ R m×m : X = X ⊤ ≻ 0}. If we equip M with the Riemannian metric M, N X := Tr X −1 M X −1 N , M, N ∈ T X M, for every X ∈ M, then the SPD manifold is a Hadamard manifold [4]. The Riemannian distance is given by where · F corresponds to the Euclidean (or Frobenius) norm, and the exponential mapping is for every X ∈ M, where exp denotes the matrix exponential. The time-varying cost function is defined as The Karcher mean for each set of measurements {A k,i } N i=1 received at time k is the unique minimiser of f k (X), i.e. x ⋆ k := arg min X∈M f k (X), for all k ∈ N 0 . The cost function (15) is known to be geodesically strongly convex with σ = 1 and geodesically ζ-smooth (i.e. (5) holds with L = ζ), see e.g. [67] and [66]. We consider δ = 0.001, N = 10, two problem sizes m ∈ {3, 9}, and the manifold dimension is d = m(m + 1)/2. For this example, we estimated 2 V = 0.5 and ζ = 1.5. The step size and oracle's precision for each m ∈ {3, 9} are chosen as per Theorem 2, which givesᾱ 3 = 0.0074, α 9 = 0.0015,η 3 = 0.0089, andη 9 = 0.005. The matrices {A k,i } N i=1 were randomly generated using the Manopt toolbox in Matlab, see [12]. Figure 1 (left) depicts the average tracking error E [e k ] after implementing the zeroth-order algorithm (3) for 100 random runs and different values of problem size. We can see the error converges, and its asymptotic value is upper bounded by some ∆. That is, in expectation, the algorithm can track the optimisers up to an asymptotic error. The theoretical asymptotic bounds from Corollary 1 are ∆ 3 = 543.73 and ∆ 9 = 2666, for m ∈ {3, 9}. Comparing with the asymptotic values in Figure 1, we can conclude that these bounds can become conservative depending on the application. Obtaining less conservative bounds is an open question for future research. It can also be seen that the larger the problem size, the longer E [e k ] takes to converge.
Lastly, we have also included a comparison with the first-order version of (3), in which the gradient of the function is fully available at every iteration k. We note that this is for illustration purposes only, since the results of this paper assume only cost function evaluations are available. We note a difference in rate of convergence, as expected, but also in terms of the asymptotic ∆ they achieve. We can see that if we have access to the gradient, then ∆ is smaller. This can actually be interpreted from the theory developed above. Particularly, from Proposition 1, we know the oracle is a biased estimator of the gradient of the cost function, and this difference can be upper bounded. These error terms wound not appear if the gradient is available at every k, and thus ∆ in Corollary 1 would be larger for the zeroth-order case. This, in link with Remark 1, can lead to interesting future work related to the choice of oracles and the way u is sampled when computing them.

Remark 4
In the context of our example, it is worth mentioning that for applications such as diffusion tensor imaging, the Euclidean averaging of SPD matrices often leads to a "swelling effect", artificial extra diffusion introduced in computation, see e.g. [3]. Particularly, it means the determinant of the Euclidean mean can be strictly larger than the original determinants. In diffusion tensor imaging, diffusion tensors correspond to covariance matrices of the local Brownian motion of water molecules. Introducing more diffusion is physically unacceptable in this context. Therefore, Riemannian approaches such as the one presented in this paper would be preferred.

Conclusions
A gradient-free algorithm for the minimisation of time-varying cost functions on Hadamard manifolds was proposed. Bounds on the expectation of the tracking error and on dynamic regret were derived, and choices for algorithm parameters such that the asymptotic tracking error bound is minimised were provided. Finally, the theoretical results were validated via numerical experiments.
Future work includes the extension to a more general class of Riemannian manifolds by trying to relax the way we sample the random vector u in (4), and also the extension to the non-convex time-varying case. In addition, looking at the problem in a different set of coordinates at each step is also an interesting future direction. Lastly, exploring different gradient approximations in the oracle for this time-varying context is also of interest.

A Auxiliary lemmas
Lemma 2 ( [4]) Let M be a Hadamard manifold and X ⊂ M a closed convex set. Then, the mapping P X (x) := {y ∈ X : dist(x, z) = inf z∈X dist(x, z)} is single-valued and nonexpansive, that is, we have dist(P X (x), P X (y)) ≤ dist(x, y) for every x, y ∈ M.

Lemma 4 ( [39]
) Suppose X is a d-dimensional subspace of R n , with orthogonal projection matrix P ∈ R n×n , u 0 ∼ N (0, I n ), and u = P u 0 is the orthogonal projection of u 0 onto X . Then,

Lemma 5 ( [67])
If a, b, c are the sides (i.e. lengths) of a geodesic triangle in an Alexandrov space with curvature lower bounded by κ, and A is the angle between sides b and c, then B Proof of Lemma 1.

C Proof of Proposition 1.
(a) We complete the proof with the following steps. Essentially we want to quantify how well the expectation of the oracle approximates the real gradient.
To compute a bound on E g η,k + (x, u) 2 , we proceed by definition of g η,k + (x, u), and thus first we bound f k + (Exp x (ηu)) − f k (x). We start by adding a convenient zero, that is, we add and subtract both f k + (x) and gradf k + (x), ηu , and then we use Assumptions 1(a) and 1(c).
That is, Consequently, The proof is thus complete by means of Lemma 4(b),(c).

D Proof of Theorem 1.
In this proof our main goal is to obtain (9). To achieve this, we note that (8) gives us a relation between E [ e k+1 | x k ] and E [ē k | x k ]. Therefore, we first compute the following. We take conditional expectations in (7) and obtain where the last inequality follows from Proposition 1(b). Below, we focus on computing the term (17). We add a convenient zero, in this case we add and subtract gradf k + (x k ) Now, note that since f is geodesically strongly convex as per Assumption 1(a), it is immediate to show that each f satisfies for all x ∈ M, where x ⋆ := arg min x∈M f (x). Inequality (18) is the Riemannian counterpart of the restricted secant inequality in R n . Then, by the Cauchy-Scharwz inequality and (18), where the last inequality follows from Proposition 1(a).
There is now only one term in (17) that remains to be bounded, which is gradf k + (x k ) 2 . To that end, note that from (5) and the reverse triangle inequality, gradf (x) − Γ x y gradf (y) ≤ Ldist(x, y), which, in turn, implies that gradf k + (x k ) ≤ Le k . Therefore, by using the latter together with (19) into (17), we obtain By Jensen's inequality we get E [ē k | x k ] 2 ≤ E ē 2 k x k , and the proof is thus complete from applying (8).
E Proof of Corollary 1.
The first part of the proof consists in simplifying the expression for ψ(e k ) that comes from Theorem 1 given the choice of constant step-size and also the bound on ζ. Note that Assumption 1(e) implies that ζ(κ, e k ) ≤ ζ(κ, R). Then, for α k = α, The second part of the proof boils down to using the above simplified expression for ψ(e k ) to obtain a recursive expression for E [e k ] so that we can iterate it and compute (10). Recall from the proof of Theorem 1 that E [ē k | x k ] 2 ≤ ψ(e k ) = ae 2 k + be k + c. Therefore, Lemma 3 implies E [ē k | x k ] ≤ ρe k + D. We then use the latter inequality in (8) to get E [ e k+1 | x k ] ≤ ρe k + D + 2V . Applying expectation then leads to the following recursive equation for E [e k ], Since ρ < 1 for 0 < α < σ 2L 2 (d+4)ζ(κ,R) , we can iterate (20) and obtain (10) in the limit k → ∞, concluding the proof.
F Proof of Theorem 2.
Next, we obtainᾱ that minimises ∆. To that end, we compute the derivative of ∆ evaluated at η with respect to α and set it to be equal to zero. Note that and then Lastly, we simply group terms in (22) to write Aα 2 + Bα + C = 0 with A, B and C as per the theorem statement, completing the proof.

G Proof of Theorem 3.
To prove Theorem 3 we need the following two intermediate lemmas.
Lemma 6 Let ρ k := 2(d + 4)L 2 ζ(κ, R)α 2 k − σα k + 1, and suppose α k ∈ (0, σ/(2L 2 (d+4)ζ(κ, R))) for all k ∈ N 0 . Then, Proof We know that e k+1 ≤ē k + dist(x ⋆ k + , x ⋆ k + +1 ), and thus Note that this is a relaxed version of (8) since we do not use the change on minimiser bound from Assumption 1(b). Therefore, from the proof of Corollary 1, it is easy to see that for any α k > 0, Summing both sides of (23) over and adding E [e 0 ] to both sides, Lemma 7 Let ρ k := 2(d + 4)L 2 ζ(κ, R)α 2 k − σα k + 1, and suppose α k ∈ (0, σ/(2L 2 (d+4)ζ(κ, R))) for all k ∈ N 0 . Then, Proof From the proof of Corollary 1, we can conclude that E [ē k ] ≤ ρ k E [e k ] + D k . Then, by the triangle inequality of the Riemannian distance, we can write We proceed similarly to the proof of Lemma 6, that is, completing the proof. Now we can proceed with the proof of Theorem 3. Note that the geodesic strong convexity of f k + , the Cauchy-Schwarz inequality, and Assumption 2 imply Similarly, it is immediate to show that Therefore, by means of Lemmas 6, 7, and the regret definitions in (12), we can immediately obtain the following upper bounds on Reg Track.
T and Reg Est.
T , Reg Track.
Reg Est.
The next step of the proof follows by using the so-called doubling-trick from [55, Section 2.3.1], which divides the algorithm rounds into periods of increasing size, specifically, in periods of 2 m rounds, m ∈ N 0 . Formally, this allows us to write, for any T ≥ 1, Then, if for each period of 2 m rounds we can find a step size and oracle's precision such that D k ≤c/ √ 2 m for somec > 0, k ∈ [2 m − 1, 2 m+1 − 2], m ∈ N 0 , then we would have which would prove the regret bound for Reg Track.
T in (14a). We can also bound T k=1 D k exactly as above to prove Reg Est.
T in (14b). Therefore, to conclude the proof, we have to show that the choices of α k and η k in the theorem statement indeed imply D k ≤c/ √ 2 m in each period of length 2 m , which is what we do below.
Since A ≥ 0 and C > 0, a necessary condition such that x 1 , x 2 > 0 is that B k < 0, which implies that α k < c 2 2Lδ(d+4) 2 ζ(κ,R)T k . We also need that ∆ k := B 2 k − 4AC ≥ 0, which imposes and extra condition on the step size α k , that is, .
The solutions toĀα 4 k +Bα 2 k +C = 0 are We analyse the determinant∆ k :=B 2 k − 4ĀC k , that is, which is always non-negative. Additionally, note thatĀ ≥ 0 for d < 4, andĀ ≤ 0 for d ≥ 5. We will use this fact to write a closed-form choice for the step size α k in each period T k of length 2 m .