{'fraction_non_alphanumeric': 0.03146743804237354, 'fraction_numerical': 0.013056055038699086, 'mean_word_length': 4.5290685109141995, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated Algorithms for Convex and Non-Convex Optimization on Manifolds

We propose a general scheme for solving convex and non-convex optimization problems on manifolds. The central idea is that, by adding a multiple of the squared retraction distance to the objective function in question, we"convexify"the objective function and solve a series of convex sub-problems in the optimization procedure. One of the key challenges for optimization on manifolds is the difficulty of verifying the complexity of the objective function, e.g., whether the objective function is convex or non-convex, and the degree of non-convexity. Our proposed algorithm adapts to the level of complexity in the objective function. We show that when the objective function is convex, the algorithm provably converges to the optimum and leads to accelerated convergence. When the objective function is non-convex, the algorithm will converge to a stationary point. Our proposed method unifies insights from Nesterov's original idea for accelerating gradient descent algorithms with recent developments in optimization algorithms in Euclidean space. We demonstrate the utility of our algorithms on several manifold optimization tasks such as estimating intrinsic and extrinsic Fr\'echet means on spheres and low-rank matrix factorization with Grassmann manifolds applied to the Netflix rating data set.


Introduction
Optimization is a near ubiquitous tool used in a wide-range of disciplines including the physical sciences, applied mathematics, engineering and the social sciences. Formally, it aims to maximize or minimize some quantitative criteria, namely, the objective function with respect to some parameters of interest. In many broad, complex learning in modern data science, the parameters are naturally defined over to be on a manifold. The emerging field of statistics on manifolds based on Fréchet means (Bhattacharya and Bhattacharya, 2012;Bhattacharya and Lin, 2017) can be viewed as one of the notable examples of optimization on general manifolds.
Another example can be found in building scalable recommender systems where extracting a low-rank matrix involves an optimization problem over a Grassmann manifold (Boumal et al., 2019). Recent development in geometric deep learning, where the input or output layer constrained to be on a Riemannian manifold (Lohit and Turaga, 2017;, constitutes another important class of applications. Other applications arise in diverse areas ranging from medical imaging analysis, Procrustes shape matching, dimension reduction, dynamic subspace tracking, and cases involving ranking and orthogonality constraints-among many others. This proliferation of manifold-valued applications demands fundamental development of models, algorithms and theory for solving optimization problems over non-Euclidean spaces. The current literature on optimization over manifolds mainly focuses on extending existing Euclidean space algorithms, such as Newton's method (Smith, 2014;Ring and Wirth, 2012), conjugate gradient descent (Edelman et al., 1998;Nishimori et al., 2008), steepest descent (Absil et al., 2010), trust-region methods (Absil et al., 2007;Boumal and Absil, 2011) and others. Many of the objective functions in manifold optimization problems are very complex. One of key challenges for solving such problems lies in the difficulty in verifying the convexity and the degree of convexity of the objective function. Current approaches cannot adapt to the complexity of the problem at hand in manifold spaces.
We take a major step to address these issues by proposing a general scheme to solve convex and non-convex optimization problems on manifolds using gradient-based algorithms originally designed for convex functions. The key idea is to "convexify" the objective function by adding a multiple of the squared retraction distance. The proposed algorithm does not require knowledge of whether the objective function is convex but will automatically converges to an optima if the function is strongly convex. When the objective is non-convex, it achieves rapid convergence to a stationary point. The proposed algorithm is a generalization of Nesterov acceleration (Nesterov, 2004), which improves the convergence rate of gradient descent algorithms. Our algorithm (which we call A 2 ) takes any general existing optimization method (which we call A ), originally designed for convex functions, and converts it into a method applicable for non-convex functions.
Similar schemes have been explored for optimization problems in Euclidean space (Paquette et al., 2018). Generalizations to arbitrary manifolds, however, require fundamentally novel theoretical development. In the Euclidean case, the gradient steps are taken towards lines, whereas for the manifold we use the retraction curves which crucially affects the result and raises the difficulty in proving convergence. Also for manifolds, it is not trivial to correctly convexify the 'weakly-convex function', a broad class of non-convex functions on manifolds we consider which account for most of the interesting examples of non-convex functions in machine learning. We propose a novel idea to convexify the objective locally with the help of the retraction. Key features of our algorithm include adaptation to the unknown weak convexity of the objective function and automatic Nesterov acceleration. The proposed algorithm can be used to accelerate a broad class of A algorithms including gradient descent as well as parallel optimization approaches (see Saparbayeva et al., 2018).
Our paper is organized as follows: In Section 2, we introduce related work on accelerated optimization algorithms. Next, we present our proposed acceleration algorithm on manifolds in Section 3 and present theoretical convergence results. In Section 4, we consider a simulation study of estimating Fréchet means and a real data example using the Netflix prize data set in a matrix completion problem. Liu et al. (2017) propose accelerated first-order methods for geodesically convex optimization on Riemannian manifolds. This is a direct generalization of Nesterov's original linear extrapolation mechanism to general Riemannian manifolds via a non-linear operator. One drawback of Liu et al. (2017) is that the accelerated step of their algorithm involves exact solving of non-trivial implicit equations.

Related work
Zhang and Sra (2018a) later proposed a computationally tractable accelerated gradient algorithm and a novel estimation sequence for convergence analysis. Our approach is fundamentally different from theirs. We regularize an objective function with a squared retraction distance (see Proposition 1), solve a sequence of convex subproblems, adapt to the degree of weak convexity of the objective function, and produce accelerated rates for convex objectives. Even in the convex case, our approach can deal with a much broader class of retraction-based convex functions. Paquette et al. (2018) proposes a general scheme called "Catalyst acceleration" for solving general optimizations in Euclidean space, which has inspired development of some ideas for our work. Similar ideas have been explored for convex functions in Euclidean space in both theory and practice . However, optimization problems on manifolds are of fundamentally different nature and require development of substantially new tools and theory.
There is an interesting line of work on proposing fast algorithms for stochastic optimization on manifolds (see Zhang et al., 2018;Zhou et al., 2019;Bonnabel, 2013) which employ very different techniques such as minibatching, variance reduction and utilizing the uncertainty of inputs. Methods like Zhang and Sra (2018b) propose optimization methods that are analogous to Nesterov-type algorithms for manifold spaces. 3 Accelerated algorithms for optimization on manifolds 3.1 Weakly convex functions on manifolds with respect to retraction mapping We first define general retraction-based, weakly convex, convex, and strongly convex functions by generalizing from their geodesic-based counterparts . We then prove an important proposition that can transform a non-convex function into a convex one simply by adding a multiple of the squared retraction-distance to the objective function.
Definition 1. A retraction on a manifold M is a smooth mapping from its tangent bundle R : T M → M with the following properties: where 0 θ denotes the zero vector on the tangent space 2. For any point θ ∈ R the differential d (R θ ) of the retraction mapping at the zero vector 0 θ ∈ T θ M has to be equal to the identity mapping on T θ M , that is The exponential map on a Riemannian manifold can be viewed as a special case of the retraction map, and the inverse-exponential map is a special case of the inverse-retraction map. A good choice of retraction map can lead to substantial reduction in computation burden compared to the exponential map. We see an example in Section 4.2 on the choice of a retraction map for Grassmannian; Figure 1 provides a visualization of a retraction map.
We first define the retraction distance function on M Since at zero the differential of the retraction map is the identity, there is a small enough neighborhood D of the point θ where the inverse retraction map R −1 θ is bi- where ϑ 1 , ϑ 2 ∈ D, and K 1 ≥ 1, and K 2 ≥ 1.
In addition, we also require the squared retraction distance function to be 2R 1strongly retraction convex around ϑ-that is, for some δ > 0 and constant 0 ≤ R 1 ≤ 1 the following inequality holds: Due to the fact that at the zero vector 0 ϑ ∈ T ϑ M the differential of R ϑ is equal to identity mapping, we can see that in a small neighborhood of ϑ, the square retraction distance function behaves like the square normal function which is strongly convex.
We now define the notion of convex functions on manifolds with respect to the retraction map.
Definition 3. A function f is convex with respect to the retraction R if for any points θ 1 , θ 2 ∈ M the inequality holds (2) Now we are ready to define one of the most important classes of non-convex functions called weakly-convex functions which constitute many interesting applications of non-convex functions in machine learning.
Definition 4. A function f is ρ-weakly convex with respect to the retraction R if for any points θ 1 , θ 2 ∈ M the inequality holds Given the strong retraction convexity of the squared retraction distance (see (1)), we can regularize the weakly convex function f by adding the term κ 2 d 2 R (θ, ϑ) and turn it into a convex function through the following proposition.
Proposition 1. Let d R be a retraction distance that is strong-retraction convex or satisfies the inequality (1) in the subset D ⊂ M . Then the function f is R 1 κ-weakly convex in D if and only if the function Proof. Let f be ρ-weakly convex. Then for any θ 1 , θ 2 ∈ D and any λ ∈ [0, 1] For functions defined on an Euclidean space we have a definition of a weakly convex function that is equivalent to (3): Over the manifold, however, there is no such straightforward equivalence. This is due to the distance function d 2 R (ϑ, θ) which does not satisfy the following equality: Nevertheless in some neighborhood of ϑ,, for some δ > 0, the following inequality holds where d R (θ 1 , ϑ) < δ and d R (θ 2 , ϑ) < δ.
Therefore the function f is ρ-weakly convex with respect to the retraction R if for any points θ, ϑ ∈ M such that λ ∈ [0, 1], the approximate secant inequality holds

The acceleration algorithm on manifolds
In this section, we propose our acceleration algorithms for convex and non-convex functions on manifolds. We first minimize the convex subproblem of an objective function f for some existing approach A (such as a gradient descent algorithm) where the objective function is written as with a positive regularization parameter κ. Proposition 1 ensures the convexity of the subproblem for an appropriate level of regularization. Therefore, with an existing approach A , we define the proximal operator where ϑ is a prox-center.
To consider optimizing p(ϑ), we focus on A having linear convergence rates. Specifically, a minimization algorithm A , generating the sequence of iterates (θ k ) k≥0 , has a linear convergence rate if there exists τ A , f ∈ (0, 1) and a constant C A , f ∈ R such that There are multiple optimization algorithms on manifolds with linear convergence rates for strongly-convex functions on manifold. These include gradient descent, conjugate gradient descent, MASAGA (Babanezhad et al., 2018), RSVRG , and many others.
For a proximal center ϑ and a smoothing parameter κ, we let At the k-th iteration, given a previous iterate θ k−1 and the extrapolation termθ k−1 , we perform the following steps: 2. Accelerated proximal point step.
One needs a stopping criterion, since we cannot use the functional gap as a stopping criterion here as in the convex case. A stationarity stopping criterion is adopted which consists of two conditions: Here, dist(·, ·) denotes the standard Euclidean distance on the tangent space.
Recall that a quadratic of the retraction distance is added to f to make the subproblem convex. So if the weak-convexity parameter ρ is known, then one should set κ > ρ to make the problem convex. In this case, it is proven that the number of inner calls to A for the subproblems can be bounded by proper initialization point ϑ 0 : • if f is smooth, then set ϑ 0 = θ; However, in general one does not have knowledge of ρ. Thus we propose a method that allows algorithm A 1 (Algorithm 1) to handle the convexity problem adaptively.
Our idea is to let A run on the subproblem for T predefined iterations, output the pointθ T , and check if a sufficient decrease occurs. If the subproblem is convex, then the aforementioned descent and adaptive stationary conditions are guaranteed. If either of the conditions are violated, then the subproblem is deemed non-convex. In this case, we double the value κ and repeat the previous steps The tuning parameter κ should be chosen big enough to ensure the convexity of the subproblems and simultaneously small enough to obtain the optimal complexity by not letting the subproblem deviate too far away from the original objective function. Thus we introduce κ cv x as an Adependent smoothing parameter. Notice that the linear convergence rate τ A ,h κ of A is independent of the prox-center and varies with κ. We define κ cv x as by running T iterations of A , using the initialization strategy described below Equation (7).
Finally, for an initial estimate θ 0 ∈ M , smoothing parameters κ 0 , κ c v x , an optimization algorithm A , and a stopping criterion based on a fixed budget T and S, we have the following acceleration algorithm, A 2 , for the manifold (Algorithm 2).

Choose θ k to be any point satisfying
until the stopping criterion is dist(∂ f (θ k ), 0θ) < ε; Remark 1. Note that there are two sequences θ k and θ k in Algorithm A 2 . Since the extrapolation step is designed for the convex case, the second sequence {θ k } approximates the optimal point with accelerated rate which means that it approaches the optimal point faster than the first sequence θ k above. Intuitively, when the first sequence is chosen it uses the initial algorithm A and adapts the smoothing parameter to our objective-implying that the Nesterov step failed to accelerate convergence.
In the adaptation method A 1 (θ k−1 , κ k−1 , T ), the resultingθ k and κ k have to satisfy the following inequalities dist 0θ The resultingθ k , needs to satisfy the condition that if the function f is convex, then We then have the following lemma: then the inequality holds: Proof. We can find v ∈ ∂h κ (θ, ϑ) with v ≤ ε. Taking into account ∂h κ (θ, ϑ) = ∂ f (θ) + κ∇d 2 R (θ, ϑ) the result follows.
Since we assume retraction distance function d R is continuous, we can deduce that the vector field ∇d 2 R (θ, ϑ) is continuous, so the conditions of Lemma 1 are very mild. Also, as mentioned previously, the square retraction distance function d 2 R (·, ϑ) acts like a square normal function in a small neighborhood of ϑ.
We define the following retraction-based strongly convex function: Definition 5. A function f is µ-strongly convex with respect to the retraction R if for any points θ 1 , θ 2 ∈ M and µ > 0 the inequality holds Then we have the following convergence analysis for the acceleration algorithm A 2 : Theorem 1. Fix real-valued constants κ 0 , κ c v x > 0 and the point θ 0 ∈ M . Set κ max = max k≥1 κ k . Suppose that the number of iterations T is such thatθ k satisfies (17), and Then for any N ≥ 1, the iterated sequence generated by the acceleration algorithm satisfies If in addition the function f is κ c v x (K 4 1 K 4 2 − R 1 )-strongly convex and S k is chosen so thatθ k satisfies (20), then where θ * is any minimizer of the function f .
The detailed proof of this theorem can be found in the Appendix.
Remark 2. If the original method A has a linear rate of convergence then our method A 2 also converges to the local minimum for the strongly convex case. If the knowledge of the strong-convexity is given, then some existing method can achieve optimal linear rate for smooth and convex functions , however, it is overall an extremely difficult to verify convexity of a function on a manifold, and our method adapts to that without requiring the knowledge of the convexity. Note that our algorithm also applies to the subgradient descent method, where instead of gradient of the function one takes the subdifferential, for non-smooth functions. In this case, for the stronglyconvex objective, the subgradient method converges to the optimum with O(1/N ) rate of convergence (see ). Thus our accelerated rate O(1/N 2 ) can be considered optimal for strongly-convex functions on the manifold.

Simulation study and data analysis
To examine the convergence and acceleration rates of our proposed algorithm, we first apply our method to the estimation of both intrinsic and extrinsic Fréchet means on spheres, in which one has the exact optima for comparison in the case of extrinsic mean. We also apply our algorithm to the Netflix movie-ranking data set as an example of optimization over Grassmannian manifolds in the low-rank matrix completion problem.

Estimation of intrinsic Fréchet means on manifolds
We first consider the estimation problem of Fréchet means on manifolds (Fréchet, 1948). In this simple example, we have observations {x 1 , . . . , x N } that lie on a sphere S d and our goal is to estimate the sample mean: If ρ is the embedded distance metric in the Euclidean space, then there exists a closed which is the projection of the Euclidean mean x onto the sphere (Bhattacharya and Bhattacharya, 2012). This is called the extrinsic mean. When ρ is taken to be the geodesic or intrinsic distance,θ is called the intrinsic mean. We will consider estimation of both extrinsic and intrinsic means using our method compared to other optimization techniques.
One simple examples of a retraction map for S d is where | · | is the Euclidean norm in R d +1 . Therefore the inverse retraction has the following expression We first compare our accelerated method against gradient descent optimization and a Newton-type optimization scheme, DANE (Shamir et al., 2014), and a Nesterov method, RAGD (Zhang and Sra, 2018b), adapted for manifolds. For all the experiments in this section, we optimized the step size of the optimizer using an Armijo condition backtracing line seach (Armijo, 1966) where we reduce the step size by a factor of .95 until the difference between the old loss function evaluation and the new one is 10 −5 × .95. For our Catalyst algorithm manifold we set the A2 budget to S = 10, the A1 number of iterations to T = 5, and cutoff parameter for A1 is initialized at .1. For the DANE results, we set the regularization term to 1, For RAGD we set the shrinkage parameter to 1. Our synthetic data set is 10,000 observations generated i.i.d from a 100 dimensional N(0, I ) distribution projected onto S 99 .
We run each optimization routine for 100 iterations. Figure 2 and Figure 3 shows that our novel accelerated method converges, for an intrinsic mean as well as an extrinsic mean example, to an optima in fewer iterations than the other competing methods, both in terms of the loss function value and the norm of the loss function gradient. Moreover, we can see in the intrinsic mean example, our method is able to obtain a smaller loss function and gradient norm than the competing methods. In the extrinsic mean example, our method obtains a comparable loss function value and MSE between the learned parameter and the closed-form expression of the sample mean with other methods in fewer iterations and obtains a smaller gradient norm than the competing methods.
By explicit calculation we show the objective functions are strongly convex over a neighborhood of any point on the manifold (see the Appendixx for a proof). This is a highly non-trivial task for general objective functions, hence necessitating an adaptive method such as ours. Moreover, in the extrinsic mean example, since we have a closed form expression of the Fréchet mean we also show that our optimization approach converges to the true extrinsic mean in terms of mean squared error faster than the other optimization methods.

Real data analysis: the Netflix example
Next, we consider an application of our algorithm to the Netflix movie rating dataset. This dataset of over a million entries, X ∈ R M ×N , consists of M = 17770 movies and N = 480189 users, in which only a sparse subset of the users and movies have ratings. In order to build a better recommendation systems to users, we can frame the problem of predicting users' ratings for movies as a low-rank matrix completion problem by learning the rank-r Grassmannian manifold U ∈ Gr(M, r) which optimizes for the set of observed entries (i , j ) ∈ Ω the loss function where W is an r -by-N matrix. Each user k has the loss function L (U , k) = 1 This results in the following gradient For this problem on Grassman manifolds, we have the retraction map: Figure 4: Results for the parallel (right) and reduced (left) Netflix example. . and the inverse retraction map: We look at a comparison of our method against a standard gradient descent method on a subset of the data where we only observe a million ratings (≈ 1.5% of the full data set). In this setting we fix the matrix rank r = 5 and the regularization parameter λ = .01. Figure 4 shows that our accelerated method obtains a smaller loss function value, a smaller identical test set MSE, and nearly identical loss gradient norm faster than RAGD, DANE, or a typical gradient descent approach.
On a large scale, we apply a parallelized version of our accelerated method and a communication-efficient parallel algorithm on manifolds proposed in (Saparbayeva et al., 2018, ILEA) on the full Netflix dataset. We randomly distribute the data across 64 processors and run the optimization routine for 200 iterations. In Figure 4, again we can see steady acceleration that our method provides in terms of the loss function value across iterations and the loss of gradient norm though ILEA obtains slightly better test set MSE than our method.

Conclusion and Discussion
We propose a general scheme for solving non-convex optimization on manifolds which yields theoretical guarantees of convergence to a stationary point when the objective function is non-convex. When the objective function is convex, it leads to accelerated convergence rates for a large class of first order methods, which we show in our numerical examples. One of the interesting future directions we want to pursue is proposing accelerated algorithms on statistical manifolds (manifolds of densities or distributions) by employing information-geometric techniques, and applying the algorithms to accelerate convergence and mixing MCMC algorithms.
6 Appendix 6.1 Proof to Theorem 2 We first introduce a simple lemma.
Lemma 2. Suppose the sequence {α k } k≥1 is produced by A 2 . Then, the following bounds hold for all k ≥ 1 Proof of Theorem 2. The descent condition in dist 0θ implies { f (θ k )} k≥0 are monotonically decreasing. From this Using condition (17), we apply Lemma 2 with ϑ = θ k−1 , θ =θ k and ε = κK d R (θ k , θ k−1 ); hence Combining the above inequality with (18), one has Summing j = 1 to N , we can conclude Set µ k = 1 k+1 . Completing the square yields thereby with inequality from Lemma 1 we conclude

Strong convexity of the objective function in estimating the intrinsic Fréchet means on the sphere
We provide a proof that the objective functions in estimating both the intrinsic and extrinsic Fréchet means on the sphere in Section 4 is strongly convex.
Proof. In order to prove the strong-convexity of the intrinsic mean on the sphere S n , we will prove the strong-convexity of the square intrinsic distance function from the point x 0 ∈ S n d 2 g (x 0 , x) = arccos 2 (x T 0 x).
So for the geodesic from the point x 1 ∈ S n to the point x 2 ∈ S 2 γ(λ) = exp x 1 λ log x 1 x 2 = cos λ arccos(x T 1 x 2 ) x 1 + sin λ arccos(x T 1 x 2 ) we need to show following inequality For the sake of briefness let's use the following notations d 1 = arccos(x T 0 x 1 ), d 2 = arccos(x T 0 x 2 ), d 2 = arccos(x T 1 x 2 ).
In order to proof the strong-convexity of Fréchet function in estimating extrinsic mean on the sphere S n , we will prove the strong-convexity of the square extrinsic distance function from the point x 0 ∈ S n d 2 e (x 0 , x) = 2(1 − x T 0 x).