{'fraction_non_alphanumeric': 0.04938153957574949, 'fraction_numerical': 0.01501839181230822, 'mean_word_length': 4.15017667844523, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Riemannian accelerated gradient methods via extrapolation

In this paper, we propose a simple acceleration scheme for Riemannian gradient methods by extrapolating iterates on manifolds. We show when the iterates are generated from Riemannian gradient descent method, the accelerated scheme achieves the optimal convergence rate asymptotically and is computationally more favorable than the recently proposed Riemannian Nesterov accelerated gradient methods. Our experiments verify the practical benefit of the novel acceleration strategy.


Introduction
In this paper, we consider the following optimization problem where M is a Riemannian manifold and f : M − → R is a smooth, real-valued function. Optimization on a Riemannian manifold naturally appears in various fields of applications, including principal component analysis [22,61], matrix completion and factorization [35, 56,13], dictionary learning [17,27], optimal transport [49,40,26], to name a few. Riemannian optimization [2,12] provides a universal and efficient framework for problem (1) that respects the intrinsic geometry of the constraint set. In addition, many non-convex problems turns out to be geodesic convex (a generalized notion of convexity) on the manifold, which yields better convergence guarantees for Riemannian optimization methods.
Nevertheless, it remains unclear whether there exists a simple strategy to accelerate firstorder algorithms on Riemannian manifolds. Existing research on accelerated gradient methods focus primarily on generalizing Nesterov acceleration [42] to Riemannian manifolds, including [37,4,63,6,31,36]. However, most of the algorithms are theoretic constructs and are usually less favourable in practice. This is due to the excessive use of exponential map, inverse exponential map, and parallel transport (formally defined in Section 3) each iteration. In addition, the Nesterov acceleration based methods require the knowledge of smoothness and strong convexity constants, which are often unknown. Furthermore, recent studies [23,19] show global acceleration cannot be achieved on manifolds in general, thus further questioning the utility of Nesterov type of acceleration for Riemannian optimization.
In this paper, we undertake a different approach for Riemannian acceleration compared to the aforementioned studies. The main idea is based on extrapolation, which aims to combine a (converging) sequence of iterates on manifolds to produce an accelerated sequence. We highlight that the iterate sequence can be given by any first-order solver, not necessarily the Riemannian gradient descent. In the Euclidean space, such an idea dates back to Aitken's ∆ 2 process [5] and Shanks sequence transformation [48,16], which are collectively referred to as nonlinear acceleration. Similar methods are also known as -algorithm [60], minimal polynomial extrapolation [50] and Anderson acceleration [58], just to name a few. The core idea of these methods is to compute extrapolation as a linear combination of the iterates where the weights depend nonlinearly on the iterates.
Particularly, we focus on the formulation in [46], which introduces regularization to nonlinear acceleration that shows better stability when optimizing non-quadratic functions. Given a sequence of iterates x 0 , ..., x k , [46] produces extrapolation k i=0 c i x i where the weights {c i } minimize the regularized combination of residuals, i.e., k i=0 c i (x i+1 − x i ) 2 2 + λ k i=0 c 2 i . It has been shown such a scheme achieves optimal convergence rate asymptotically without knowing the function-specific constants.
A natural question is whether such simple idea of extrapolation can be generalized to Riemannian optimization so that we achieve acceleration on manifolds both theoretically and practically. In this paper, we answer this question affirmatively. In particular, our contributions are as follows.
• We propose a simple acceleration strategy for Riemannian optimization based on the idea of extrapolation, which we call the Riemannian nonlinear acceleration (RiemNA) strategy. To extrapolate on manifolds, we develop several averaging schemes, which generalize weighted averaging in the Euclidean space from various perspectives.
• When the iterates are generated by the Riemannian gradient descent method, we show RiemNA leads to the same optimal convergence rate as the Riemannian Nesterov gradient methods [37, 63,4,36] asymptotically, while being more computationally efficient. We also show the convergence is robust to the choice of different averaging schemes on manifolds.
• Different to most analysis for Riemannian accelerated gradient methods [37, 63,4,36] that focus on exponential map and parallel transport, we show the use of general retraction and vector transport achieves the same convergence guarantees.
• Various experiments demonstrate the superiority of RiemNA compared to both Riemannian gradient descent and Riemannian Nesterov accelerated gradient methods.

Related work on Riemannian acceleration
The first paper that generalizes the Nesterov acceleration strategy [42] from the Euclidean space to the Riemannian manifold space is [37]. The algorithm however requires to compute the exact solution of a nonlinear equation each iteration and is usually infeasible in practice. In [63], a tractable algorithm for Riemannian Nesterov acceleration, called RAGD, is proposed and is proved to achieve local acceleration around the optimal solution. Furthermore, [4] shows global convergence and eventual acceleration of RAGD with carefully chosen parameters. Recently, [36] develops practical algorithms that achieve acceleration for both geodesic convex and strongly convex functions by assuming a bounded domain and curvature. Finally, [31] proposes a general framework for analyzing acceleration on manifolds where RAGD is a special case. Another direction of acceleration research on manifolds is inspired by the continuous dynamics formulation of Nesterov acceleration in the Euclidean space [52,59]. In [6], an ordinary differential equation (ODE) is introduced to model Nesterov acceleration on manifolds, and later [20] provides a variational framework for understanding Riemannian acceleration. Although the ODEs achieve acceleration, discretization of the continuous processes does not necessarily preserves the same accelerated rate of convergence. Lastly, there also exist acceleration results for specific manifolds, including sphere and hyperbolic manifolds [39], Stiefel manifold [51] and for nonconvex optimization [18]. For more related works on Riemannian acceleration, we refer to [19, Appendix A].

Preliminaries
Basic Riemannian geometry. A Riemannian manifold M is a smooth manifold endowed with a smooth inner product structure, namely a Riemannian metric on the tangent space T x M, for all x ∈ M. The Riemannian inner product between any u, v ∈ T x M is written as u, v x and the induced norm of a tangent vector u is u x = u, u x . A 'straight' line on manifold is called a geodesic γ : [0, 1] − → M, which is a locally distance minimizing curve with zero acceleration. Riemannian distance between x, y ∈ M is d(x, y) = inf γ 1 0 γ (t) γ(t) dt where γ(0) = x, γ(1) = y. Exponential map Exp x : T x M − → M maps a tangent vector u ∈ T x M to γ(1) with γ(0) = x, γ (0) = u. If between x, y ∈ M, there exists a unique geodesic connecting them, the exponential map has an inverse Exp −1 x (y) and the distance can be computed as d(x, y) = Exp −1 x (y) x = Exp −1 y (x) y . In this work, we only consider a unique-geodesic subset X of the manifold, which we explicitly assume in Section 5. Parallel transport Γ y x : T x M − → T y M allows tangent vectors to be transported along a geodesic that connects x to y such that the induced vector fields in parallel. It is known that parallel transport is isometric, i.e., Γ y x u, Γ y x v y = u, v x for any u, v ∈ T x M. In this paper, we also consider the more general retraction and vector transport that include exponential map and parallel transport as special cases. A retraction, Retr x : T x M − → M is a first-order approximation to the exponential map and a vector transport T y x : T x M − → T y M is a linear map between tangent spaces that approximates the parallel transport.
Function classes on Riemannian manifolds. For a differentiable, real-valued function f : represents the directional derivative of f along u and ∇f (x) is the Euclidean gradient and ·, · 2 represents the Euclidean inner product. The Riemannian Hessian at x, Hessf (x) : T x M − → T x M is a self-adjoint operator, defined as the covariant derivative of the Riemannian gradient. On Riemannian manifolds, one can extend the notion of gradient and Hessian Lipschitzness in the Euclidean space to geodesic Lipschitz gradient and Hessian associated with the exponential map and parallel transport. Some equivalent characterizations, including function smoothness and bounded Hessian norm also exist for the Riemannian counterparts. Furthermore, the notion of convexity can be similarly generalized to geodesic convexity where the convex combination is defined along the geodesics. Similar notions are also properly defined with respect to general retractions. The formal definitions are deferred to Appendix E. See also [12] for a more thorough treatment.
Metric distortion. Due to the curved geometry of Riemannian manifolds, many of the metric properties in the linear space are lost. To perform convergence analysis, we require the following geometric lemmas on manifolds that provide bounds on the metric distortion. Lemma 1 ([4,53]). Consider a compact subset X ⊆ M with unique geodesic. Let x, y = Exp x (u) ∈ X for some u ∈ T x M. Then for any , where X has curvature upper bounded by κ in magnitude and C κ (r) : . 4,32,38,53]). For a compact subset X ⊆ M with unique geodesic, there exists constants C 0 > 0, C 1 , C 2 ≥ 1 that depend on the curvature and diameter of X such that for all x, y, z ∈ X , u ∈ T x M we have

Riemannian nonlinear acceleration
We consider a generalization of nonlinear acceleration for Riemannian optimization via a weighted Riemannian average on the manifold. Specifically, for a set of weights {c i } k i=0 and a set of points {x i } k i=0 , we define the weighted Riemannian averagex c,x as computed according to the following recursion:x In the Euclidean space, which is a special case of Riemannian manifold, (Avg.1) recovers the weighted mean asx c,x = k i=0 c i x i (see Lemma 11 in Appendix B for more details). The coefficients {c i } k i=0 are determined by minimizing a weighted combination of the residuals Exp −1 x i (x i+1 ) ∈ T x i M, i = 0, ..., k. Given the residuals lie on different tangent spaces, we use parallel transport Γ x k x i to map them to a fixed tangent space, i.e., Following the work [46], we also consider a regularization in the coefficients, which leads to the following optimization problem: We show in Proposition 1 that the optimal coefficients c * has a closed-form solution.
The Riemannian nonlinear acceleration (RiemNA) strategy is presented in Algorithm 1, which takes a sequence of non-diverging iterates from any solver as input and constructs an extrapolation using weights that solve (2). The extrapolation is performed in parallel to the update of the iterate sequence. It is worth mentioning that when the manifold is the Euclidean space, Algorithm 1 exactly recovers the nonlinear acceleration algorithm in [46].
Computational comparison to Riemannian Nesterov acceleration. Based on Algorithm 1, each extrapolation step requires evaluation of 2k times the inverse exponential map, k times the parallel transport and exponential map as well as solving a linear system in dimension k + 1. Nevertheless, in practice, we often consider a limited-memory version of Algorithm 1 where only the most recent m iterates are used and the extrapolation only happens every m iterations (see more details in Section 8). In addition, if the iterates are given by Riemannian gradient descent, we only require evaluation of parallel transport, exponential map and inverse exponential map once per iteration on average. This is more efficient than even the most practical implementation of the Riemannian Nesterov accelerated gradient methods [63,36] (with algorithms given in Appendix D), which require at least twice the number of exponential and inverse exponential map per iteration. Such claims are also verified empirically in Section 8.

Convergence acceleration for Riemannian gradient descent
This section analyzes the convergence acceleration of RiemNA (Algorithm 1) for the iterates generated by the Riemannian gradient descent (RGD) method [2,12]. In particular, we show that the extrapolated point (output of Algorithm 1) eventually becomes a good estimate of the optimal point with more iterations. We start by making the following assumption throughout the paper. Assumption 1. Let x * ∈ M be a local minimizer of f . The iterates generated, i.e., x 0 , x 1 , ... stay within a neighbourhood X around x * with unique geodesic. Furthermore, the sequence of iterates is The former condition in Assumption 1 ensures the exponential map is invertible and is standard for analyzing accelerated gradient methods on manifolds [4,31,36]. This condition is satisfied for any non-positively curved manifolds, such as symmetric positive definite (SPD) manifold with the affine-invariant metric [9]. In addition, this also holds true for any sufficiently small subset X of any manifold.
Linear iterates and error decomposition in the Euclidean space. First, we recall that the convergence analysis for the Euclidean nonlinear acceleration [46] relies critically on a sequence of linear fixed-point iterates that satisfyx i − x * = G(x i−1 − x * ) for some positive semi-definite and contractive matrix G. The main idea is to show the algorithm converges optimally onx i and then bound the deviation arising from the nonlinearity. Particularly, let {x i } k i=0 be the given iterates and {x i } k i=0 be the linear iterates. Consider c * ,ĉ * as the coefficients solving (2) in the Euclidean setup respectively. The convergence analysis in [46] aims to bound each term from the following error decomposition: From linearized iterates to iterates on manifolds. On general Riemannian manifolds, due to the curved geometry, it becomes nontrivial to generalize such error decomposition. To mitigate this, we start with identification of linearized iterates on manifolds in the tangent space of x * . For notational convenience, we denote ∆ x := Exp −1 x * (x) for any x ∈ X . We now consider the linearized iteratesx i are produced by the following progression as for some G : T x * M − → T x * M as a self-adjoint, positive semi-definite operator with G x * ≤ σ < 1, where we denote A x * as the operator norm for any linear operator A on the tangent space T x * M. In fact, we show in the Lemma 3 that the progression of iterates from the Riemannian gradient descent method is locally linear on the tangent space of the local minimizer x * . This requires the following regularity assumption on the objective function f .
where id denotes the identity operator and ε i x * = O(d 2 (x i , x * )) and ε 0 = 0. Lemma 3 suggests that it is reasonable to consider the linearized iterates {x k } defined in (3) where G = id − η Hessf (x * ). It is clear that for a local minimizer x * , there exists µ, L > 0 such that µ id Hessf (x * ) L id. This is irrespective of whether the function f is geodesic strongly convex or has geodesic Lipschitz gradient. Thus, for proper choices of η, we can always ensure G is positive semi-definite and contractive.
In this paper, the convergence analysis focus on the case when G = id − η Hessf (x * ) and {x i } are given by Riemannian gradient descent to simplify the bounds. However, we highlight that most of the analysis holds for more general and symmetric G. Now, we approach the convergence ofx c * ,x computed from manifold weighted average recursion in (Avg.1) by considering the following error bound decomposition: where we denoteĉ * as the coefficients solving (2) with the residualsr i = ∆x i+1 − ∆x i from the linearized iterates {x i } in (3) andxĉ * ,x ,x c * ,x as weighted average computed using pairs respectively. Before we bound each of the error term, we first present a lemma relating the averaging on manifolds to averaging on the tangent space.
Remark 2. Lemma 4 shows that the error between the averaging on the manifold and averaging on the tangent space is on the order of O(d 3 (x 0 , x * )). This relies heavily on the metric distortion bound given in Lemma 1, 2, which only holds for the case of exponential map and parallel transport. Nevertheless, we highlight that when the general retraction and vector transport are used, we can follow the idea of [54,Lemma 12] to show the error is on the order of O(d 2 (x 0 , x * )). See Proposition 3 and Section 7 for more details where we discuss the convergence under the more general setup.
Error bound from the linear term. As in the Euclidean space, we see that the extrapolation using the linearized iterates converges in a near-optimal rate, via the regularized Chebyshev polynomial.
Definition 1 (Regularized Chebyshev polynomial [46]). The regularized Chebyshev polynomial of degree k, in the range of [0, σ] with a regularization parameter α, denoted as C : deg(p) = k, p(1) = 1} as the set of polynomials of degree k with coefficients summing to 1 and p 2 is the Euclidean norm of the coefficients of the polynomial p. We write the maximum valued as Next, we present the error bound coming from the linear term in Lemma 5, which follows from the definition of regularized Chebyshev polynomial and Lemma 4. We notice that due to the curvature of the manifold, we see an additional error term 1 compared to the Euclidean counterpart.
Lemma 5 (Error from the linear term). Under Assumption 1, letxĉ * ,x be computed from (Avg.1) Error bound from coefficient stability. Now we come to bound the deviation between the optimal coefficients computed via the Riemannian gradient descent iterates {x i } and the linearized iterates {x i }. To this end, we require the following results on the coefficients.
It should be noted that in the Euclidean space, ψ = k i=0 ∆ x i − ∆x i 2 = x i −x i 2 and also can be shown to have an order of O(d 2 (x 0 , x * )) under some Lipschitz conditions on the function (See [46,Proposition 3.8]). On manifolds, such error also suffers from additional distortion, which is on the order of O(d 2 (x 0 , x * )).
Next, based on Lemma 6, we show that the error from coefficient stability can be bounded as follows. The proof follows from linearizing the weighted average on the tangent space T x * M where we bound the deviation arising from the coefficients. Hence, an extra error 2 appears in the bound.
Lemma 7 (Error from coefficient estimation). Under the same settings as in Lemma 6, letxĉ * ,x , Error bound from nonlinearity. Next, we show that the nonlinearity term can be bounded in Lemma 8, which follows a similar idea of linearization on a fixed tangent space. Additional error 3 is again due to the curvature of the manifold, which vanishes in the case M is the Euclidean space.
Lemma 8 (Error from the nonlinearity). Under the same settings as in Lemma 6, we have Finally, we combine Lemma 5, 7 and 8 to obtain the following convergence result for Algorithm 1 applying on the iterates generated from the Riemannian gradient descent (RGD).
Theorem 1 (Convergence of RiemNA with RGD iterates). Under Assumptions 1, 2, let {x i } k i=0 be given by the Riemannian gradient descent method, i.e., We next show that even with additional distortion from the curved geometry of the manifold, the asymptotic optimal convergence is still guaranteed. This is mainly due to the fact that all errors incurred by the metric distortion, i.e., 1 , 2 , 3 are on the order of at least O(d 2 (x 0 , x * )), which is primarily attributed to Lemma 4.
Proposition 2 (Asymptotic optimal convergence rate of RiemNA with RGD iterates). Under the same settings as in Theorem 1, Remark 3. The asymptotic optimal convergence rate holds as long as 1 , 2 , 3 are on the order of Suppose at the local minimizer, we have 0 ≺ µ id Hessf (x * ) L id, which is satisfied beyond geodesic µ-strongly convex and geodesic L-smooth function. Then, we see that by choosing

Alternative averaging schemes on manifolds
In this section, we propose alternative averaging schemes on manifolds used for extrapolation. For the iterates obtained from the Riemannian gradient descent method, we show the schemes ensure the same asymptotically optimal convergence rate obtained in Proposition 2.
The first scheme we consider is based on the following equality in the Euclidean space for the weighted mean, i.e., We define an alternative weighted averaging as Based on the earlier analysis, to show the convergence ofx c,x defined in (Avg.2), we only require to show that Lemma 4 holds for the new scheme, with an error of order at least O(d 2 (x 0 , x * )). We formalize this claim in the next lemma and show the error is on the order of O(d 3 (x 0 , x * )).
In addition, we may consider the weighted Fréchet mean, which can also be used in place of the two aforementioned averaging schemes. The Fréchet mean is defined as the solution to the following optimization problemx Nevertheless, for general manifolds, it is not guaranteed the existence and uniqueness of the solution.
In fact, one can ensure the uniqueness of the solution when the function 1 2 d 2 (x, x ) is geodesic τ -strongly convex in x. From [6, Lemma 2], we see that the geodesic strong convexity of problem (Avg.3) holds for sufficiently small X on any manifold as well as for any non-positively curved manifold. Specifically, when M is non-positively curved, we have τ = 1. While for other manifolds, let D be the diameter of X and κ + > 0 be the upper curvature bound. Then, geodesic strong convexity is satisfied with τ < 1 when D < π 2 √ κ + .
. Under the additional assumption of geodesic strong convexity, Lemma 10 shows an extra tighter bound on d( Thus, we see the error from the linear term does not suffer from metric distortion ( 1 = 0). The error bound from coefficient stability and nonlinearity terms however, still incur additional errors as the previous two averaging schemes.
Lemma 9 and 10 allows convergence under the two averaging schemes to be established by exactly following the same steps as before. This is sufficient to show the same convergence bound holds (i.e., Theorem 1 and Proposition 2).

Convergence under general retraction and vector transport
It is noticed that nearly all studies on Riemannian accelerated gradient methods analyze convergence only under the exponential map and parallel transport. In this section, for the proposed Riemannian nonlinear acceleration on iterates from Riemannian gradient descent, we highlight that similar analysis and the same convergence bounds can be derived for general retraction and vector transport under additional assumptions. The main assumptions include bounding the deviation between retraction and exponential map as well as between vector transport and parallel transport. In addition, the Lipschitz gradient and Hessian need to be compatible with retraction and vector transport (defined in Appendix E.2). Assumption 3. The neighbourhood X is totally retractive where retraction has a smooth inverse. Function f has retraction Lipschitz gradient and Lipschitz Hessian. Assumption 4. There exists constants a 0 , a 1 , a 2 , δ a 0 ,a 1 > 0 such that for all x, y, z ∈ X , Retr −1 Assumption 5. The vector transport T y x is isometric and there exists a constant a 3 > 0 such that for all x, y ∈ X , T y Remark 5. Assumption 4 is required to bound the deviation from the retraction to the exponential map, which can be considered natural given retraction approximates the exponential map to the first-order. In fact, Assumption 4 has been commonly used in [45,34,24] for analyzing Riemannian first-order algorithms using retraction and can be satisfied for a sufficiently small neighbourhood (see for example [44,29]). Similarly, Assumption 5 is used to bound the deviation between the vector transport to parallel transport, which is also standard in [30, 34, 24]. One can follow the procedures in [30] to construct isometric vector transport that satisfies such condition for common manifolds like SPD manifold [30], Stiefel and Grassmann manifold [28].
In this section, we only show convergence under the recursive weighted average computation for extrapolation, i.e.,x Nevertheless, similar analysis can be also performed on the alternative two averaging schemes discussed in Section 6. The next theorem shows that using retraction and vector transport can also achieve the asymptotic optimal convergence rate. This proof follows similarly as the case for exponential map and parallel transport while using the Assumptions 4, 5. In particular, both these two assumptions ensure the deviations between retraction and exponential map, vector transport and parallel transport are on the order of O(d 2 (x 0 , x * )). Thus, the additional error terms 1 Theorem 2 (Convergence of RiemNA with RGD iterates under general retraction and vector transport). Under Assumption 1, 3, 4 and 5, let {x i } k i=0 be given by Riemannian gradient descent via retraction, i.e., Then, using retraction and vector transport in Algorithm 1 and lettingx c,x be computed from (4), it satisfies that 3 ), the same asymptotic optimal convergence rate (Proposition 2) holds.

Experiments
In this section, we test the Riemannian nonlinear acceleration (RiemNA) strategy for several applications on various manifolds, including the sphere, Stiefel, Grassmann, and symmetric positive definite manifolds. To this end, we apply the proposed strategy RiemNA (Algorithm 1) for accelerating iterates from the Riemannian gradient descent method with fixed stepsize (RGD), which we call RGD+RiemNA. The concrete steps are shown in Algorithm 2. Specifically, we run the Riemannian gradient descent to produce x 0 , . . . , x m−1 , where m is the memory depth. Then, we producex c * ,x with these iterates by Algorithm 1. We only consider the recursive weighted average in (Avg.1) for the experiments. We then restart Riemannian gradient descent with x 0 =x c * ,x for the next epoch.
Such an implementation has also been considered in [46] for the Euclidean case. Based on Section 7, we can apply general retraction and vector transport to implement Algorithm 2. for i = 1, ..., m do 5: end for 8: 10:

11:
Restart with x 0 =x c * ,x . 12: end while We compare the proposed RGD+RiemNA (Algorithm 2) with the recent development of the Riemannian Nesterov accelerated gradient (RNAG) methods [36], which correspond to the current state-of-the-art. We also include RAGD (a variant of Nesterov acceleration on manifolds proposed in [63]) and RGD as baselines. In particular, we compare with RNAG-C [36] (designed for geodesic convex functions) and RNAG-SC [36] and RAGD [63] (designed for geodesic strongly convex functions) regardless of whether the objective is of the particular class. We include more details of the algorithms in Appendix D.
According to [36], RNAG-C, RNAG-SC, and RAGD require the knowledge of geodesic Lipschitz constant L. Further, RNAG-SC and RAGD require the geodesic strong convexity parameter µ. In particular, the stepsize of RNAG-C, RNAG-SC and RAGD should be set as 1/L. If such constants are available, we set them accordingly. Otherwise, we tune over the parameters L, µ for RNAG-C, RNAG-SC to obtain the best results and set the same parameters for RAGD for comparability. Following [36], the additional parameters ξ, ζ are fixed to be 1 for RNAG-C, RNAG-SC and β = µ/L/5 for RAGD. We set stepsize of RGD to be 1/L if available and tune the stepsize otherwise. For the proposed RiemNA, we fix λ = 10 −8 and choose memory depth over {5, 10}.
For fair comparisons, we use exponential map, inverse exponential map, and parallel transport for all the algorithms whenever such operations are properly defined. For other cases, we use retraction, inverse retraction, and vector transport. It should be noted that we maintain consistency of the use of such operations across all the algorithms.
All the algorithms are stopped once gradient norm reaches below 10 −6 . All experiments are coded in Matlab using the Manopt toolbox [15]. The codes can be found on https://github.com/andyjm3.
Leading eigenvector computation. First, we consider the problem of computing the leading eigenvector of a symmetric matrix A of size d × d, by solving For the experiment, we generate a positive definite matrix A with condition number 10 3 and exponentially decaying eigenvalues in dimension d = 10 3 . As shown in [36, Proposition 7.1], the problem has geodesic L-Lipschitz gradient with L to be the eigengap of matrix A, i.e., the difference between maximum and minimum eigenvalues of A.
The stepsize is thus set as 1/L for all methods. For RNAG-SC and RAGD, we set µ = 10. For  RiemNA, we set memory depth to be m = 10. We use exponential and inverse exponential map as well as projection-type vector transport for all algorithms including RGD+RiemNA. The operations are given in Appendix C.
In Figure 1a, we compare the algorithms in optimality gap measured as f (x t ) − f (x * ) against both iteration number and runtime. We notice that RiemNA performs competitively with Nesterov acceleration in terms of iterations and significantly outperforms all baselines in runtime. This is due to the light computational cost compared to Nesterov type of acceleration on manifolds (according to the discussion in Section 4). Finally we remark that in the initial phase, RiemNA does not necessarily ensure descent in the objective, only in the later phase where acceleration takes place. This is in accordance with our local convergence analysis.
In Figure 1b, we test the sensitivity of the choices of regularization parameter λ (on the left with m = 10 fixed) and memory depth m (on the right with λ = 10 −8 fixed). The results show robustness of RiemNA under various choices of regularization parameter λ and memory depth m.
Fréchet mean of SPD matrices. The second example is to compute the Fréchet mean of several symmetric positive definite (SPD) matrices under the affine-invariant metric [9]. Specifically, given a set of SPD matrices ++ as the set of SPD matrices of size d × d, · F as the Frobenius norm, and logm(·) is the matrix logarithm. The geometry of the SPD manifold is given in Appendix C. For the experiments, we use exponential map and its inverse as well as the parallel transport for all the algorithms. As commented previously, the geometry is negatively curved, and thus, the Fréchet mean problem is geodesic 1-strongly convex. Hence µ = 1. For this problem, we generate random N = 100 SPD matrices of dimension d = 10. The stepsize for all methods are tuned and set to be 0.5. For RiemNA, we set memory depth m = 5. In Figure 1c, we observe superiority of RiemNA in both iterations and runtime. We particularly notice in this case, the extrapolation step leads to a significant acceleration of convergence.
Orthogonal procrustes problem. We also consider the orthogonal procrustes problem on Stiefel manifold. Suppose we are given A ∈ R r×r , B ∈ R p×r , the objective is min X∈St(p,r) XA − B 2 F where St(p, r) := {X ∈ R p×r : X X = I} is the set of column orthonormal matrices, which forms the so-called Stiefel manifold with the canonical metric. To implement the algorithms, we use QR-based retraction and inverse retraction as well as projection-type vector transport. The definitions are given in Appendix C. We generate random matrices A, B where the entries are normal distributed. We set p = 100, r = 5. For this problem, both L and µ are unknown. Hence, we tune and set stepsize to be 1 for all methods. For RNAG-SC and RAGD, we select µ = 0.005 and for RiemNA we set memory depth m = 5. Figure 1d shows similar patterns as previous two experiments where RiemNA consistently converges faster than all the baselines.
Nonlinear eigenspace problem. Finally, the problem of computing nonlinear eigenspace arises as the total energy minimization on Grassmann manifold [64], i.e., min X∈Gr(p,r) 1 2 tr(X LX) + 1 4 ρ(X) L −1 ρ(X) where ρ(X) := diag(XX ) and L is a discrete Laplacian operator. The Grassmann manifold Gr(p, r) is the set of all r-dimensional subspaces in R p , which can be identified as an equivalence class of Stiefel manifold. For experiment, we implement the algorithms with QR-based retraction and inverse retraction as well as projection-based vector transport similar to Stiefel manifold (details given in Appendix C). We generate L as a tridiagonal matrix with main diagonal entries to be 2 and sub-and super-diagonal entries to be −1. The stepsize is tuned and set to be 0.1 for all methods and for RNAG-SC, RAGD, µ = 5 and for RiemNA, m = 5. We show the convergence results in Figure 1e where we also observe that RiemNA significantly outperforms other methods.

Concluding remarks and perspectives
In this paper, we introduce a simple strategy for accelerating first-order Riemannian optimization algorithms on manifolds, which is based on the idea of extrapolation and nonlinear acceleration. The extrapolation step is performed via novel intrinsic weighted averaging schemes on manifolds. By carefully bounding the metric distortion, we show that the Riemannian acceleration scheme achieves the same optimal convergence rate as the Euclidean counterparts. When applied to the iterates from the Riemannian gradient descent method, the new acceleration scheme asymptotically matches the convergence guarantee obtained by the Nesterov type of acceleration algorithms on manifolds, while being more practically favourable. We verify such claims for a variety of applications.
For simplicity, this work only analyze convergence of the acceleration with the Riemannian gradient descent method. This is because locally on the tangent space at optimality, Riemannian gradient descent leads to a linear progression given by a symmetric operator G as shown in Lemma 3. However, for general solvers, including conjugate gradient and even Nesterov accelerated gradient methods on manifolds, Lemma 3 may not hold. It would be interesting to explore nonlinear acceleration for more advanced solvers and also to consider stochastic settings.

A Extensions
In this section, we consider various extensions to the proposed nonlinear acceleration on manifolds.

A.1 Online Riemannian nonlinear acceleration
Following [47,10], we can extend Algorithm 1 to the online setting, where the extrapolated point x c,x is used to update the iterate sequence. The idea is to add a mixing step by updatingx c,x in the direction of the weighted average of the gradients, i.e., gradf (x c,x ) = k i=0 c i Γx c,x x i gradf (x i ). For the averaging schemes (Avg.1), (Avg.3), the next iteration starts with Expx c,x (−δ gradf (x c,x )) for some mixing parameter δ > 0. Particularly for the tangent space averaging scheme (Avg.2), we show a more efficient strategy of mixing, which we focus in this paper. The averaging and mixing steps are both performed on the same tangent space. Specifically, let x −1 = x 0 , we define the following progression of the online nonlinear acceleration on manifolds.
The complete procedures are presented in Algorithm 3. Compute Solve c * = arg min c∈R k+1 :c 1=1 k i=0 c i r i 2 x k + λ c 2 2 .

A.2 Practical considerations
Here are some practical considerations to use nonlinear acceleration on manifolds.
Iterates from Riemannian gradient descent with line-search. Suppose the iterates {x i } k i=0 are generated from x i = Exp x i (−η i gradf (x i−1 )) where the stepsize is determined from a line-search procedure (such as backtracking line-search [14]) and thus varies across iterations. Nevertheless, Lemma 3 still holds with G i = id − η i Hessf (x * ). Suppose the stepsize is chosen such that G i ≤ σ < 1. Then the analysis still holds under this setting.
Safeguarding decrease. Due to the curved geometry of the manifold and nonlinearity of the objective function, it is not guaranteed that f (x c,x ) will decrease. In the main text, we only show local convergence of the acceleration strategy. A typical globalization technique is to only keep the extrapolated point if it shows sufficient decrease compared to previous iterates, i.e., f (x c,x ) ≤ τ min i=0,...,k f (x i ) for some τ < 1. In [46], an adaptive regularization strategy has been proposed to select regularization parameter λ. Here we adapt the same strategy on manifolds, which we show in Algorithm 4. As noticed in [46], a higher value of λ pushes the weights close to uniform and thus stays closer to x 0 . Thus the line-search over t tries to enhance the progress compared to the initialization. In addition, for online Riemannian nonlinear acceleration specifically, we may consider performing a line-search over the parameter δ to ensure a sufficient descent condition is met.
Limited-memory and extrapolation frequency. Rather than keeping all the previous iterates for extrapolation, we can set a memory depth of m and using only the most recent m iterates to compute the extrapolated point. In practice, m is usually set to be less than 10. In addition, we notice that compared to the Euclidean version, the computational cost for the Riemannian nonlinear acceleration can be high due to the use of parallel transport. Hence to mitigate this issue, we may only compute the extrapolated point every m iteration.
Efficient update of the residual matrix R. Recall for each application of Riemannian nonlinear acceleration, we need to compute where we write using (isometric) vector transport and general retraction. This includes parallel transport and exponential map as special cases. By isometry, in the next iteration when we receive r k+1 , the update of R only requires computing Γ x k+1 x k r i , r k+1 x k+1 , i = 0, ..., k + 1. Denote the vector x k r i , r k+1 x k+1 ] 0≤i≤k . Then the updated residual matrix is . Algorithm 4 Adaptive regularized Riemannian nonlinear acceleration (AdaRiemNA) 1: Input: A sequence of iterates x 0 , ..., x k+1 . Tentative regularization parameters {λ j } k j=1 .

B From Euclidean averaging to Riemannian averaging
To extend the idea of weighted average to manifolds, we first rewrite the weighted average on the Euclidean space as follows.
Lemma 11 (Weighted average recursion). Given a set of coefficients {c i } k i=0 with k i=0 c i = 1 and a set of iterates {x i } k i=0 . Let the streaming weighted average be defined asx Proof. For some γ 1 , ..., γ k , the streaming weighted average is defined asx . We first show the streaming weighted average has the form We prove such argument by induction. For i = 1, it is clear thatx 1 = (1 − γ 1 )x 0 + γ 1 x 1 and satisfies the form. Suppose at i = k , the equality is satisfied, then for i = k + 1, we havẽ which satisfies the equality. Hence this argument holds for all i ∈ [k]. Finally, at i = k, we see that the choice that γ i = c i i j=0 c i leads to the matching coefficients.

C Geometry of specific Riemannian manifolds
Sphere manifold. The sphere manifold S d−1 is an embedded submanifold of R d with the tangent space identified as T x S d−1 = {u ∈ R d : x u = 0}. The Riemannian metric is given by u, v = u, v 2 for u, v ∈ T x S d−1 . We use the exponential map derived as Exp x (u) = cos( u 2 )x+sin( u 2 ) u u 2 and the inverse exponential map as Exp −1 x is the orthogonal projection of any v ∈ R d to the tangent space T x S d−1 . The vector transport is given by the projection operation, i.e., T y x u = Proj y (u).
Symmetric positive definite (SPD) manifold. The SPD manifold of dimension d is denoted as S d ++ := {X ∈ R d×d : X = X, X 0}. The tangent space T X M is the set of symmetric matrices. The affine-invariant Riemannian metric is given by U, V X = tr(X −1 U X −1 V ) for any U, V ∈ T X S d ++ . We make use of the exponential map, which is Exp X (U ) = Xexpm(X −1 U ) where expm(·) is the matrix exponential. The inverse exponential map is derived as We consider the parallel transport given by Stiefel manifold. The Stiefel manifold of dimension p × r is written as St(p, r) := {X ∈ R p×r : X X = I}. The Riemannian metric is the Euclidean inner product defined as U, V X = U, V 2 . We consider the QR-based retraction Retr X (U ) = qf(X + U ) where qf(·) returns the Q-factor from the QR decomposition. The inverse retraction is derived as for X, where R is solved from the system X Y R + R Y X = 2I. The vector transport is given by the orthogonal projection, which is Grassmann manifold. The Grassmann manifold of dimension p × r, denoted as Gr(p, r), is the set of all r dimensional subspaces in R p (p ≥ r). Each point on the Grassmann manifold can be identified as a column orthonormal matrices X ∈ R p×r , X X = I and two points X, Y ∈ Gr(p, r) are equivalent if X = Y O for some O ∈ O(r), the r × r orthogonal matrix. Hence Grassmann manifold is a quotient manifold of the Stiefel manifold. We consider the popular QR-based retraction, i.e. R X (U ) = qf(X + U ) where for simplicity, we let X to represent the equivalence class and U represents the horizontal lift of the tangent vector. The inverse retraction is also based on QR factorization, i.e. R −1

D Riemannian Nesterov accelerated gradient methods
Here we include the practical algorithms of Riemannian Nesterov accelerated gradient methods presented in [63,36]. It is worth noticing that all the algorithms are analyzed under the exponential map and inverse exponential map and parallel transport. In contrast, the proposed RiemNA applies to general retraction and vector transport. We first present the (constant-stepsize) RAGD method in [63,Algorithm 2], which is included in Algorithm 5. We see the algorithm requires three times evaluation of the exponential map and two times the inverse exponential map.

E Function classes on Riemannian manifolds
This section briefly reviews the various functions classes on Riemannian manifolds.

E.1 Geodesic gradient Lipschitzness and Hessian Lipschitzness
First we provide several equivalent characterizations for the gradient and Hessian Lipschitzness. For proof and more detailed discussions, see [12,Section 10.4].
y k (x k+1 ) 9: end for 10: Output: x K Lemma 12 (Geodesic gradient Lipschitzness and function smoothness). A function f : M − → R has geodesic L-Lipschitz gradient in X ⊆ if for all x, y = Exp x (u) ∈ X in the domain of the exponential map, we have Γ x γ(t) gradf (γ(t)) − gradf (x) x ≤ L tu x , for all t ∈ [0, 1] and γ(t) := Exp x (tu). This is equivalent to function having bounded Hessian as Hessf (x) x := max u∈TxM: u x=1 Hessf (x) [u] x ≤ L, where Hessf (x) x denotes the operator norm of Riemannian Hessian. If function f has geodesic L-Lipschitz gradient, then function f is geodesic L-smooth, which satisfies Lemma 13 (Geodesic Hessian Lipschitzness). A function f has geodesic ρ-Lipschitz Hessian in X ⊆ M if for all x, y = Exp x (u) ∈ X in the domain of the exponential map, we have for all t ∈ [0, 1] and γ(t) := Exp x (tu). If function f has geodesic ρ-Lipschitz Hessian, then function f satisfies

E.2 Retraction gradient Lipschitzness and Hessian Lipschitzness
In this section, we define the gradient and Hessian Lipschitzness with respect to a retraction, which generalizes the definitions in Section E.1.

Definition 3 (Retraction Hessian Lipschitzness)
. A function f : M − → R has retraction ρ-Lipschitz Hessian in X ⊆ M if for all x, y ∈ R x (u) ∈ X in the domain of the retraction, we have where we denote c(t) = R x (tu).

E.3 Geodesic convexity and strong convexity
We start with the notion of geodesic convex set. A subset X ⊆ M is called geodesic convex if for any two points in the set, there exists a geodesic joining them that lies entirely within the set.
A similar notion of convexity with respect to retraction also exists by replacing the geodesic curve γ(t) with retraction curve c(t) = Retr x (tu). See [30] for more details.

F Main proofs
Before we proceed with the proofs of the results in the paper, we introduce a lemma that will be used often in the course of the proof.
Lemma 14. Under Assumption 1, for any w, x, y, z ∈ X , we have Γ Proof of Lemma 14.
where we apply Lemma 1 and 2.

F.1 Proof of Proposition 1
Proof of Proposition 1. Let µ ∈ R be the dual variable. Then we have c * , µ * satisfy the KKT system: Solving the system yields the desired result.

F.2 Proof of Lemma 3
Proof of Lemma 3. First, we consider the pushforward operator Exp y x : T x M − → T y M for any x, y ∈ M, defined as Exp y x (u) := Exp −1 y (Exp x (u)) for any u ∈ T x M. The differential of Exp y x at 0 along u ∈ T x M is derived as where we denote T y The second equality is due to Exp x (0) = 0, DExp x (0) = id and the third equality follows from the inverse function theorem. Then by Taylor's theorem for Exp x * x i around 0, we have x i . Then by Hessian Lipschitzness (Lemma 13), we have around x * Combining (5) with (6) yields To show the desired result, it remains to show the operator (Γ x * x i T x i x * ) −1 is locally identity. This is verified in [54,Lemma 6] for general retraction. We restate here and adapt to the case of exponential map.
Consider the function H(u) := (Γ where the second equality is due to the property of parallel transport (see for example [12,Proposition 10.37]). In addition, from [57,Theorem A.2.9], we see the second order derivative of H is given by where we denote Riem x as the Riemann curvature tensor evaluated at x. We notice that Riem x (u, ·)u : T x M − → T x M is symmetric with respect to the Riemannian metric (see for example [7]). For ∈ T x M, we apply the Taylor's theorem for H up to second order, which yields . Substituting this result into (7), we obtain .
. We can bound the error term as follows.
where we use the bounds on e i x * , i x * as well as Hessf (x * )[∆ x i ]+e i = Γ x * x i gradf (x i ) and geodesic gradient Lipschitzness (Lemma 12) such that gradf (x i ) 2 ≤ L ∆ i

F.3 Proof of Lemma 4
Proof of Lemma 4. The proof is by induction. Let γ i = c i i j=0 c j and first we rewrite the averaging on tangent space as following the recursion defined as ∆ To show the difference between ∆x c,x , based on Lemma 2, it suffices to show the distance betweenx k =x c,x and Exp x * ( ∆ x k ) is bounded.
To this end, we first notice thatx 0 = x 0 = Exp x * ( ∆ x 0 ) and we consider bounding the difference betweenx 1 and Exp x * ( ∆ x 1 ). To derive the bound, we first observe that by Lemma 1, where we see where the last inequality is from the proof of Lemma 14. Thus combining (9), (10) leads to .
The first term on the right of (11) can be bounded as where we again use the result from the proof of Lemma 14. To see (12) is on the order of O(d 3 (x 0 , x * )), we only need to show ∆x i 2 , which can be seen by a simple induction argument. First, it is clear that ∆x 0 2 x * = d 2 (x 0 , x * ). Then suppose for any i < k, we have d(x i , x * ) = O(d(x 0 , x * )). Then from Lemma 2, we have Thus, using d( (12) is on the order of O(d 3 (x 0 , x * )). Now we bound the second term on the right of (11). Particularly, where we apply Lemma 2 multiple times. Combining the previous results, we see Now applying Lemma 2, we obtain . Thus the proof is complete.

F.4 Proof of Lemma 5
Proof of Lemma 5. Directly combining Lemma 15 and Lemma 4 gives the result.

Lemma 15 (Convergence of the linearized iterates). Consider the linearized iterates {x
Proof of Lemma 15.
x * where p ∈ P 1 k and {c i } k i=0 are the corresponding coefficients. Then we obtain whereλ = λ/d 2 (x 0 , x * ) and we use the fact that G − id x * ≤ 1. Then where we see that

F.5 Proof of Lemma 6
Proof of Lemma 6. From Proposition 1 and following [46, Proposition 3.2], we obtain Now we bound R 2 . First we see R can be rewritten as R G x k R, where G x k ∈ R r×r is the positive definite metric tensor at x k and R = [ r i ] ∈ R r×k is the collection of tangent vector in an orthonormal basis and r is the intrinsic dimension of the manifold. Thus we can write Riemannian inner product as r i , r j x k = r i G x k r j and On the other hand, denote the perturbation matrix P = R −R. Then from Proposition 1 and following [46, Proposition 3.2], we have where we use Lemma 2. Now we respectively bound each of the two terms on the right. First we see from Lemma 14, Combining (15), (14), (13) leads to Finally, recall we can write R = R G x k R and similarly forR =R G x * R whereR = [ r i ]. By isometry of parallel transport, we have Then Notice that Also

F.6 Proof of Lemma 7
Proof of Lemma 7. From Lemma 2, we first observe that d(xĉ * ,x ,x c * ,x ) ≤ C 1 ∆xĉ * ,x − ∆x c * ,x x * . Now we derive a bound on the term ∆xĉ * ,x − ∆x c * ,x x * . Notice that from Lemma 4, we have ) and we denote δ c = c * −ĉ * . The bound on δ c 2 is from Lemma 6.

F.7 Proof of Lemma 8
Proof of Lemma 8. Similarly to Lemma 7, we first see d(x c * ,x ,x c * ,x ) ≤ C 1 ∆x c * ,x − ∆x c * ,x x * due to Lemma 2. Again using Lemma 4, we see Thus we can bound Then using Lemma 6 to bound c * 2 , we obtain

F.8 Proof of Theorem 1
Proof of Theorem 1. Following the decomposition of error, we show Now we maximize the bound over ĉ * . From [46, Proposition A.1], we see the maximum of which completes the proof.

F.10 Proof of Lemma 9
Proof of Lemma 9. First, we write By Lemma 1, we obtain where we use the fact that C κ (x) is increasing for x > 0. In addition, from Lemma 2, Using Lemma 14, we obtain Let e = ∆x c,x − k i=0 c i ∆ x i . Now combining the above results gives Under Assumption 1 and C κ (x) = O(x 2 ), we see e x * = O(d 3 (x 0 , x * )).

F.11 Proof of Lemma 10
Proof Under the assumption and definition of geodesic µ-strongly convex, we see g(t) is µ-strongly convex in t. Thus, we have g (t 0 ) − g (0) ≥ µt 0 for any t 0 . Let t 0 = ∆x c,x x * and denote the geodesic γ(t) := Exp x * (tη). We derive that g (t) = gradD(Exp x * (tη)), γ (t) by chain rule. Then we have g (t 0 ) = gradD(x c,x ), γ (t 0 ) xc,x = 0 and g (0) = gradD(x * ), η . Finally, we see where the first inequality is due to Cauchy-Schwarz inequality. The first claim is proved by noticing gradD(x * ) x * = k i=0 c i ∆ x i x * and ∆x c,x x * = d(x c,x , x * ). For the second claim, we first observe from the proof of Lemma 14 that where the order can be seen due to that d(x c,x , x * ) ≤ 1 µ k i=0 c i d(x i , x * ) = O(x 0 , x * ) from the first claim. Thus letε := Exp −1 xc,x (x i ) − Γx c,x x * Exp −1 x * (x i ) − Exp −1 x * (x c,x ) , we have ε xc,x = O(d 3 (x 0 , x * )). From the first order stationarity, we see Taking the norm and using the isometry of parallel transport, we obtain the desired result.

G Proofs under general retraction and vector transport
Here we show that when we use general retraction Retr in place of the exponential map Exp, thus avoiding the lemma on metric distortion (Lemma 1, 2), we can still show a similar result as Lemma 4 but with an error on the order of O(d 2 (x 0 , x * )) instead of O(d 3 (x 0 , x * )) as for the case of exponential map. The main idea of proof follows from [54]. The next proposition formalizes such claim. For this section, we denote ∆ x = Retr −1 x * (x) for any x ∈ X where the retraction has a smooth inverse. For general retraction, the deviation is on the order of O( ∆ x 0 Proposition 3. Suppose all iterates x i ∈ X , a neighbourhood where retraction has a smooth inverse. Consider the weighted averagex c,x =x k given by (Avg.1) with retraction. Assume the sequence of iterates is non-divergent in retraction, i.e. ∆ x i x * , ∆x i x * = O( ∆ x 0 x * ). Then we have ∆x c,x = k i=0 c i ∆ x i + e, with e x * = O( ∆ x 0 2 x * ), Proof. The proof generalize the proof of [54, Lemma 12]. First denote Retr y x := Retr −1 y • Retr x and we notice that where we denote γ i = c i i j=0 c j and F : In addition, it can be verified that F (∆x i ) = ∆x i . Now by chain rule, we have where the third inequality uses the inverse function theorem. Hence the Taylor expansion of F at ∆x i up to second order gives where the second equality follows by expanding the first equality. Let i = k − 1, this leads to where we let e = k−1 x * ) due to the non-divergent assumption. The proof is complete.