Nonsmooth Riemannian Optimization

• Nonsmooth Riemannian optimization is the study of optimization on Riemannian manifolds with objectives that lack traditional smoothness, integrating intrinsic geometric constraints.
• It employs methodologies such as tangent-space proximal subproblems, Moreau-envelope smoothing, and operator-splitting techniques to address gradient and feasibility challenges.
• Applications span sparse PCA, robust matrix completion, and graph matching, highlighting scalable strategies for complex geometric problems.

Searching arXiv for papers on nonsmooth Riemannian optimization to ground the article in the current literature. Nonsmooth Riemannian optimization studies optimization problems on Riemannian manifolds when the objective is not everywhere differentiable, or when nonsmooth regularization is coupled with geometric constraints such as orthogonality, fixed rank, or sphere constraints. In this setting, the manifold geometry is intrinsic to feasibility, while the nonsmooth term is often interpreted in the ambient Euclidean space. The area includes composite models of the form $f+h$, split formulations such as $f(x)+g(Ax)$, difference-of-convex objectives, geodesically convex problems on Hadamard manifolds, and stochastic or multiobjective variants. A central difficulty is that many standard Riemannian methods rely on smoothness to define and compute Riemannian gradients, retractions, and manifold subproblems, whereas nonsmoothness introduces generalized subgradients, smoothing homotopies, tangent-space proximal subproblems, augmented Lagrangians, and bundle models (Peng et al., 2022).

1. Problem classes and geometric settings

A large part of the literature focuses on embedded submanifolds of Euclidean space, especially the Stiefel manifold, spheres, fixed-rank manifolds, and compact submanifolds. Typical composite problems take forms such as

$$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$

with $f$ smooth and $g,h$ convex, weakly convex, locally Lipschitz, or lower semicontinuous, depending on the model. Orthogonality-constrained formulations are especially prominent; for example,

$$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$

where the feasible set is the Stiefel manifold $St(m,n)$, and $g$ may be a proximal-friendly weakly convex function (Zhu et al., 5 Apr 2026).

The subject is not confined to compact embedded manifolds. On Hadamard manifolds, proximal point and bundle methods exploit the fact that every pair of points is joined by a unique minimizing geodesic and that the exponential map is globally defined. In that setting, nonsmooth geodesically convex optimization becomes the natural analogue of Euclidean convex optimization, and the squared distance is itself geodesically convex (Amaral et al., 18 Nov 2025). At the opposite extreme, locally Lipschitz multiobjective problems have been studied on complete Riemannian manifolds without convexity or differentiability assumptions on the objectives (Tang et al., 2023).

Application structure is a major driver of model design. The cited literature repeatedly uses sparse PCA, compressed modes, graph matching with penalties, robust low-rank matrix completion, robust subspace recovery, DPCP, generalized eigenvalue problems, orthogonal dictionary learning, geometric median-type problems on spheres, and regularized classification on the sphere as representative nonsmooth manifold problems (Kangkang et al., 2019).

2. Stationarity, subdifferentials, and manifold-specific variational notions

Because ordinary gradients may not exist, nonsmooth Riemannian optimization relies on generalized first-order conditions. Several papers use the Riemannian Clarke subdifferential, defined through pullbacks by a retraction; in one common formulation, a point is stationary when

$$0\in \operatorname{grad} f(X)+\operatorname{Proj}_{T_X\mathcal{M}}\partial h(X),$$

or, more generally,

$0_x\in \partial f(x).$

For split constrained models, $f(x)+g(Ax)$0-KKT notions are also standard; for example, for orthogonality constraints one uses

$f(x)+g(Ax)$1

with approximate versions based on $f(x)+g(Ax)$2 and $f(x)+g(Ax)$3 (Wang et al., 2020).

A notable feature of the area is the coexistence of several stationarity concepts, each matched to a different analytical framework. Moreau-envelope smoothing methods define generalized $f(x)+g(Ax)$4-stationarity through a nearby prox-corrected point and an ambient-space subgradient of the original nonsmooth term. Riemannian ADMM methods use approximate KKT conditions involving tangent projections, subdifferential residuals, and feasibility residuals. Goldstein-type analyses define $f(x)+g(Ax)$5-stationarity through the minimum norm of a transported $f(x)+g(Ax)$6-subdifferential. Multiobjective methods replace scalar stationarity by necessary conditions for Pareto optimality such as

$f(x)+g(Ax)$7

and by $f(x)+g(Ax)$8-criticality based on the convex hull of Riemannian $f(x)+g(Ax)$9-subgradients (Sahinoglu et al., 24 Oct 2025).

The literature also emphasizes that subgradients live in different tangent spaces and must be compared by parallel transport, vector transport, or projection. This point is central in stochastic and multiobjective analyses, where one repeatedly constructs transported gradients or $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$0-subgradients before taking convex hulls, recursive averages, or benchmark comparisons. In embedded settings, many first-order conditions reduce to ambient Euclidean gradients projected onto tangent spaces; in intrinsic settings, retraction pullbacks and tangent-space localization are essential (Tang et al., 2023).

3. Retraction-based, proximal, smoothing, and Bregman methods

The earliest widely used algorithmic pattern in the area is to solve a tangent-space proximal subproblem and then return to the manifold by a retraction. For the Stiefel manifold, Riemannian stochastic proximal gradient methods solve

$$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$1

and combine this with variance reduction in the finite-sum and online settings (Wang et al., 2020). Related deterministic methods on structured nonsmooth problems over embedded submanifolds use proximal-gradient-subgradient steps on the tangent space followed by retraction and backtracking, with subsequential criticality and KL-based full-sequence convergence (Li et al., 2024).

A second major pattern is smoothing. Moreau-envelope homotopies replace a weakly convex or Lipschitz nonsmooth term by a smooth envelope with a decreasing smoothing parameter. For compact embedded submanifolds, this yields Riemannian smoothing gradient and stochastic gradient methods for

$$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$2

with deterministic $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$3 complexity and stochastic $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$4 complexity for generalized $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$5-stationarity (Peng et al., 2022). Later single-loop stochastic smoothing methods combine dynamic smoothing with momentum-based variance reduction and attain $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$6 when $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$7 is Lipschitz continuous, while a distinct indicator-function regime yields $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$8 under an error bound condition (Deng et al., 14 May 2025).

A third line extends proximal and Bregman-gradient ideas beyond globally Lipschitz gradients. Relative smoothness methods on embedded submanifolds use a reference function $$\min_{x\in M} f(x)+g(x), \qquad \min_{x\in M} f(x)+h(Ax), \qquad \min_{x\in M} f(x)+h(x)-g(x),$$9 and the Bregman distance

$f$0

The retraction-based method computes

$f$1

and handles nonsmooth convex $f$2 directly. When the reference function is $f$3, the constraint subproblem admits a closed-form solution, and both deterministic methods achieve $f$4 iteration complexity (He et al., 5 Aug 2025).

4. Augmented Lagrangian, ADMM, and primal-dual frameworks

Operator-splitting methods form another major branch of nonsmooth Riemannian optimization. Early manifold augmented Lagrangian methods split the nonsmooth term by introducing an auxiliary variable $f$5, smooth the inner subproblem by the Moreau envelope, and then solve a sequence of smooth manifold subproblems. This strategy underlies the manifold inexact augmented Lagrangian method for

$f$6

with convergence to a critical point under suitable assumptions, and with convergence to a global minimizer if the approximate global minimizer of the iteration subproblem could be obtained (Kangkang et al., 2019).

Riemannian ADMM adapted this splitting principle to updates that are easy to compute on manifolds. For problems of the form

$f$7

a smoothed augmented Lagrangian, a Riemannian gradient step for the manifold block, and a proximal step for the nonsmooth block yield an $f$8 complexity bound for obtaining an $f$9-stationary point (Li et al., 2022). Subsequent work removed the need for smoothing. The adaptive Riemannian ADMM for

$g,h$0

introduces an auxiliary splitting variable and an adaptive coordination of the stepsizes and penalty parameters, uses only one Riemannian gradient evaluation and one proximal update per iteration, and achieves $g,h$1 for an $g,h$2-approximate KKT point without any smoothing approximation (Deng et al., 21 Oct 2025).

Retraction-free primal-dual approaches go one step further for orthogonality constraints. A linearized smoothing augmented Lagrangian method uses a Moreau-envelope-style smoothing and a strictly convex surrogate around the current iterate, but avoids retractions, manifold subproblem solving, semismooth Newton inner loops, matrix inversion, and the exponential map. In the Stiefel setting, this single-loop primal-dual framework attains $g,h$3 iterations for an $g,h$4-KKT point and establishes asymptotic sequential convergence under the standard Kurdyka–Łojasiewicz property (Zhu et al., 5 Apr 2026).

Stochastic splitting methods have also matured rapidly. The first stochastic Riemannian ADMM with provable near-optimal complexity, MARS-ADMM, treats

$g,h$5

by combining one proximal $g,h$6-update, one retraction step, a momentum-based variance-reduced estimator for the smooth component, and adaptive penalty and dual stepsize updates. Under mild assumptions, it achieves $g,h$7, improving upon the previously best-known bound of $g,h$8 for stochastic Riemannian operator-splitting methods (Jin et al., 28 Dec 2025).

5. Bundle methods, conjugate-subgradient methods, proximal point methods, and inexact primitives

Not all nonsmooth manifold algorithms are proximal-gradient or splitting methods. Bundle methods transfer classical Euclidean cutting-plane stabilization to tangent spaces. A restricted memory quasi-Newton bundle method minimizes a locally Lipschitz function over a complete Riemannian manifold, uses subgradient aggregation to avoid solving the time-consuming quadratic programming subproblem when calculating the candidate descent direction, incorporates Riemannian SR1 and BFGS updates, and proves that if the serious iteration steps are finite, then the last serious iterate is stationary; otherwise, every accumulation point of the serious iteration sequence is stationary (Tang et al., 2024). A related conjugate-subgradient line uses a convex combination of the negative selected subgradient and the previous search direction transported to the current tangent space, together with an interval-reduction line search, and is presented as the first CG-type method developed for Riemannian nonsmooth optimization (Tang et al., 7 Sep 2025).

Hadamard-manifold methods provide a different geometry. An adaptive proximal point framework for

$g,h$9

works on complete, simply connected Riemannian manifolds with nonpositive sectional curvature, exploits the global exponential map, and obtains $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$0 iteration complexity both when the Lipschitz constant of the smooth part is known and when it is learned adaptively by backtracking (Amaral et al., 18 Nov 2025). For geodesically convex objectives, bundle methods have recently been extended to inexact manifold primitives: a proximal bundle method on Hadamard manifolds establishes the first oracle-complexity bounds that rely only on subgradients and inexact primitives, obtains sublinear rates for general objectives, and optimal linear convergence under sharp function growth (Díaz et al., 29 Apr 2026).

A further development is the introduction of nonsmooth Riemannian difference-of-convex optimization. For

$$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$1

with $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$2 and $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$3 convex and possibly nonsmooth, the inexact Riemannian proximal DC framework returns an $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$4-Riemannian critical point within $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$5 outer iterations, while one practical variant attains an overall $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$6 complexity. The same work also establishes, for the first time in the manifold setting, equivalence between suitably chosen nonsmooth DC terms and $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$7-regularized or $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$8-constrained counterparts on the sphere (Jiang et al., 10 Sep 2025).

6. Stochastic finite-time theory, local convergence, and geometric interpretations

Stochastic finite-time theory for fully nonsmooth manifold objectives has developed along two complementary directions. One direction studies composite problems with explicit smooth-plus-nonsmooth structure, as in stochastic proximal-gradient, stochastic smoothing, and stochastic ADMM methods. The other direction treats fully nonsmooth nonconvex objectives through Goldstein-type stationarity. In the latter setting, RO2NC and its zeroth-order variant extend online-to-nonconvex conversion to complete embedded Riemannian manifolds, define the Riemannian $$\min_{X\in\mathbb{R}^{m\times n}} \ell(X)+g(X)\quad \text{s.t.}\quad X^\top X=I_n,$$9-subdifferential by transported Clarke subgradients, and prove the first-ever finite-time guarantee for fully nonsmooth, nonconvex optimization on manifolds: $St(m,n)$0 to find a $St(m,n)$1-Goldstein stationary point (Sahinoglu et al., 24 Oct 2025).

Alongside nonasymptotic theory, local and asymptotic analyses clarify why some augmented Lagrangian and subgradient methods work well in practice. A local variational theory for nonsmooth optimization on manifolds introduces the manifold variational sufficient condition, shows that its strong version is equivalent to the manifold strong second-order sufficient condition in certain circumstances, constructs a local dual problem, and proves linear convergence of the RALM without imposing uniqueness of the multiplier (Zhou et al., 2023). At a more geometric level, generic nonsmooth subgradient dynamics can be interpreted as an approximate Riemannian gradient method on an active manifold that captures the function’s nonsmooth activity. In that picture, off the manifold, subgradients point back toward it; on the manifold, the function is smooth; and strict saddles are avoided by perturbed dynamics (Davis et al., 2021).

This suggests a unifying interpretation of the field. Some methods preserve manifold geometry by tangent-space subproblems and retractions; some avoid retractions altogether by primal-dual ambient-space formulations; some replace nonsmoothness by smoothing or by bundle models; and some show that, locally, a nonsmooth problem behaves like smooth optimization on a distinguished submanifold. The recurring technical obstacles are likewise stable across the literature: tangent spaces vary with the iterate, geodesics and parallel transport may be expensive, retraction errors can accumulate, inner subproblems can dominate cost, and feasibility control is delicate under nonsmoothness (Aspman et al., 2023).

7. Applications, empirical patterns, and recurrent themes

Sparse PCA is the most common benchmark across the literature. It appears in proximal-gradient, smoothing, ADMM, primal-dual, and DC formulations on the Stiefel manifold and sphere, typically with $St(m,n)$2, capped-$St(m,n)$3, or $St(m,n)$4 regularization. Robust matrix completion, compressed modes, DPCP, robust subspace recovery, graph matching, geometric median-type models, regularized classification on the sphere, and low-rank quadratic sensing play a similar role in other geometric settings (Jiang et al., 10 Sep 2025).

Across these experiments, several empirical regularities recur. Single-loop methods are repeatedly reported to be attractive because they avoid nested inner solvers. Explicit proximal or linearized updates often reduce per-iteration cost. Retraction-free or smoothing-free formulations are motivated by feasibility drift, retraction errors, or the overhead of tangent-space manifold subproblems. In large-scale Stiefel problems, methods such as LSALM, ARADMM, and stochastic smoothing schemes are reported to scale better than earlier state-of-the-art algorithms, while bundle and conjugate-subgradient methods derive their efficiency from avoiding expensive quadratic subproblems (Zhu et al., 5 Apr 2026).

A recurrent misconception is that nonsmooth Riemannian optimization is merely Euclidean nonsmooth optimization with a manifold projection appended at the end. The literature consistently argues otherwise. Tangent spaces vary from point to point; generalized subgradients must be transported or localized; retractions are first-order approximations whose errors enter the analysis; curvature modifies model validity and oracle complexity; and even the notion of stationarity is not unique across problem classes. Another recurrent misconception is that Riemannian methods are inherently restricted to smooth objectives. The current body of work shows instead that nonsmoothness can be handled through tangent-space proximal models, Moreau-envelope smoothing, ADMM and augmented Lagrangian splitting, bundle stabilization, Goldstein stationarity, and DC decompositions, depending on the geometry and the objective structure (Díaz et al., 29 Apr 2026).

Taken together, these developments define nonsmooth Riemannian optimization as a mature intersection of manifold optimization, variational analysis, nonsmooth analysis, and operator splitting. The field now contains deterministic, stochastic, zeroth-order, multiobjective, convex, weakly convex, DC, and geodesically convex theories, as well as local second-order analyses and finite-time guarantees. A plausible implication is that future progress will continue to be organized less by a single canonical algorithm than by the interaction between geometry, nonsmooth structure, and the computational cost of manifold primitives.