Adaptive Single-Loop Methods for Stochastic Minimax Optimization on Riemannian Manifolds

Abstract: Stochastic minimax optimization on Riemannian manifolds has recently attracted significant attention due to its broad range of applications, such as robust training of neural networks and robust maximum likelihood estimation. Existing optimization methods for these problems typically require selecting stepsizes based on prior knowledge of specific problem parameters, such as Lipschitz-type constants and (geodesic) strong concavity constants. Unfortunately, these parameters are often unknown in practice. To overcome this issue, we develop single-loop adaptive methods that automatically adjust stepsizes using cumulative Riemannian (stochastic) gradient norms. We first propose a deterministic single-loop Riemannian adaptive gradient descent ascent method and show that it attains an $ε$-stationary point within $O(ε<sup>{-2})$ iterations. This deterministic method is of independent interest and lays the foundation for our subsequent stochastic method. In particular, we propose the Riemannian stochastic adaptive gradient descent ascent method, which finds an $ε$-stationary point in $O(ε<sup>{-6})$ iterations. Under additional second-order smoothness, this iteration complexity is further improved to $O(ε<sup>{-4})$, which even outperforms the corresponding complexity result in Euclidean space. Some numerical experiments on real-world applications are conducted, including the regularized robust maximum likelihood estimation problem, and the robust training of neural networks with orthonormal weights. The results are encouraging and demonstrate the effectiveness of adaptivity in practice.

• This paper introduces two single-loop adaptive algorithms, RAGDA and RSAGDA, for nonconvex-geodesically-strongly-concave minimax optimization on Riemannian manifolds.
• RAGDA and RSAGDA achieve an $\epsilon$-stationary point in $\mathcal{O}(\epsilon^{-2})$ iterations for deterministic settings and $\mathcal{O}(\epsilon^{-6})$ in stochastic settings, with an option for $\mathcal{O}(\epsilon^{-4})$ under additional second-order smoothness.
• The pivotal factor is the coupling of primal and dual stepsize via $\max\{v^x,v^y\}$. Two state-of-the-art algorithms are proposed: RAGDA for deterministic settings, and RSAGDA for stochastic settings.

Problem setting and motivation

The paper studies stochastic minimax optimization on Riemannian manifolds:

$$\min_{x \in \mathcal{M}_x}\max_{y\in\mathcal{M}_y} f(x,y) = \mathbb{E}_{\xi \sim \mathcal{D}}[f(x,y;\xi)],$$

where $f$ is nonconvex in $x$ and geodesically strongly concave in $y$, and both variables lie on complete Riemannian manifolds. Two applications motivate the formulation: adversarially robust training of neural networks with orthonormal (Stiefel-constrained) weights against universal perturbations on a sphere, and robust maximum likelihood estimation over the sphere and the positive-definite cone (2602.09840).

The central gap addressed is adaptivity. Existing Riemannian minimax methods—Riemannian Hamiltonian methods, corrected extragradient, gradient descent ascent variants—require prior knowledge of Lipschitz constants and strong concavity parameters. The only adaptive method in this space, AdaRHD, is double-loop and analyzed only for deterministic bilevel problems. The paper proposes the first single-loop adaptive algorithms for Riemannian minimax optimization on manifolds, in both deterministic and stochastic settings.

Algorithm design

Both proposed methods follow a cumulative-gradient-norm adaptation scheme inspired by AdaGrad-style methods and TiADA. The deterministic method, RAGDA, maintains accumulators $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$ and $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$, then updates via retractions:

$$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$

Two design choices are theoretically consequential. First, the primal stepsize uses $\max\{v^x,v^y\}$ rather than $v^x$ alone; the analysis shows dual dynamics are controlled through primal gradients, so coupling the stepsizes is what enables a single-loop scheme with guarantees. Second, distinct decay exponents $\alpha,\beta$ generalize the standard choice of $f$0 and allow tuning of the time-scale separation between primal and dual updates. The stochastic variant, RSAGDA, replaces exact gradients with unbiased bounded stochastic oracles while retaining the same structure.

The analysis rests on standard assumptions: joint Lipschitz smoothness under parallel transport, geodesic strong concavity in $f$1, boundedness of the value function $f$2 and its gradient, Hadamard structure on $f$3 with a curvature-dependent trigonometric-distance constant $f$4, and second-order retraction accuracy.

Deterministic convergence

The proof strategy bounds the accumulated squared gradient norms by controlling the auxiliary sequences $f$5 and $f$6. A key lemma shows that after an initial phase ($f$7 iterations until $f$8 exceeds a threshold), accumulated primal gradients control accumulated dual gradients, with the interaction made explicit through $f$9. Because the primal stepsize depends on $x$0, the analysis splits into the regimes $x$1 (self-consistent bound on $x$2) and $x$3 (indirect bound via $x$4).

Theorem 1 establishes that for any $x$5, RAGDA attains $x$6, i.e., an $x$7-stationary point in $x$8 iterations—matching the best known rates for nonadaptive Riemannian GDA-type methods without requiring any problem parameters.

Stochastic convergence

Under unbiasedness and uniform boundedness of stochastic gradients, the analysis proceeds in expectation. Lemma 5 bounds expected cumulative dual suboptimality $x$9 in terms of expected accumulated primal gradients and moments of $y$0; the bound is converted to dual gradient norms via smoothness-plus-concavity, and primal norms are controlled through a descent argument on $y$1 with a carefully chosen constant $y$2 that cancels cross terms.

Theorem 2 gives, for $y$3:

$y$4

optimized at $y$5 to yield $y$6, hence $y$7 complexity. This matches the rate TiADA achieves in Euclidean space under comparable assumptions.

Improved rate under second-order smoothness

Adding Lipschitz continuity of the mixed derivative $y$8 and Hessian $y$9 yields Lipschitz smoothness of the solution map $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$0 (with constant $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$1). This permits a sharper local approximation of $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$2 around $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$3: the error term becomes quadratic in the primal stepsize rather than linear, which relaxes the exponent constraint from $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$4 to $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$5 (via a slack parameter $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$6).

Theorem 3 then gives $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$7 at $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$8, i.e., $v^x_{t+1}=v^x_t+\|\nabla_x f(x_t,y_t)\|^2$9 complexity. Notably, the authors point out this sharpens even the Euclidean TiADA result: TiADA's rate was stated as $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$0 for some small $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$1, whereas the refined hyperparameter analysis here achieves exactly $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$2 per squared-gradient-norm iteration. This is a bold claim relative to prior art, resting on the more flexible choice of decay exponents enabled by the relaxed constraint.

Numerical experiments

Experiments compare RAGDA (MATLAB/Manopt) against RHM, RCON, RGDA, RCEG, TSRGDA, and a modified AdaRHD on regularized robust maximum likelihood estimation ($v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$3, $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$4, $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$5), where RAGDA converges faster and to higher accuracy in both iteration count and wall-clock time. RSAGDA (PyTorch/Geoopt) is compared against RSGDA and Acc-RSGDA on robust training of a three-layer network with orthonormal hidden weights on MNIST, FashionMNIST, and CIFAR-10. Representative results:

Dataset	Method	Natural acc.	PGD$v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$6 (largest $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$7)	FGSM (largest $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$8)
MNIST	RSGDA / Acc-RSGDA / RSAGDA	97.13 / 96.77 / 97.66	89.49 / 89.25 / 90.96	89.97 / 89.66 / 91.37
FashionMNIST	RSGDA / Acc-RSGDA / RSAGDA	87.16 / 87.71 / 88.78	77.81 / 77.90 / 78.29	77.93 / 78.10 / 78.58
CIFAR-10	RSGDA / Acc-RSGDA / RSAGDA	48.50 / 47.31 / 49.30	41.75 / 41.64 / 42.86	41.81 / 41.67 / 42.91

RSAGDA consistently achieves higher training/test accuracy, lower loss, faster convergence per gradient call, and better accuracy under all attack radii. One caveat: the baselines from Huang & Gao were originally designed for convex dual domains, and the authors substitute projections with retractions/exponential maps when updating the maximization variable—an adaptation not covered by those baselines' original analyses.

Limitations and open questions

Several assumptions constrain applicability. The dual manifold must be Hadamard with sectional curvature lower bounded by $v^y_{t+1}=v^y_t+\|\nabla_y f(x_t,y_t)\|^2$9, and the curvature correction constant $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$0 must be uniformly bounded along iterates—a condition that may fail if iterates drift far from $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$1. Stochastic gradients are assumed uniformly bounded, which excludes heavy-tailed noise. The improved $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$2 rate requires second-order smoothness, and its advantage over the general result depends on that assumption holding. The theory guarantees only first-order stationarity of the saddle-point residuals, not global optimality of the nonconvex primal problem, and no lower bounds are given to certify that $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$3 (or $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$4) is optimal for this adaptive single-loop class. Experiments are limited to small-scale networks and moderate-dimensional problems; scalability to large models remains untested. Open questions raised by the authors include variance-reduced adaptive extensions, nonsmooth objectives, constrained product-manifold settings, and integration into federated learning systems.

Conclusion

This paper closes a gap in Riemannian optimization by providing the first single-loop, parameter-adaptive methods for nonconvex–geodesically-strongly-concave minimax problems. RAGDA achieves $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$5 deterministic complexity and RSAGDA achieves $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$6, improvable to $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$7 under second-order smoothness—the latter matching and slightly sharpening the best known Euclidean adaptive rate. The coupled $$x_{t+1} = \mathrm{Retr}_{x_t}\!\left(-\frac{\eta^x}{\max\{v^x_{t+1},v^y_{t+1}\}^{\alpha}}\nabla_x f(x_t,y_t)\right), \quad y_{t+1} = \mathrm{Retr}_{y_t}\!\left(\frac{\eta^y}{(v^y_{t+1})^{\beta}}\nabla_y f(x_t,y_t)\right).$$8 stepsize mechanism is identified as the structural ingredient enabling single-loop guarantees, and empirical results on robust likelihood estimation and adversarially robust network training support the practical value of the adaptive mechanism.