{'fraction_non_alphanumeric': 0.05527223765042498, 'fraction_numerical': 0.015147689977278465, 'mean_word_length': 4.210101718695195, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 15, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Riemannian Hamiltonian methods for min-max optimization on manifolds

In this paper, we study the min-max optimization problems on Riemannian manifolds. We introduce a Riemannian Hamiltonian function, minimization of which serves as a proxy for solving the original min-max problems. Under the Riemannian Polyak--{\L}ojasiewicz (PL) condition on the Hamiltonian function, its minimizer corresponds to the desired min-max saddle point. We also provide cases where this condition is satisfied. To minimize the Hamiltonian function, we propose Riemannian Hamiltonian methods (RHM) and present their convergence analysis. We extend RHM to include a consensus regularization and to the stochastic setting. We illustrate the efficacy of the proposed RHM in applications such as subspace robust Wasserstein distance, robust training of neural networks, and generative adversarial networks.


Introduction
In this paper, we consider the Riemannian manifold constrained min-max problem where M x , M y are Riemannian manifolds and f : M x × M y − → R is a jointly smooth real-valued function. The aim is to find a global saddle point (x * , y * ) that satisfies for all (x, y) ∈ M x × M y , f (x * , y) ≤ f (x * , y * ) ≤ f (x, y * ). ( Examples of Riemannian manifolds of interest include the sphere manifold, the Stiefel manifold, the manifold of orthogonal matrices, the manifold of doubly stochastic matrices, and the symmetric positive definite manifold, to name a few [3,12,67,10]. When both M x , M y are the Euclidean space, problem (1) reduces to the classical min-max problem, which has been widely studied for applications including adversarial training [40], robust learning [19], generative adversarial networks [22,5], constrained optimization [9], and fair statistical inference [39], among others. When f is convex in x and concave in y (convex-concave), the existence of a global saddle point is guaranteed by the well-established minimax theorem [53,69]. Algorithms converging to such saddle points include the optimistic gradient descent ascent (OGDA) algorithm [59] and the extra-gradient algorithm (EG) [21], which have been analyzed in [52,48,49,47]. For

Contributions
In this paper, we propose a class of methods for solving the min-max problem (1) on Riemannian manifolds, which we call Riemannian Hamiltonian methods (RHM). The idea is to minimize the squared norm of the Riemannian gradient of (1), known as the Riemannian Hamiltonian. Minimizing the Hamiltonian function serves as a good proxy for solving problem (1). Under the Riemannian Polyak-Łojasiewicz (PL) condition [78] on the Hamiltonian function, its minimizer recovers the desired saddle point. In the Euclidean space, the Hamiltonian gradient methods have been considered for solving min-max problems [1,6,44,38], which show great promise in accelerating and stabilizing the convergence to saddle points. This paper generalizes many of those analysis to Riemannian manifolds.
It should be emphasized that the proposed generalization to manifolds is nontrivial as the analysis for the Euclidean counterparts, such as in [1], rely heavily on the matrix properties of the Jacobian. Generalization to Riemannian manifolds require adherence to Riemannian operations independent of the matrix structure. Another challenge is to deal with the varying inner product (Riemannian metric) structure on manifolds. We handle the above by devising novel proof strategies and proposing a metric-aware Riemannian Hamiltonian function that respects the manifold geometry.
In particular, we show global linear convergence of any Riemannian solver to saddle points of problem (1) as long as the Riemannian Hamiltonian of f satisfies the Riemannian PL condition [78]. We show this occurs when f is geodesic strongly convex geodesic strongly concave, and also for some nonconvex functions with sufficient geodesic linearity. We additionally extend the proposed RHM to incorporate a consensus regularization and to the stochastic setting, and prove their convergence. Existing Riemannian algorithms for solving (1) such as [79] make use of the exponential map to update the iterates on the manifolds. In this work, we discuss convergence results with exponential as well as general retraction maps on manifolds.
We empirically show the convergence of our proposed RHM algorithms for different min-max functions and compare them with existing baselines. We further demonstrate the usefulness of RHM algorithms in various applications such as learning subspace robust Wasserstein distance, robust training of neural networks and training of generative adversarial networks.

Organizations
The rest of the paper is organized as follows. Section 2 reviews the preliminary knowledge on Riemannian geometry and Riemannian optimization as well as introduces various functions classes on Riemannian manifolds. We also briefly discuss the existing literature on mix-max optimization in the Euclidean space and on Riemannian manifolds. In Section 3, we propose the Riemannian Hamiltonian function and RHM algorithms, as well as analyze their convergence under the Riemannian PL condition. We provide three cases when such condition is satisfied. Section 4 introduces and analyzes the Riemannian Hamiltonian consensus method. Sections 5 and 6 extend the proposed methods to stochastic settings and to the case of retraction. In Section 7, we empirically compare our algorithms with different baselines on various applications. Section 8 concludes the paper.

Preliminaries
In this section, we give a brief overview of Riemannian geometry and relevant ingredients required for Riemannian optimization. For a more complete treatment of the topic, see [3,12]. We also briefly discuss some of the existing works on min-max optimization.

Basic Riemannian geometry
Riemannian manifold M is a manifold with a Riemannian metric, which is a smooth, symmetric positive definite function g : T p M × T p M − → R on every tangent space T p M, with p ∈ M. It is usually written as an inner product ·, · p . The metric structure induces a norm for any tangent vector ξ ∈ T p M, which is ξ p := ξ, ξ p . For a linear operator on the tangent space H : T p M − → T p M, its operator norm is defined as H p := max ξ∈TxM: ξ p=1 H[ξ] p .
A geodesic on the manifold γ : [0, 1] − → M is the locally shortest curve with zero acceleration. The exponential map at p, Exp p : T p M − → M is defined as the end point of a geodesic along the initial velocity. That is, Exp p (ξ) = γ(1) where γ (0) = ξ, γ(0) = p for any ξ ∈ T p M. Riemannian distance is computed as d(p, q) = 1 0 γ (t) γ(t) dt where γ(t) is the distance minimizing geodesic connecting p, q ∈ M. In a totally normal neighbourhood Ω where there exists a unique geodesic between any p, q ∈ Ω, the exponential map has a well-defined inverse Exp −1 p : M − → T p M and the Riemannian distance can be written as d(p, q) = Exp −1 p (q) p = Exp −1 q (p) q . Parallel transport Γ q p : T p M − → T q M transports tangent vector along the geodesic while being isometric, i.e., ξ, ζ p = Γ q p ξ, Γ q p ζ q for any ξ, ζ ∈ T p M.

Riemannian product manifolds
The product of Riemannian manifolds M = M x ×M y is a Riemannian manifold with the Riemannian metric defined as, for any p = (x, y) ∈ M, and (u, u ), My y , where ·, · Mx , ·, · My are Riemannian metrics on M x , M y respectively. From the metric, one can derive the geodesic, the exponential map, parallel transport, Riemannian distance, which also admit a product structure. See more details in [12].

Riemannian optimization ingredients
Riemannian optimization treats the constrained problem as an unconstrained problem on manifold by generalizing the notions of gradient and Hessian. For a differentiable function h : M − → R, the Riemannian gradient at p, gradh(p) is a tangent vector that satisfies gradh(p), ξ p = Dh(p) [ξ] for any ξ ∈ T p M where Dh(p)[ξ] is the directional derivative of h along ξ. The Riemannian Hessian of h, Hessh(p) : T p M − → T p M is a symmetric linear operator, defined as the covariant derivative of the Riemannian gradient. For a bi-function f : M x × M y − → R, we can similarly define Riemannian partial gradient grad x f (x, y), grad y f (x, y) as Riemannian gradient for x, y, holding the other variable constant. The Riemannian cross derivative grad 2 and similarly for grad 2 yx f (x, y).

Riemannian geodesic convex optimization
Geodesic convexity [76,12] generalizes the notion of convexity to Riemannian manifold. A geodesic convex set Ω ⊆ M requires for any two points in the set, there exist a geodesic (on M) connecting them that lies entirely in the set. From this definition, any connected, complete Riemannian manifold is geodesic convex itself. A function h : Ω − → R is geodesic convex if for any p, q ∈ Ω, it satisfies that h(γ(t)) ≤ (1 − t)h(p) + th(q) for t ∈ [0, 1] and γ is a geodesic connecting p, q. A function is geodesic linear if it is both geodesic convex and geodesic concave.
We call a function h(p) g-(strongly)-convex if it is geodesic (strongly) convex. Similarly, we call a function f (x, y) g-(strongly)-convex-concave if it is geodesic (strongly) convex in x and geodesic (strongly) concave in y.
Next, we define the spectrum of a linear operator on the tangent space, which is used to analyze the Riemannian Hessian as well as the Riemannian cross derivatives in the subsequent sections.
Definition 1 (Spectrum of a linear operator). Consider a linear operator T : V − → W where V, W are two inner product spaces. If V = W , and T is symmetric, i.e., T = T * , where T * is the adjoint operator of T , then we say (λ, v) is an eigenpair of T if T [v] = λv. In general, when V = W , the singular value σ of T is the square root of the eigenvalues of T * • T .
We use λ min /λ max and σ min /σ max to represent the smallest/largest eigenvalues and singular values, respectively. We also use λ |min| to denote the minimum eigenvalue in magnitude. Below, we introduce several function classes on manifolds, generalizing the Lipschitz continuity as well as the Polyak-Łojasiewicz condition from the Euclidean space [62,58] Definition 2 (Lipschitz continuity [12]). Let L 0 , L 1 , L 2 > 0. [78,33,23] The following lemma shows the connection between smoothness of a function on manifold and its Lipschitz Riemannian gradient, which is fundamental for convergence analysis.
Lemma 1 (Lipschitz Riemannian gradient and smoothness [12]). For a function h : M − → R, its Riemannian gradient is L 1 -Lipschitz continuous if and only if Hessh(p) p ≤ L 1 for all p ∈ M. Suppose h has L 1 -Lipschitz Riemannian gradient, then h is L 1 -smooth on M with |h(q) − h(p) − gradh(p), s p | ≤ L 1 2 s 2 p , for all q = Exp p (s) ∈ Ω and p ∈ M.

Min-max optimization
Here we discuss related works on min-max optimization both in the Euclidean space and on Riemannian manifolds.

In Euclidean space
In the Euclidean space (i.e., R n ), the standard gradient descent ascent (GDA) that follows the minmax gradient is known to cycle or diverge for simple convex-concave objectives [43]. To address the cycling issue, the optimistic gradient descent ascent algorithm (OGDA) [59] modifies the GDA update to include an additional gradient momentum. On the other hand, the extra-gradient algorithm (EG) [21] employs an additional min-max gradient step at every iteration. As shown in [46,47], both OGDA and EG methods approximate the proximal point method [61] and converge sublinearly under convex-concave settings [52,48] and linearly under strongly-convex-strongly-concave settings [74,46]. However, for the more general nonconvex-nonconcave settings, finding a global saddle point satisfying (2) is difficult and several existing works [16,4,41] aim to find a local saddle point that satisfies (2) in a local neighbourhood. It should be noted that when the function is convex-concave, all local saddle points are global.
A necessary set of conditions for the saddle points is that they satisfy the first-order stationarity, i.e., the gradients with respect to x and y vanish. This motivates the Euclidean Hamiltonian gradient descent (EHGD) [44,6,1,38] approach for solving the min-max problem, which minimizes the sum of the squares of the gradient norms with respect to x and y. It should be noted that EHGD works under the assumption that all such stationary points are global min-max saddle points [1,38]. Cases are discussed where this assumption is satisfied, which allows EHGD to converge to a global min-max saddle point of the original min-max problem [1,38]. Further, studies [44,6,1,38] demonstrate good empirical performance of EHGD in a variety of applications.
It should be noted that EHGD approaches have only been studied for unconstrained problems in the Euclidean space. Challenges in the constrained settings appear with definition of the Hamiltonian and subsequent analysis.

On Riemannian manifolds
There is a growing theoretical and empirical interest in solving min-max problems under Riemannian optimization framework [37,28,26,79]. An extension of the GDA algorithm to manifolds, named RGDA, has been proposed in [26]. However, [26] considers a min-max setting in which the minimization problem (in x) is on a manifold, but the maximization problem (in y) is on a convex set. In addition, it analyzes the convergence when the maximization problem over y is strongly concave. Hence, [26] does not study the general Riemannian min-max problem (1). It discusses the convergence of their algorithm to first-order stationary points of the min-max problem. Additionally, they propose different stochastic extensions of their algorithm and analyze their convergence.
Recently, [79] has proposed a Riemannian corrected extra-gradient algorithm (RCEG) for the Riemannian min-max problems (1), which contains two steps. First, RCEG takes a step similar to the RGDA update. Then, starting from the newly obtained point, RCEG combines the RGDA direction with the direction of the first step. In the g-convex-concave settings, this correction allows [79] to prove convergence of RCEG to global min-max saddle points of (1).
More details on the RGDA and RCEG algorithms are in Appendix C.

Riemannian Hamiltonian gradient methods
As mentioned earlier, the Euclidean Hamiltonian approach [44,6,1,38] is a popular approach to tackle the min-max problem (1) when M x and M y are restricted to the Euclidean space. Specifically, the Euclidean Hamiltonian function E is defined as, where ∇ x f (x, y) and ∇ x f (x, y) are the partial derivatives of f with respect to x and y, respectively.
Here, · 2 denotes the Frobenius norm. The global minimum of the function E is attained when E(x, y) = 0, i.e., ∇ x f (x, y) = 0 and ∇ y f (x, y) = 0. This corresponds to a first-order stationary point of the function f . Hence, minimization of E in (3), becomes a good proxy to solve the original min-max problem. Building on the Euclidean Hamiltonian approach, a generalization to the Riemannian min-max problem (1) requires understanding of first-order stationary points on manifolds M x and M y . These are necessarily identified with the points where the Riemannian gradient of f vanishes. This leads to our proposed Riemannian Hamiltonian function as where grad x f (x, y) and grad y f (x, y) are the Riemannian partial gradients of f with respect to x and y respectively. Here, grad x f (x, y) 2 x is the square of the gradient norm in the Riemannian metric sense on M x . Similarly, grad y f (x, y) 2 y = grad y f (x, y), grad y f (x, y) 2) Equation (4) generalizes the Euclidean metric considered in (3) to a Riemannian metric. This generalization allows to use other varying metrics for min-max problems in the Euclidean space, e.g., the Fisher information metric [18] or real-projective space metrics [3, Chapter 2].
It should be noted that the Riemannian Hamiltonian (4) can be viewed on the product manifold M = M x × M y , i.e., for p = (x, y) ∈ M, the Riemannian gradient is grad p f (p) = (grad x f (x, y), grad y f (x, y)), and therefore, H(x, y) = grad p f (p) 2 p . Hence, we propose to solve the following problem on the product manifold as Similar to the EHGD approaches [1,38], we work with the following assumption. It is worth noticing that under Assumption 1, solving (5) is equivalent to solving (1). On Riemannian manifolds, Assumption 1 holds when f is g-convex-concave.
We now show that the Riemannian gradient of the Riemannian Hamiltonian H(p) admits a simple expression.
Proof. First, we see that H is a smooth function on the manifold due to the smoothness of f and its Riemannian gradient (formally characterized later in Proposition 2). For any smooth vector field U : M − → T M, denoted as U ∈ X(M), we have U H = gradH, U , where ·, · is the Riemannian metric (on any tangent space). Let ∇ be the Riemannian connection (or the Levi-Civita connection) of M, which provides a way to differentiate vector fields on manifolds. By definition, the Riemannian connection satisfies the metric compatibility property [3,12]  Remark 3. It should be noted that for the Euclidean case when x ∈ R m , y ∈ R n , existing works [6,1,38] analyze the Hamiltonian methods in the form of J v, where J is an asymmetric Jacobian matrix and v is the min-max gradient (∇ x f (x, y), −∇ y f (x, y)). For the same setting, however, Proposition 1 obtains the Hamiltonian gradient as H∇f , where H and ∇f are the (Euclidean) Hessian matrix and gradient vector ∇f = (∇ x f (x, y), ∇ y f (x, y)), respectively. This is not surprising as J v = H∇f . Proposition 1 allows to analyze the performance of the Riemannian Hamiltonian approach in terms of the symmetric Riemannian Hessian operator. The analysis in [1,38] heavily rely on the matrix structure of J and makes use of the linear algebraic properties of the Jacobian. Our approach, thanks to Proposition 1, adheres to general Riemannian manifolds as we directly deal with the operator, which is independent of the matrix structure. Hence, many of the subsequent analysis in this paper differ from [1,38]. Compute the step ξ(p t ) from the gradient gradH(p t ) = Hessf (p t )[gradf (p t )].
To minimize the Riemannian Hamiltonian (5), one can apply first-order Riemannian solvers including Riemannian steepest descent [75], Riemannian conjugate gradient [60], or second-order solvers, such as Riemannian trust-regions [2,11], provided the Hessian (or approximated Hessian) of the Hamiltonian is available. We refer to such class of methods for solving min-max problems on manifolds collectively as Riemannian Hamiltonian methods (RHM). Its procedures are outlined in Algorithm 1, where the step ξ(p t ) is computed depending on the selected solver. Now we analyze the performance of the proposed RHM. In particular, we aim to obtain the global minimizer p * of H, which satisfies H(p * ) = 0 with RHM. However, this may not always be numerically tractable without additional structures on the Riemannian Hamiltonian. One such structure is assuming the Riemannian Hamiltonian is g-convex, for which RHM converges to the optimal p * (g-convexity guarantees convergence to global optimality). This, however, may not lead to interesting problem classes for f . Moreover, there is no guarantee that H is a g-convex even when f is g-convex-concave.
Another interesting structure is the Polyak-Łojasiewicz (PL) condition. The PL condition [58] amounts to a sufficient condition to establish linear convergence for gradient-based methods to global optimality [32]. The Riemannian version of the PL condition (Definition 3) has been studied in [78,33,80,23]. In Section 3.1, we impose the Riemannian PL condition on the Hamiltonian H as it allows convergence of RHM to global optimality. It should be noted that functions satisfying the Riemannian PL condition subsume g-(strongly)-convex functions. In Section 3.2, we discuss many interesting function classes of f that allow the Hamiltonian H to satisfy the condition.

Convergence analysis
To analyze the convergence of RHM, we focus on the Riemannian steepest descent direction in the main text, i.e., ξ(p t ) = −η t gradH(p t ) with either fixed stepsize or variable stepsize computed from backtracking line-search [13,12]. We include the details of implementing the Riemannian conjugate gradient and Riemannian trust-region methods together with their convergence analysis in Appendix F.
We make the following standard assumption [3,12,78,66] throughout the rest of the paper.
In the next proposition, we show that the Riemannian Hamiltonian H is L-smooth.
Proof. According to Lemma 1, it is sufficient to show that the Riemannian gradient of H is L-Lipschitz. From Proposition 1 and Assumption 2, we have for any p ∈ M, q = Exp p (s) ∈ Ω, the domain of exponential map around p, where we apply the triangle inequality and the isometry property of parallel transport.
If the Hamiltonian H satisfies the Riemannian PL condition, then we show that Algorithm 1 with the steepest descent update (RHM-SD) converges linearly to the global minimizer of H. We begin with the convergence result for RHM-SD with fixed stepsize.
Theorem 1 (Linear convergence of RHM-SD with fixed stepsize). Let f satisfy Assumption 2 and H satisfy the Riemannian PL condition, i.e., 1 2 From the smoothness of the Riemannian Hamiltonian H (Proposition 2, Lemma 1) and the gradient update in Algorithm 1, we have where the last inequality employs the Riemannian PL condition. This leads to H(p t+1 ) ≤ (1− δ L )H(p t ). Applying this result recursively completes the proof.
Line-search methods are practically favourable because they adapt the stepsize without requiring the knowledge of the Lipschitz constant L. Here, we consider the backtracking line-search for choosing stepsize η t for Riemannian steepest descent, which is commonly used in practice. Given an initial stepsizeη, the backtracking line-search iteratively decreases the stepsize by a factor of ∈ (0, 1) until the Armijo-type sufficient decrease condition is satisfied, i.e., for some update direction ζ(p t ). The complete procedure is included in Appendix B. We next present the convergence for RHM-SD with backtracking linesearch.
Theorem 2 (Linear convergence of RHM-SD with backtracking line-search). Under the same setting as in Theorem 1, consider Algorithm 1 using the steepest descent direction with backtracking line-search, parameters r 1 , ρ ∈ (0, 1), and an initial stepsizeη. Then, the iterates p t satisfy gradf (p t ) 2 pt ≤ 1 − 2 min ηr 1 , Proof. Given H is L-smooth, the proof follows from [12,Lemma 4.12] and the Riemannian PL condition.

Important problem classes for RHM
We now discuss the instances of f where the Riemannian Hamiltonian satisfies the Riemannian PL condition (Definition 3). This allows RHM (Algorithm 1) to converge to global min-max saddle points of (1). From the expression of gradH(p) in Proposition 1, we observe that if all eigenvalues of Hessf (p) are lower bounded in magnitude ( i.e., |λ| ≥ α > 0), then the Riemannian Hamiltonian H satisfies the Riemannian PL condition with δ = α 2 . This is because Our aim, therefore, is to identify classes of f that satisfy the right hand side of (7). We provide three cases where the Riemannian PL condition is naturally satisfied on the Riemannian Hamiltonian H, which generalize the results in [1] to Riemannian manifolds. These include the cases when the objective f is g-strongly-convex-concave and when f is smooth with sufficient geodesic linearity.
In order to analyze function classes of f that lead to (7), we require the following results on the Riemannian Hessian Hessf (p) of the product manifold M (which are of independent interest as well).

1) Decomposition of the Riemannian Hessian Hessf (p) and adjoint property of the cross derivatives.
This is shown in Appendix D.
2) We establish general lower bounds on the eigenvalue magnitude of the Riemannian Hessian, which we include in Appendix E.
The above results help to bound the eigenvalues of Hessf (p) in terms of the spectrum of Hess x f (x, y), Hess y f (x, y), and the cross derivatives grad 2 xy f (x, y), grad 2 yx f (x, y). We now present the main results below.
Proposition 3 (Geodesic strongly convex strongly concave). Let f (x, y) be geodesic strongly convex in x and geodesic strongly concave in y with parameter µ > 0. Then, H satisfies the Riemannian PL condition (7) with δ = µ 2 .
Proof. We show that if there exists an eigenpair (λ, ξ) of Hessf (p) such that |λ| < µ with p = (x, y), ξ = (u, v), then it leads to a contradiction. From the expression of the Riemannian Hessian in Proposition 9, we have This can be equivalently written as From (8), we obtain where id is the identity operator. From the symmetry of the Riemannian cross derivatives (Proposition 10), we can substitute (10) into (9), which gives The geodesic strong convexity in x and geodesic strong concavity in y leads to Hess x f (x, y) µ id and Hess y f (x, y) −µ id respectively. Thus the LHS of (11) is smaller than −µ, which contradicts |λ| < µ. Thus, all eigenvalues of Hessf (p) satisfies |λ| ≥ µ.

Proposition 4 (Smooth and geodesic linear).
Let σ min (grad 2 xy f (x, y)) ≥ τ > 0 and let f (x, y) be geodesic linear in one variable and has L 1 -Lipschitz Riemannian gradient in another variable. Then, H satisfies the Riemannian PL condition (7) with δ = τ 4 Proof. Without loss of generality, we assume f (x, y) has L 1 -Lipschitz gradient in x and geodesic linear in y. The geodesic linearity in y implies that Hess y f (x, y) = 0, and therefore, we can apply Lemma 2, which shows Finally, from the assumption σ min (grad 2 xy f (x, y)) ≥ τ , the proof is complete.
It is worth noticing that the sufficient geodesic bilinearity condition in Proposition 5 can be interpreted as requiring a sufficiently large weight on the geodesic bilinear component in the objective where f 0 is geodesic linear in each x and y (i.e. bilinear) with the weight c l > 0 and f 1 , f 2 have L 1 -Lipschitz Riemannian gradient. Because by definition, Riemannian Hessian of a geodesic linear function is zero, f has 2L 1 -Lipschitz Riemannian gradient (by Lemma 1). Let τ 0 , Υ 0 be the minimum and maximum singular values of grad 2 xy f 0 (x, y). Then, τ = c l τ 0 , Υ = c l Υ 0 . The sufficient geodesic bilinearity condition is satisfied for

Riemannian Hamiltonian consensus method
In the Euclidean setting, [44] proposes the consensus method for solving min-max problems in the Euclidean space. The consensus method has also been viewed as a perturbation of the Euclidean Hamiltonian method [1]. In this section, we propose an extension of RHM with steepest descent update, namely the Riemannian Hamiltonian consensus method (RHM-CON), by combining the Hamiltonian gradient direction with the min-max gradient direction. The update of RHM-CON is given by is the min-max gradient When γ = 0, this reduces to RHM-SD. The RHM-CON method is formalized in Algorithm 2. Below, we provide the convergence result for RHM-CON.
Theorem 3 (Linear convergence of RHM-CON). Under Assumption 2 with L = L 0 L 2 + L 2 1 , suppose that the Riemannian Hamiltonian H satisfies the PL condition. Let c > 0 such that Proof. First, we highlight that From the smoothness of Riemannian Hamiltonian (Proposition 2, Lemma 1) and the update in where the second inequality follows from η ≤ 1 L (which gives η 2 − η 2 L 2 ≥ 0) and the lower bound on ζ(p t ) 2 pt . The last inequality uses the PL condition and η ≤ 1 L < 1+c Lc , which ensures − η 2 − ηc 2 + η 2 Lc 2 < 0. From the choice of η and γ as well as the definition of ν, we have ν > 0. This is because Applying this result recursively completes the proof.
From Theorem 3, we see that linear convergence is achieved provided that the weight γ on min-max gradient direction is sufficiently small. Also, we highlight that a uniform parameter c > 0 always exists in a compact set as long as γv(p t ) = −gradH(p t ). This can be ensured by choosing a small value for γ.

Stochastic min-max optimization
Applications such as domain generalization, robust training, and generative adversarial networks yield a min-max problem with a stochastic function f , e.g., with a finite sum structure of the function [38]. Under the stochastic setting, the objective function in (1) can be expressed as an expectation, i.e., where ω ∈ D is a random variable following a certain distribution D. This implies an expectation structure on the Riemannian Hamiltonian as for ω, ϕ ∈ D. Modifying Proposition 1 for the stochastic setting leads to where S = {ω 1 , ..., ω |S| }, S = {ϕ 1 , ..., ϕ |S | } are randomly selected subsets with ω i , ϕ j ∈ D.
The stochastic Hamiltonian gradient provides an unbiased estimate of the full gradient, i.e., E S,S [gradH S,S (p)] = gradH(p). We now show the convergence result of RHM-SGD.
Theorem 4 (Convergence of RHM-SGD with fixed and decaying stepsize). Let Assumption 2 hold with L = L 0 L 2 + L 2 1 , and let the Riemannian Hamiltonian H satisfy the PL condition with parameter δ. Assume also that the variance of the stochastic gradient is bounded, i.e., E ω,ϕ gradH ω,ϕ (p t ) 2 pt ≤ G. Then, RHM-SGD with fixed stepsize The proof follows from [32,Theorem 4] and can be easily adapted to the Riemannian manifold setting, and therefore, is omitted.
We can similarly consider the stochastic version of RHM-CON, which we denote as RHM-SCON, with the update step as where v S (p t ) is the stochastic min-max gradient on sample set S. Theorem 4 can be adapted to prove the convergence of RHM-SCON following similar assumptions and analysis in Theorem 3.

Convergence under retraction
Existing algorithms for solving (1), such as RCEG [79], employs the exponential map to update iterates on the manifolds. However, in many cases, the computational cost of implementing the exponential map for many Riemannian manifolds is prohibitive. An alternative is to consider the more general retraction operation [3,Chapter 4]. In this section, we show that the use of retraction (instead of the exponential map) in RHM algorithms guarantees similar convergence under an additional mild assumption.
Retraction R p : T p M − → M is a map that satisfies for all p ∈ M, (1) R p (0) = p (2) DR p (0)[ξ] = ξ for all ξ ∈ T p M. From the definition, we observe that the exponential map is a special case of retraction. In practice, when an efficient retraction is available, the Hamiltonian gradient update can be performed via retraction, i.e., p t+1 = R pt (−η gradH(p t )). To analyze the convergence, we make the following additional assumption that bound the differential operator of the retraction map.
Assumption 3. There exists constants θ 1 , θ 2 > 0 such that the retraction curve c(t) := R p (tξ) with This assumption is always satisfied for a compact manifold M. The compactness appears to be necessary for retraction-based analysis for first-order algorithms [23,66,33,12]. We remark that for the case of the exponential map, the retraction curve coincides with the geodesic curve. Then, p ξ c(t) = 1 by isometric property of parallel transport. Also, θ 2 = 0 from the definition of the geodesic. Proposition 6. Under Assumptions 2 and 3, the Riemannian Hamiltonian H is retraction L Rsmooth with L R = θ 2 1 L + θ 2 L 1 L 0 , i.e., for any p ∈ M, q = R p (ξ) ∈ U, we have H(q) ≤ H(p) + gradH(p), ξ p + L R 2 ξ 2 p .
Proof. For any retraction curve c(t) = R p (tξ) with ξ p = 1 and t ≥ 0 such that c(t) ∈ U, we obtain where the second inequality applies the gradient of Hamiltonian is L-Lipschitz (Proposition 2, Lemma 1) and Assumption 3. The last inequality follows from Assumption 2. The proof from (13) to L R -smoothness of H is due to [29,Lemma 3.2], which we include here for completeness. For any ξ ∈ T p M such that R p (ξ) ∈ U, let α = ξ p , ζ = ξ/ ξ p and hence ξ = αζ with ζ p = 1. Applying Taylor's Theorem on H • R p gives wheret ∈ [0, α]. Thus, the proof is complete.
Using Proposition 6, we show below that RHM-SD attains a linear convergence rate with retraction.
Theorem 5 (Linear convergence of RHM-SD under retraction). Under same settings as in Theorem 1, suppose Assumption 3 holds, and the iterates stay in the compact set U. Then, RHM-SD with retraction and η = 1/L R converges with gradf (p t ) 2 Similar analysis with retraction operation can be performed for other variants of RHM including RHM-CON, RHM-SGD, and RHM-SCON.

Experiments
In this section, we discuss empirical performance of the proposed Riemannian Hamiltonian methods for various min-max optimization problems on manifolds. The algorithms are implemented in Matlab using the Manopt package [14] except for Section 7.4, 7.5 where we use Pytorch with the Geoopt package [34]. We highlight that there exist many other manifold optimization packages, such as ROPTLIB [30], Manopt.jl [8], Pymanopt [73], McTorch [42], and RiemOpt [70], where RHM can also be implemented efficiently. We use the following acronyms for the various RHM algorithms considered in this section.
• RHM-SD: RHM with steepest descent direction with backtracking line search.
• RHM-TR: RHM with the trust-region method where we use Hessian approximation with finite differentiation [11].
• RHM-SGD: RHM with stochastic gradient (Section 5) • RHM-SCON: RHM with stochastic consensus method (Section 5) We compare the proposed Riemannian Hamiltonian methods with the Riemannian gradient descent ascent (RGDA) [26] and the Riemannian corrected extra-gradient (RCEG) [79]. As discussed previously, RGDA has not been studied and analyzed for solving the general min-max problem (1), but when M y is a convex subset of the Euclidean space [26]. In our experiments, however, we extend RGDA to solve (1).
For all the experiments, we implement the algorithms with exponential map for comparability with RCEG, except for the applications of subspace robust Wasserstein distance (Section 7.3), robust training (Section 7.4) and generative adversarial networks (Section 7.5) where we implement with retraction map because the manifolds considered do not have a well-defined logarithm map. Hence, for these applications, RCEG is excluded for comparison. In robust training and generative adversarial network experiments, we also test stochastic algorithms for RGDA and RHM. The codes are available at https://github.com/andyjm3.

Geodesic quadratic bilinear optimization
The first example we consider is where X, Y ∈ S d ++ , the set of d × d symmetric positive definite (SPD) matrices. The weights c l , c q ≥ 0 control the balance between the linear and quadratic terms.
For X ∈ S d ++ , the tangent space T X S d ++ is the set of symmetric matrices. When endowed with the affine-invariant (AI) metric, i.e., U, V X = tr(X −1 UX −1 V), for any U, V ∈ T X S d ++ , one can derive the geodesic, exponential map, and other Riemannian optimization ingredients [24,10,56]. We include the expressions in Appendix A. Here, we use M SPD to represent the SPD manifold with the AI metric. It is worth noticing that the function (14) is nonconvex-nonconcave in the Euclidean space (with details included in Appendix G).
However, the log-det function is geodesic linear on SPD manifold with the AI metric [71] and we show in the following proposition that f (X, Y) is g-convex-concave, although not necessarily g-strongly-convex-concave. (14) is g-convex-concave on M SPD but not g-strongly-convex-concave.

Proposition 7. The function
We next prove that the Riemannian Hamiltonian H of the objective (14) satisfies the PL condition, which allows linear convergence of the proposed RHM algorithms. In Proposition 8, we see that there exist a continuum of global saddle points. Consequently, we define an optimality gap criterion as | det(X) − 1| + | det(Y) − 1| for a candidate point (X, Y).

Experiment settings and results
We consider d = 30 and discuss results on various combinations of c q , c l . We compare our RHM with RGDA [26] and RCEG [79]. All the choices of stepsize are tunned to reflect the best performance except for RHM-SD, RHM-CG, RHM-TR where the stepsizes are selected adaptively by the algorithms. For RHM-CON, we set γ = 0.5. Convergence of an algorithm is measured in terms of gradf (p t ) pt , which is equivalent to 2H(p t ). This measure of convergence has also been considered in [79] for min-max problems on manifolds. Algorithms are stopped either when gradient norm falls below 10 −10 or the max iteration has been reached. Results are reported in Fig. 1.
From Fig. 1, we observe rapid convergence of RHM algorithms in all the settings. The convergence for RGDA varies across different choices of c q , c l where it converges faster when the weight on the quadratic term (c q ) is relatively higher and is not able to converge when c l increases. We also observe convergence for RCEG in all cases but the rate is slower compared to RHM algorithms. In Fig. 1f, we further compare the optimality gap where we observe all the proposed RHM algorithms reach below 10 −10 at a faster rate than the baselines. The slopes of RHM-SD-F and RHM-CON are steeper than that of RCEG (indicating better theoretical rates for RHM). Additional results on optimality gap comparisons are in Fig. 4 in Appendix H. Finally, Fig. 1g shows the runtime performance of various algorithms, with the markers indicating the progress of respective algorithms per iteration. We observe that the per-iteration computational cost of RHM is higher than RGDA. This is because RHM exploits second-order information of f to compute the gradient of H. Also, we see that RCEG can be costly because it requires evaluation of the exponential map twice and the logarithm map once per iteration.

Robust Geometry-aware PCA
Geometry-aware principal component analysis (PCA) on M SPD [25] concerns dimensionality reduction for SPD matrices while preserving geometric structures on the manifold. The robust PCA (or robust Fréchet mean) on SPD manifolds has been considered in [79]. For a set of SPD matrices M i ∈ S d ++ , i = 1, ..., n, the aim is to find the Fréchet mean M ∈ S d ++ that is bounded away from zero, i.e., where α > 0 and S d−1 := {x ∈ R d : x 2 = 1} denotes the sphere manifold and dist : S d ++ × S d ++ is the Riemannian distance on M SPD .
We first note that the function in (15) is geodesic strongly convex in M and geodesic nonconcave in x. Also, it is difficult to verify the Riemannian PL condition on the Hamiltonian of (15). Hence, this is a challenging problem instance as it does not fall into the studied settings of the existing works [26,79] including ours.

Experiment settings and results
For this problem, we follow the same settings as discussed in [79] for generating the SPD matrices M i with the eigenvalues bounded in [µ 0 , µ 1 ]. Following [79], we choose d = 50, n = 40, µ 0 = 0.2, and µ 1 = 4.5. The convergence results are presented in Figs. 2a and 2b, where we only include RHM-SD-F, RHM-CON (γ = 0.5), and RHM-CG for clarity (RHM-TR performs similar to RHM-CG). We observe that although RGDA and RCEG converge faster than RHM when α = 3, they fail  Figure 1: Experiments on the geodesic quadratic bilinear problem (14) with d = 30, under varying weights c q , c l . We observe that our RHM algorithms converge quickly in all settings while baselines such as RGDA [26] and RCEG [79]. The performance of RGDA varies greatly with the settings where it converges only for a few settings and for the others RGDA fails to converge. RCEG presents a relatively more stable convergence behavior than RGDA but with a rate that is slower than our proposed RHM algorithms.
to converge when α = 0.1. The latter finding is not surprising as both RGDA and RCEG seem to perform poorly on approximately bilinear problems (as also observed in Section 7.1). In contrast, we observe that RHM algorithms converge in both the settings, which is also validated by our analysis in Section 3.2. It is known that the conjugate gradient based methods outperforms steepest descent methods on more challenging optimization problems. This explains the faster convergence of RHM-CG over RHM-SD-F and RHM-CON. Overall, the results in Fig. 2 show the benefit of the Riemannian Hamiltonian modeling in non standard settings.

Subspace robust Wasserstein distance
We next consider the problem of learning subspace robust Wasserstein distance [55,37,28], where the aim is to compute the Wasserstein distance over the worst-case optimal transport cost on a low- where δ x is the Dirac at location x. The weights a i , b j belong to the probability simplex, i.e., i a i = j b j = 1. The objective (with entropy regularization) is then given as where St(d, r) := {U ∈ R d×r : U U = I} is the set of column orthonormal matrices (d ≥ r), known as the Stiefel manifold.
is the set of couplings, which forms the so-called doubly stochastic manifold (or coupling manifold) [18,67,45].

Experiment settings and results
We follow the same experiment settings as in [37,28] and consider a uniform distribution over hypercube [−1, 1] d and a pushforward map defined as T (x) = x + 2 sign(x) ( k i=1 e i ), where sign(x) extracts the sign of x elementwise and {e i } d i=1 are the canonical basis of R d . We choose d = 30, r = 5, k = 2, n = 100, = 0.2 and compare the proposed RHM-SD-F, RHM-CON (γ = 0.5), RHM-CG with RGDA in Fig. 2c. RCEG cannot be implemented to solve (16) because the doubly stochastic manifold does not have a well-defined logarithm map. From the results, we see similar convergence speed of all methods while due to the inbuilt line-search algorithm of RHM-CG, it converges to a point with a smaller gradient norm.

Robust training of neural networks with orthonormal weights
We next consider adversarial robust training of deep neural networks with orthonormal weights [26]. Adversarial training of neural networks provide robust prediction against small data perturbations. Orthonormality on parameters has shown to improve generalization accuracy as well as accelerate and stabilize convergence of neural network models [7,17,77,27]. This corresponds to optimization over the Stiefel manifold.
In particular, we consider the adversarial training to defend against a universal perturbation p proposed in [50]. The perturbation set we consider is the sphere manifold S d−1 (r) := {p ∈ R d : p 2 = r} with radius r. This requires the perturbed samples to stay a certain distance away from  the original ones, a strategy also applied in [36]. Given a set of data-target pairs where L(·, ·) is a loss function and h(·) represents the forward function of a neural network.

Experiment settings and results
The adversarial training is implemented for classification tasks on MNIST images [35] where we include two hidden layers of size 16 with the orthonormality constraint. We compare the proposed stochastic version of RHM (RHM-SGD), detailed in Section 5, with Riemannian stochastic gradient descent ascent (RSGDA) algorithm [26]. We highlight that RHM-SCON performs similarly to RHM-SGD, and thus, we exclude its result for clarity. Because we require dual sampling per-iteration to compute the stochastic Hamiltonian gradient gradH S,S (p t ) = 1

|S||S |
ω∈S,ϕ∈S gradH ω,ϕ (p), we choose the batch size to be 32 for both S, S and 64 for RSGDA. Hence, the per-iteration sampling cost is identical. We measure convergence in terms of the relative Hamiltonian H(p t )/H(p 0 ), where the Hamiltonian is evaluated on the full training set. The stepsize is fixed for both the algorithms.
We plot the convergence results (with the best tuned stepsize) in Fig. 3a, which are averaged over five different runs. We see a clear advantage of RHM-SGD compared to RSGDA with faster and more stable convergence.

Orthonormal generative adversarial networks
Generative adversarial networks (GAN) [22,6] are popular in generating synthetic samples by optimizing a min-max game between a generator and a discriminator. The orthonormality constraint on weight parameters of the discriminator has shown to benefit the training of GANs [15,51]. In particular, given samples {x i } n i=1 we consider the following min-max problem where D(·), G(·) represent the discriminator and generator with {W D }, {W G } denoting their network weight parameters respectively. Here, σ(·) is the sigmoid function and the prior z i is sampled from the standard normal distribution.

Experiment settings and results
Following [6], we train the GAN model on 2d samples from a multimodal mixture of Gaussian distribution. The ground truth is shown in Fig. 3c. Both the generator and discriminator have 5 hidden layers with 128 units and ReLU activation. The dimension of the prior z i is 64. For simplicity, we add the orthonormal constraint only for the penultimate layer of the discriminator model. For this experiment, we apply RHM-SCON with γ = 0.5 and compare against RSGDA, both with fixed stepsize. The batch size is chosen to be 128 for RHM-SCON and 256 for RSGDA. Similarly, the best choices of stepsize are reported, and the results are averaged over five different runs. The convergence in terms of the relative Hamiltonian are shown in Fig. 3b, where we see RSGDA diverges while RHM-SCON is more stable. We also examine the solution quality by providing the generated samples from both algorithms at iteration 10 4 , 2 × 10 4 , and 3 × 10 4 in Figs 3d and 3e respectively. We note that RSGDA results in undesired mode collapse, an observation also made in [6] for training SGDA on the Euclidean space. In contrast, RHM-SCON quickly converges and recovers the ground truth distribution. Even though RHM-SGD converges to a lower Hamiltonian value, its performance in recovery of the ground truth is poor, as shown in Fig. 5 in Appendix H where the generated samples collapse to a single point. It indicates that RHM-SGD converges to a stationary point which is not a saddle point (not surprising as Assumption 1 may not be satisfied). This also highlights the practical benefit of consensus regularization for RHM (Section 4), as evidenced in the good performance of RHM-SCON.

Concluding remarks
Building on the success of the Hamiltonian methods for solving min-max problems in the Euclidean space, we have considered a more general problem on manifolds, and proposed a Riemannian Hamiltonian function H that respects the manifold geometry. This leads to a gradient expression (in Proposition 1) that allows simple analysis for the resulting optimization methods. Adapting the proofs from the Euclidean space to Riemannian manifolds requires to forgo the matrix structure of the ingredients, which includes addressing a varying inner product (Riemannian metric). The proposed Riemannian Hamiltonian methods (RHM) come with convergence guarantees and various extensions. The experiments validate the good performance of RHM in different applications. As future work, one direction is to develop variance reduction stochastic extension to RHM. In addition, the current convergence analysis is measured in the Riemannian Hamiltonian, which is the gradient norm squared of the original objective f . It remains a question whether linear convergence can be maintained in terms of the optimality gap on function value of f .

A Riemannian geometries of the considered manifolds
In this section, we review the Riemannian optimization-related ingredients of several manifolds that are considered in the experiments section. The expressions are from the works [3,12,71,67,18,45].

A.1 Symmetric positive definite manifold
Consider the set of the symmetric positive definite matrices of size d × d, S d ++ := {X : R d×d : X = X, X 0}, equipped with the affine-invariant Riemannian metric. The geodesic from X to Y is given by γ(t) = X 1/2 (X −1/2 YX −1/2 ) t X 1/2 . At X ∈ S d ++ , the exponential map is derived as The Riemannian gradient of a function f : S d ++ :− → R is given by gradf (X) = X∇f (X)X, where ∇f (X) is the Euclidean partial derivative of f at X.

A.2 Sphere manifold
It is defined as It can be endowed with the standard inner product at the Riemannian metric, i.e., u, v The Riemannian gradient of f is Proj x (∇f (x)), where ∇f (x) is the Euclidean partial derivative of f at x.

A.3 Stiefel manifold
It is the set St(d, r) = {X ∈ R d×r : X X = I}. It is a generalization of the sphere manifold to higher dimensions and can be similarly endowed with the standard inner product as metric U, V X = U, V 2 For the experiments, we consider the popular QR-based retraction for approximating the exponential map, i.e., R X (U) = qf(X + U), where qf(·) returns the Q-factor from the QR decomposition for any tangent vector U.

A.4 Doubly stochastic manifold
The doubly stochastic manifold (or coupling manifold) between two discrete probability measures j} endowed with the Fisher information Riemnnanian metric. The geometry has been developed in [18,67,45].
approach for balancing non-negative matrices to satisfy the row-sum and column sum constraint and later adapted to solve the optimal transport problem efficiently [57]. Let A ∈ R m×n , A i,j > 0, and denote Sinkhorn(A) as the output of applying the Sinkhorn-Knopp algorithm on A with constraint defined by Π(µ, ν), i.e., Sinkhorn(A) ∈ Π(µ, ν). Subsequently, the retraction is given by R Γ (U) = Sinkhorn(Γ exp(U Γ)), where exp, , and are elementwise exponential, product, and division operations, respectively.

B Line-search methods and Wolfe conditions on Riemannian manifolds
In this section, we provide details of implementing line-search methods on Riemannian manifolds. We first present the Riemannian versions of the Armijo, Wolfe, and strong Wolfe conditions [65].
Definition 4. Consider an iterative algorithm for minimizing h : M − → R, producing p t+1 = Exp pt (ϑ t ξ t ) for some direction ξ t ∈ T pt M and stepsize ϑ t ∈ R. The Armijo condition is h(p t ) − h(p t+1 ) ≥ r 1 ϑ t −gradh(p t ), ξ t , for some r 1 ∈ (0, 1). The (weak) Wolfe condition is the Armijo condition together with (17) and the strong Wolfe condition is the Armijo condition with (18), where for some r 2 ∈ (r 1 , 1). Here, DExp is the differential of the exponential operation.
One can generalize the analysis from the Euclidean space to show that there exists a stepsize that satisfy the three conditions for arbitrary direction ξ t . The backtracking line-search for satisfying the Armijo condition is in Algorithm 3. This has been used in Riemannian steepest descent method [13]. The procedures that return stepsizes satisfying the Wolfe conditions are in [63,54].

C Review of RGDA and RCEG
In this section, we provide the details of the Riemannian gradient descent ascent [26] and Riemannian corrected extra-gradient [79] algorithms for min-max optimization on manifolds.
RGDA simultaneously updates the variables in the direction of the min-max Riemannian gradient, i.e., x t+1 = Exp xt (−η t grad x f (x t , y t )), y t+1 = Exp yt (η t grad y f (x t , y t )).
Given that Hess x and Hess y satisfy the self-adjoint property, we obtain We can see (20) holds for any choice of (u, v), (w, z) and this only happens when grad 2 , v y holds for any (u, v). To see this, consider the vectorization of the tangent vectors as u, v, w, z. We also denote B xy , B yx as the matrix representation of the linear operators grad 2 xy f (x, y), grad 2 yx f (x, y) at (x, y) respectively. Then (20) can be rewritten as where G x , G y are the (symmetric positive definite) metric tensors at x, y. This is equivalent to which is satisfied for any u, v, w, z and any G x , G y as metric tensors. Hence, G y B xy = B yx G x and the proof is complete.
Remark 5. Proposition 10 shows that the Riemannian cross derivatives are symmetric with respect to Riemannian metric on respective manifolds. When M x , M y are the Euclidean spaces, then Proposition 10 is equivalent to the Schwarz's theorem of symmetric second-order derivatives.

E Essential lemmas
The following lemmas generalize [1,Lemmas 17,28] to linear operators specifically in terms of the Riemannian Hessian operator. We first highlight that for two operators T , T * that are adjoint, we have λ(T • T * ) = λ(T * • T ) = σ 2 (T ) = σ 2 (T * ). Proof. We consider the operator Hessf (p) • Hessf (p) and study its eigenvalue. First, we see that for any p = (x, y) ∈ M x × M y and ξ = (u, v) ∈ T x M x × T y M y , we have , and therefore, Suppose (δ, ξ) is an eigenpair of the operator Hessf (p) • Hessf (p), which gives Let B xy = grad 2 xy f (x, y), B yx = grad 2 yx f (x, y), and H x = Hess x f (x, y). Suppose δ < σ 4 min (Bxy) 2σ 2 min (Bxy)+ Hx 2 x < σ 2 min (B xy ). Then, we have B xy • B yx − δ id is invertible where we use the fact that B xy and B yx are adjoint. Hence, from (22) . Substituting the expression of v into (21) yields We next show that when then (23) does not have a nontrivial solution in u (i.e., u = 0), which leads to a contradiction that ξ is an eigenvector. It suffices to show that for any δ satisfying the condition (24), the following inequality −δ H x Nevertheless, we can verify when δ < b a , (29) does not have any nontrivial solution for u, which gives a contradiction. Specifically, we show the following inequality is always satisfied under the condition on δ, which violates (29) for any u = 0, because (30) would imply that subsequently (29) implies u = 0, hence, ξ = 0, a contradiction. It remains to show that under δ < b a , (30) is satisfied. That is, the roots of (30) are given by 1 2 (a ± √ a 2 − 4b). We have shown This implies (30) is always satisfied and results in a contradiction. Hence, δ ≥ b a , which completes the proof.

F Analysis of RHM with conjugate gradient and trust-region update steps
We provide details of the implementation and convergence analysis of minimizing the Riemannian Hamiltonian with the Riemannian conjugate gradient and trust-region methods, i.e., we consider Algorithm 1 with the update step ξ(p t ) computed as conjugate gradient direction and trust-region step.

F.1 RHM with conjugate gradient (RHM-CG)
For the Riemannian conjugate gradient methods (CG), the step ξ(p t ) is given by a conjugate direction of the Riemannian gradient as ξ(p t ) = η t ζ(p t ) = η t − gradH(p t ) + β t Γ pt p t−1 ξ(p t−1 ) , where β t is the weight on the previous direction. There exists many choices of β t , which leads to different types of conjugate gradient direction. See [65,54] for a review. Given the conjugate gradient direction ξ(p t ) is a modified gradient direction, we can always ensure a linear convergence under the Riemannian PL condition provided the direction stays 'close' to the negative gradient direction, i.e., ξ(p t ) is a descent direction with ξ(p t ), −gradH(p t ) pt > 0. For practical implementation, such condition can always be safeguarded by choosing β t = 0 when ξ(p t ), −gradH(p t ) ≤ 0.
Theorem 6 (Linear convergence of RHM-CG). Under the same settings as in Theorem 1, consider Algorithm 1 with conjugate gradient direction ξ(p t ) where β t and η t are chosen such that ξ(p t ), −gradH(p t ) ≥ c gradH(p t ) 2 pt for some c > 0 and the Armijo condition (Definition 4) is satisfied. Letη = min i=0,...,t η i . Then, the iterates p t satisfies gradf (p t ) 2 pt ≤ (1−2r 1η cδ) t gradf (p 0 ) 2 p 0 . Proof. From the Armijo condition, we have for the stepsize η t , where the last inequality follows from the definition ofη and H(p t ) ≥ 0 for all p t . Applying the result recursively completes the proof.
We notice that the bound only requires a descent direction and a sufficient function decrease. Hence, we suspect a tighter bound exists when analyzing specific types of conjugate gradient.

F.2 RHM with trust-region (RHM-TR)
For the Riemannian trust-region (TR) method, the update step ξ(p t ) is computed by (approximately) solving the trust-region subproblem on the tangent space, i.e., ξ(p t ) = arg min where H t : T pt M − → T pt M is a self-adjoint linear operator that approximates the Hessian HessH(p t ). The model m pt (ξ) provides a second-order approximation to the pull-back function H pt := H • Exp pt . The step computed from solving (31) is tentative and accepted when the model provides a sufficient approximation accuracy. Let For some threshold ρ ∈ [0, 1/4), the step is only accepted when ρ t > ρ . In addition, the trust-region radius δ t is changed according to ρ t . To solve (31), truncated conjugate gradient is applied and with specific stopping condition, we can ensure local superlinear convergence. See [2,3,13] for more details.
Next, we show that RHM-TR achieves global linear convergence and local superlinear convergence to the global optimality. The former is guaranteed by a sufficient decrease in the model value and the latter is a corollary of the standard analysis in [2,3]. For simplicity of the analysis, we consider only a subsequence of iterates T where for all t ∈ T , ρ t > ρ and the trust-region step is accepted.
Proof. First from Assumption 2, gradH(p t ) pt = Hessf (p t )[gradf (p t )] pt ≤ L 1 L 0 and the operator norm of H t is bounded as H t pt ≤ H t − HessH(p t ) pt + HessH(p t ) pt ≤ L H L 0 L 1 + L.
Also, the trust-region direction ξ(p t ) returned by the truncated conjugate gradient method satisfies a so-called Cauchy decrease inequality [3, eq. (7.14)], which gives m pt (0) − m pt (ξ(p t )) ≥ 1 2 gradH(p t ) pt min ∆ t , gradH(p t ) pt H t pt ≥ 1 2 gradH(p t ) pt min c gradH(p t ) pt , gradH(p t ) pt H t pt where the second inequality follows from the definition of c and Assumption 2 where Furthermore, from the acceptance rule, Hence, the linear convergence is proved by recursively applying the result. The superlinear convergence simply follows from [3, Theorem 7.4.11] around any local minima.
Proof. It is sufficient to verify whether g(t) := c q (log det(Z + tV)) 2 + c l log det(Z + tV) log det(Y) where Z ∈ S d ++ , V is symmetric, and Z + tV ∈ S d ++ is convex in t. For this, we derive g (t) = 2c q log det(Z + tV) + c l log det(Y) where λ i , i = 1, ..., d are the eigenvalues of Z −1/2 VZ −1/2 . We see g (t) is not necessarily non-negative, and therefore, f is not convex in X in the Euclidean sense. Applying the same argument to Y completes the proof.
Proof of Proposition 8. The Riemannian gradient of f is derived as grad X f (X, Y) = c l log det(Y) + 2c q log det(X) X grad Y f (X, Y) = c l log det(X) − 2c q log det(Y) Y.

H.1 Optimality gap for geodesic quadratic bilinear optimization
We include additional convergence results in Fig. 4 on the optimality gap for the geodesic quadratic bilinear optimization problem in Section 7.1.

H.2 Results of RHM-SGD for orthonormal GAN
We show the sample collapse of RHM-SGD in Fig. 5.

H.3 trace-logarithm bilinear optimization
We consider the 'bilinear' example of [79] on the symmetric positive definite (SPD) manifold (endowed with the affine-invariant metric), i.e.,   (14) with d = 30, under different weights c q , c l . We observe that the RHM algorithms show a good rate of convergence in all the settings. In particular, RHM-SD-F and RHM-CON significantly outperforms RCEG in all the settings indicating better theoretical rates. For the experiment setting, we consider γ = 0.2 for RHM-CON and X 0 = Y 0 = I. The convergence results are shown in Fig. 6, where we notice that both RGDA and RCEG oscillate while all the RHM algorithms are convergent. RHM-CON and RHM-SD-F converge rapidly initially but subsequently have a slow rate of convergence due to the hardness of the problem. RHM-CG, on the other hand, has a faster rate of convergence.