Constrained optimization problems arise throughout machine learning, in classical settings such as dimension reduction [Boumal and Absil, 2011], dictionary learning [Sun et al., 2016a,b], and deep neural networks [Huang et al., 2018], but also in emerging problems involving decision-making and multi-agent interactions. While simple convex constraints (such as norm constraints) can be easily incorporated in standard optimization formulations, notably (proximal) gradient descent [Raskutti and Mukherjee, 2015, Giannou et al., 2021b,a, Antonakopoulos et al., 2020, Vlatakis-Gkaragkounis et al., 2020], in a range of other applications such as matrix recovery [Fornasier et al., 2011, Candes et al., 2008], lowrank matrix factorization [Han et al., 2021] and generative adversarial nets [Goodfellow et al., 2014], the constraints are fundamentally nonconvex and are often treated via special heuristics.
Thus, a general goal is to design algorithms that systematically take account of special geometric structure of the feasible set [Mei et al., 2021, Lojasiewicz, 1963, Polyak, 1963]. A long line of work in the machine learning (ML) community has focused on understanding the geometric properties of commonly used constraints and how they affect optimization; [see, e.g., Ge et al., 2015, Anandkumar and Ge, 2016, Sra and Hosseini, 2016, Jin et al., 2017, Ge et al., 2017, Du et al., 2017, Reddi et al., 2018, Criscitiello and Boumal, 2019, Jin et al., 2021]. A prominent aspect of this agenda has been the re-expression of these constraints through the lens of Riemannian manifolds. This has given rise to new algorithms [Sra andHosseini, 2015, Hosseini andSra, 2015] with a wide range of ML applications, inclduing online principal component analysis (PCA), the computation of Mahalanobis distance from noisy measurements [Bonnabel, 2013], consensus distributed algorithms for aggregation in ad-hoc wireless networks [Tron et al., 2012] and maximum likelihood estimation for certain non-Gaussian (heavyor light-tailed) distributions [Wiesel, 2012].
Going beyond simple minimization problems, the robustification of many ML tasks can be formulated as min-max optimization problems. Well-known examples in this domain include adversarial machine learning [Kumar et al., 2017, Chen et al., 2018], optimal transport [Lin et al., 2020a], and online learning [Mertikopoulos and Sandholm, 2018, Bomze et al., 2019, Antonakopoulos et al., 2020]. Similar to their minimization counterparts, non-convex constraints have been widely applicable to the minmax optimization as well [Heusel et al., 2017, Daskalakis and Panageas, 2018, Balduzzi et al., 2018, Mertikopoulos et al., 2019, Jin et al., 2020]. Recently there has been significant effort in proving tighter results either under more structured assumptions [Thekumprampil et al., 2019, Nouiehed et al., 2019, Lu et al., 2020, Azizian et al., 2020, Diakonikolas, 2020, Golowich et al., 2020, Lin et al., 2020c,b, Liu et al., 2021, Ostrovskii et al., 2021, Kong and Monteiro, 2021], and/or obtaining last-iterate convergence guarantees [Daskalakis and Panageas, 2018, 2019, Mertikopoulos et al., 2019, Adolphs et al., 2019, Liang and Stokes, 2019, Gidel et al., 2019, Mazumdar et al., 2020, Liu et al., 2020, Mokhtari et al., 2020, Lin et al., 2020c, Hamedani and Aybat, 2021, Abernethy et al., 2021, Cai et al., 2022] for computing min-max solutions in convex-concave settings. Nonetheless, the analysis of the iteration complexity in the general non-convex non-concave setting is still in its infancy [Vlatakis-Gkaragkounis et al., 2019, 2021]. In response, the optimization community has recently studied how to extend standard minmax optimization algorithms such as gradient descent ascent (GDA) and extragradient (EG) to the Riemannian setting. In mathematical terms, given two Riemannian manifolds M, N and a function f : M × N → R, the Riemannian min-max optimization (RMMO) problem becomes min x∈M max y∈N f (x, y).
The change of geometry from Euclidean to Riemannian poses several difficulties. Indeed, a fundamental stumbling block has been that this problem may not even have theoretically meaningful solutions. In contrast with minimization where an optimal solution in a bounded domain is always guaranteed [Fearnley et al., 2021], existence of such saddle points necessitates typically the application of topological fixed point theorems [Brouwer, 1911, Kakutani, 1941], KKM Theory [Knaster et al., 1929]). For the case of convex-concave f with compact sets X and Y, Sion [1958] generalized the celebrated theorem [Neumann, 1928] and guaranteed that a solution (x , y ) with the following property exists min x∈X f (x, y ) = f (x , y ) = max y∈Y f (x , y).
However, at the core of the proof of this result is an ingenuous application of Helly's lemma [Helly, 1923] for the sublevel sets of f , and, until the work of Ivanov [2014], it has been unclear how to formulate an analogous lemma for the Riemannian geometry. As a result, until recently have extensions of the min-max theorem been established, and only for restricted manifold families [Komiya, 1988, Kristály, 2014, Park, 2019]. Zhang et al. [2022] was the first to establish a min-max theorem for a flurry of Riemannian manifolds equipped with unique geodesics. Notice that this family is not a mathematical artifact since it encompasses many practical applications of RMMO, including Hadamard and Stiefel ones used in PCA [Lee et al., 2022]. Intuitively, the unique geodesic between two points of a manifold is the analogue of the a linear segment between two points in convex set: For any two points x 1 , x 2 ∈ X , their connecting geodesic is the unique shortest path contained in X that connects them.
Even when the RMMO is well defined, transferring the guarantees of traditional min-max optimization algorithms like Gradient Ascent Descent (GDA) and Extra-Gradient (EG) to the Riemannian case is non-trivial. Intuitively speaking, in the Euclidean realm the main leitmotif of the last-iterate analyses the aforementioned algorithms is a proof that δ t = x t -x * 2 is decreasing over time. To achieve this, typically the proof correlates δ t and δ t-1 via a "square expansion," namely:
.
(1.1) Notice, however that the above expression relies strongly on properties of Euclidean geometry (and the flatness of the corresponding line), namely that the the lines connecting the three points x t , x t-1 and x * form a triangle; indeed, it is the generalization of the Pythagorean theorem, known also as the law of cosines, for the induced triangle (ABC) := {(x t , x t-1 , x * )}. In a uniquely geodesic manifold such triangle may not belong to the manifold as discussed above. As a result, the difference of distances to the equilibrium using the geodesic paths d 2 M (x t , x * ) -d 2 M (x t-1 , x * ) generally cannot be given in a closed form. The manifold's curvature controls how close these paths are to forming a Euclidean triangle. In fact, the phenomenon of distance distortion, as it is typically called, was hypothesised by Zhang et al. [2022, Section 4.2] to be the cause of exponential slowdowns when applying EG to RMMO problems when compared to their Euclidean counterparts.
Multiple attempts have been made to bypass this hurdle. Huang et al. [2020] analyzed the Riemannian GDA (RGDA) for the non-convex non-concave setting. However, they do not present any last-iterate convergence results and, even in the average/best iterate setting, they only derive suboptimal rates for the geodesic convex-concave setting due to the lack of the machinery that convex analysis and optimization offers they derive sub-optimal rates for the geodesic convex-concave case, which is the problem of our interest. The analysis of Han et al. [2022] for Riemannian Hamiltonian Method (RHM), matches the rate of second-order methods in the Euclidean case. Although theoretically faster in terms of iterations, second-order methods are not preferred in practice since evaluating second order derivatives for optimization problems of thousands to millions of parameters quickly becomes prohibitive. Finally, Zhang et al. [2022] leveraged the standard averaging output trick in EG to derive a sublinear convergence rate of O(1/ ) for the general geodesically convex-concave Riemannian framework. In addition, they conjectured that the use of a different method could close the exponential gap for the geodesically strongly-convex-strongly-convex scenario and its Euclidean counterpart.
Given this background, a crucial question underlying the potential for successful application of first-order algorithms to Riemannian settings is the following:
Is a performance gap necessary between Riemannian and Euclidean optimal convex-concave algorithms in terms of accuracy and the condition number?
Our aim in this paper is to provide an extensive analysis of the Riemannian counterparts of Euclidean optimal first-order methods adapted to the manifold-constrained setting. For the case of the smooth objectives, we consider the Riemannian corrected extragradient (RCEG) method while for non-smooth cases, we analyze the textbook Riemannian gradient descent ascent (RGDA) method. Our main results are summarized in the following table.
For the definition of the acronyms, Det and Stoc stand for deterministic and stochastic, respectively. GSCSC and GCC stand for geodesically strongly-convex-strongly-concave (cf. Assumption 3.1 or Assumption D.1) and geodesically convex-concave (cf. Assumption 3.2 or Assumption D.2). Here ∈ (0, 1) is the accuracy, L, the Lipschitzness of the objective and its gradient, κ = /µ is the condition number of the function, where µ is the strong convexity parameter, (τ 0 , ξ 0 , ξ 0 ) are curvature parameters (cf. Assumption 2.1), and σ 2 is the variance of a Riemannian gradient estimator.
Our first main contribution is the derivation of a linear convergence rate for RCEG, answering the open conjecture of Zhang et al. [2022] about the performance gap of single-loop extragradient methods. Indeed, while a direct comparison between d 2 M (x t , x * ) and d 2 M (x t-1 , x * ) is infeasible, we are able to establish a relationship between the iterates via appeal to the duality gap function and obtain a contraction in terms of d 2 M (x t , x * ). In other words, the effect of Riemannian distance distortion is quantitative (the contraction ratio will depend on it) rather than qualitative (the geometric contraction still remains under a proper choice of constant stepsize). More specifically, we use d 2 M (x t , x )+d 2 N (y t , y ) and d 2 M (x t+1 , x ) + d 2 N (y t+1 , y ) to bound a gap function defined by f (x t , y ) -f (x , ŷt ). Since the objective function is geodesically strongly-convex-strongly-concave, we have f
). Then, using the relationship between (x t , y t ) and (x t , ŷt ), we conclude the desired results in Theorem 3.1. Notably, our approach is not affected by the nonlinear geometry of the manifold.
Secondly, we endeavor to give a systematic analysis of aspects of the objective function, including its smoothness, its convexity and oracle access. As we shall see, similar to the Euclidean case, better finite-time convergence guarantees are connected with a geodesic smoothness condition. For the sake of completeness, in the paper's supplement we present the performance of Riemannian GDA for the full spectrum of stochasticity for the non-smooth case. More specifically, for the stochastic setting, the key ingredient to get the optimal convergence rate is to carefully select the step size such that the noise of the gradient estimator will not affect the final convergence rate significantly. As a highlight, such technique has been used for analyzing stochastic RCEG in the Euclidean setting [Kotsalis et al., 2022] and our analysis can be seen as the extension to the Riemannian setting. For the nonsmooth setting, the analysis is relatively simpler compared to smooth settings but we still need to deal with the issue caused by the nonlinear geometry of manifolds and the interplay between the distortion of Riemannian metrics, the gap function and the bounds of Lipschitzness of our bi-objective. Interestingly, the rates we derive are near optimal in terms of accuracy and condition number of the objective, and analogous to their Euclidean counterparts.
We present the basic setup and optimality conditions for Riemannian min-max optimization. Indeed, we focus on some of key concepts that we need from Riemannian geometry, deferring a fuller presentation, including motivating examples and further discussion of related work, to Appendix A-C.
Riemannian geometry. An n-dimensional manifold M is a topological space where any point has a neighborhood that is homeomorphic to the n-dimensional Euclidean space. For each x ∈ M, each tangent vector is tangent to all parametrized curves passing through x and the tangent space T x M of a manifold M at this point is defined as the set of all tangent vectors. A Riemannian manifold M is a smooth manifold that is endowed with a smooth ("Riemannian") metric •, • x on the tangent space T x M for each point x ∈ M. The inner metric induces a norm • x on the tangent spaces.
A geodesic can be seen as the generalization of an Euclidean linear segment and is modeled as a smooth curve (map), γ : [0, 1] → M, which is locally a distance minimizer. Additionally, because of the non-flatness of a manifold a different relation between the angles and the lengths of an arbitrary geodesic triangle is induced. This distortion can be quantified via the sectional curvature parameter κ M thanks to Toponogov's theorem [Cheeger andEbin, 1975, Burago et al., 1992]. A constructive consequence of this definition are the trigonometric comparison inequalities (TCIs) that will be essential in our proofs; see Alimisis et al. [2020, Corollary 2.1] and Zhang and Sra [2016, Lemma 5] for detailed derivations. Assuming bounded sectional curvature, TCIs provide a tool for bounding Riemannian "inner products" that are more troublesome than classical Euclidean inner products.
The following proposition summarizes the TCIs that we will need; note that if κ min = κ max = 0 (i.e., Euclidean spaces), then the proposition reduces to the law of cosines.
Proposition 2.1 Suppose that M is a Riemannian manifold and let ∆ be a geodesic triangle in M with the side length a, b, c and let A be the angle between b and c. Then, we have 1. If κ M that is upper bounded by κ max > 0 and the diameter of M is bounded by π √ κmax , then
where ξ(κ, c) := 1 for κ ≤ 0 and ξ(κ, c)
Also, in contrast to the Euclidean case, x and v = grad x f (x) do not lie in the same space, since M and T x M respectively are distinct entities. The interplay between these dual spaces typically is carried out via the exponential maps. An exponential map at a point x ∈ M is a mapping from the tangent space T x M to M. In particular, y := Exp x (v) ∈ M is defined such that there exists a geodesic γ : [0, 1] → M satisfying γ(0) = x, γ(1) = y and γ (0) = v. The inverse map exists since the manifold has a unique geodesic between any two points, which we denote as Exp -foot_0
x : M → T x M. Accordingly, we have d M (x, y) = Exp -1
x (y) x is the Riemannian distance induced by the exponential map. Finally, in contrast again to Euclidean spaces, we cannot compare the tangent vectors at different points x, y ∈ M since these vectors lie in different tangent spaces. To resolve this issue, it suffices to define a transport mapping that moves a tangent vector along the geodesics and also preserves the length and Riemannian metric •, • x ; indeed, we can define a parallel transport Γ y
x : T x M → T y M such that the inner product between any u, v ∈ T x M is preserved; i.e., u, v x = Γ y x (u), Γ y x (v) y .
Riemannian min-max optimization and function classes. We let M and N be Riemannian manifolds with unique geodesic and bounded sectional curvature and assume that the function f : M × N → R is defined on the product of these manifolds. The regularity conditions that we impose on the function f are as follows.
Definition 2.1 A function f : M×N → R is geodesically L-Lipschitz if for ∀x, x ∈ M and ∀y, y ∈ N , the following statement holds true:
Additionally, if function f is also differentiable, it is called geodesically -smooth if for ∀x, x ∈ M and ∀y, y ∈ N , the following statement holds true,
where (grad x f (x , y ), grad y f (x , y )) ∈ T x M × T y N is the Riemannian gradient of f at (x , y ), Γ x x is the parallel transport of M from x to x, and Γ y y is the parallel transport of N from y to y.
Definition 2.2 A function f : M × N → R is geodesically strongly-convex-strongly-concave with the modulus µ > 0 if the following statement holds true,
where (subgrad x f (x , y ), subgrad y f (x , y )) ∈ T x M × T y N is a Riemannian subgradient of f at a point (x , y ). A function f is geodesically convex-concave if the above holds true with µ = 0.
2. M, N admit unique geodesic paths for any (x, y), (x , y ) ∈ M × N .
3. The sectional curvatures of M and N are both bounded in the range [κ min , κ max ] with κ min ≤ 0.
If κ max > 0, we assume that the diameter of manifolds is bounded by π √ κmax .
Under these conditions, Zhang et al. [2022] proved an analog of Sion's minimax theorem [Sion, 1958] in geodesic metric spaces. Formally, we have
which guarantees that there exists at least one global saddle point (x , y ) ∈ M × N such that min x∈M f (x, y ) = f (x , y ) = max y∈Y f (x , y). Note that the unicity of geodesics assumption is algorithm-independent and is imposed for guaranteeing that a saddle-point solution always exist. Even though this rules out many manifolds of interest, there are still many manifolds that satisfy such conditions. More specifically, the Hadamard manifold (manifolds with non-positive curvature, κ max = 0) has a unique geodesic between any two points. This also becomes a common regularity condition in Riemannian optimization [Zhang andSra, 2016, Alimisis et al., 2020]. For any point (x, ŷ) ∈ M × N , the duality gap f (x, y ) -f (x , ŷ) thus gives an optimality criterion.
In the setting where f is geodesically strongly-convex-strongly-concave with µ > 0, it is not difficult to verify the uniqueness of a global saddle point (x , y ) ∈ M × N . Then, we can consider the distance gap (d(x, x )) 2 + (d(ŷ, y )) 2 as an optimality criterion for any point (x, ŷ) ∈ M × N .
If f is also geodesically -smooth, we denote κ = µ as the condition number.
Given the above definitions, we can ask whether it is possible to find an -saddle point efficiently or not. In this context, Zhang et al. [2022] have answered this question in the affirmative for the setting where f is geodesically -smooth and geodesically convex-concave; indeed, they derive the convergence rate of Riemannian corrected extragradient (RCEG) method in terms of time-average iterates and also conjecture that RCEG does not guarantee convergence at a linear rate in terms of last iterates when f is geodesically -smooth and geodesically strongly-convex-strongly-concave, due to the existence of distance distortion; see Zhang et al. [2022, Section 4.2]. Surprisingly, we show in Section 3 that RCEG with constant stepsize can achieve last-iterate convergence at a linear rate. Moreover, we establish the optimal convergence rates of stochastic RCEG for certain choices of stepsize for both geodesically convex-concave and geodesically strongly-convex-strongly-concave settings.
In this section, we revisit the scheme of Riemannian corrected extragradient (RCEG) method proposed by Zhang et al. [2022] and extend it to a stochastic algorithm that we refer to as stochastic RCEG. We present our main results on an optimal last-iterate convergence guarantee for the geodesically stronglyconvex-strongly-concave setting (both deterministic and stochastic) and a time-average convergence guarantee for the geodesically convex-concave setting (stochastic). This complements the time-average convergence guarantee for geodesically convex-concave setting (deterministic) [Zhang et al., 2022, Theorem 4.1] and resolves an open problem posted in Zhang et al. [2022, Section 4.2].
Input: initial points (x 0 , y 0 ) and stepsizes η > 0.
Algorithm 2 SRCEG Input: initial points (x 0 , y 0 ) and stepsizes η > 0. for t = 0, 1, 2, . . . , T -1 do Query (g t x , g t y ) as a noisy estimator of Riemannian gradient of f at a point (x t , y t ).
xt
x , ĝt y ) as a noisy estimator of Riemannian gradient of f at a point (x t , ŷt ).
The recently proposed Riemannian corrected extragradient (RCEG) method [Zhang et al., 2022] is a natural extension of the celebrated extragradient (EG) method to the Riemannian setting. Its scheme resembles that of EG in Euclidean spaces but employs a simple modification in the extrapolation step to accommodate the nonlinear geometry of Riemannian manifolds. Let us provide some intuition how such modifications work.
We start with a basic version of EG as follows, where M and N are classically restricted to be convex constraint sets in Euclidean spaces:
Turning to the setting where M and N are Riemannian manifolds, the rather straightforward way to do the generalization is to replace the projection operator by the corresponding exponential map and the gradient by the corresponding Riemannian gradient. For the first line of Eq. (3.1), this approach works and leads to the following updates:
However, we encounter some issues for the second line of Eq. (3.1): The aforementioned approach leads to some problematic updates,
indeed, the exponential maps Exp xt (•) and Exp yt (•) are defined from T xt M to M and from T yt N to N respectively. However, we have -grad x f (x t , ŷt ) ∈ T xt M and grad y f (x t , ŷt ) ∈ T ŷt N . This motivates us to reformulate the second line of Eq. (3.1) as follows:
In the general setting of Riemannian manifolds, the terms x t -xt and y t -ŷt become Exp -1 xt (x t ) ∈ T xt M and Exp -1 ŷt (y t ) ∈ T ŷt N . This observation yields the following updates:
We summarize the resulting RCEG method in Algorithm 1 and present the stochastic extension with noisy estimators of Riemannian gradients of f in Algorithm 2.
We present our main results on global convergence for Algorithms 1 and 2. To simplify the presentation, we treat separately the following two cases:
Assumption 3.1 The objective function f is geodesically -smooth and geodesically strongly-convexstrongly-concave with µ > 0.
Assumption 3.2 The objective function f is geodesically -smooth and geodesically convex-concave.
Letting (x , y ) ∈ M × N be a global saddle point of f (which exists under either Assumption 3.1 or 3.2), we let D 0 = (d M (x 0 , x )) 2 + (d N (y 0 , y )) 2 > 0 and κ = /µ for geodesically strongly-convexstrongly-concave setting. For simplicity of presentation, we also define a ratio τ (•, •) that measures how non-flatness changes in the spaces: τ ([κ min , κ max ], c) = ξ(κ min ,c) ξ(κmax,c) ≥ 1. We summarize our results for Algorithm 1 in the following theorem.
Theorem 3.1 Given Assumptions 2.1 and 3.1, and letting η = min{1/(2 √ τ 0 ), ξ 0 /(2µ)}, there exists some T > 0 such that the output of Algorithm 1 satisfies that
an -saddle point of f in Definition 2.4) and the total number of Riemannian gradient evaluations is bounded by
where τ 0 = τ ([κ min , κ max ], D) ≥ 1 measures how non-flatness changes in M and N and ξ 0 = ξ(κ max , D) ≤ 1 is properly defined in Proposition 2.1.
Remark 3.1 Theorem 3.1 illustrates the last-iterate convergence of Algorithm 1 for solving geodesically strongly-convex-strongly-concave problems, thereby resolving an open problem delineated by Zhang et al. [2022]. Further, the dependence on κ and 1/ cannot be improved since it matches the lower bound established for min-max optimization problems in Euclidean spaces [Zhang et al., 2021]. However, we believe that the dependence on τ 0 and ξ 0 is not tight, and it is of interest to either improve the rate or establish a lower bound for general Riemannian min-max optimization.
Remark 3.2 The current theoretical analysis covers local geodesic strong-convex-strong-concave settings. The key ingredient is how to define the local region; indeed, if we say the set of {(x, y) :
} must be contained in the above local region and the objective function is also geodesic strong
, our theoretical analysis guarantees the last-iterate linear convergence rate. Such argument and definition of local region were standard for min-max optimization in the Euclidean setting; see Liang and Stokes [2019, Assumption 2.1]. For an important optimization problem that is globally geodesically strongly-convex-strongly-concave, we refer to Appendix B where Robust matrix Karcher mean problem is indeed the desired one.
In the scheme of SRECG, we highlight that (g t x , g t y ) and (ĝ t x , ĝt y ) are noisy estimators of Riemannian gradients of f at (x t , y t ) and (x t , ŷt ). It is necessary to impose the conditions such that these estimators are unbiased and has bounded variance. By abuse of notation, we assume that
where the noises (ξ t x , ξ t y ) and ( ξt x , ξt y ) are independent and satisfy that
We are ready to summarize our results for Algorithm 2 in the following theorems.
Theorem 3.2 Given Assumptions 2.1 and 3.1, letting Eq. (3.2) and Eq. (3.3) hold with σ > 0 and letting η > 0 satisfy η = min{ 1 24
µT }, there exists some T > 0 such that the output of Algorithm 2 satisfies that E[(d(x T , x )) 2 + (d(y T , y )) 2 ] ≤ and the total number of noisy Riemannian gradient evaluations is bounded by
where τ 0 = τ ([κ min , κ max ], D) ≥ 1 measures how non-flatness changes in M and N and ξ 0 = ξ(κ max , D) ≤ 1 is properly defined in Proposition 2.1.
Theorem 3.3 Given Assumptions 2.1 and 3.2 and assume that Eq. (3.2) and Eq. (3.3) hold with σ > 0 and let η > 0 satisfies that η = min{ 1 4 √ τ 0 , 1 σ D 0 ξ 0 T }, there exists some T > 0 such that the output of Algorithm 2 satisfies that E[f (x T , y ) -f (x , ȳT )] ≤ and the total number of noisy Riemannian gradient evaluations is bounded by
, where τ 0 = τ ([κ min , κ max ], D) measures how non-flatness changes in M and N and ξ 0 = ξ(κ min , D) ≥ 1 is properly defined in Proposition 2.1. The time-average iterates (x T , ȳT ) ∈ M × N can be computed by (x 0 , ȳ0 ) = (0, 0) and the inductive formula:
Remark 3.3 Theorem 3.2 presents the last-iterate convergence rate of Algorithm 2 for solving geodesically strongly-convex-strongly-concave problems while Theorem 3.3 gives the time-average convergence rate when the function f is only assumed to be geodesically convex-concave. Note that we carefully choose the stepsizes such that our upper bounds match the lower bounds established for stochastic minmax optimization problems in Euclidean spaces [Juditsky et al., 2011, Fallah et al., 2020, Kotsalis et al., 2022], in terms of the dependence on κ, 1/ and σ 2 , up to log factors.
The last-iterate linear convergence rate in terms of Riemannian metrics is only limited to geodesically strongly convex-concave cases but other results, e.g., the average-iterate sublinear convergence rate, are derived under more mild conditions. This is consistent with classical results in the Euclidean setting where geodesic convexity reduces to convexity; indeed, the last-iterate linear convergence rate in terms of squared Euclidean norm is only obtained for strongly convex-concave cases. As such, our setting is not restrictive. Moreover, Zhang et al. [2022] showed that the existence of a global saddle point is only guaranteed under the geodesically convex-concave assumption.
We present numerical experiments on the task of robust principal component analysis (RPCA) for symmetric positive definite (SPD) matrices. In particular, we compare the performance of Algorithm 1 and 2 with different outputs, i.e., the last iterate (x T , y T ) versus the time-average iterate (x T , ȳT ) (see the precise definition in Theorem 3.3). Note that our implementations of both algorithms are based on the manopt package [Boumal et al., 2014]. All the experiments were implemented in MATLAB R2021b on a workstation with a 2.6 GHz Intel Core i7 and 16GB of memory. Due to space limitations, some additional experimental results are deferred to Appendix G.
Experimental setup. The problem of RPCA [Candès et al., 2011, Harandi et al., 2017] can be formulated as the Riemannian min-max optimization problem with an SPD manifold and a sphere manifold. Formally, we have max
In this formulation, α > 0 denotes the penalty parameter,
As demonstrated by Zhang et al. [2022], the problem of RPCA is nonconvex-nonconcave from a Euclidean perspective but is locally geodesically strongly-convexstrongly-concave and satisfies most of the assumptions that we make in this paper. In particular, the SPD manifold is complete with sectional curvature in [-1 2 , 1] [ Criscitiello and Boumal, 2022] and the sphere manifold is complete with sectional curvature of 1. Other reasons why we use such example are: (i) it is a classical one in ML; (ii) Zhang et al. [2022] also uses this example and observes the linear convergence behavior; (iii) the numerical results show that the unicity of geodesics assumption may not be necessary in practice; and (iv) this is an application where both min and max sides are done on Riemannian manifolds.
Following the previous works of Zhang et al. [2022] and Han et al. [2022], we generate a sequence of data matrices M i satisfying that their eigenvalues are in the range of [0.2, 4.5]. In our experiment, we fix α = 1.0 and also vary the problem dimension d ∈ {25, 50, 100}. The evaluation metric is set as gradient norm. We set n = 40 and n = 200 in Figure 1 and 2. For RCEG, we set η = 1 2 where > 0 is selected via grid search. For SRCEG, we set η t = min{ 1 2 , a t } where , a > 0 are selected via grid search. Additional results on the effect of stepsize are summarized in Appendix G.
Experimental results. Figure 1 summarizes the effects of different outputs for RCEG; indeed, RCEG-last and RCEG-avg refer to Algorithm 1 with last iterate and time-average iterate respectively. It is clear that the last iterate of RCEG consistently exhibits linear convergence to an optimal solution in all the settings, verifying our theoretical results in Theorem 3.1. In contrast, the average iterate of RCEG converges much slower than the last iterate of RCEG. The possible reason is that the problem of RPCA is only locally geodesically strongly-convex-strongly-concave and averaging with the iterates generated during early stage will significantly slow down the convergence of RCEG.
Figure 2 presents the comparison between SRCEG (with either last iterate or time-average iterate) and RCEG with last-iterate; here, SRCEG-last and SRCEG-avg refer to Algorithm 2 with last iterate and time-average iterate respectively. We observe that SRCEG with either last iterate or average iterate converge faster than RCEG at the early stage and all of them finally converge to an optimal solution. This demonstrates the effectiveness and efficiency of SRCEG in practice. It is also worth mentioning that the difference between last-iterate convergence and time-average-iterate convergence is not as significant as in the deterministic setting. This is possibly because the technique of averaging help cancels the negative effect of imperfect information [Kingma andBa, 2015, Yazıcı et al., 2019].
Inspired broadly by the structure of the complex competition that arises in many applications of robust optimization in ML, we focus on the problem of min-max optimization in the pure Riemannian setting (where both min and max player are constrained in a smooth manifold). Answering the open question of Zhang et al. [2022] for the geodesically (strongly) convex-concave case, we showed that the Riemannian correction technique for EG matches the linear last-iterate complexity of their Euclidean counterparts in terms of accuracy and conditional number of objective for both deterministic and stochastic case. Additionally, we provide near-optimal guarantees for both smooth and non-smooth min-max optimization via Riemannian EG and GDA for the simple convex-concave case.
As a consequence of this work numerous open problems emerge; one immediate open question for future work is to explore whether the dependence on the curvature constant is also tight. Additionally, another generalization of interest would be to consider the performance of RCEG in the case of Riemannian Monotone Variational inequalities (RMVI) and examine the generalization of Zhang et al. [2022] existence proof. Finally, there has been recent work in proving last-iterate convergence in the convexconcave setting via Sum-Of-Squares techniques [Cai et al., 2022]. It would be interesting to examine how one could leverage this machinery in a non-Euclidean but geodesic-metric-friendly framework.
The literature for the geometric properties of Riemannian Manifolds is immense and hence we cannot hope to survey them here; for an appetizer, we refer the reader to Burago et al. [2001] and Lee [2012] and references therein. On the other hand, as stated, it is not until recently that the long-run nonasymptotic behavior of optimization algorithms in Riemannian manifolds (even the smooth ones) has encountered a lot of interest. For concision, we have deferred here a detailed exposition of the rest of recent results to Appendix A of the paper's supplement. Additionally, in Appendix B we also give a bunch of motivating examples which can be solved by Riemannian min-max optimization.
Minimization on Riemannian manifolds. Many application problems can be formulated as the minimization or maximization of a smooth function over Riemannian manifold and has triggered a line of research on the extension of the classical first-order and second-order methods to Riemannian setting with asymptotic convergence to first-order stationary points in general [Absil et al., 2009]. Recent years have witnessed the renewed interests on nonasymptotic convergence analysis of solution methods. In particular, Boumal et al. [2019] proved the global sublinear convergence results for Riemannian gradient descent method and Riemannian trust region method, and further demonstrated that the Riemannian trust region method converges to a second-order stationary point in polynomial time; see also similar results in some other works [Kasai and Mishra, 2018, Hu et al., 2018, 2019]. We are also aware of recent works on problem-specific methods [Wen and Yin, 2013, Gao et al., 2018, Liu et al., 2019] and primal-dual methods [Zhang et al., 2020].
Compared to the smooth counterpart, Riemannian nonsmooth optimization is harder and relatively less explored [Absil and Hosseini, 2019]. A few existing works focus on optimizing geodesically convex functions over Riemannian manifold with subgradient methods [Ferreira and Oliveira, 1998, Zhang and Sra, 2016, Bento et al., 2017]. In particular, Ferreira and Oliveira [1998] provided the first asymptotic convergence result while Zhang and Sra [2016] and [Bento et al., 2017] proved an nonasymptotic global convergence rate of O( -foot_2 ) for Riemannian subgradient methods. Further, Ferreira and Oliveira [2002] assumed that the proximal mapping over Riemannian manifold is computationally tractable and proved the global sublinear convergence of Riemannian proximal point method. Focusing on optimization over Stiefel manifold, Chen et al. [2020] studied the composite objective function and proposed Riemannian proximal gradient method which only needs to compute the proximal mapping of nonsmooth component function over the tangent space of Stiefel manifold. Li et al. [2021] consider optimizing a weakly convex function over Stiefel manifold and proposed Riemannian subgradient methods that drive a near-optimal stationarity measure below within the number of iterations bounded by O( -4 ).
There are some results on stochastic optimization over Riemannian manifold. In particular, Bonnabel [2013] proved the first asymptotic convergence result for Riemannian stochastic gradient descent, which is extended by a line of subsequent works [Zhang et al., 2016, Tripuraneni et al., 2018, Becigneul and Ganea, 2019, Kasai et al., 2019]. If the Riemannian Hessian is not positive definite, some recent works have suggested frameworks to escape saddle points [Sun et al., 2019, Criscitiello andBoumal, 2019].
Min-Max optimization in Euclidean spaces. Focusing on solving specifically min-max problems, the algorithms under euclidean geometry have a very rich history in optimization that goes back at least to the original proximal point algorithms [Martinet, 1970, Rockafellar, 1976] for variational inequality (VI) problems; At a high level, if the objective function is Lipschitz and strictly convex-concave, the simple forward-backward schemes are known to converge -and if combined with a Polyak-Ruppert averaging scheme [Ruppert, 1988, Polyak and Juditsky, 1992, Nemirovski et al., 2009], they achieve an O(1/ 2 ) complexity 2 without the caveat of strictness [Bauschke and Combettes, 2011]. If, in addition, the objective admits Lipschitz continuous gradients, then the extragradient (EG) algorithm [Korpelevich, 1976] achieves trajectory convergence without strict monotonicity requirements, while the time-average iterate converges at O(1/ ) steps [Nemirovski, 2004]. Finally, if the problem is strongly convex-concave, forward-backward methods computes an -saddle point at O(1/ ) steps; and if the operator is also Lipschitz continuous, classical results in operator theory show that simple forward-backward methods suffice to achieve a linear convergence rate [Facchinei andPang, 2007, Bauschke andCombettes, 2011].
Min-Max optimization on Riemannian manifolds. In the case of nonlinear geometry, the literature has been devoted on two different orthogonal axes: a) the existence of saddle point for min-max objective bi-functions and b) the design of algorithms for the computation of such points. For the existence of saddle point, a long line of recent work tried to generalize the seminal minima theorem for quasi-convex-quasi-concave problems of Sion [1958]. The crucial bottleneck of this generalization to Riemannian smooth manifolds had been the application of both Knaster-Kuratowski-Mazurkiewicz (KKM) theorem and Helly's theorem in non-flat spaces. Before Zhang et al. [2022], the existence of saddle points had been identified for the special case of Hadamard manifolds [Komiya, 1988, Kristály, 2014, Bento et al., 2017, Park, 2019].
Similar with the existence results, initially the developed methods referred to the computation of singularities in monotone variational operators typically in hyperbolic Hadamard manifolds with negative curvature [Li et al., 2009]. More recently, Huang et al. [2020] proposed a Riemannian gradient descent ascent method (RGDA), yet the analysis is restricted to N being a convex subset of the Euclidean space and f (x, y) being strongly concave in y. It is worth mentioning that for the case Hadamard and generally hyperbolic manifolds, extra-gradient style algorithms have been proposed [Wang et al., 2010, Ferreira et al., 2005] in the literature, establishing mainly their asymptotic convergence. However it was not until recent Zhang et al. [2022] that the riemannian correction trick has been analyzed for the case of the extra-gradient algorithm. Bearing in our mind the higher-order methods, Han et al. [2022] has recently proposed the Riemannian Hamiltonian Descent and versions of Newton's method for for geodesic convex geodesic concave functions. Since in this work, we focus only on first-order methods, we don't compare with the aforementioned Hamiltonian alternative since it incorporates always the extra computational burden of second-derivatives and hessian over a manifold.
We provide some examples of Riemannian min-max optimization to give a sense of their expressivity. Two of the examples are the generic models from the optimization literature [Ben-Tal et al., 2009, Absil et al., 2009, Hu et al., 2020] and the two others are the formulations of application problems arising from machine learning and data analytics [Pennec et al., 2006, Fletcher and Joshi, 2007, Lin et al., 2020a].
We can consider a rather straightforward generalization of constrained optimization problem from Euclidean spaces to Riemannian manifolds [Bergmann and Herzog, 2019]. This formulation finds a wide range of real-world applications, e.g., non-negative principle component analysis, weighted max-cut and so on. Letting M be a finitedimensional Riemannian manifold with unique geodesic, we focus on the following problem:
where g := (g 1 , g 2 , . . . , g m ) : M → R m and h := (h 1 , h 2 , . . . , h n ) : M → R n are two mappings. Then, we can introduce the dual variables λ and µ and reformulate the aforementioned constrained optimization problem as follows, min
Suppose that f and all of g i and h i are geodesically convex and smooth, the above problem is a geodesicconvex-Euclidean-concave min-max optimization problem.
Example B.2 (Distributionally robust Riemannian optimization) Distributionally robust optimization (DRO) is an effective method to deal with the noisy data, adversarial data, and imbalanced data. We consider the problem of DRO over Riemannian manifold; indeed, given a set of data samples {ξ i } N i=1 , the problem of DRO over Riemannian manifold M can be written in the form of
where p = (p 1 , p 2 , . . . , p N ) and S = {p ∈ R N : N i=1 p i = 1, p i ≥ 0}. In general, (x; ξ i ) denotes the loss function over Riemannian manifold M. If is geodesically convex and smooth, the above problem is a geodesic-convex-Euclidean-concave min-max optimization problem.
We consider a robust version of classical matrix Karcher mean problem. More specifically, the Karcher mean of N symmetric positive definite matrices {A i } N i=1 is defined as the matrix X ∈ M = {X ∈ R n×n : X 0, X = X } that minimizes the sum of squared distance induced by the Riemannian metric:
The loss function is thus defined by
which is known to be nonconvex in Euclidean spaces but geodesically strongly convex. Then, the robust version of classical matrix Karcher mean problem is aiming at solving the following problem:
where γ > 0 stands for the trade-off between the computation of Karcher mean over a set of {Y i } N i=1 and the difference between the observed samples {A i } N i=1 and {Y i } N i=1 . It is clear that the above problem is a geodesically strongly-convex-strongly-concave min-max optimization problem.
Example B.4 (Projection robust optimal transport problem) We consider the projection robust optimal transport (OT) problem -a robust variant of the OT problem -that achieves superior sample complexity bound [Lin et al., 2021]. Let {x 1 , x 2 , . . . , x n } ⊆ R d and {y 1 , y 2 , . . . , y n } ⊆ R d denote sets of n atoms, and let (r 1 , r 2 , . . . , r n ) and (c 1 , c 2 , . . . , c n ) denote weight vectors. We define discrete probability measures µ = n i=1 r i δ x i and ν = n j=1 c j δ y j . In this setting, the computation of the kdimensional projection robust OT distance between µ and ν resorts to solving the following problem:
It is worth mentioning that the above problem is a geodesically-nonconvex-Euclidean-concave min-max optimization problem with special structures, making the computation of stationary points tractable. While the global convergence guarantee for our algorithm does not apply, the above problem might be locally geodesically-convex-Euclidean-concave such that our algorithm with sufficiently good initialization works here.
In addition to these examples, it is worth mentioning that Riemannian min-max optimization problems contain all general min-max optimization problems in Euclidean spaces and all Riemannian minimization or maximization optimization problems. It is also an abstraction of many machine learning problems, e.g,. principle component analysis [Boumal and Absil, 2011], dictionary learning [Sun et al., 2016a,b], deep neural networks (DNNs) [Huang et al., 2018] and low-rank matrix learning [Vandereycken, 2013, Jawanpuria andMishra, 2018]; indeed, the problem of principle component analysis resorts to optimization problems on Grassmann manifolds for example.
To generalize the first-order methods in Euclidean setting, we introduce several basic concepts in metric geometry [Burago et al., 2001], which are known to include both Euclidean spaces and Riemannian manifolds as special cases. Formally, we have Definition C.1 (Metric Space) A metric space (X, d) is a pair of a set X and a distance function d(•, •) satisfying: (i) d(x, x ) ≥ 0 for any x, x ∈ X; (ii) d(x, x ) = d(x , x) for any x, x ∈ X; and (iii) d(x, x ) ≤ d(x, x ) + d(x , x ) for any x, x , x ∈ X. In other words, the distance function d(•, •) is non-negative, symmetrical and satisfies the triangle inequality.
A path γ : [0, 1] → X is a continuous mapping from the interval [0, 1] to X and the length of γ is defined as length(γ) := lim n→+∞ sup 0=t 0 <...<tn=1 n i=1 d(γ(t i-1 ), γ(t i )). Note that the triangle inequality implies that sup 0=t 0 <...<tn=1 n i=1 d(γ(t i-1 ), γ(t i )) is nondecreasing. Then, the length of a path γ is well defined since the limit is either +∞ or a finite scalar. Moreover, for ∀ > 0, there exists n ∈ N and the partition 0 = t 0 < . . . < t n = 1 of the interval [0, 1] such that length(γ) ≤ n i=1 d(γ(t i-1 ), γ(t i )) + .
Definition C.2 (Length Space) A metric space (X, d) is a length space if, for any x, x ∈ X and > 0, there exists a path γ : [0, 1] → X connecting x and x such that length(γ) ≤ d(x, x ) + .
We can see from Definition C.2 that a set of length spaces is strict subclass of metric spaces; indeed, for some x, x ∈ X, there does not exist a path γ such that its length can be approximated by d(x, x ) for some tolerance > 0. In metric geometry, a geodesic is a path which is locally a distance minimizer everywhere. More precisely, a path γ is a geodesic if there is a constant ν > 0 such that for any t ∈ [0, 1] there is a neighborhood I of [0, 1] such that, d(γ(t 1 ), γ(t 2 )) = ν|t 1 -t 2 |, for any t 1 , t 2 ∈ I.
Note that the above generalizes the notion of geodesic for Riemannian manifolds. Then, we are ready to introduce the geodesic space and uniquely geodesic space [Bacak, 2014]. d) is a geodesic space if, for any x, x ∈ X, there exists a geodesic γ : [0, 1] → X connecting x and x . Furthermore, it is called uniquely geodesic if the geodesic connecting x and x is unique for any x, x ∈ X.
Trigonometric geometry in nonlinear spaces is intrinsically different from Euclidean space. In particular, we remark that the law of cosines in Euclidean space (with • as 2 -norm) is crucial for analyzing the convergence property of optimization algorithms, e.g.,
where a, b, c are sides of a geodesic triangle in Euclidean space and A is the angle between b and c. However, such nice property does not hold for nonlinear spaces due to the lack of flat geometry, further motivating us to extend the law of cosines under nonlinear trigonometric geometry. That is to say, given a geodesic triangle in X with sides a, b, c where A is the angle between b and c, we hope to establish the relationship between a 2 , b 2 , c 2 and 2bc cos(A) in nonlinear spaces; see the main context for the comparing inequalities. Finally, we specify the definition of section curvature of Riemannian manifolds and clarify how such quantity affects the trigonometric comparison inequalities. More specifically, the sectional curvature is defined as the Gauss curvature of a 2-dimensional sub-manifold that are obtained from the image of a two-dimensional subspace of a tangent space after exponential mapping. It is worth mentioning that the above 2-dimensional sub-manifold is locally isometric to a 2-dimensional sphere, a Euclidean plane, and a hyperbolic plane with the same Gauss curvature if its sectional curvature is positive, zero and negative respectively. Then we are ready to summarize the existing trigonometric comparison inequalities for Riemannian manifold with bounded sectional curvatures. Note that the following two propositions are the full version of Proposition 2.1 and will be used in our subsequent proofs.
Proposition C.1 Suppose that M is a Riemannian manifold with sectional curvature that is upper bounded by κ max and let ∆ be a geodesic triangle in M with the side length a, b, c and A which is the angle between b and c. If κ max > 0, we assume the diameter of M is bounded by π √ κmax . Then, we have
where ξ(κ, c) := 1 for κ ≤ 0 and ξ(κ, c)
Proposition C.2 Suppose that M is a Riemannian manifold with sectional curvature that is lower bounded by κ min and let ∆ be a geodesic triangle in M with the side length a, b, c and A which is the angle between b and c. Then, we have
where ξ(κ, c)
Remark C.1 Proposition C.1 and C.2 are simply the restatement of Alimisis et al. [2020, Corollary 2.1] and Zhang and Sra [2016, Lemma 5]. The former inequality is obtained when the sectional curvature is bounded from above while the latter inequality characterizes the relationship between the trigonometric lengths when the sectional curvature is bounded from below. If κ min = κ max = 0 (i.e., Euclidean spaces), we have ξ(κ min , c) = ξ(κ max , c) = 1. The proof is based on Toponogov's theorem and Riccati comparison estimate [Petersen, 2006, Proposition 25] and we refer the interested readers to Zhang and Sra [2016] and Alimisis et al. [2020] for the details.
In this section, we propose and analyze Riemannian gradient descent ascent (RGDA) method for nonsmooth Riemannian min-max optimization and extend it to stochastic RGDA. We present our results on the optimal last-iterate convergence guarantee for geodesically strongly-convex-strongly-concave setting (both deterministic and stochastic) and time-average convergence guarantee for geodesically convexconcave setting (both deterministic and stochastic).
Compared to Riemannian corrected extragradient (RCEG) method, our Riemannian gradient descent ascent (RGDA) method is a relatively straightforward generalization of GDA in Euclidean spaces. More specifically, we start with the scheme of GDA as follows (just consider M and N as convex constraint sets in Euclidean spaces),
Algorithm 3 RGDA Input: initial points (x 0 , y 0 ) and stepsizes η t > 0. for t = 0, 1, 2, . . . , T -1 do Query (g t x , g t y ) ← (subgrad x f (x t , y t ), subgrad y f (x t , y t )) as Riemannian subgradient of f at a point (x t , y t ).
x t+1 ← Exp xt (-η t • g t x ). y t+1 ← Exp yt (η t • g t y ). end for Algorithm 4 SRGDA Input: initial points (x 0 , y 0 ) and stepsizes η t > 0. for t = 0, 1, 2, . . . , T -1 do Query (g t x , g t y ) as a noisy estimator of Riemannian subgradient of f at a point (x t , y t ).
) is one subgradient of f . By replacing the projection operator by the corresponding exponential map and the gradient by the corresponding Riemannian gradient, we have
. where (g t x , g t y ) ← (subgrad x f (x t , y t ), subgrad y f (x t , y t )) is one Riemannian subgradient of f . Then, we summarize the resulting scheme of RGDA method in Algorithm 3 and its stochastic extension with noisy estimators of Riemannian gradients of f in Algorithm 4.
We present our main results on the global convergence rate estimation for Algorithm 3 and 4 in terms of Riemannian gradient and noisy Riemannian gradient evaluations. The following assumptions are made throughout for geodesically strongly-convex-strongly-concave and geodesically convex-concave settings. Imposing the geodesically Lipschitzness condition is crucial to achieve finite-time convergence guarantee if we do not assume the geodesically smoothness condition. Note that we only require the lower bound for the sectional curvatures of manifolds and this is weaker than that presented in the main context. Letting (x , y ) ∈ M × N be a global saddle point of f (it exists under either Assumption D.1 or D.2), we let D 0 = (d M (x 0 , x )) 2 + (d N (y 0 , y )) 2 > 0 and summarize our results for Algorithm 3 in the following theorems.
Theorem D.1 Under Assumption D.1 and let η t > 0 satisfies that η t = 1 µ min{1, 2 t }. There exists some T > 0 such that the output of Algorithm 3 satisfies that (d(x T , x )) 2 + (d(y T , y )) 2 ≤ and the total number of Riemannian subgradient evaluations is bounded by
where ξ 0 = ξ(κ min , D) measures the lower bound for the change of non-flatness in M and N .
Theorem D.2 Under Assumption D.2 and let η t > 0 satisfies that η t = 1 L D 0 2ξ 0 T . There exists some T > 0 such that the output of Algorithm 3 satisfies that f (x T , y ) -f (x , ȳT ) ≤ and the total number of Riemannian subgradient evaluations is bounded by
, where ξ 0 = ξ(κ min , D) measures the lower bound for the change of non-flatness in M and N , and the time-average iterates (x T , ȳT ) ∈ M × N can be computed by (x 0 , ȳ0 ) = (0, 0) and the inductive formula:
Remark D.1 Theorem D.1 and D.2 establish the last-iterate and time-average rates of convergence of Algorithm 3 for solving Riemannian min-max optimization problems under Assumption D.1 and D.2 respectively. Further, the dependence on L and 1/ can not be improved since it has matched the lower bound established for the nonsmooth min-max optimization problems in Euclidean spaces.
In the scheme of SRGDA, we highlight that (g t x , g t y ) is a noisy estimators of Riemannian subgradient of f at (x t , y t ). It is necessary to impose the conditions such that these estimators are unbiased and has bounded variance. By abuse of notation, we assume that
where the noises (ξ t x , ξ t y ) satisfy that
We are ready to summarize our results for Algorithm 4 in the following theorems.
Theorem D.3 Under Assumption D.1 and let Eq. (D.2) and Eq. (D.3) hold with σ > 0 and let η t > 0 satisfies that η t = 1 µ min{1, 2 t }. There exists some T > 0 such that the output of Algorithm 4 satisfies that E[(d(x T , x )) 2 + (d(y T , y )) 2 ] ≤ and the total number of noisy Riemannian gradient evaluations is bounded by
where ξ 0 = ξ(κ min , D) measures the lower bound for the change of non-flatness in M and N .
Theorem D.4 Under Assumption D.2 and let Eq. (D.2) and Eq. (D.3) hold with σ > 0 and let η t > 0 satisfies that η t = 1 2 D 0 ξ 0 (L 2 +σ 2 )T . There exists some T > 0 such that the output of Algorithm 4 satisfies that E[f (x T , y ) -f (x , ȳT )] ≤ and the total number of noisy Riemannian gradient evaluations is bounded by
, where ξ 0 = ξ(κ min , D) measures the lower bound for the change of non-flatness in M and N , and the time-average iterates (x T , ȳT ) ∈ M × N can be computed by (x 0 , ȳ0 ) = (0, 0) and the inductive formula:
Remark D.2 Theorem D.3 and D.4 establish the last-iterate and time-average rates of convergence of Algorithm 4 for solving Riemannian min-max optimization problems under Assumption D.1 and D.2. Moreover, the dependence on L and 1/ can not be improved since it has matched the lower bound established for nonsmooth stochastic min-max optimization problems in Euclidean spaces.
In this section, we present some technical lemmas for analyzing the convergence property of Algorithm 1 and 2. We also give the proofs of Theorem 3.1, 3.2 and 3.3.
We provide two technical lemmas for analyzing Algorithm 1 and 2 respectively. Parts of the first lemma were presented in Zhang et al. [2022, Lemma C.1]. For the completeness, we provide the proof details.
Lemma E.1 Under Assumption 3.1 and let {(x t , y t ), (x t , ŷt )} T -1 t=0 be generated by Algorithm 1 with the stepsize η > 0. Then, we have
where (x , y ) ∈ M × N is a global saddle point of f .
Proof. Since f is geodesically -smooth, we have the Riemannian gradients of f , i.e., (grad x f, grad y f ), are well defined. Since f is geodesically strongly-concave-strongly-concave with the modulus µ ≥ 0 (here µ = 0 means that f is geodesically concave-concave), we have
Since (x , y ) ∈ M × N is a global saddle point of f , we have f (x t , y ) -f (x , ŷt ) ≥ 0. Recalling also from the scheme of Algorithm 1 that we have
By the definition of an exponential map, we have
This implies that
Putting these pieces together yields that
Equivalently, we have By Proposition C.1 and using that max{d M (x t , x ), d N (ŷ t , y )} ≤ D, we have
By Proposition C.2 and using that max{d M (x t , x ), d N (ŷ t , y )} ≤ D, we have
and Exp -1 ŷt (y t+1 ), Exp -1 ŷt (y ) ≤ 1 2 ξ 0 (d N (ŷ t , y t+1 )) 2 + (d N (ŷ t , y )) 2 -(d N (y t+1 , y )) 2 . By the definition of an exponential map and Riemannian metric, we have
(E.4) Further, we see from the scheme of Algorithm 1 that we have
By the definition of an exponential map, we have
Using the definition of a parallel transport map and the above equations, we have
Since f is geodesically -smooth, we have
Plugging the above inequalities into Eq. (E.4) yields that max {d M (x t , x t+1 ), d N (ŷ t , y t+1 )} ≤ η (d M (x t , x t ) + d N (ŷ t , y t )).
Therefore, we have
Plugging the above inequalities and Eq. (E.3) into Eq. (E.2) yields the desired inequality.
The second lemma gives another key inequality that is satisfied by the iterates generated by Algorithm 2.
Lemma E.2 Under Assumption 3.1 (or Assumption 3.2) and the noisy model (cf. Eq. (3.2) and (3.3)) and let {(x t , y t ), (x t , ŷt )} T -1 t=0 be generated by Algorithm 2 with the stepsize η > 0. Then, we have
where (x , y ) ∈ M × N is a global saddle point of f .
Proof. Using the same argument, we have (µ = 0 refers to geodesically convex-concave case)
Combining the arguments used in Lemma E.1 and the scheme of Algorithm 2, we have
Putting these pieces together with Eq. (3.2) yields that
By the same argument as used in Lemma E.1, we have
Since f is geodesically -smooth and Eq. (3.2) holds, we have
Therefore, we have
Plugging the above inequalities and Eq. (E.6) into Eq. (E.5) yields that
Taking the expectation of both sides and using Eq. (3.3) yields the desired inequality.
Since Riemannian metrics satisfy the triangle inequality, we have
Plugging the above inequality into the inequality from Lemma E.1 yields that
2µ }, we have 4ξ 0 η 2 2 + µη -ξ 0 ≤ 0. By the definition, we have τ 0 ≥ 1, κ ≥ 1 and ξ 0 ≤ 1. This implies that
Putting these pieces together yields that
This completes the proof.
Since Riemannian metrics satisfy the triangle inequality, we have
Plugging the above inequality into the inequality from Lemma E.2 yields that
By the definition, we have τ 0 ≥ 1, κ ≥ 1 and ξ 0 ≤ 1. This implies that
By the inductive arguments, we have
Putting these pieces together yields that
This completes the proof.
By the inductive formulas of xt+1 = Exp xt ( 1 t+1 • Exp -1 xt (x t )) and ȳt+1 = Exp ȳt ( 1 t+1 • Exp -1 ȳt (ŷ t )) and using Zhang et al. [2022, Lemma C.2], we have
Plugging the above inequality into the inequality from Lemma E.2 yields that (recall that µ = 0 in geodesically convex-concave setting here)
T , and 3ξ 0 ησ 2 ≤ 3σ ξ 0 D 0 T . Putting these pieces together yields that
T . This completes the proof.
In this section, we present some technical lemmas for analyzing the convergence property of Algorithm 3 and 4. We also give the proofs of Theorem D.1, D.2, D.3 and D.4.
We provide two technical lemmas for analyzing Algorithm 3 and 4 respectively. The first lemma gives a key inequality that is satisfied by the iterates generated by Algorithm 3.
Lemma F.1 Under Assumption D.1 (or Assumption D.2) and let {(x t , y t )} T -1 t=0 be generated by Algorithm 3 with the stepsize η t > 0. Then, we have
Proof. Since f is geodesically strongly-concave-strongly-concave with the modulus µ ≥ 0 (here µ = 0 means that f is geodesically concave-concave), we have
Recalling also from the scheme of Algorithm 3 that we have
By the definition of an exponential map, we have
This implies that By Proposition C.2 and using that max{d M (x t , x ), d N (y t , y )} ≤ D, we have Letting {b t } t≥1 be a nonnegative sequence such that a t+1 ≤ (1 -P t )a t + Q t 2 where P > 1 and Q > 0. Then, Chung [1954] proved that a t ≤ Q P -1 1 t . Therefore, we have (d M (x t , x )) 2 + (d N (y t , y )) 2 ≤ 8ξ 0 L 2 µ 2 t , for all t ≥ 2.
This completes the proof. Since (x , y ) ∈ M × N is a global saddle point of f , we have E[f (x t , y ) -f (x , y t )] ≥ 0. Plugging this inequality into the inequality from Lemma F.2 yields that E (d M (x t+1 , x )) 2 + (d N (y t+1 , y )) 2 ≤ (1 -µη t )E (d M (x t , x )) 2 + (d N (y t , y )) 2 + 4ξ 0 η 2 t (L 2 + σ 2 ).
Since η t = 1 µ min{1, 2 t }, we have E (d M (x t+1 , x )) 2 + (d N (y t+1 , y )) 2 ≤ (1-2 t )E (d M (x t , x )) 2 + (d N (y t , y )) 2 + 16ξ 0 (L 2 +σ 2 )
, for all t ≥ 2.
Applying the same argument as used in Theorem D.1, we have (d M (x t , x )) 2 + (d N (y t , y )) 2 ≤ 16ξ 0 (L 2 +σ 2 )
, for all t ≥ 2.
This completes the proof.
F.5 Proof of Theorem D.4
Using the same argument, we have
Plugging the above inequality into the inequality from Lemma F.2 yields that (recall that µ = 0 in geodesically convex-concave setting and
This together with (d M (x 0 , x )) 2 + (d N (y 0 , y )) 2 ≤ D 0 implies that E[f (x T , y ) -f (x , ȳT )] ≤ D 0 2ηT + 2ξ 0 η(L 2 + σ 2 ).
Since η = 1 2 D 0 ξ 0 (L 2 +σ 2 )T , we have f (x T , y ) -f (x , ȳT ) ≤ 2 ξ 0 (L 2 +σ 2 )D 0 T . This completes the proof.
We present some additional experimental results for the effect of different choices of α as well the effect of different choices of η for for RCEG. In our experiment here, we set n = 40 consistently. Figure 3 presents the performance of RCEG when α = 2.0. We observe that the results are similar to that summarized in Figure 1. In particular, the last iterate of RCEG consistently achieves the linearly convergence to an optimal solution in all the settings. In contrast, the average iterate of RCEG converges much slower than the last iterate of RCEG. Figure 4 summarizes the effect of different choices of η in RCEG. We observe that setting η as a relatively larger value will speed up the convergence to an optimal solution while all of the choices here lead to the linear convergence. This suggests that the choice of stepsize η in RCEG can be aggressive in practice.
In particular, our assumed upper and lower bounds κmin, κmax guarantee that TCIs in Proposition
2.1 can be used in our analysis for proving finite-time convergence.
For the rest of the presentation, we adopt the convention of presenting the fine-grained complexity performance measure for computing an O( )-close solution instead of the convergence rate of a method. Thus a rate of the form xtx * ≤ O(1/t 1/p ) typically corresponds to O(1/ p ) gradient computations and the geometric rate xtx * ≤ O(exp(-µt)) matches usually up with the O(ln(1/ )) computational complexity.