Zeroth-Order Riemannian Gradient Estimator

• Zeroth-order Riemannian gradient estimator is a surrogate that approximates the true Riemannian gradient using only function evaluations and intrinsic geometric constructions.
• It leverages mechanisms like tangent-space finite differences, ambient perturbation with projection, and pullback-based coordinate differences to respect manifold geometry.
• The estimator’s bias, variance, and moment bounds depend on intrinsic dimension and curvature, guiding its use in stochastic and online Riemannian optimization.

A zeroth-order Riemannian gradient estimator is a tangent-valued surrogate for $\operatorname{grad} f(x)$ on a manifold $\mathcal M$ constructed from function evaluations alone. In the literature, this notion covers several distinct mechanisms: tangent-space finite differences pushed back to the manifold by $\operatorname{Exp}_x$ or $\operatorname{Retr}_x$, projection-based ambient perturbations for embedded manifolds, pullback-based coordinate-wise differences in a fixed tangent space, diffusion-map-style kernel estimators from manifold samples, and structure-exploiting algebraic surrogates for quotient-of-quadratics problems on spheres and generalized spheres. What unifies these constructions is that the estimated object is intrinsically Riemannian: a vector in $T_x\mathcal M$ intended to approximate, or in special cases recover up to a known factor, the Riemannian gradient governing first-order stationarity (Li et al., 2020, Li et al., 2023, Bresch, 29 Apr 2026).

1. Geometric target and formal role

In stochastic Riemannian optimization, the basic problem is

$$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$

where $\mathcal M$ is a $d$-dimensional complete Riemannian manifold and only noisy function values are available. The target stationarity notion is intrinsically first-order: a point $\bar x$ is an $\epsilon$-approximate first-order stationary point if

$\mathcal M$0

The Riemannian gradient is defined by

$\mathcal M$1

for every curve $\mathcal M$2 with $\mathcal M$3 and $\mathcal M$4 (Li et al., 2023).

For embedded submanifolds $\mathcal M$5 with the induced Euclidean metric, the target vector can be written as the tangent projection of the ambient Euclidean gradient: $\mathcal M$6 This identity is central in embedded-manifold constructions, because it makes clear that zeroth-order estimation on manifolds is not merely estimation of an ambient derivative; it is estimation of a tangent vector whose codomain changes with the iterate (Li et al., 2020, Maass et al., 2020).

The manifold geometry determines how perturbations are generated and how function evaluations are made. Classical constructions use tangent perturbations $\mathcal M$7, then map $\mathcal M$8 back to the manifold through $\mathcal M$9 or a retraction. Other constructions replace tangent perturbation by ambient perturbation followed by projection, or bypass explicit tangent geometry entirely by using local kernel averages on sampled manifold points. These distinctions are not cosmetic: they determine whether the estimator is biased, whether its variance depends on intrinsic or ambient dimension, and which geometric assumptions are needed for analysis (Li et al., 2023, Wang et al., 30 Jul 2025, Gomez et al., 2021).

2. Principal estimator constructions

The foundational intrinsic construction is the tangent-space finite-difference estimator. In the embedded-submanifold framework of stochastic zeroth-order Riemannian optimization, one samples $\operatorname{Exp}_x$0, projects it to the tangent space as $\operatorname{Exp}_x$1, and defines

$\operatorname{Exp}_x$2

Its stochastic minibatch version replaces $\operatorname{Exp}_x$3 by $\operatorname{Exp}_x$4 and averages over $\operatorname{Exp}_x$5 samples. This is a one-sided Gaussian finite-difference estimator on the pullback $\operatorname{Exp}_x$6, and its moment bounds depend on the intrinsic dimension $\operatorname{Exp}_x$7, not the ambient dimension $\operatorname{Exp}_x$8 (Li et al., 2020).

A closely related construction uses geodesics or retractions directly on abstract manifolds. With $\operatorname{Exp}_x$9 in $\operatorname{Retr}_x$0,

$\operatorname{Retr}_x$1

and the practical retraction analogue is

$\operatorname{Retr}_x$2

Because each summand is a scalar multiple of $\operatorname{Retr}_x$3, the estimator lies in the correct tangent space by construction. This class underlies the Zo-RASA framework and related manifold zeroth-order stochastic approximation methods (Li et al., 2023).

For online optimization on Hadamard manifolds, the same intrinsic logic appears in a two-point asynchronous form: $\operatorname{Retr}_x$4 The two function values may be taken at different times, which introduces an online bias term absent in time-invariant settings. Geodesic perturbation remains intrinsic, but the estimator is explicitly designed for changing objectives rather than static stochastic losses (Maass et al., 2020).

Several later works depart from tangent-Gaussian smoothing. One line fixes a base point $\operatorname{Retr}_x$5, works in the Euclidean tangent model $\operatorname{Retr}_x$6, and applies coordinate-wise central differences to the pullback

$\operatorname{Retr}_x$7

The resulting estimator is

$\operatorname{Retr}_x$8

which estimates $\operatorname{Retr}_x$9, and at $T_x\mathcal M$0 estimates $T_x\mathcal M$1 directly (He et al., 2024).

Another line avoids tangent sampling altogether on embedded manifolds with tractable projection. The projection-based estimator is

$T_x\mathcal M$2

with $T_x\mathcal M$3 uniform on the ambient unit sphere in $T_x\mathcal M$4. Here the perturbation is purely Euclidean, and feasibility is restored by the manifold projection $T_x\mathcal M$5. This eliminates tangent-basis construction and tangent-space random vector generation, which the paper identifies as a computational bottleneck on manifolds such as $T_x\mathcal M$6 (Wang et al., 30 Jul 2025).

A more data-driven construction appears in the diffusion-map literature. For an embedded manifold $T_x\mathcal M$7, the estimator

$T_x\mathcal M$8

converges to the intrinsic Riemannian gradient at interior smooth points. This method does not require explicit tangent spaces, charts, exponential maps, or retractions for gradient estimation itself; instead it uses ambient displacements and local kernel averages over sampled manifold points (Gomez et al., 2021).

Finally, certain quotient-of-quadratics problems admit algebraic surrogates that are not generic smoothing estimators. For generalized Rayleigh quotient maximization on the sphere, full-sphere random directions can be used to recover the Riemannian gradient in expectation without sampling in the tangent space. This is a structure-specific departure from standard manifold zeroth-order methodology (Bresch, 29 Apr 2026).

3. Bias, approximation, and moment identities

Most generic zeroth-order Riemannian estimators are biased at fixed smoothing radius. In the tangent-Gaussian retraction framework, the bias and second-moment bounds are explicit: $T_x\mathcal M$9 and

$$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$0

For the stochastic minibatch estimator, the mean-squared error includes a smoothing term and a variance term: $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$1 These formulas are the basic nonasymptotic justification for tangent-space Gaussian smoothing on manifolds (Li et al., 2020).

The Zo-RASA analysis adopts the same intrinsic estimator family and controls its bias relative to the true Riemannian gradient through

$$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$2

It also bounds the second moment and MSE by quantities scaling with intrinsic dimension $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$3, batch size $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$4, and the stochastic first-order variance $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$5. In this framework, $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$6 controls smoothing bias, while the moving-average recursion controls variance accumulation across iterates (Li et al., 2023).

In online Hadamard-manifold optimization, the bias has two sources. For

$$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$7

the bound is

$$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$8

The term proportional to $$\min_{x\in \mathcal M} f(x) := \mathbb{E}_{\xi}[F(x,\xi)],$$9 is the usual smoothing/discretization error; the term proportional to $\mathcal M$0 is the asynchrony penalty from objective variation between the two evaluations. This separation is one of the paper’s main conceptual points (Maass et al., 2020).

Projection-based ambient perturbation produces the same qualitative structure. Under proximal smoothness and Euclidean smoothness assumptions, the projection-based estimator satisfies

$\mathcal M$1

and

$\mathcal M$2

Thus the estimator is asymptotically unbiased as $\mathcal M$3, with mean-square error decomposing into an $\mathcal M$4 approximation term and an $\mathcal M$5 sampling term (Wang et al., 30 Jul 2025).

The intrinsic two-point estimator on geodesically complete manifolds is analyzed directly in mean square rather than through an ambient embedding. For

$\mathcal M$6

with $\mathcal M$7 uniform on the Riemannian unit sphere, the intrinsic MSE bound is

$\mathcal M$8

Here curvature enters through the factor $\mathcal M$9, making the dependence on sectional curvature explicit rather than embedding-dependent (Ma et al., 12 Jan 2026).

The main exception to the generic “biased unless $d$0” pattern arises in structure-specific quotient problems. For generalized operator norm and generalized Rayleigh quotient optimization on $d$1, the estimator

$d$2

satisfies

$d$3

with a matching second-moment identity

$d$4

The estimator is therefore unbiased up to a known positive scalar factor, even though the sampled directions are drawn from the full ambient sphere rather than the tangent space (Bresch, 29 Apr 2026).

4. Second-order surrogates and structure-aware variants

Zeroth-order Riemannian estimation is not restricted to first-order information. In the tangent-Gaussian framework, a Riemannian Hessian estimator is built from symmetric function differences: $d$5 This estimator underlies zeroth-order cubic-regularized Newton methods and comes with operator-norm approximation guarantees under Hessian-Lipschitz assumptions formulated via parallel transport (Li et al., 2020).

The quotient-of-quadratics setting yields even sharper structure-aware second-order formulas. For $d$6, the matrix

$d$7

is unbiased for $d$8 under full-sphere sampling, while a tangent-sphere variant is unbiased for a multiple of $d$9. The same paper then combines empirical gradient and Hessian surrogates in a zeroth-order quasi-Newton or Newton-style update, with the step-size subproblem still solvable in closed form because the directional reduction is an exact rational quadratic (Bresch, 29 Apr 2026).

Pullback-based acceleration uses a different second-order principle. The estimator

$\bar x$0

approximates $\bar x$1 with deterministic error

$\bar x$2

This estimator is then evaluated at extrapolated tangent-space points inside an accelerated inner loop, so the main novelty is not stochastic smoothing but the integration of a deterministic zeroth-order pullback gradient into an accelerated tangent-space dynamics (He et al., 2024).

A related structure-aware construction appears for generalized Rayleigh quotients on the $\bar x$3-sphere $\bar x$4. With a random tangent unit direction $\bar x$5, the scalar

$\bar x$6

produces the tangent estimator $\bar x$7, and the paper proves

$\bar x$8

together with

$\bar x$9

Here the estimator is explicitly tied to an exact line-search update and avoids both $\epsilon$0 and $\epsilon$1 (Bresch et al., 5 Dec 2025).

5. Algorithmic frameworks and complexity results

The estimator enters optimization in several distinct ways. In the basic retraction-Gaussian framework, smooth nonconvex Riemannian gradient descent uses

$\epsilon$2

or the stochastic minibatch analogue with $\epsilon$3. The resulting oracle complexities depend on the intrinsic dimension $\epsilon$4: $\epsilon$5 in deterministic smooth nonconvex optimization, $\epsilon$6 in stochastic smooth nonconvex optimization, and $\epsilon$7 in the geodesically convex setting (Li et al., 2020).

Zo-RASA embeds the intrinsic tangent-Gaussian estimator inside a transported moving average,

$\epsilon$8

or its practical vector-transport variant. The main result is that one-sample or constant-order batches suffice to achieve the same $\epsilon$9 oracle complexity for $\mathcal M$00-first-order stationarity, while retractions and vector transport replace exponential maps and parallel transport in the practical algorithm (Li et al., 2023).

The pullback-based accelerated framework improves the $\mathcal M$01-dependence. RAZGD attains $\mathcal M$02 function-query complexity for an $\mathcal M$03-approximate first-order stationary point, and $\mathcal M$04 for an $\mathcal M$05-approximate second-order stationary point with high probability. It also permits larger smoothing parameters, including the second-order choice

$\mathcal M$06

which the paper interprets as improved robustness of finite-difference estimation (He et al., 2024).

Projection-based zeroth-order federated learning uses the ambient-perturbation estimator inside local projected updates and proves sublinear convergence in a projected-gradient-type stationarity metric. Under the stated conditions,

$\mathcal M$07

after choosing $\mathcal M$08 sufficiently small. The paper emphasizes that this matches the rate of the corresponding first-order federated method up to controllable zeroth-order estimation terms (Wang et al., 30 Jul 2025).

For geodesically incomplete manifolds, the algorithmic role of the estimator changes because $\mathcal M$09 may be undefined under the original metric. The proposed remedy is to replace the original metric $\mathcal M$10 by a structure-preserving complete metric $\mathcal M$11, then run retraction-based SGD with the intrinsic symmetric two-point estimator under $\mathcal M$12. Under suitable conditions, the resulting method again achieves the standard nonconvex zeroth-order complexity order $\mathcal M$13, and $\mathcal M$14-stationarity can be related back to the original metric (Ma et al., 12 Jan 2026).

6. Scope, limitations, and related directions

The term “zeroth-order Riemannian gradient estimator” does not denote a single canonical formula. Some estimators are intrinsic tangent-space smoothers; some are projection-based embedded-manifold constructions; some are deterministic pullback coordinate differences; some are data-driven kernel operators; and some are exact algebraic surrogates available only for quotient-of-quadratics geometry. A common misconception is that tangent-space sampling is mandatory. The literature contains counterexamples: ambient projection-based perturbations on embedded manifolds, full-sphere sampling for generalized Rayleigh quotient problems, and diffusion-map estimators that never construct tangent bases explicitly (Wang et al., 30 Jul 2025, Bresch, 29 Apr 2026, Gomez et al., 2021).

Another misconception is that zeroth-order Riemannian estimators are generically unbiased. The opposite is typical: in tangent-Gaussian smoothing, online two-point estimation, projection-based perturbation, and pullback finite differences, fixed-$\mathcal M$15 bias is part of the basic theory. Exact or scalar-unbiased identities arise only in specialized settings, such as the quotient-of-quadratics estimators on spheres or $\mathcal M$16-spheres (Li et al., 2020, Maass et al., 2020, Bresch, 29 Apr 2026, Bresch et al., 5 Dec 2025).

The literature also shows that geometry matters at several different levels. Incomplete metrics can invalidate perturbation queries unless the metric is modified; non-Euclidean metric spheres require careful direction sampling, since naive Gaussian-rescaling is biased; and curvature can enter the MSE bound explicitly, not only through global convergence constants. These facts make manifold zeroth-order estimation qualitatively different from simply projecting a Euclidean black-box method onto a constraint set (Ma et al., 12 Jan 2026).

A final related direction comes from Euclidean unbiased zeroth-order estimation. The telescoping constructions of “On the Optimal Construction of Unbiased Gradient Estimators for Zeroth-Order Optimization” develop exact unbiasedness from randomized finite-difference corrections in $\mathcal M$17. This suggests that an exactly unbiased Riemannian analogue might be built by combining tangent-space isotropy with geodesic or retraction-based telescoping differences, but that extension is not developed there and remains a plausible implication rather than a stated theorem (Ma et al., 22 Oct 2025).

Taken together, the current literature defines the zeroth-order Riemannian gradient estimator as a family of intrinsically motivated tangent surrogates whose design is governed by manifold geometry, oracle model, and problem structure. The classical form is a tangent-space finite difference composed with $\mathcal M$18 or $\mathcal M$19; the modern variants relax tangent sampling, exploit pullback or projection geometry, or leverage exact algebraic identities in structured objectives. The central technical question across all variants is the same: how to recover enough of $\mathcal M$20 from function evaluations to support first-order or second-order Riemannian optimization with provable complexity and feasible computation.