Riemannian Zeroth-Order Gradient Estimation with Structure-Preserving Metrics for Geodesically Incomplete Manifolds
Authors: Shaocong Ma
Heng Huang
Venue: ICLR 2026 Poster (2026)
Keywords: zeroth-order optimization
riemannian optimization
stochastic gradient descent
Primary area: optimization

Abstract: In this paper, we study Riemannian zeroth-order optimization in settings where the underlying Riemannian metric $g$ is geodesically incomplete, and the goal is to approximate stationary points with respect to this incomplete metric. To address this challenge, we construct structure-preserving metrics that are geodesically complete while ensuring that every stationary point under the new metric remains stationary under the original one. Building on this foundation, we revisit the classical symmetric two-point zeroth-order estimator and analyze its mean-squared error from an intrinsic perspective, depending only on the manifold’s geometry rather than any ambient embedding. Leveraging this intrinsic analysis, we establish convergence guarantees for stochastic gradient descent (SGD) with this intrinsic estimator. Under additional suitable conditions, an $\epsilon$-stationary point under the constructed metric $g'$ also corresponds to an $\epsilon$-stationary point under the original metric $g$, thereby matching the best-known complexity in the geodesically complete setting.  Empirical studies on synthetic problems confirm our theoretical findings, and experiments on a practical mesh optimization task demonstrate that our framework maintains stable convergence even in the absence of geodesic completeness.

Published as a conference paper at ICLR 2026
RIEMANNIAN ZEROTH-ORDER GRADIENT ESTIMA- TION WITH PRREERIGE-PRESSEVVING METRICS FOR GEODESICALLY INCOMPLETE MANIFOLDS
Shaocong Ma, Heng Huang" Department of Computer Science University of Maryland College Park, MD 20742, USA (scma 0908, heng) ) @umd. edu
ABSTRACT
In this paper, we study Riemannian zeroth-order optimization in settings where the underlying Riemannian metric g is geodesically incomplete, and the goal is to ap- proximate stationary points with respect to this incomplete metrie. To address this challenge, we construct structure-prestt prreserig- metries that are geodesically complete while ensuring that every stationary point under the metric remains station- ary under the original one. Building on this foundation, we revisit the classical symmetric two-point zeroth- estimator and analyze its mean-squared error from a purely intrinsic perspective, depending only on the manifold's ' geometry rather than any ambient embedding. Leveraging this intrinsic analysis, we establish convergence guarantees for stochastic gradient descent with this intrinsic estimator. Under additional suitable conditions, an e-stationary point under the constructed metric g' also corresponds to an e-stationary point under the original metric g, thereby matching the best-known complexity in the geodesically complete setting. Empirical studies on synthetic problems confirm our theoretical findings, and ex- periments on practical mesh optimization task demonstrate that our framework maintains stable convergence even in the absence of geodesic completeness.
INTRODUCTION
In this work, we consider the stochastic optimization problem on the smooth manifold M equipped with Riemannian metric g:
min f(p) (g==lf(p:E))
(I)
PEM
where (M, g) forms d-dimensional Riemannian manifold, the individual loss f(::)) - is smooth function depending on a random data point & drawn from a distribution E. The Riemannian metric g allows for defining the first-order gradient f in the tangent space at each E M, leading to the standard first-order Riemannian stochastic gradient method (Ring & Wirth, 2012: Bonnabel, 2013; Smith, 2014; Sato, 2021). In many practical scenarios, especially when the system incorporates non non-differentiable external solvers or black-boy objective functions especially when dealing with non-differentiable modules black-box objective functions, the explicit gradient of the objective function is either unavailable or prohibitively expensive to compute. This practical challenge necessitates the use of zeroth-order optimization technique to approximate the gradient direction solely using the function evaluation (Nesterov & Spokoiny, 2017; Li et al., 2023b), given by
{p(EE) ((x)))( - (=v)), U,
2/1
where v is random vector sampled from distribution over the tangent space T,MM and μ 0 the perturbation stepsize. The exponential map exp, : C T.M + M sends a tangent vector
"This work was partially supported by NSF IIS 2347592, 2348169, DBI 2405416, CCF 2348306, CNS 2347617, RISE 2536663.

---

Published as a conference paper at ICLR 2026
E T.MM to the manifold M along the geodesic starting at p, with B denoting an open ball centered at the origin in Tp.M In practice, the exponential map is often replaced by a first-order approximation known retraction (Definition B.3).
1.1 CHALLENGES IN RIEMANNIAN ZEROTH- ZEROTH-ORDER OPTIMIZATION
While existing analyses of Riemannian zeroth-order optimization establish convergence guarantees under various algorithms and assumptions (Chattopadhyay et al., 2015; Fong & Tiño, 2019; Wang et al., a, Wang & Feng, 2022; Maass et al., 2022; Nguyen & Balasubramanian, 2023; Li et al., 2023b;a; Wang. 2023; He al.. 2024; Wang et al., 2023; Goyens al., 2024; Zhou et al., 2025; Ochoa & Poveda, 2025), fundamental yet often overlooked issue arises from the local nature of the exponential map (or, more generally, retractions). In practice, Riemannian zeroth-order methods often endow M with an Eaclidean metric IE by viewing it as a submanifold of an ambient Euclidean space IR" and inheriting the metric from the embedding. This setting helps simplify numerical computations, but it has fundamental limitation: the inherited Euclidean metric gE may not be geodesically complete. Specifically, for a point p E M,, the exponential map exp, is not necessary globally defined over the entire tangent space T M.MM Consequently, randomly sampled tangent vector v E TpM may fall outside the domain of exp,, making (xp,(e) undefined. Theoretically. one could instead begin with a geodesically complete metric, under which the exponential map exp TM + .M is globally defined on the full tangent bundle T.M The Nomizu-Ozek theorem (Nomizu & Ozeki, 1961; Lee, 2018) guarantees the existence of such complete metric on any smooth manifold without boundary. Then by applying the Nash embedding theorem (Nash, 1956), one could, in principle, obtain an equivalent geodesically complete Euclidean metric, allowing direct application of existing convergence analyses. However, the constructive proof of Nash's 's theorem is numerically nontrivial, making it infeasible for practical optimization algorithms.
This challenge motivates us to consider the following natural question:
Q: How can we perform Riemannian zeroth-order optimization when the canonical Euclidean metric is geodesically incomplete?
To answer this question, we need to develop a Riemannian zeroth-order optimization algorithm for a given metric g that may not be geodesically complete, yet remains capable of finding a stationary point. Our contributions are outlined in the following subsection.
1.2 CONTRIBUTIONS
Contribution SSrracture-Preserving Metric Construction): To address the potential geodesic incompleteness of the given metric g. we construct the sructure-presserving metrics g' (Definition 2.5) in Theorem 2.6 that: (i) is geodesically complete, (ii) is conformally equivalent to the original metric g. and (iii) ensures any e-stationary point under g is also an e-stationary point under g'. These properties allow us to work with the new metric g while maintaining the desired property as the original metric g. However, adopting the structure- structure-preservee preserving metric raises fundamental challenge:
the geometry induced by g' generally differs from that of g. In particular, g' is typically no longer an Euclidean metric inherited from the original ambient Euclidean space, which precludes the direct use of standard Riemannian zeroth-order gradient estimators (Li et al., 2023a;b). Overcoming this mismatch between estimator design and underlying geometry leads to second contribution.
Contribution 2 (Intrinsic Zeroth-Order Gradient Estimation): Rather than finding a new ambient Euclidean space for the structure-pre preettrvip metric g', develop an intrinsic framework for zeroth- order optimization under non-Euclidean Riemannian metrics that relies solely on the manifold structure itself, and not on any embedding or representation in a larger ambient space. Under this intrinsic framework, we further analyze the meca--squaeed error (MSE) of the classical symmetric two-point zeroth-order gradient estimator (Fquation (2)) under an arbitrary geodesically complete metric g in Theorem 2.7, revealing the fundamental connection between the approximation error of gradient estimator and the curvature of the underlying manifold:
Euw-Uni( fp: d
1
+ nr?r?
+
where v ~ Und is uniformly drawn from the unit sphere gd-1 c TpM induced by g', S v) is the gradient estimator given by Equation (2), and is uniform upper bound on the
2

---

Published as a conference paper at ICLR 2026
absolute sectional curvature of (M.g'). In the flat approximation error for zeroth-order gradient estimation in Euclidean spaces. Building on this result, Theorem 2.9 establishes the convergence of SGD under general Riemannian metric g.
O, the bound reduces to the classical
Contribution 3 (Efficient Sampling under General Metrics): Moreover, sampling uniformly from the unit sphere S 1 C M. with respect to general Riemannian metrie g nontrivial. We show that the commonly used rescaling approach (i.e. drawing a Gaussian vector and normalizing it to g-unit length) introduces an inherent bias under non-Buclidean metrics. To overcome this issue, we apply the rejection sampling method (Devroye, 2006) to Algorithm 1, an unbiased sampling procedure for generating g g-unit-length tangent vectors. Proposition 2.8, we prove that the output distribution of our method is exactly uniform over Sd-1
Contribution 4 (Empirical Validation): Lastly, validate our theoretical results and demonstrate the empirical effectiveness of the proposed framework, we conduct extensive experiments on both synthetic and the practical experiments. Synthetic experiments examine: (i) the impact of sampling bias arising from rescaling sampling, and (ii) the influence of geometric curvature on estimation accuracy. the mesh optimization task, our method further shows practical effectiveness scenarios where geodesic completeness is absent.
1.3 APPLICATIONS
In this section, we highlight several applications of Riemannian zeroth-order optimization where the underlying manifold is geodesically incomplete.
Mesh Optimization In physical simulations, mesh optimization is essential for improving dis- cretized surface quality. Modern neural physical models, such as CFD-GCN (Belbute- Peres et al., 2020), adjust vertex positions by optimizing a quality metric, usually involving an external PDE solver. A major bottleneck is the requirement to implement auto-differentiaatio through this solver to obtain gradients, which is fundamentally difficult. Riemannian zeroth-order optimization offers compelling alternative by avoiding this gradient calculation. In this setting, the manifold consists of the valid configuration space of vertex positions. This manifold, however, is geodesically incomplete under the Euclidean metric, because configurations on the boundary (e.g., a vertex on an edge) are excluded to prevent numerical instability.
Irrigation System Layout Design This application seeks to optimize the physical coordinates of sprinklers to maximize water coverage. The coverage objective function is often a complex, non-differentiable simulation (e.g., modeling spray overlap, pressure, and wind), making it difficult to compute gradients. Riemannian zeroth-order optimization provides a gradient-free solution. The underlying manifold is the configuration space of valid sprinkler positions, defined by the open set within the field's ' boundaries. This manifold is geodesically incomplete, as typically we cannot directly put the sprinklers on the boundary of the field.
Covariance Matrix Estimation This fundamental problem in multivariate statistics and machine learning, essential for tasks like PCA and Gaussian modeling. The goal is to find matrix that best represents the data's covariance, often by minimizing a loss function (e-g.. maximizing likelihood). The underlying manifold is the set of all d x d positive definite matrices, denoted S A matrix C' is in this manifold if it is symmetric and x Ca > 0 for all non-zero vectors x E Rd This manifold geodesically incomplete because it is an open convex cone.
In summary, the incompleteness in these examples poses fundamental challenge, existing literature typically requires geodesic completeness for gradient estimation. This limitation motivates our work to develop framework that can perform Riemannian zeroth-order optimization without geodesic completeness.
2
MAIN RESULTS
In this section, present the main results of this paper: (i) We propose the concept of structure- preserving metric (Definition 2.5) and provide its construction based on an arbitrary given metric (Theorem 2.6). (ii) Then we derive the approximation upper bound of the two-point zeroth-order
3

---

Published as a conference paper at ICLR 2026
gradient estimator intrinsically, that is, it does not rely how the manifold is embedded into the ambient space (Theorem 2.7). (ii)) To numerically obtain the gradient estimator under a general Riemannian metric 9. adopt the rejection sampling algorithm (Algorithm 1) to sample from the g-unit sphere. Later, Proposition 2.8 guarantees that the sampled vector satisfies the desired property. (iv) In Theorem 2.9, we establish the convergence of SGD under general Riemannian metrie g-
Here we summarize the assumptions used our theoretical analysis. A brief manifold preliminary is included in Appendix B. Detailed discussions of each assumption are provided in Appendix C.I.
Assumption 2.1. In the optimization problem given by Equation (1), the individual loss function f(=E)= : M - IR satisfies the following two properties: (a) L-Bounded Hessian; for all p E M,, the Hessian matrix the point p is bounded by L. (b) Lower boundedness; the infimum f exists almost surely with E.
The following assumption imposes a regularization condition on the retraction used in Theorem 2.9. While is always possible to construct pathological retraction that deviates substantially from the exponential map, such choices may still scale with ||m||p but would negatively affect the final convergence rate.
Assumption 2.2. Let f M R be a smooth function. There exists constant CRet 2 0 such that
(S(Ret,(v) [(exp,(u)) CReellulle-
Assumption 2.3. There exist constants O and M3, M1 > 0 such that || f(2)( Iis S M3 and
V+f(2)) I M4 for all q E Bp(p, p), where Bp(p, p) denotes the geodesic ball of radius p and Iss is the Hilbert-Schmidi norm.
Assumption 2.4. There exists constant 2 0 such that the sectional curvature of the Riemannian
manifold (M.g) g) satisfies [Ko(o)] I K, for every point p E M and every 2-plane a C TpM.
Equivalently, -k Kp(o) for all p and
2.1 STRUCTURE- PURUCTURE-PREESERVING PRESERVING METRIC
We begin with the definition of stracture-preserving metric associated with given metric g. Since the exponential map of an arbitrary Riemannian metric y not necessarily globally defined on the entire tangent bundle TM (Proposition B.2), we seek an alternative metric g that is geodesi- cally complete while preserving the essential geometric behavior the original metric g'. This consideration motivates the following definition:
Definition 2.5. Let (M,g) be Riemannian manifold. A Riemannian metric g' is called structure- preserving with respect to if it satisfies:
(a) (Geodesic completeness) There exists p > 0 such that for any p E M,, the domain of the exponential map exp, Ty.M- + M contains the ball B,(p) == {o E TrMM ||l||a S p}.
(b) (Conformal equivalence) There exists a positive smooth function h M - R such that (v,,, w) h(p)g (v, w) for all p E M and all v, w E T,M.M
(c) (e-stationarity preservation) For any smooth function f . - IR and e > 0, every e-stationary point of f under g is also e-stationary point of f under g':
Here, we include brief discussion on the motivation for introducing each condition: (i) The first condition (geodesic completeness) ensures that if we set the perturbation stepsize μ < p and fix the random vector the g-unit sphere gd-1 C T,M, the perturbed point μv E T.M will always be within the domain the exponential map. (îi) The conformal equivalence condition preserves the set of stationary points; that is, for any smooth function f .M- R. if p is stationary point under g. then it also a stationary point under g', and vice versa. (iii) The e-stationarity preservation condition gives rise to the name ceationary-preeerving metric". It states that any e-stationary point under remains an e-stationary point under g', ensuring that the transformation leaves the original set of e-stationary points unchanged. We emphasize, however, that the converse need not hold: an
'A point p E M called an e-stationary point of the smooth function f under the Riemannian metric g the length of its gradient p less than e; that is, (V)(), f(()) <.
4

---

Published as a conference paper at ICLR 2026
e-stationary point under g' is generally not e-stationary point under g. Nevertheless, under suitable conditions, this asymmetry does not affect the overall complexity guarantees as we will discuss it in Corollary 2.10.
In the following theorem, we demonstrate that given a metric g, it is always possible to construct a metric g' which is structure- preservving with respect g.
Theorem 2.6. Let M be smooth manifold (possibly non-compact), and let be any Riemannian metric M. Then there exists Riemannian metric g' on M which is brructue-preserviing with respect to
Proof. The proof follows the classical construction presented by Nomizu & Ozeki (1961) with modifying the conformal coefficient h M -> (0, +o0) to ensure the e-stationarity preservation condition presented in Definition 2.5. The full proof provided in Appendix C.3.
As illustrated in Figure 1b, the metries constructed in this theorem ensure that geodesics remain within the manifold for all directions and lengths, eliminating concerns that random perturbations in zeroth-order gradient estimation could map outside the domain of the exponential map. Moreover, the conformal equivalence condition given by Definition 2.5 preserves the set of stationary points; therefore, in Riemannian zeroth-order optimization, it suffices work with the new metric g'.
(a) Euclidean Metric
(b) Structure Preserving Metrics
Geodesic contours centered at = (0.2, 0.2, 0.6) under the Euclidean metric (Figure 1a) and three pruuctur-ppreeevvin metrics (Figure 1b). Radii range from 0.1 to 0.9 in steps of 0.15. Under each structure- preserring metric, geodesics from p never exit the probability simplex, regardless of direction or length.
Figure 1:
Challenges Arising from the Setucture-Preservinn Metric Although Theorem 2.6 ensures that the constructed metric g' satisfies the desired properties, existing results in Riemannian zeroth-order optimization cannot be applied directly to establish convergence guarantees under g'. This limitation arises because much of the current literature assumes Euclidean setting, where M is embedded in Euclidean space and the gradient estimation determined by that embedding. In contrast, the new metric of is generally Eunculiida with respect to the original ambient Euclidean space of g. To address this obstacle, we are motivated to develop an intrinsic zeroth-order optimization framework that operates solely on the manifold's geometry, without requiring M be viewed subset of any Euclidean space.
2.2 INTRINSIC ZEROTH-ORDER GRADIENT ESTIMATION UNDER NON-EUCLIDEAN METRIC
In this section, we introduce the intrinsic approach to estimate the gradient of the function M R without relying the ambient space. We take geodesically complete metric and consider the classical symmetric estimator
V for=ffeeppp (we)), -f(eex,
24x
where exp, T..M + M is the exponential map. As noted by Bonnabel (2013), common to replace the exponential map with the retraction (Definition B.3).
The following theorem characterizes the mean mean-squared error (MSE) this zeroth-order gradient estimator, establishing connection between its approximation error and the intrinsic geometric
5

---

Published as a conference paper at ICLR 2026
(a) Rescaling Sampling
(b) Rejection Sampling (Algorithm 1)
Figure Illustration of sampling on the unit sphere induced by the non- non-Euclidean Riemannian metric g. The na've rescaling sampler (Left Panel) produces visibly non-uniform distribution, leading to biased estimator. Our rejection sampler (Right Panel) presented in Algorithm eliminates the bias and yields an even, truly uniform distribution.
properties of the underlying Riemannian manifold. The result is derived under the assumptions of bounded third- and fourth-order derivatives (Assumption 2.3) and globally bounded sectional curvature (Assumption 2.4). The full upper bound and the proof is deferred Appendix C.4.
Theorem 2.7. Let (M,9) g) be complete d-dimensional Riemannian manifold and E M. Let J:M R be smooth function and suppose that Assumptions 2.3 and 2.4 hold. Fix a perlurbation stepsize μ > 0 satisfying p² S min --1 + 6
d + J ), and for any unit vector E TpM define the symmetric zeroth-order estimator as in Equation (3). Then, for v Unif(Sd-1) uniformly sampled from the 9y-umil sphere in M,M,
1 122
|V
The bound Theorem 2.7 reveals how the estimation error connects the intrinsic geometry of the manifold. In particular, the sectional curvature term K quantifies the influence of local geometry on the estimator's variance. When & the curvature contribution disappears, and the bound reduces to the standard Euclidean variance expression.
2.3 SAMPLING FROM THE NON-EUCLIDEAN UNIT SPHERE
As the Riemannian metric defines a bilinear form on the tangent space T,M, uniformly sampling the g-unit sphere B := {0 E Tp.M (p v) 1} is equivalent to uniformly sample from the following compact set C {v E gd: v AA Ae 1} for some positive definite matrix A pdxd The matrix > 0 determined by the Riemannian metric and the choice of local coordinates: in practice, commonly use the local coordinate system spanned by the basis { a In this basis, the entries of are given by A.j == Sp 3 82s
Challenges in Sampling from the g-Unit Sphere In Euclidean space, sampling from the unit sphere is relatively straightforward: one can sample from the standard Gaussian distribution and rescale the vector to have unit length. However, this method does not extend to the g-unit sphere. As illustrated in Figure 2, rescaling- baased based sampling results in points being predominantly concentrated along the minor axes. To achieve a truly uniform distribution over the g-unit sphere, we adopt the rejection sampling method in Algorithm I that generates random vectors uniformly distributed over the compact set
The following proposition confirms that sampling strategy yields the desired properties; that is, the resulting output v exactly follows the uniform distribution over the unit sphere determined by the Riemannian metric g. The detailed proof is deferred Appendix C.5.
Proposition 2.8. Let the vector v be generated by Algorithm Then it follows the uniform distribution the compact set C = E IRd: o Au 1).

---

Published as a conference paper at ICLR 2026
Algorithm 1: Uniform Sampling on Ellipsoid {ar E Rdd A'Ar 1)
Input: positive definite matrix A Rdxd
Output: v E Rd
Q. + eig(A)
1/ Eigenvalue Decomposition: A=QAQ' A diag(Ai....... Aa)
L t -AQA- 11/2, 1-1-1122 Amax max(diag(A)
while True do
Draw z ~ [00, Ia); set +- 2/||=|1:0 +- Ls
// proposal point on C
Draw u ((0 1)
u < Vu" A Amax then
return
2.4 CONVERGENCE OF ZEROTH-ORDER SGD UNDER NON-EUCLIDEAN METRIC
In the previous section, have shown that the accuracy of Riemannian zeroth-order gradient estimator is improved as the underlying geometric structure selected to flatter. However, there is no free lunch in simply flattening the manifold. As have commented in Definition 2.5, the e-stationary point under the new metric g' may not be the c-stationary point under the original metric gi so one must balance estimator accuracy with optimization dynamics.
To solve the optimization problem in Equation (I), we employ the SGD algorithm. Starting from an initial parameter pt, the updates given by
Pe+1 Retp, (n V (P):
for / = 1,2,..., where Ret TM M is the retraction (Definition B.3), E R is the learning rate, (6)) is the stochastic data sample accessed at the t-th update, and (p;;)) is the Riemannian zeroth-order gradient estimation f(: at the point pt defined as
V5(p) ((Retp(w)) [(Retp(- V,
244
Now we build the convergence analysis of Riemannian SGD algorithm. We write S b if there exists a constant >> such that S Cb. The constant C may depend only on fixed problem parameters. Besides the boundedness assumption made in Theorem 2.7, we additionally require the L-smoothness (Assumption 2.1) and the regularization condition on the retraction (Assumption 2.2).
Theorem 2.9. Let (M.g) be geodesically complete d-dimensional Riemannian manifold. Let : M + R be a smooth function and suppose that Assumptiios 2.1 to 2.4 hold. Define the symmetric zeroth-order estimator as in Equation (5). Let {p.y-1 be the SGD dynamic finding the stationary point of Equation (1) generated by the update rule Equation (4) with requiring ni S v
and 112 S V (specified in Equation (21)), then there exists constants C1 C2, C3 >> such that
d
min 1S1T
[P)()
C1 np +2 + C3 pp?
In particular, choosing S + / yields min /V < V
Proof. The proof directly follows the standard convergence analysis of SGD in Euclidean space (Mishchenko et al., 2020). We may further relax the L-smoothness assumption to the expected smooth- ness condition proposed by (Khaled Richtárik, 2022). The zeroth-order gradient approximation eттor term is bounded using Theorem 2.7. See Appendix C.6 for the full proof.
Importantly, the upper bound in Theorem 2.9 not our final goal. We typically begin with a canonical Euclidean metric 9æ, which may fail to be geodesically complete. To overcome this issue, we construct a new metric g == hgg via Theorem 2.6 and then apply the convergence analysis under this metric (using Theorem 2.7 and Theorem 2.9). However, e-stationary point delivered by

---

Published as a conference paper at ICLR 2026
SGD under often is not e-stationary point under ge, unless the additional condition stated in the following corollary is imposed:
Corollary 2.10. Let gE the Euclidean metric, and let g be a structure-prexerving metric with respect to gE- Under the sume assumptions Theorem 2.9, suppose that either of the following conditions holds:
(a) 9E geodesically complete; or
(b) the set of e-statonary points under 9E, K { E E::M: [ S Seh,
is compuct.
Then it requires at most T < O (1) iterations to achieve minigist E [ S
Proof. Under either condition, the conformal coefficient h constructed in Theorem 2.6 admits uniform upper bound. Consequently, an e-stationary point with respect to the new metric g := hgg is also an e-stationary point with respect to the original metrie gE. up to constant scaling factor. This structure allows the complexity bound established Theorem 2.9 to transfer directly to the metric JE- See Appendix C.7 for the full proof.
Item (a) corresponds to the classical setting in which the original metric is geodesically complete. Item (b), on the other hand, specifies conditions under which an e-stationary point under the new metric is also an e-stationary point under the original metric. We emphasize that Theorem 2.9 establishes convergence even in more general scenarios, though with potentially worse complexity bounds than the geodesically complete case. This phenomenon highlights key distinction between the framework studied in our work and the traditional geodesically complete setting. Building on this result, we extend the best-known complexity bound for Riemannian zeroth-order SGD on smooth objectives from the special case of manifolds equipped with a Euclidean metric to a much broader class of manifolds endowed with general Riemannian metrics.
3 EXPERIMENTS
In the experimental section, we aim to validate the theoretical findings presented in Section 2. The two synthetic experiments are designed investigate the following questions: (i) How does sampling bias influence the convergence behavior of Riemannian zeroth-order SGD? (ii) How does the curvature of the underlying manifold affect the accuracy of gradient estimation? In addition, we conduct real-world experiment on mesh optimization (Hoppe et al., 1993; Belbute-Peree Pertee et al.. 2020; Ma et al., 2025), practical application in which the positions of nodes are naturally represented as points on the probability simplex.
3.1 SYNTHETIC EXPERIMENT: IMPACT SAMPLING BIAS
In this experiment, we investigate the impact of sampling bias in zeroth-order Riemannian optimiza- tion. Specifically, consider objective functions defined on the Euclidean space 12 equipped with non- Con-Euclidean Riemannian metric given by (u, v) u." Av:
fquadratic (2) = = 2B x T T (B+ E)x,
flogistic (z) = E(c.v) log(1 + exp(-y( x)) + 2 Bx,
where each entry of & is independently drawn from (000, 1), and (C,v) v) is sampled from fixed categorical data distribution. The matrix B E Iddd is a pre- rre-generated positive definite matrix. We compare two sampling strategies for Riemannian gradient estimation in the zeroth-order setting: (i) Rejection sampling (Algorithm 1), which produces uniform samples from the Riemannian unit sphere and is unbiased shown in Proposition 2.8. (ii) Rescaling sampling, which samples Gaussian vector then normalizes it to the unit sphere with respect to the Riemannian metric J.A.
Experimental Implications For each configuration, we report the average objective value over 16 independent runs using the same hyperparameter settings for the SGD optimizer. As shown in Figure 3, the rejection sampling method (Algorithm !) consistently outperforms the traditional rescaling approach; the rescaling method even leads to divergence under the same hypereparameter
8

---

Published as a conference paper at ICLR 2026
setting for the logistic loss objective (right panel of Figure 3). These results highlight the importance of using Algorithm to ensure an unbiased uniform distribution over the Riemannian g g-unit sphere, which is critical for stable and effective training. The complete experimental details are included in Appendix D.I.
0000
105
$730
5000
3000
23328
13"4
2000
0.724
Sectiondl canyeture Maa)
Sectionll conetuue Wed Lnifarm Dieribllinn
Sampling
Samsiing
Figure 3:
The impact of sampling bias on the Figure 4: The impact of sectional curvatures on convergence of Riemannian zeroth-order SGD,
the gradient estimation accuracy.
3.2 SYNTHETIC EXPERIMENT: IMPACT OF SECTIONAL CURVATURE
In this experiment, we investigate the impact of sectional curvature on the accuracy of zeroth-order gradient estimation. Specifically, evaluate gradient estimation errors at fixed point po under four conformally equivalent Riemannian metrics with different curvatures, We consider two objective functions commonly used in the optimization problem on probability simplex:
d
(p) = (L((pll)) = -En Pi log (dp:),
/Euctittean (p) = [p q112
D
2
2
where y L denotes the centroid of the simplex. We measure the accuracy of gradient estimation using the -quared error (MSE) under its own Riemannian metric, computed over 50,000 independent trials of zeroth-order gradient estimation (Equation (3)). The complete experimental details are included in Appendix D.2.
Experimental Implications As depicted in Figure 4, the Riemannian MSE of zeroth-order gradient estimation decreases as the sectional curvature K (po) decreases. This empirical finding aligns with our theoretical upper bound presented in Theorem 2.7, illustrating a clear connection between gradient estimation accuracy and the intrinsic geometric properties of the underlying manifold. particular, higher curvature consistently results in larger estimation errors for both objective functions.
3.3 GRADIENT-BASER MESH OPTIMIZATION
In modern physical simulation, solving PDEs often relies on finite-volume methods with spatial discretizations and external solvers that lack automatic differentiation support (Belbute-Peres et al., 2020; Ma et al., 2025), making the zeroth-order approach an ideal tool for optimizing mesh positions.
In this experiment, we consider the gradient- based mesh optimization problem for solving the Helmholtz. equation (Goodman, 2017; Engquist & Zhao, 2018),
v²f= = -
where 2 denotes the Laplace operator, k = 10 is the wave number, and f is the eigenfunction. The ground-truft solution computed on fine mesh with resolution 200 x 200. Our goal is to optimize the node positions of regular coarse mesh with resolution 20 x 20 so that its performance approximates that of the ground-truth solution. The mesh node (in our setting, boundary nodes are fixed and excluded from optimization) is represented using simplex formulation: each trainable node p = (x. y) is expressed as a convex combination of its six neighbors under the regular triangular initialization. This parameterization naturally leads to manifold optimization problem. However, the coordinate simplex, under its canonical embedding, is geodesically incomplete. To ensure the exponential map remains well-defined and prevent perturbed nodes from crossing mesh edges,

---

Published as a conference paper at ICLR 2026
JYv
Unconstrained
(a) Soft Projection
(b) Reversion
(c) Pructure-Preserving
Figure 5: The leftmost panel illustrates an invalid optimization step on mesh node; it crosses the edge, causing potential error in the external PDE solver. Figure (a) illustrates the Soft Projection approach, which resolves the issue by repeatedly reducing the perturbation stepsize μ along the perturbation direction until the movement becomes valid. Figure (b) shows the Reversion approach, which instead handles invalid steps by reverting to the original position. Figure (c) takes the advantage of the structure- structure-preservvie preserving metric, which twists the underlying Riemannian structure ensuring that the perturbation won't move the point out of the domain.
1e.1
Damerion
Starra heconura
Opeende
50ap
Figure 6: The left panel shows the ground-truth prediction (background), the initial mesh (blue), and the optimized mesh (red) using our proposed method. The nodes adaptively concentrate around the critical region while preserving the overall mesh structure. The right panel presents the loss curves for different approaches. Our method achieves both stable and efficient convergence.
adopt our proposed crtucture-preserviin approach and compare it against several natural baselines, as illustrated Figure 5.
Figure 6 presents the loss curves of the up-sampled prediction over 20,000 optimization steps. The unconstrained method often violates mesh validity, leading to unstable fluctuations, most notably around the 16,000th step. The reversion prevents invalid updates but quickly stalls after 8,000 steps; similarly, the soft projection stabilizes training but progresses slowly, showing little improvement beyond 14,000 steps. In contrast,
approach consistently reduces the error throughout training, achieving the lowest final MSE without instability. These findings highlight that structure-preserving approaches only maintain feasibility but also enable effective convergence.
truccure-presenvee
CONCLUSION
In this work, consider the zeroth-order optimization problem Riemannian manifolds when the underlying metric might be geodesically incomplete. We propose the structure-preservvigg metric that is geodesically complete, while preserving the original of stationary points (Theorem 2.6). Building on this foundation, we intrinsically derive the accuracy upper bound of the classical two-point gradient estimator and reveal the role of manifold curvature (Theorem 2.7). We further propose an unbiased rejection sampling scheme for generating perturbation directions under general Riemannian metrics (Proposition 2.8). Our theoretical analysis establishes convergence guarantees that extend the best-known complexity results beyond the Euclidean setting to broader class of Riemannian manifolds (Theorem 2.9). Empirical studies, including synthetic experiments and a mesh optimization task, demonstrate that structure- prw.uurrigreeeeigggggg approaches enable stable and effective convergence. These findings extend the theoretical understanding of zeroth-order optimization methods Riemannian manifolds and provide practical tools for Riemannian black-box optimization.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
P.-. Absil, C. G. Baker, and K. A. Gallivan. Trust-region methods on riemannian manifolds. Foun- dations of Computational Mathematics, 7(3):303-330, 2007. doi: 0000://1000000000077990
P.-A. Absil, Mahony, and R. Sepulchre. Optimization Algorithms Matrix Manifolds. Princeton University Press, Princeton, NJ, 2008.
Naman Agarwal, Nicolas Boumal, Brian Bullins, and Coralia Cartis. Adaptive regularization with cubics on manifolds. Mathematical Programming. 188:85-134, 2021.
Kwangjun Ahn and Suvrit Sra. From nesterov's estimate sequence to riemannian acceleration. In Jacob Abernethy and Shivani Agarwal (eds.), Proceedings of the 33rd Conference on Learning Theory, volume 125 of Proceedings of Machine Learning Research, pp. 1-35. PMLR, 2020.
Foivos Alimisis, Antonio Orvieto, Gary Becigneul, and Aurelien Lucchi. Momentum improves optimization on riemannian manifolds. In Proceedings of the 24th International Conference on Artificial Intelliencee and Statistics, volume 130 of Proceedings of Machine Learning Research, pp. 1979-1987. PMLR, 2021.
Gary Bécigneul and Sctavian-Eugee Ganea. Riemannian adaptive optimization methods. In Interna- tional Conference on Learning Representations (ICLR), 2019.
Filipe De Avila Belbute-Peree Thomas Economon, and Zico Kolter. Combining differentiable pde solvers and graph neural networks for fluid flow prediction. In international conference on machine learning, pp. 2402-2411. PMLR, 2020.
Silvère Bonnabel. Stochastic gradient descent on riemannian manifolds. IEEE Transactions on Automatic Control, 58(9):2217-222)) 2013. doi: 10. 20.10997AC.2211.22556111
Amit Chattopadhyayy. Suviseshamuthu E. Selvan, and Umberto Amato. A derivative-free riemannian powell's ' method, minimizing hartley-entropy-baaee ica contrast. IEEE Trans. Neural Networks and Learning Systems, 1993):98831195, 2015.
Christopher Criscitiello and Nicolas Boumal. Negative curvature obstructs acceleration for strongly geodesically convex optimization, even with exact first-order oracles. Proceedings of the 35th Conference Learning Theory, volume 178 of Proceedings of Machine Learning Research, pp. 1-47, 2022.
Lue Devroye.
Nonuniform random variate generation. Handbooks operations research and management science, 13:83-121, 2006.
Manfredo Perdigão do Carmo. Riemannian Geometry. Mathematics: Theory & Applications. Birkhäuser Boston, Inc., Boston, MA, 1992. ISBN 0-8176-3490-8, Translated from the second Portuguese edition by Francis Flaherty.
James Dugundji. Topology. Allyn and Bacon, Boston, 1966. ISBN 978-0-697-06889-7.
Björn Engquist and Hongkai Zhao. Approximate separability of the green's 's function of the helmholtz equation in the high frequency limit. Communications on Pure and Applied Mathematics, 71(11): 2220-2274 2018.
Herbert Federer. Geometric Measure Theory. Classics in Mathematics. Springer Berlin, Heidelberg, 1996. ISBN 000-3-540-6666--5 doi: 10.1007/978- 3-642-62010-3 org/10. Grundlehren der mathematischen Wissenschaften.
URL https: : /doi. 522-55222220000 2. Originally published volume 153 in the series:
Robert Simon Fong and Peter Tiño. Stochastic derivative-free optimization on riemannian manifolds.
arXiv preprint arXiv:1908.068 1908.0678, 2019.
Joseph W. Goodman. Introduction Fourier Optics. W. H. Freeman Macmillan Higher Education, New York, 4 edition, 2017. ISBN 978-1319119164,
Florentin Goyens and Clément W. Royer. Riemannian trust-region methods for strict saddle functions with complexity guarantees. Mathematical Programming, pp. 1-43, 2024.

---

Published as a conference paper at ICLR 2026
Florentin Goyens, Coralia Cartis, and Stéphane Chrétien. Registration of algebraic varieties using riemannian optimization. arXiv preprint arXiv:2401.08562, 2024.
Linus Hamilton and Ankur Moitra. A no-go theorem for robust acceleration in the hyperbolic plane. In Advances in Neural Information Processing Systems, volume 34, pp. 3330-334444 1230-3342,-0 2021.
Chang He, Zhaoye Pan, Xiao Wang, and Bo Jiang. Riemannian accelerated zeroth-order algorithm: Improved robustness and lower query complexity. In Proc. 41st Int. Conf. Machine Learning (ICML), 2024.
Heinz Hopf and Willi Rinow. Ueber den begriff der vollständigen differrentiaggoometrishhen fläche. Commentarti Mathematici Helvetici, 3(1):209-225, 1931, doi: 10.1007//F001001811
Hugues Hoppe, Tony DeRose, Tom Duchamp, John McDonald, and Werner Stuetzle. Mesh opti- mization. In Proceedings of the 20th annual conference Computer graphics and interactive techniques, pp. 19-26, 1993.
Wen Huang and Ke Wei. Riemannian proximal gradient methods. Mathematical Programming, 194 (1-2):371-413, 2022.
Wen Huang, P.-A. Absil, and Kyle A. Gallivan. A riemannian bfgs method without differentiated retraction for nonconvex optimization problems. SIAM Journal on Optimization, 28(1) ::4704499, 2018. doi: 10.1137/17MI 127582.
Leon Isserlis. On certain probable errors and correlation coefficients of multiple frequency dis- tributions with skew regression. Biometrika, 11(3):185-190, 1916. ISSN 0006-3444 doi: 10.10930 wt.1099/bomee///// 183.185. URL
185. JS-
TOR 2331846.
Leon Isserlis. On formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables. Biometrika, 12(1-2):134 (1111:3434139 1918. doi: 10. 10.10///biomet/ 12.1-2.134.
Ahmed Khaled and Peter Richtárik. Better theory for sgd in the nonconvex world. Transactions on Machine Learning Research, 2022.
Jungbin Kim and Insoon Yang. Accelerated gradient methods for geodesically convex optimiza- tion: Tractable algorithms and convergence analysis. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvári, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th In- ternational Conference on Machine Learning, volume 162 of Proceedings of Machine Learn- ing Research, pp. 11255-11282, Baltimore, Maryland, USA, 17-23 Jul 2022. PMLR. URL https: httpp//pceeeenn.ceeeeeee/oe mlr press /vlsss/vvle http:/vv22kkki/ kim22k html.
Lambert H. Koopmans. The Spectral Analysis of Time Series, volume 22 Probability and Mathe- matical Statistics. Academic Press, New York, 1974. ISBN 0124192505.
Vyacheslav Kungurtsev, Francesco Rinaldi, and Damiano Zeffiro. Retraction- direct search methods for derivative free riemannian optimization. Journal of Optimization Theory and Applica- tions, 2023. to appear.
Zhijian Lai and Akiko Yoshise. Riemannian interior point methods for constrained optimization manifolds. Journal of Optimization Theory and Applications, 201:433-469, 2024.
John M. Lee. Smooth Manifolds, pp. 1-29. Springer New York. New York, NY, 2003. ISBN 978- 0-387-21752-9, doi: 10. 1111111/97770--87--1117711
URL
John M. Lee. Introduction to Riemannian Manifolds, volume 176 of Graduate Texts in Mathematics.
Springer Cham, edition, 2018. ISBN 978-3-319-91754-2 doi: 10.1077/978-33319-915555
Jiaxiang Li. Krishnakumar Balasubramanian, and Shiqian Ma. Zeroth-order riemannian averaging
stochastic approximation algorithms. arXiv preprint arXiv:2309. 14906, 2023a.
12

---

Published as a conference paper at ICLR 2026
Jiaxiang Li, Krishnakumar Balasubramanian, and Shiqian Ma. Stochastic zeroth-order riemannian derivative estimation and optimization. Mathematicc of Operations Research, (832):113-112111 2023b.
Yuanyuan Liu, Fanhua Shang, James Cheng. Hong Cheng, and Licheng Jiao. Accelerated first-order methods for geodesically convex optimization on riemannian manifolds. Advances in Neural Information Processing Systems, 30, 2017.
Shaocong Ma and Heng Huang. On the optimal construction of unbiased gradient estimators for zeroth-order optimization. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025a. URL https: weppereeview.nee//oouu/a222
6ONt.
Shaocong Ma and Heng Huang. Revisiting zeroth-order optimization: Minimum-variance two- point estimators and directionally aligned perturbations. In The Thirteenth International Confer- ence on Learning Representations, 2025b. URL https:/ /oopnreeeiiieeeeeeee eeeieewweee nww.//e/// forum? ? yWFOSIT9ik.
Shaocong Ma, James Diffenderfer, Bhavya Kailkhura, and Yi Zhou. Deep learning of pde correction and mesh adaption without automatic differentiation. Machine Learning, 114(3):1-25, 2025.
Alejandro I. Maass, Chris Manzie, Dragan Nešić, Jonathan H. Manton, and Iman Shames. Tracking and regret bounds for online zeroth-order euclidean and riemannian optimization. SIAM Journal on Optimization, 32(2):445-44) 469, 2022.
Kanti V. Mardia and Peter E. Jupp. Directional Statisties. Wiley Series in Probability and Statistics. doi: 10.1002/9780470311697)
John Wiley & Sons, Chichester, UK, 1999 ISBN 9780471953333,
URL https: /oii www.ooog///..o 1002/ 97199977007777777777
Konstantin Mishchenko, Ahmed Khaled, and Peter Richtárik. Random reshuffling: Simple analysis with vast improvements. Advances in Neural Information Processing Systems, 33:17309-17320, 2020.
John Nash. The imbedding problem for riemannian manifolds. Annals of Mathematics, Second Series, 63(1):20-63, 1956. doi: 10.2307/1969989.
Yu Nesterov. Gradient methods for minimizing composite functions. Mathematical programming, 140(11) 2013a.
Yu.6. Nesterov. A method of solving a convex programming problem with convergence rate of F Dokl, Akad. Nauk SSSR, 264(3): 1983.
Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013b.
Yurii Nesterov and Vladimir Spokoiny. Random gradient-fre free minimization of convex functions. Foundations of Computational Mathematics, 17:527-566 2017.
Anthony Nguyen and Krishnakumar Balasubramanian. Stochastic zeroth-order functional constrained optimization: Oracle complexity and applications. INFORMS Journal on Optimization, 5(3):266- 272, 2023.
Katsumi Nomizu and Hideki Ozcki. The existence of complete ricmannian metrics. Proceedings of the American Mathematical Society, 12(6):889 891, 1961.
Daniel F. Ochoa and Jorge Poveda. Robust global optimization on smooth compact manifolds via hybrid gradient-free dynamics. Automatica, 171:11 11916, 2025.
Wolfgang Ring and Benedikt Wirth. Optimization methods riemannian manifolds and their
application to shape space. SIAM Journal on Optimization, 22(2):596-627, 2012.
Hiroyuki Sato. Riemannian Optinization and Its Applications. SpringerBriefs in Electrical and Computer Engineering. Springer Cham, 2021. ISBN 978-34030-66388-00 doi: 10.1007/ 978-3-030-62291-3 URL https: //tes:///oo...oo org /tt.111.100000 0000000777788-00-2533300
13

---

Published as a conference paper at ICLR 2026
Hiroyuki Sato. Riemannian conjugate gradient methods: General framework and specific algorithms with convergence analyses. SIAM Journal on Optimization, 32(4):2690-23) 2717, 2022. doi: 10.1137/ 21M1464178.
Hiroyuki Sato and Toshihiro Iwai. new, globally convergent riemannian conjugate gradient method. arXiv preprint, 2013.
Steven Thomas Smith. Optimization techniques on riemannian manifolds. arXiv:1407.5965, 2014.
arXiv preprint
Michael Spivak. Calculus. Publish or Perish, Houston, TX, edition, 1994, ISBN 978-0-914098-89- 8.
Tianyu Wang- On sharp stochastic zeroth-order hessian estimators over riemannian manifolds. Information and Inference: A Journal of the IMA, 12(2):787-813, 2023.
Tianyu Wang and Yasong Feng. Convergence rates of stochastic zeroth-order gradient descent for') ojasiewiez functions. arXiv preprint urXiv:2210. 16997, 2022.
Tianyu Wang, Yifeng Huang, and Didong Li. From the greene-wu convolution to gradient estimation riemannian manifolds. arXiv preprint arXiv:2108.07406, 2021a.
Xi Wang, Zhipeng Tu, Yiguang Hong, Yingyi Wu, and Guodong Shi. No-regret online learning riemannian manifolds. Advances in Neural Information Processing Systems, 34:28323-28335. 2021b.
Xi Wang. Zhipeng Tu, Yiguang Hong, Yingyi Wu, and Guodong Shi. Online optimization over riemannian manifolds. Journal of Machine Learning Research, 24(130):1-67, 2023.
Melanie Weber and Suvrit Sra. Riemannian optimization via frank-wolft methods. Mathematical Programming, 199:525 556, 2023.
Teng-Teng Yao, Zhi Zhao, Zheng-Jian Bai, and Xiao-Qing Jin. A riemannian derivative-free polak- ribière-polyak method for tangent vector field. Numerical Algorithms, 86(1):325-355, 2021.
Hongyi Zhang and Suvrit Sra. First-order methods for geodesically convex optimization. In Proceedings of the 29th Anmual Conference on Learning Theory, volume 49 of Proceedings of Machine Learning Research, pp. 1617-1638, 2016.
Hongyi Zhang and Suvrit Sra. An estimate sequence for geodesically convex optimization. In Proceedings of the 31st Conference Learning Theory, volume 75 of Proceedings of Machine Learning Research, pp. 1703-1723, 2018.
Hongyi Zhang, Sashank J. Reddi, and Suvrit Sra. Riemannian svrg: Fast stochastic optimization on riemannian manifolds. In Advances in Neural Information Processing Systems, volume 29, pp. 4592 4502-4600, 2016.
Juan Zhou, Kangkang Deng, Hongxia Wang, and Zheng Peng. Inexact riemannian gradient descent method for nonconvex optimization with strong convergence. Journal of Scientific Computing, 103(3):1-24, 2025.
14

---

Published as a conference paper at ICLR 2026
Appendix
Table of Contents
Related Literature
Optimization on Riemannian Manifolds
A.2 Riemannian Zeroth-Order Gradient Estimators
16
B Preliminaries
17
C Main Results
19
C.I Assumptions
19
C.2 Supporting Lemmas
20
C,3 Proof of Theorem 2.6
27
C.4 Proof of Theorem 2.7
28
C.5 Proof of Proposition 2.8
32
C.6 Proof of Theorem 2.9
33
C.7 Proof of Corollary 2.10
34
D Experimental Details
35
D.1 Synthetic Experiment: Impact of Sampling Bias
35 
D.2 Synthetic Experiment: MSE Curvature
(...............................
36
D.3 Gradient-Based Mesh Optimization
37
D.4 Additional Discussions Stability
38
E Limitations and Future Work
38
F The Use Large Language Models (LLMs)
38
A RELATED LITERATURE
A.1 OPTIMIZATION ON RIEMANNIAN MANIFOLDS
First-Order Methods Riemannian first-order optimization adapts gradient-based methods to Riemannian manifolds. For geodesically convex functions, Riemannian gradient descent enjoys convergence guarantees akin to Euclidean GD, with complexity O(L/c) for L-smooth objectives. Zhang & Sra (2016) established global complexity bounds Hadamard manifolds with curvature- dependent rates. Stochastic Riemannian gradient descent converges almost surely under standard assumptions (Bonnabel, 2013), while sariance-reducee variants such as R-SVRG (Zhang et al., 2016) and R R-SRG/SPIDER improve convergence for finite-sum problems. Adapting acceleration (Nesterov, 2013a; 1983; 2013b) to manifolds proved challenging due the absence of global linearity. Early methods (Liu al., 2017) were shown computationally impractical; Zhang & Sra (2018) and Ahn & Sra (2020) addressed this issue by controlling metric distortion, achieving accelerated rates under bounded curvature. Alimisis et al. (2021) proposed momentum-baan based RAGDsDR, while Kim & Yang (2022) achieved optimal accelerated rates with RNAG, matching the O O(VL/e) Euclidean complexity. There are still some fundamental limits remained: Hamilton & Moitra (2021) and Criscitiello & Boumal (2022) showed that curvature may prevent acceleration entirely negatively curved manifolds. These negative impacts would be eliminated using the second-order methods.
Second-Order Methods Riemannian second-order methods utilize curvature via Hessians and connections. Newton-type methods achieve quadratic local convergence using the Riemannian
15

---

Published as a conference paper at ICLR 2026
Hessian (Absil et al., 2008), though global convergence requires safeguards like line search or trust- region strategies. Trust-region methods (Absil et al., 2007) solve quadratic models in the tangent space and retract back, ensuring convergence second-order points. Recent improvements analyze their behavior near strict saddles (Goyens & Royer, 2024). Alternatively, Riemannian ARC (Agarwal et al., 2021) uses cubic regularization to achieve optimal O(e D(2=)/2) complexity. Quasi-Newton methods generalize BFGS to manifolds via vector transports. Ring & Wirth (2012) initiated this line, and Huang et al. (2018) showed global convergence (and superlinear rates) under mild assumptions. Limited-memory variants (R-LBFGS) scale better large problems. Overall, second-order methods offer faster local convergence but require careful geometric handling of Hessians and transports.
Zeroth-Order Methods When gradients are unavailable, zeroth-order methods estimate descent directions via sampling. Li et al. (2023b;a) applied Gaussian smoothing in tangent spaces using exponential maps to construct unbiased gradient estimators with variance bounds that depend on curvature and dimension. A stochastic zeroth-order Riemannian gradient descent achieves O(n/ e) convergence for smooth nonconvex functions. Wang al. (2023) proposed two-point bandit methods (2-2-BAN) for online geodesically convex optimization, showing regret bounds matching Euclidean rates up to curvature factors. Other derivative-fres free approaches include retraction-baser direct search methods, as in Kungurtsev et al. (2023), with convergence guarantees for smooth and nonsmooth objectives. Yao et al. (2021) developed a only function values and nonmonotone line search, achieving global convergence and hybridizing with Newton steps for improved performance.
conjugate gradient method using
Hybrid and Other Emerging Directions Several novel methods extend optimization frameworks to the Riemannian setting. Adaptive methods such as Riemannian Adagrad and Adam (Bécigneul & Ganea, 2019) address the challenge of accumulating gradients across varying tangent spaces by working on product manifolds, yielding convergence results for geodesically convex problems. Riemannian conjugate gradient (CG) methods, which define conjugacy across tangent spaces via vector transport, have been shown to converge globally under standard line-search assumptions (Sato & Iwai, 2013; Sato, 2022; Kim & Yang. 2022). Projection-free methods like Riemannian Frank-Wolfe avoid expensive retractions by solving a linear oracle at each step. Weber & Sra (2023) showed that Riemannian Frank- Wolfe method converges sublinearly in general and linearly under geodesic strong convexity. For composite objectives with nonsmooth regularizers, Riemannian proximal gradient methods offer convergence guarantees; Huang & Wei (2022) proved an (11/k) rate under retraction-basee baaed convexity. Finally, primal-dual interior-poin methods have also been adapted: Lai & Yoshise (2024) introduced Riemannian interior point algorithm with local superlinear convergence and global guarantees, mirroring the classical barrier method behavior in curved spaces.
A.2 RIEMANNIAN LEROTH-ORDER GRADIENT ESTIMATORS
In this section, we discuss several widely used gradient estimators in Riemannian optimization and highlight their connections to our work. Importantly, all of these estimators are developed under the assumption of a complete Riemannian manifold. In contrast, our setting differs from this convention by considering optimization over possibly geodesically incomplete Riemannian manifolds.
Wang et al. (2021b) This paper extends the one-point bandit estimator to homogeneous Hadamard manifolds. At the point x E M and given y uniformly sampled from the geodesic sphere centered at I with the radius ó, by using the gradient estimator
exp= (v)
f())) = f(y)
eexp (v)
this work established the best-possible regret rate O(73/1) for g-convex losses in the online regret optimization problem.
Wang et al. (2023) This journal version further develops a two-point bandit estimator on symmetric Hadamard manifolds. Uniformly draw y from the geodesic sphere centered at with the radius 8 and defined -y as the antipodal point of y. The gradient estimator is given by
f( := f(y) f(-w)
v) exp." (v)
2
expx (v)
16

---

Published as a conference paper at ICLR 2026
The regret improves ((VT) for g-convex and (log T) for strongly g-convex losses.
Li et al. (2023b) Maass Riemannian zeroth-order oracle for the online setting (Maass et al., 2022) and the expected loss setting (Li et al., 2023b). With a tangent perturbation E T.M (obtained by projecting an ambient Gaussian onto M), the gradient estimator is
(2022) These papers introduce con-symmetric two-point
f(()) = fo exp, (10) f(x) V.
Here have adjusted the estimator from Maass et al. (2022) to the time-invariant expected objective function setting to align with our problem setup. This estimator is the direct generalization of the one-side Gaussian smoothing estimator widely used in Euclidean zeroth-order optimization.
He et al. (2024) This work extends coordinate-wiss finite differences to manifolds. Using an orthonormal basis {ei} of TyM,, the deterministic coordinate-wil zeroth-order estimator is
(pe:) - f exp, Ex.
f((x)
2/u
In summary, compared to approaches that rely on projecting from the ambient Euclidean space, our analysis is purely intrinsic, that is, the gradient estimator depends only on Riemannian structure and is independent of any particular embedding. In contrast prior intrinsic estimators, which primarily focus on geodesically convex problems, our work addresses the non-convex setting. As result, our contributions extend the scope of existing research Riemannian zeroth-order optimization.
B PRELIMINARIES
In this section, review some basic definitions and results from Riemannian geometry that are used our analysis. For a full review, we refer the reader to some classical textbook (Lee, 2003; 2018). For convenience, we summarize our notations in Table 2.
Smooth Manifolds d-dimensional smooth manifold M is a second-countable Hausdorff topo- logical space such that at any point p M, there exists Up M,, a neighborhood of p. such that Up diffeomorphism to the Euclidean space IRd Let C (U) be all smooth functions over U C M. A deviation at E M is a linear mapping v - (U,) + IR satisfying
v(fg) vff) g(p) v(g) fp)
for all f,g E C (Up). Then the langent space at p. denoted by T,., is the real vector space of all deviation at p. The tangent bundle the disjoint union of all tangent spaces
TM {p. v) p M, T.M)..
smooth map f : - R" called an immersion if its differential df lplp: TpM + Tym R", defined by df lp (v) == o(f) for each e T,M, is an injective function at every E M; is called an embedding if is immersion and is also homeomorphic onto its image f(M) == (f(p) p M).
Riemannian Manifolds A d-dimensional Riemannian manifold (M,9) is d-dimensional smooth manifold equipped with Riemannian metric g, which assigns each point p E M an inner product
9p .T.M x M
where T denotes the tangent space p M. We also write ( represent Op and || l|p for the norm it induces. Let ó M + 1. be an embedding from the smooth manifold M to the Euclidean space R" Then M inherits Riemannian metric from the ambient Euclidean structure via the pullback metric
90 (v, u) = (d)),(v), (d)()) (u), = (o(v), ((u)),
where ( denotes the Euclidean inner product on R In this case, we say the metric g induced by the embedding d. and refer to I the ambient Euclidean space. To distinguish between Riemannian metrics that may be induced by embeddings into different ambient spaces, we introduce the following definition:
17

---

Published as a conference paper at ICLR 2026
Definition B.1 (n-Euclidean metric). A Riemannian metric g is called n-Euclidean if there exists smooth embedding Φ : M + IR" such that g is induced by Ø.
Notably, given an arbitrary d-dimensional Riemannian manifold (M. g). the Nash embedding theorem (Nash, 1956; Lee, 2018) states that there always exists n E N such that the Riemannian metric g n-Fuclidean. However, if we consider different Riemannian metric g' on the same manifold M,, there is no guarantee that g' can also be realized as an n-Euclidean metric for the same n. This observation motivates us to develop an intrinsic analysis framework that does not depend on any specific embedding.
Geodesic A vector field on MM is a smooth section X : + TM of the canonical Lungent-bundle projection æ TM - M: equivalently, it smooth map satisfying w X idy. Let (MM) be the space of all vector fields on Riemannian manifold (M.g). The Levi-Civita connection is the unique affine connection
V (MM) X(M) + X(M), (
satisfying torsion-free and metric-comppaible." Let I C R be an open interval containing 0. A smooth curve : + M is called a geodesic over I if its velocity vector 71(1) == dr la (xr E M satisfies the geodesic equation
(1)) =0
for all t 1. Given a point p E MM and an initial velocity v E T,MM, there always exists a unique (Theorem 4.10, Lee (2018)). The exponential map
geodesic 7 such that 7(0) p and 70(0) at p. denoted exp, T.M + M,, is defined by cxp,,() == =(1). Importantly, the existence of geodesic does not guarantee that ry can be defined over an open interval containing [0, 1]; that is, the exponential map be undefined for some (p, v) TM. We summarize this observation in the
following proposition:
Proposition B.2 (Proposition 5.7, Lee (2018)). The exponential map exp, TyM + M is locally defined an open neighbor of E T.MMM
Remark. This proposition reveals a fundamental difference between Riemannian and Fuclidean zeroth-order optimization: in the Riemannian setting, one cannot simply apply small perturbation in the direction v at the point E M, since the exponential map (xp) may be undefined. Developing a zeroth-order gradient estimator that operates within this local geometric structure is one of the central goals of our work.
Computing exp, (v) involves solving a differentiable equation, which often costly or intractable; hence, existing Riemannian optimization literature typically uses the first-order approximation called the retraction to approximate the exponential map.
Definition B.3 (Retraction). A retraction manifold M is a smooth map Ret TM .M such that for all E M:
1. (0) = p, where Tp.M the zero vector;
2. The differential dRetplo Tp.M- > Tp.M satisfies d'Rtplo = idr,MM-
Here, Ret, T,MM + MM denotes the restriction of to the tangent space at p. Intuitively, a retraction approximates exp,(v) by preserving the first-order geometry of geodesics while being easier to compute.
The following lemma further characterizes the relation between the exponential map and the retraction. We present it here without providing the proof.
Lemma B.4 (Theorem 2, Bonnabel (2013)). Let (M,g) be a smooth Riemannian manifold.
"W call an affine connection torsion-free if X V [X,Y], where the Lie bracket [X, [X is defined by [X, YJ( (((()) Y()()) for any C (M)), and metric-compattble if X (g((, 2))
9(VxY, 2) + s(Y, V Z) for all X,Y.Z.
MMore explicitly, we choose extension vector field X E X satisfying ((((1)) 71(()))))))))))" for all I. Then we define V / (t) == V X l1t)- Here directly use V ( for our convenience, as this definition does rely the choice extension (sec Lemma 4.9, Lee (2018)).
18

---

Published as a conference paper at ICLR 2026
(i) The exponential map exp TM + M is retraction.
(ii) For every p E M, the geodesic distance d(:::) M M 10, +00) berween exp,(1) and Ret,(U) is upper bounded as
d (xp, Relp(()) 2 $C11lo1
for any u and any retraction Ret.
Let the cotangent space T 7,M M be the dual space of M,M that is, the space of all linear mappings & T + Riemannian metric g:
Gradient
There a natural isomorphism between M and Tp M induced by the
bp Tp. + TTM M, e 4 gp(1, :);
1, M
w w wt satisfying 9p(wF, (0= w(v).
for all u M. Let f: - R be a smooth real-value functiio. The differential of f p e M. given by df (v) == vf. naturally defines covector in the cotangent space; that is, df T.M. The gradient of f. denoted by Vf E X(M). is vector field given by
p p Vf(p) = (dflp))
In this paper, investigate the approach of estimating V ( given only the access to the function evaluation. There have been a rich body of literature in this direction, which we summarize in Appendix A. In contrast, our approach is purely intrinsic, making result distinct from existing literature.
C MAIN RESULTS
C.I ASSUMPTIONS
The following assumption is standard in stochastic optimization literature (Mishchenko al., 2020; Khaled & Richtárik, 2022). In the context of Riemannian optimization, it is often coupled with Assumption 2.2 to define the L-smoothness of the pullback function (Bonnabel, 2013; Li al., 2023b; He et al., 2024). In contrast, we decouple these two assumptions to make their respective roles and dependencies more transparent.
Assumption C.1. In the optimization problem given by Equation (1), the individual loss function
[(:E) : M R
satisfies the following two properties:
(a) I.-Bounded Hessian; for all p E M, the Hessian matrix at the point p is bounded by L..
(b) Lower boundedness; the infimum ft == exists almost surely with 5 ~~ E.
The following assumption imposes regularization condition the retraction used in Theorem 2.9. While it is always possible to construct pathological retraction that deviates substantially from the exponential map, such choices may still scale with ||u||p but would negatively affect the final convergence rate.
Assumption C.2. Let f M + P. be a smooth function. There exists constant CRet 2 such that
(f(Ret,(v) f(exp, (u)) < CRei
Remark. This assumption can indeed be replaced with stronger but more widely used boundedness assumption (e.g. the bounded gradient assumption). Bonnabel (2013) has shown that the geodesic distance between the (first-order) retraction Ret,(e) and the exponential map exp,(u) is of the order o[[uu]]) (see Theorem 2, Bonnabel (2013)). In Lemma C,11, we show that given appropriate smoothness and boundedness conditions, the gap between f (Ret,, (v)) and (())p,(v)) is also the order of) ()). which implics Assumption 2.2. Here we present this weaker assumption to avoid introducing the bounded gradient assumption.
19

---

Published as a conference paper at ICLR 2026
Assumption C.3. There exist constants 0 and M3, M4 O such that
Hs S Ma. | S M4,
for all q E Bp(p, p). where (p p) denotes the geodesic ball of radius p and || ||Hs is the Hilbert- Schmidt norm.
Remark. Unlike the Euclidean setting, optimization Riemannian manifolds often relies on ad- ditional boundedness assumptions. For example, He et al. (2024) and Li et (2023b) impose Lipschitz continuity condition on the Hessian of the pullback objective (Assumption 4.2 in He et al. (2024), Assumption 2.2 in et al. (2023b)), which can be viewed as a variant of Assumption 2.3. The assumption of bounded fourth-order derivatives in Assumption 2.3 is less common in the literature. However, we emphasize that it plays crucial role our analysis: enables us to capture the depen- dence sectional curvature in the accuracy of zeroth-order gradient estimation (see Theorem 2.7). From our perspective, introducing this assumption leads to a novel and more refined result that has yet been explored in existing work.
Assumption C.4. There exists a constant K. 2 0 such that the sectional curvature of the Riemannian manifold gM.g) satisfies
[Kp(a) SK, for every point p E M and every 2-plane TpM.
Equivalently, -K S Ky(a) S x for all p and a.
Remark. Many existing literature (Wang et al., 2021b; 2023) also made assumptions on the sectional curvature (lower) boundedness. Here we present a slightly stronger assumption: we assume the sectional curvature uniformly bounded (i.e. both upper and lower boundedness). We note that this assumption has also been used in existing literature (see Assumption 1, Alimisis et al. (2021)).
C.2 SUPPORTING LEMMAS
The following lemma generalizes the expected smoothness widely used in non-convex optimization (Mishchenko et al., 2020; Khaled & Richtárik, 2022; Ma & Huang, 2025b;a).
Lemma C.5. Let fē = inf f(p:&) and f := inf, Be~af(:;E). Suppose that Assumption 2.1 is satisfied and l < +00. Then there exists A, B 2 0 such that for any E M.
E 13 < A((()) + B.
Proof. Let : F be the geodesic with 7() = p and 7/ = u (with the unit length). Then the composite of smooth mappings f E): : + R is smooth function. By applying the coercive inequality, i.e. Proposition (Mishchenko al. (2020)), to <)) (it is equivalent to apply the Taylor formula expanding fp: {) at the minima point P arg arg,in, f f(: E)):
(6(0) S(EE + 2T.
| [V(((())
where the inequality we apply the L-bounded Hessian assumption (Assumption 2.1). Let f ( (i.e. ft = infpem_ f (p: {)) be the minima of the individual loss f(;:);; then we obtain
E fx) < 2Lf(p) 2LEx-=ff
2L(f(p) f*1 2L Be-afel,
where f the minima of the objective loss function. Typically, we have
inf Eq-af(::) 2 Bç-EʃE EçL= inf f(-15) p
by Jensen's 's inequality using the convexity of the inf operator. Therefore, B 2 0. The proof is completed by defining
2L and B 2 Ec-
Lemma C.6. Let M be a smooth manifold. Then there exists smooth function p M + [0, +00) that is proper; iee. for every compact set C R, p (C') is compact in M.
20

---

Published as a conference paper at ICLR 2026
Proof. This result directly comes from Proposition 2.28 (Lee, 2003) and it can be directly generalized for arbitrary Hausdorff paracompact topological space, as for a Hausdorff space, the paracompactness is equivalent the existence of partitions of unity (Dugundji, 1966). Here present proof without using the partitions of unity.
By Proposition A.60 (Lee, 2003), the smooth manifold M admits exhaustion by compact sets*; that is, sequence of compact sets {K, ) in M, such that
K2 CJCK for all ;;
Uß=1 K2
For each j. we can always have a smooth function b : M - [0. such that 1b on Kj and
supp( C K J+1 This existence is guaranteed by Proposition 2.25 (Lee, 2003). Define a smooth function p M + [0, +00) by
p()) = Σ (1- (1-b)(p)). J=1
For any fixed there exists j with E Ky; as the result, there is at most finite entries in this series non-zero. The finite-sum of smooth functions is also smooth. Moreover, p ppppp-~((paaae ((-o0, c)) C K[a [o]++
which is compact. Since is always non-negative, it implies that p is proper.
Remark. If the manifold M compact (e.g., sphere), then every continuous function serves as an exhaustion function. This offers an alternative perspective on the aracture-preserving metric: for compact manifold, we do not need to worry about the exponential map sending points outside the manifold, as all metrics constructed in Theorem C.16 are automatically geodesically complete.
Lemma C.7 (Isserlis). Let (M,g) g) be d-dimensional smooth Riemannian manifold, E M,, and f: M > TR be a smooth functioon Denote by
aB = {u TyM: ||uس| = =1
the unit sphere in the tangent space. Write UnB) for the corresponding uniform probability measures. If v= (V1, U2...., vd) Unif(@B) then
n,
Evi Vig ... Ui
ad+22id (0422-2)))10012222 Σ Epaire
2a. (rrr)) 8ir 2
where P2 represents the set of all pairings of {1, 2,..,,..
2k} (i.e. all distinct ways of partitioning is the Kronecker delta.
0 i i j, {1. n. into pairs s)), and ō.j = 1 i= j.
Proof. This result is known as the generalization of Isserlis's theorem (Isserlis, 1916; 1918). Our presented version is taken from Wikipedia, which refers to Koopmans (1974); Mardia & Jupp (1999).
Lemma C.S. Let (M,g) be a d-dimensional Riemannian manifold. Assume there exists a constant 2 0 such that the sectional curvature satisfies
[Kp(a) S
for every point p E M and every 2-plane C. T.M.
Then, for every p M the Ricci curvature obeys the operator-norm bound
[R, v)
[Ricpllop
sup
|
< (d 1)
T U
'We always require the manifold to be second-coumtable and Hausdorff; and all topological spaces locally homomorphism the Fuclidean space are locally compact.
21

---

Published as a conference paper at ICLR 2026
Proof. Fix point p and non-zero vector v E T.M. Extend to an orthonormal basis {u/lollps e2,..., ea) of MpM. By the classical formula relating Ricci and sectional curvature,
(U, =)= = Σ (M)spanfu,e:))lpp
1=2
Taking absolute values and using |K| S K gives
(Ri v) & 1) () |el..
Dividing by ||ll and taking the supremum over all non-zero e yields Ricpllop $ (d - 1)H, as claimed.
Lemma C.9. Let L E Rdxd be invertible diffeomorphism defined
L: L: C := {v E R² v Av = 1},
= Ls,
where LL = Ia. Denote by dga-1 and de the (d - 1) dimienesil Hausdorff measures on Sd- and C, respectively. Then aga-1 I- is absolutely continuous w.r.f. ae. and
d(ogu-1 1-1
((s) = det Li -
doc
J(L1o)'
Proof. The result immediately follows Theorem 3.2.3 (Federer, 1996). Here the linear map J the (d - dimdinensil Jacobian of I. defined as
d-1 ((s) = Jo- = A dL(s)]lo.
where d((s) T,5á-1 - T.C is the differential of L, A is the wedge product, and || Ilm denotes standard operator norm || == sup | [)x) As a linear map, the wedge product gives
A Ad-1 dL(s) (det D)- Taking the norm yields
((s) = det LII((L --1sl.
Then it completes the proof.
Lemma C.10. Let y be a geodesic defined over the open interval I t O satisfying (i) 700) = p and (ii) (0) Let F I - R. be scalar function over defined
F(t) ::ex=,pppp.
exp, (7()).
Then the following relations hold:
(1) Ft( V[(()))) (V (VS(p),)), r)p
(2) (t) V=1(H(0))0(0))((0)))000000) v22)()))) (p)())
(3) (t) ("(((()1) F( v v3
Proof:
(1) As I- R, the chain rule gives
dF, = df o dry: T.R + TT(ty + T Ro)
We take 0 or E T.R. Then
F((t) = dF
df (()))
01
[v((((t)) (7())
where applies the isomorphism between T and TPM given by 2. When treating ((((t)) as element in TØ M through this isomorphism, we also write:
V/(7()))) (0) (V)((()) (+(())-
Here, ff(p) to represent that the gradient Vf(p) is understood as 1-form mapping from M, M to R. When t = 0, we immediately obtain F'(0) = (Vf(p), v), by 70) = p and 7((0)
22

---

Published as a conference paper at ICLR 2026
(2) The chain rule gives
TR P = (ro) + df(()(dP))) : x TRR + To
We take of T.P. Then
θ 0
=t) Fel
Ot' 0t )t== dl fx() (1)))(())) + (Vran
As (f M + fo R = R linear function, it always maps 0. By the property of geodesic, V y (t) 0, leading to
= (o
= v (6(0)))00)))00
Here, we directly take dL (T(0) (()((() as it has been a 2-form in M @ M. To align the same notation used in V, still use [= When O, we immediately obtain F*(0) (0) V f(p)(e, u| by using (() = p and 7(0)) = v.
(3) The chain rule gives
F²F. = dL Jo dr, dr) + 3d (dyr, drr)) + d5 d'
We take ou T.R. As Asy + M geodesic. the last two terms are zeros. Then
Ft) = d3 fx
= v (0,0))
Now the proof is completed.
Lemma C.11. Let f M - f is uniformly bounded by a constant G > for all p E M, then there exists a constant
R be smooth function. Suppose that Assumption 2.1 holds.
CRet 2 O such that
(f(ett(()) f(exp, (v)) S CRee|lullp
Proof. It suffices to apply the standard Taylor formula (Spivak, 1994) to both functions
f o Retp TrM = R + R and f foeepp RR + IR,
then evaluate their difference. We set (t) == exp, (tv) as the geodesic and TRee (t) := Retp(tx) (tv) as the first-order approximation of the geodesic. The Taylor formula gives
fo (x)xp(() f()) + (VI(p),)) + S (1 (exp))
I Relpel( = (()) + (V/(p), [ (((1 (1-1) 1-= ((Ret, (lw)) (1), Thee (ke (1)el + i,
where is the correction term reflecting the curvature from the approximated geodesic Niee, given by
-1- (L-t)(V((mell (f)), V the teet (f)) (m))
( (1
When Ret = exp. the Levi-Civita connection V X(M) x X(M) + x(M) automatically gives (tota(tt>> (tv) 0, which recovers the zero approximation error.
When considering general first-order retraction, we can further upper bound it using the bounded gradient assumption. Since the gradient (TRee (f)) uniformly bounded, we can also upper bound its directional derivative
V Thmm (1) Keet (4) ) here
this uniform upper bound as t.
14
(1- 1) (VV((neee (1), V Then (0) Thee (1)) (m)))mnaeeee ce
10
(1 (1-1) [V]((R((()) (Taa))))))))
V Thee maale) dt
23

---

Published as a conference paper at ICLR 2026
(>)
Ge 2 |101
-
where (i) applies the Cauchy-Schwand inequality, and (ii) applies the uniformly bounded gradient of f M + R and the uniformly bounded Hessian of Ret M + R.
We take the difference of the above two equations. The bounded Hessian assumption implies
|f (et,(u) - f exp
-
T.11011 2
2 1
2
Ge
$ (L + 2
|
we obtain the final upper bound by setting CRet L
2
Lemma C.12. Let (M. g) be smooth. d-dimensional, geodesically complete Riemannian manifold and let f: :M- + R. be a smooth function. Suppose that Assumption 2.1 and Assumption 2.2 hold. Given unit-length vector M and the perturbation stepsize μ >> 0, define
fo (wx) - expp(-wo)
ʃo Reep)(wo) I Reep(-uuo)
h
U,
2/1
2/1
Then
11 - h11p S CREEM.
Proof. We directly take the difference bewtween vectors:
f o exexp(uu)) f expp(-wv) fo Ret,(wu) f Retp(-HU0)
lh i
24
241
If exp)(1r) expp(-Hu)) - fo (et,(uw) + Retp(-ww)) /(24)
If fooep (AV) f Retp(ue) ((2p) + expp(--oo - f o Retp((uu)) (2p)
(10)
CREIM
where (i) applies that E MpMM is the unit-length, (ii) applies the triangle inequality, and (iii) applies Assumption 2.2.
Lemma C.13. Let (M,g) be a smooth, d-dimensional, geodesically complete Riemannian manifold and let f: M + IR be smooth function. Suppose that Assumptions 2.1 to 2.4 hold. Suppose that there exists constant CRet 2 such that
(f(Ret,(u) f(exp,(u)) CReillu1I?.
(6)
Let {pe} be the SGD dynamic solving Equation (1) generated by the update rule Equation (4). Then
[(pp) [1+ (6CRet + L 1/2 d
2 x²x²
)n²
Lud'
x2
[Ef(p))) -s*)
(d+
2 2 122
(Ef(p(+)) - -f¹) + (ORee + [ (3B
36 -3 )
d
bd 5 +
3 &dmn Ceet
where & and F are given by Equation (7) and Equation (10), respectively.
Proof. Let h (p(itt)) [(Retp (wv)) RR(tt (-av)) Tp.. (also defined in Equa- tion (3)), hr (expp, (ww)) ( (-aw)) and h4 [V(pprEE) E Tr.MM At the t-th update,
24
2u
the SGD update rule (Equation (4)) gives
Pt+1 Ret,, (-njha)-
24

---

Published as a conference paper at ICLR 2026
Let - .M be the geodesic over I D [0, 1] that satisfies 7(0) pa with the initial velocity 70(0) -nnh. The Taylor formula of the scalar function y gives
f(p(+1) = fRet,. (-nhe)) - - fexxp, (-nh+)) - (expp, (-nh+))
(1)
x "nean" |l 1? + f((:)) (VV)()))) hn he + he)ps
+ S." (1- (1-1))2)(00)))
(0)
x Cne 1lleell + f(p)) - -VV - -hx)p ((De))))) + Inf² 11
Eptt() S (CRet + I ")n²Te 1l + f(p)) - (
(iit)
(V{(p))) E,..i hu - he)ms + 1/ //l Byn l1, - hellpe
< (CRet + 2 |n heel. + +((p) 2 (ptt
L
1. + 2 || he - hell?
ad
6d
3 + dnnp CRet
where (i) applies Equation (6) and the Taylor formula, (ii) applies the bounded Hessian assumption (Assumption 2.1), and (iii) takes the expectation conditional on p: on both sides; here we use Ep, [-] to represent E pr] for convenience. The last step applies 2(au, 2 I a² [|]| + 2 || for >> 0 and Lemma C.12. Then it suffices to upper bound the variance term Ep, [i ll and the bias term [,, h. hell?.
Bounding pp [l I| First, we split it following the standard routine,
Ep. helll. Ep. lh he hr he + hell?
I 3Ep Il, - helli, + 3Ep. 11b - hell + 3Bp. I
The first term 3ppp //, he is given by Lemma C.12:
3Ep li hell?, S 3CRialli
9 The second term 3E,, ll, hel is given by Theorem C.16:
1 122²
3Ep, lh, he s 3
Ep, Vp.(( || 36
33 1 + 1rn²
(f(p)) - f*) + B] 38
d
where the second inequality applies Lemma C.5 and
[ M
M3
e =
11²
(7)
3 d3
288
The last term is upper bounded by Lemma C.5:
3
3E hhell
Te LE [
61
3B
5 d2
(f(pt) - fr1
d²
Putting all together, we obtain
Ey, |//
< 3 + k22x
6L
3B
[2L (fL()) - f*) + B] +35 + 302m//2
(f(p)) - f-++
d
d
25

---

Published as a conference paper at ICLR 2026
n2262 = =6L(
1 1+ p2
1
+ -
+ (x)
) +6
36 + 3C Rtt//
d
d
d
As d 2 1, we obtain
Ey, Whell, S 66 p2
2 nr
d [[(())
- - 3B
36 3C Reell
(8)
d
where & is given by Equation (7) and B given by Lemma C.5.
Bounding Ep - h4 l Following the same proof as Theorem C.16, we obtain the expansion of the zeroth-order gradient estimator given by Equation (14). We multiply on both sides and take the expectation:
122
11² E E[ (P())(p))) + 3Ric(:,)) +
Ep, h4 - h4 6d(d + 2)
I-
12
Then we take the squared norm to obtain the bias upper bound:
14t
11 144
[B. hr hellp -
|| f(pu: Ee)l + 99|Ricc(, (fFon
+
M2
(d + 2)
(9)
MM
11
11a
M2 +
(d (d+ [V
9d2 + 2)
144
pt MM
11
+
M2 +
(d + 2)2 2/² n²
[2L (J(pt) - -5 B)
(d 2)
144
2Lp
M
1B
, (d 2)
(f(p)) f*) 9d2 (d 2)
M2
k2
B,
(d 2)²
144
(9)
where (i) applics Lemma C.8 to upper bound the Ricci curvature by the sectional curvature. For convenience, set
1 MM
ld + 2)2
+ 1r M2
144
+
14" (d+
2)z
x 2
(10)
Combine Equation (8) and Equation (9), we obtain that
[V(D)) s 6L + " 2 x2xx )np + (d+ (d+ 2)² 1221) (Bf(p))- -5²)
Ludd
1 6d
d
(Ef(p(+)) - f*) + (ORee + (3B(
2 x2
+ 36 3C )m²
d
bd +
3 ddnpi C2 Ret
where 6 and & are given by Equation (7) and Equation (10), respectively.
Lemma C.14. Suppose that S 2 0. Let three real-valued sequences (01) ( and (G.) satisfy
0, < (1+ 8,41 + Gt,
for all 1 ST... Then the iterate is bounded by
S)1+S)T
min 0, < (1+ -1 61 max G, -
81 + max G..
ISST
I8ITT
T 1StTT
Proof. We telescope the iterative relation by using
0 < (1 + S)8 8141 + G7
(1 + 5) x 0r < (1 + 2) (1 S)65) (1 + S)GT-1
26

---

Published as a conference paper at ICLR 2026
(1+ (1+ x 01 $ (1 (1+ -1+ + (1+ ()+ G2.
We sum them together and obtain
T-1
T-1 Σ 1=0
Σ
(1=S) 5)
min 02 S (1+ + 18iST
(1 + 5)
max G.. 15t5T
2
Then we re-arrange the above inequality and obtain
min 8, - Ists
81 + max G2
- 5)*
ISES
5(1 +5²
81 + max G2
(1+ S)T -1
1S1S
ST
51 + max G..
T
IStST
where (i) applies two inequalities (1+ x) < T and (1 + (1xx)² - 12 Tx.
C.3 PROOF OF THEOREM 2.6
Theorem C.15. Let M be a smooth manifold (possibly non-compact), and let g be any Riemannian metric M. Then there exists Riemannian metric g' on M which is structure-presse pprettur-presevving with respect to
Proof. In this proof, we distinguish the norms induced by different Riemannian metrics by explicitly writing || llp.g or || llp..g.. Elsewhere in the paper, we simply use llp. as no alternative metric is under consideration.
We mainly follow the construction given by Nomizu & Ozeki (1961) to obtain conformally equivalent Riemannian metric which is geodesically complete. By Lemma C.6, there exists smooth proper function p .M- + [0, +oo). Define the conformal coefficient h M (0, +oo)
h(p) = (()) + +1)
where Vp(p) E M is the gradient of at M and 9 2 1. Then we define the conformal metric g' as
( w) = h(p)g,
Now turn to prove that (M, g") complete metric space; that is, every Cauchy sequence is convergent. Let [a, b] + M be a piecewise smooth curve segment, Then the length of with respect the metric is given by
(7) = f"
()), dt
(M006))9n00(()()))()))))) dt
(( dt
((V + 1) dt
(10 dt
d (t)) dt
27

---

Published as a conference paper at ICLR 2026
J. (7 dd
p((((b) p(>(a))]:
where (i) applies the definition of h, and (ii) applies the Cauchy-Schwarz inequality. As result, for arbitrary p. q E M, we have
lp(p) p(0) S dg (p, q)-
(11)
Let {ps} C M be Cauchy sequence with respect to g'. Then Equation (11) implies that {p(p))} C must be a Cauchy sequence. We can take a finite supremum
c= s: p(p)) < +00.
Then {p&} C p ([0,c)): that is, every Cauchy sequence belongs to a compact set by our construction (Lemma C.6), which implies the completeness of (M. g').
The Hopf-Rinow theorem (Hopf & Rinow, 1931; do Carmo, 1992) states that for a connected Riemannian manifold, geodesic completeness is equivalent the metric completeness, As we have shown that the conformally equivalent metric o/ == h(p) induces complete metric space, it automatically makes (M. g') geodesically complete manifold. If ( is not connected, this argument applies each connected component, and a geodesic contained within single component. Thus, (M,g') is geodesically complete.
Lastly, we show that if the e-stationary point under y also gives an e-stationary point under g'. Recall that we always have
Ip (VVf(p)).) = = f(p).v)
for all v E M. Then
VV u) 9p ( (p))
As it holds for all and 9p a bilinear form over the linear space M, we obtain
V f(p) = vpf(p).
Suppose that |V /g((l)) S then
fgf(pl)p
0 V f(p). f())
11/llp) V 90 (V,(()) v,(())
v1/h(p) V vof(l))ln
V Vaf(p)lpp
V (IVel))) + 1) (p(()) Se.
Therefore, complete the proof.
C.4 PROOF OF THEOREM 2.7
In this subsection, we provide the proof for Theorem 2.7.
Theorem C.16. Let (M,g) be a complete d-dimensional Riemannian manifold and p E M. Let f + R be smooth smooth,,unction and suppose that Assumptions 2.3 and 2.4 hold. Fix a perturbation stepsize μ > satisfying
1 < min d d'1'2
de).
for any unit vector T define the symmetric zeroth- estimator
f(p (exp, - (exp,(-pux))
24c
Then, for v Unif(gd-1) )) uniformly sampled from the Ip gp-unit sphere in T.M,
Enn-Uiif(2 1) () 1 {(m)
< + n2
[4 M3
M2
112
3 d3
288
28

---

Published as a conference paper at ICLR 2026
Proof. Let ((t) expy(tv) be the geodesic; it satisfies (i) (0) = p and (ii) Y'(0) For the scalar function Ftt) == f()(t)), we apply the ordinary Taylor theorem (with the integral remainder) at t = 0 up to order 4 (Spivak, 1994; Bonnabel, 2013):
F(n)= = F(0) + (0) + 1 (0)
(0)
1- (1r ())))"(()) dt.
1t
1
2
0
By applying Lemma C.10, we obtain
f(r(w)) = f(p) + (Vf(p) d)pp + 12 (lf())(,, 0) 1
((p)(
(12)
2
4ra +
/ (1 -
(1'(tt), (ut), 7 (ut), (ut))
6
0
d)p + 1² (2 v) 1" ('((()
((((-1) = (()) (VV(())
(13)
2
+ rra 6
/ (1-1)
(-
where the k-th covariant derivative at E .M is a symmetric k-lineer form in T. X ... M
copics
v f(p) TTM TT.M-----MMM-
copies
and we represent the remainder term given by the Taylor theorem
I+ = f (1 (((E)))((((
(ut), ((t))
Subtracting Equation (13) from Equation (12) and dividing by 2p we obtain
((xx))(w)) ((expp(-uux) (Vf(p). + 2² ((pl(t,,))
12
(I+ (I -I-). (I+- (14)
24
12
Multiplying on both sides, we obtain
ff 1 - ((Vf(p) u)pu d vf(p)) +p² v3 f()
1² [I_ - -1- J.
d
12
=o(vv
=:Za(v)
=:((v)
By defining these shorthand notations, we have the following compact form:
(p; V(() Zo(v) 11) R(e).
We take squared-norm on both sides and treating as the uniform distribution over the g-unit sphere gd-1 1 in M. Then we obtain
E. &.I - {(
+2,0::/-8888223220000--0222-2200 R(()) [E(13)+00113--131222++++++++11+1 22) Z(()), (1 (1++))).]22(00] 12 2E, R())
s ||
29

---

Published as a conference paper at ICLR 2026
The cross term (Zo(e), R(u))p is canceled out by Lemma C.7. More explicitly, we have
E(Zo()),R(()) EE" 12 (I+ - -I-)
-E"
(VVF()),) 0)) ((()).))
12
n (I+ -I_ (0 ( = 0,
12
where (i) applies Lemma C.7. Now suffices to bound each squared term.
1. Bounding { By Assumption 2.3 and
||0|| we have
M4 / (1 (1- 1-dtt 4
I+ [+ (1, w) -
We have the similar upper bound for [T_]. Then I__ - I I+
M1
M4
Ma.
4
2
As the result,
nr² M4 Map
|
<
p
12
2
24
Therefore, we obtain
M2ne
E. [1 S
(15)
576
2. Bounding E.]lZo(e)]] Recall that Zo(v) = (VI()), u)nu (v) Then
|Zo(v)] - - d) 1 à (p)( 2 d V
(Vf(p), (0) v) 7 vf(p))
d V f (V))))))))))
2 (V((p) ())
(VI(p). 1,
dc [)){))
2 2 (V
where (i) applies ||01| = 1 and (VV(pp),) = (Vf(p), By the symmetry of the
|| llm-norm ball, we have
1
[ & v] = -9p-
where & T.M T.M + R is the tensor product of the vector with itself and @ o (V))), VS(p)) = (v, {(())) As the result,
E. (v(()() ") = (VF6)) 1
Therefore, we have
1 Eg E. 2 d )
(16)
3. Bounding E. | [2 We choose an orthonormal frame fei,..., Ca} for TpM so that every vector E T.M with ||e||p = 1 is represented as
e,
and write its coordinate (v', 02
dd ERR. AS VD TT
T M, we write the tensor representation as
Tisk == (V3 De)(p).
Therefore, we obtain
3 f(p)(,, 1 I
Zu(v) =
6
6
30

---

Published as a conference paper at ICLR 2026
where use Einstein notation to represent the sum. By the orthonormal frame, obtain
(
1
To of of v
36
Then it suffices calculate E,, of of J. By Lemma C.7, obtain
if (i,j, -
(d++2]a+T) a(a+))d44))
By Bulu'v' v of uM
if i= j, J', otherwise
As the result, we obtain
1
E, [122(11111
6T (T Tizk: + ihx Tzak
36d ( (d 4)
Recall that Tisk = V f(p)
1 HS'
We also have
Tik
((A f) + Ri
< / A /) + -(
- 2 vス f(p + +
As the result, we obtain
E. Ez <
[4]]" f(D) + +3]]
(17)
6d(d + 2) + 4)
Combining Equations (15) to (17), obtain
B. || v) - [
1 1) (1 +
2/²² + 2
(4)) +3|
d
6d () 2)(d + 4)
Mā
288
(
22 ) 1 + 1²) |[ 6d(d + (d + 4)
2μ²² n2
(4M3 + 3322d (( + Mips
288
( + (1 + n2²) + 2d
2μ + n²
24 + 41²
M2
4M² +
6d + 2)(d + 4)
6d + 2)(d + 4)
288
where (i) applies Lemma C.8 and assumptions 2.3 and 2.4. Furthermore, we set
2μ + n
x2
3x2
6d(d + 2) + 4)
d
It solves
6
8
n <
+
+
(18)
2 d
We also let
212
< d
(19)
d-1
We obtain (1 - )1+ (1++2)
3. It concludes that
E, 1 1+ p222
p2k
4 3
M3
||
+
d
d
288
Then the proof is completed. Combining Equations (18) and (19) leads to the range of p.
31

---

Published as a conference paper at ICLR 2026
C.5 PROOF OF PROPOSITION 2.8
Proposition C.17. Let the vector v be generated by Algorithm 1. Then it follows the uniform distribution over the compact set == {v R u Au 1}.
Proof. Fix a positive definite matrix A E Iddd and consider its eigenvlue decomposition
A = QAQ", A = diag(A diag(A)...... da), 0< S ... X Amax
Recall that L = QA 1/22, Then
-1/2
d
det L det Q det A-122
A
>
We observe that for every Sd-1
(Ls) ALs 1. ALs=
indicates that Ls E C (u:v Av 1}. As the result, L defines smooth bijection linear map from the sphere gd-1 to the compact set C:
Au
L L: S Ls.
Under this notation, Mprop- the distribution of the sampled vector (without rejection) in Algorithm 1 is given by the push-forward distribution of the uniform distribution via the linear map L. That is, any measurable E' C.
Perop E)) == |
pşa-1 LE),
(20)
where jia-1 is the uniform distribution over the sphere Sd-1
Denote by 0ga-1 and de the Hausdorff measures on Sd-1 and C, respectively. Then we re-write the above distribution Mprup and μgu in the density form; that is
Hprop Ppoopdac.
Aga
Pai do 1.
For arbitrary integral function g C we have
1. g(()) (()))dpppp())
o ga-1
g(L s)d (L
i
g(L s(Ls)duuge- o (
ga g(L sddpas-
where (i) applies the definition of the pull-back measure Ppop (Equation (20)). Then obtain
f (000)pppp(
f Jg ((Ls))ss (s)dogs-1
[ g(L s) PSa-1 (s)
doc ().).
j(s)
where (i) applies Lemma C.9 with J(s) = |det L [|(L - As it holds for all measurable
function g, it solves the density of Uprop
Pga-1 I I-)(
Pprop (v)
J o L-1
1
| 12
32

---

Published as a conference paper at ICLR 2026
Then consider the rejection step and the final density. Let Pout be the density of the output vector of Algorithm 1. Recall that Algorithm accepts the candidate v = L s with probability
c. A2 Iu)
c A
a(v) == P(aceept v/e) P(u < V
V
Amax
Amax
The density of the output vector is given as
Pour( x Pprop (v)o(())(())
Vrmax
As it is a constant over the compact set C, it is the uniform distribution over C. We also note that the acceptance probability is strictly positive; hence, the loop halts almost surely. This completes the proof of Proposition 2.8.
C.6 PROOF OF THEOREM 2.9
In this section, we present the proof of Theorem 2.9. We write a S if there exists a constant such that C. The hidden constant C may depend only on fixed problem parameters.
Theorem C.18. Let (M.g) g) be a complete d-dimensional Riemannian manifold. Let f : + R be a smooth function and suppose that Assumptions 2.1 to 2.4 hold. Define the symmetric zeroth- order estimator as in Equation (5). Let {pt) be the SGD dynamic finding the stationary point of Equation (1) generated by the update rule Equation (4) with requiring n S / and μ²² z V
(explicitly specified in Equation (21)), then there exists constants C1, C2. C3 > 0 such that
d nT
| S C1 SSST
min
+C2 C3 pp
In particular, choosing μ de / 4 yields
d
min IVf Z VVT
T
Proof. By Lemma C.13, we obtain that
(I + L 2 x²x² 2 )0² +
Lptdd
6d [VS()))
x2 (Bf(p)) -1*)
d
(d (d+ 2)2
(Ef(pe+)) - 5*) + (CRoz + 1 (3B(
2 12x²2
d + 38 =3C 3Cheeli))
2
3 bd I + ddmu@Reet: Ret:
It has the same structure presented in Lemma C.14, where we set
0, 2 VV (mm) S= G
2 " 2 p2:2 )m? + (d (d+2) 2)2
2
Lpdd
k 2211
6, = (==((pp)) - -
6d
2+ n2 G. (CRett + 1. / (3B((
+ 36 + 3CReell 17 nd g
3
4 dnp 2
d
Then we obtain
5
min 0 <
81 + max Gr. 1SEST
ISST
T
leads to
(3) 6c Ef(p) - -
6d 1dd
min SIST
3r
+
dnp 2 C2 ciat
rT/dd
2
4
6d + 77
(CRet
+
(3B(
2 + 1 xp
+ 38 + 3cRlll?)) n)²
d
33

---

Published as a conference paper at ICLR 2026
where (i) selects
11 T
18L(C Ret + 2)
(21)
112 S min
T
18L2 + 2)
such that TT S e², where 2 is considered as too when x: 0. Given Equation (21), we further upper bound it
min 181ST
[V=
< T/T
(B(()) - [r]] + 3d + dddu²² PC
(22)
+ (doe + (3B( 2+ 222
3B
38-
Retll?
d
( d
[6e (p1) [m]] +3d
n M3
116
M3
1¹a
x22
B]
de
9
1/T
(d 2)²
144
(d+ 2)²
2
L ) 3BC
+ x22
4 MM
M²
+ Gdot
+pp²
+3
3
d
3 d3
288
M3 n²
M2
$ T/T
[54[E((1)) - 1 +
+
n1 + +.1 312
5d 1ef/1"0 C
3 d
48
+ 54ne + 1. [)Bn 18(CRet Ret +
T,
4 M3
Mf
T. dnpp 18(Cnet 2 C
3 d3
288
(23)
-0 1/7 ) + O()
where (i) applies the formula of & and & given by Equation (7) and Equation (10), respectively.
C.7 PROOF COROLLARY 2.10
We re-state this corollary to have consistent notation previous sections.
Corollary C.19. Let g be the Euclidean metric, and let g' be structure-peeserving metric with respect to Under the same assumptions as Theorem 2.9, suppose that either of the following conditions holds:
(a) g is geodesically complete:
(b) the set of e-stationary points under g, K {pe M: [V,f(p)lnn S {e), is compact.
Then it requires at most T < 0 (g) iterations to achieve minigeer E V f(r)) E².
Proof. For the item (a), we omit its proof as is directly implied by setting h = 1. Recall that we write S b if there exists constant C> 0 such that < Cb. Now denote o w) h(p)gp(v, w). Theorem C.18 implies that
[VS())) $ $V T 15ET
min
It suffices to prove that if K is an e-stationary point under g' then it must be an e-stationary point under g (up constant scale). Note that
|| Ngff(p))leg
Vh(p) 14 vgf(p)) pp.g.
34

---

Published as a conference paper at ICLR 2026
As the result, we obtain
1
min 9p (V(()) VS(1)) - V
1515T MaX h(p)
d
min 9p ( (p). f S
max h(p) PEM
ISST
T
We restrict two sides on the compact set (given by the condition (b))
(p: (p))llps (
Because h MM- - IR is a continuous function, then must be bounded over this compact set. Let this upper bound be C. Then we obtain (with absorbing C' into 5)
d
[V & VT IStST
min
By setting V S 6,, we obtain the complexity T 2T.
D EXPERIMENTAL DETAILS
In this section, we aim include the omitted experimental details in Section 3.
Hardware and System Environment We conducted experiments the personal laptop, equipped with AMD Ryzen 9 7940HS Mobile Processor (8-core/16-three) and NVIDIA GeForce RTX 4070 Laptop GPU; however, GPUs not required in experiments. Our codes were tested using Python version 3.12.3. Additional dependencies are specified in the supplementary 'requirements.oxt file. All source codes attached.
D.1 SYNTHETIC EXPERIMENT: IMPACT OF SAMPLING BIAS
Construction of Quadratic Objective Functions We construct quadratic objective functions of the form quuadratie (x) x (B + Z) where B is a symmetric positive definite matrix that determines the landscape's curvature properties and & is the data point independently sampled from A 1) for each entry. The matrix B is generated by first creating random matrix M E Idddd with entries drawn from a standard normal distribution (0, 1), then forming
M M + dIa to ensure positive definiteness with regularization term dT, that controls the minimum eigenvalue.
Construction of Logistic Objective Functions For logistic objective functions, we construct the empirical risk minimization problem fingistic = Ei=1 log(1 + expp(xycc x)) + 5 T Ba, where B is generated as the same way the quadratic function and {(Ci, } represents the training dataset with feature vectors E Rd and binary labels y; E (-1, The feature matrix X [z]..... ........20 xx.]" T E Raxdd is generated from standard normal distribution (O, 1). A ground truth weight vector w² E Rd is generated from N0, 1) and then normalized unit length. The are generated by first computing logits x; w², then converting to probabilities p; 1/(1 + exp(-x; xx; w))), and finally sampling yi according to Bernoulli(p:) before converting to the { (-1,+1 +1} encoding. The regularization parameter À is chosen as À = 0.1.
binary labels y/s E (-1,
Construction of Riemannian Metric 9. We design Riemannian metric on the ambient Euclidean space by defining a symmetric positive definite matrix A with extreme conditioning properties. Specifically, the metric tensor is constructed by generating random orthonormal matrix Q via QR decomposition, prescribing eigenvalues that span geometrically from Amin = 1 to Amax = 10 Xmin- and forming A QAQ where A is the diagonal matrix of these eigenvalues. This construction yields a highly anisotropic Riemannian manifold with a condition number of A equal to 1 creating challenging geometric landscapes for optimization algorithms. The resulting metric induces Riemannian gradients the form A fx),, fundamentally altering the optimization dynamics compared to standard Euclidean methods.
35

---

Published as a conference paper at ICLR 2026
Hyper-Parameters Each method uses 16 random directions per iteration with perturbation stepsize ب = 10 4 for gradient estimation. The algorithms were run for 500,000 iterations with learning rates of 10-3 (quadratic) and 10 (logistic), and results were averaged over 16 independent runs to ensure statistical reliability. All curves are smoothed using a moving average with a window size of 5,000 iterations, and confidence bands represent 10th-90th percentiles across to visualize convergence variability.
D.2 SYNTHETIC EXPERIMENT: MSE VS. CURVATURE
Riemannian Metrie Construction
We work on the d-dimensional probability simplex
d+1
A = E R
Σ >==10<<<<>>
1=1
and endow its interior (identified with the first d coordinates) with structure- pressrrving Riemannian metric (Definition 2.5) conformally equivalent to the canonical Euclidean metric g
g(وو = =000(p)
where the conformal factor is
d+1 Mp) = 1+ Σ
(p) 1/2 3 log h(p),
P d+1
Varying the scalar ,3 sharpens or flattens the metric. We examine four choices 3 {0.5, 1.0, 1.5, 2.0}. At the fixed reference point Po E A (drawn once from the Dirichlet dis- tribution and held constant throughout the experiment) we measure the mean-squared error of a symmetric zeroth-order gradient estimator (Equation (3)) with using the first-order approximation of the exponential map as the retraction, where the perturbation stepsize y is set to 0.1. We note that under this approximation, the retraction degenerates to the naive Euclidean perturbation; we note that we are using a fixed point Po, it doesn't I trigger the out-of-domain issue of the incomplete Rie- mannian manifolds when μ is appropriately selected. The MSE is evaluated using the corresponding structure-prec strueture-preserving metric instead of the original Euclidean metric; the conformal scaling hp)) is applied consistently both when sampling directions (||v 1) and when converting the Euclidean gradients of the test functions (quadratic and Kullback- Leibler distance to the uniform distribution) into true Riemannian gradients.
Sectional Curvature Evaluation Instead of using 3 as the x-axis, we compute the sectional curvature Kgos (Pu) of each metric Po reflect the true relation between the intrinsic curvature and the estimation error. Let & = Φ3; then g(0) = g is a warped Euclidean metric whose curvature depends solely on 6. We draw an orthonormal pair (v, w) in the Euclidean tangent space Tp, A via Gram-Schmidt method, rescale them so that ||v|| ||w]] l 1, and evaluate
K K(oo v, w) =c 20(()) (II (I) ( 112 (Hess (pp)) u. v) - (Hess ((p)) w, w)
- (VVo(p)). w/2 (Vo(po), w)
where gradients and Hessians are taken with respect the ambient Euclidean coordinates. Because g(P) is isotropic up to the conformal factor, a single random 2-plane suffices; the resulting scalar is recorded as K (po) for that B. These four curvature values, monotonically decreasing as 3 grows, serve the horizontal tick labels in Figure lb.
Parameeerssst For each metric, we run 50,000 independent zeroth codeh-rrddrrae gradient trials, each trial drawing one random Riemannian unit direction and applying Equation (3) estimate the gradient with using the perturbation stepsize f 0.1 and using the exponential map as the retraction. The reference point Po E A sampled once and held fixed, so that changes in estimator accuracy stem solely from the chosen metric. Closed-form gradients are available for both test functions, Euclidean and KL distance to the uniform distribution. We record the mean-squared error h(po) II V ffll² for each trial. The resulting 50,000 errors per setting are summarized with log-scale box plots whose boxes span the inter-quartile range and whiskers cover the 10th-90th percentiles (outliers omitted).
36

---

Published as a conference paper at ICLR 2026
Table Hyper-paramettr settings for gradient-based mesh optimization experiments.
Fixed Hyperparameter
Symbol
Value
Fine-mesh size
Mfne
200 200
Number of nodes in the fine mesh
Nine
40,000
Coarse-mesh size
Meuurse
20 x 20
Number of nodes in the coarse mesh
N Necurse
400
Sampled nodes in each iteration All nodes positions
120
P
Total iterations
T
20,000
Tunable Hy perparameter
Random directions
Perturbation stepsize
10-
Learning rate
7
{300, 400, 500}
D.3 GRADIENT- BASED MESH OPTIMIZATION
In our work, we consider the black-box mesh optimization problem. In the well-known CFD-GCN model (Belbute-Peris support the auto-differemtiation in the SU2 PDE solver; however, in most of existing finite-volume numerical solvers, the positions of mesh nodes typically not differentiable. Therefore, need to apply the zeroth-order optimization approach.
al., 2020), additional efforts taken to allow the position of nodes to
Construction of Mesh Objective Function Let P {p:) C R22 be interior node positions of the given mesh with boundary nodes fixed. Given the coarse mesh induced by P defines PDE state üp (solved on P). Then interpolate it into the fine mesh Mine obtain the PDE state up. The objective mean-squared error (MSE) to a fixed fine-grid reference üref:
Fmeshh (P) = Nine up
2 2 Uref
where Nine denotes the number of nodes in the fine mesh. The randomness in this objective comes from the random sampling over the nodes; instead of taking all nodes be updated, each step will only sample a part nodes to be updated. In our experiments, we set the size of coarse mesh be 20 20 and the size of fine mesh be 200 200. Each time, we will randomly sample 30% x 20 x 20 120 nodes to update.
Construction of Mesh Parameterization Each interior node is updated in barycentric coordinates b E A m- with respect to its incident cell (with vertices {u, (u))y=)). i.e., p(b) = ΣØ= by This coordinate guarantees feasibility (by > 0, E3 = 1), which naturally results in probability simplex structure. Under the canonical inclusion embedding, this manifold is geodesically incomplete and hence feasible for our proposed approach.
Construction of Structure--ressevinn Metric We endow m- with the structure- prrettrrii-peeeiinnggn conformal metric ((8) as defined in Appendix D.2, and use the first-order approximation of the exponential map of g(() as the retraction (Definition B.3). We note that this approximation requires set the length of perturbation vectors to be sufficiently small to ensure the accuracy of the retraction; this requirement be satisfied by adopting the same technique as the soft projection trick used in Figure 5a. We always assume this requirement is satisfied throughout the training.
Hyper- Parame-rrssssrss Each iteration uses 4 random directions with perturbation stepsize μ 10-1 Optimization runs for T = 20,000 iterations with learning rate 77 E (300, 400, 500} (we report the best curve among these tyype--parametenn of 2,000 iterations. For all other other estimator- dependdet parammtar,m,,rs, we have included all of them in the configuration files along with source codes.
All curves smoothed with moving-average window
37

---

Published as a conference paper at ICLR 2026
D.4 ADDITIONAL DISCUSSIONS ON STABILITY
To further evaluate the robustness of the proposed method. we conducted repeated experiments using 5 independent random seeds. The results illustrated in Figure 7. We note that both our method and the Reversion method exhibit low variance and consistently show the narrowest error bands (smallest shaded areas). indicating that they are sufficiently stable. In contrast, the naive Unconstrained method suffers from the highest variance.
Figure 7: Training loss trajectories averaged over 5 independent random seeds. The solid lines represent the mean loss, while the shaded areas indicate the min-max deviation. Our structure- preserving method demonstrates superior stability and consistently lower loss.
E LIMITATIONS AND FUTURE WORK
Our work also presents several limitations which potentially point to future research directions. A potential limitation is the rejection sampling's ' reliance on the degeopscoonooiiionn of the metric matrix, which presents significant computational bottleneck when applied to high-dimensional problems. Developing more scalable algorithms (e.g. Randomized SVD or iterative solvers) to efficiently handle high-dimensional manifolds is therefore important future direction. Additionally, while our theory guarantees the existence of structure- prreterrrrres metric (Theorem 2.6), its practical construction currently relies on a case-by-case design and lacks a general construction. Moreover, our construction is based on the conformal transformation; a valuable direction for future work is to explore whether broader class of structure- pereturepresseiing transformations exists beyond the current scaling and to develop more general constructive methods that are better compatible for zeroth-orde optimization.
F THE USE OF LARGE LANGUAGE MODELS (LLMs)
In preparing this manuscript, we employed Large Language Models (LLMs) as purppse assistive tools the following ways:
Literature review support. We used the Deep Research functionality provided by existing AI platforms to help gather references and draft preliminary summaries of related work.
Language refinement. We used AI chatbots hosted on multiple platforms to generate the abstract and improve the elarity, style, and readability of the manuscript.
Proof verification. We used AI chatbots to check the logical consistency, correctness, and completeness of our formal proofs.
Codes Generation. We also applied the Al agent to generate a part of experimental codes.
All LLM-assisted outputs critically reviewed, verified, and, where necessary, revised by the authors. We take full responsibility for the content this manuscript. LLMs were not involved in generating research ideas, drawing scientific conclusions, contributing original insights.
38

---

Published as a conference paper at ICLR 2026
Table 2: A summary of notations in Riemannian manifolds.
Notations
Definition
Smooth Manifold (M)
A d-dimensional second-countable Hausdorff topological space where each point has a neighborhood Up diffeomorphic to Rd.
Deviation (v)
A linear mapping v (U,) - IR satisfying the product rule: o(fg) = v(f) g(p) v(g) f(p).
Tangent Space (T..)
The real vector space of all deviations at point M.
Cotangent Space (T
The dual space of the tangent space MrM,; the space of all linear maps + TTM + R.
Tangent Bundle (T.))
The disjoint union of all tangent spaces: TM == ((p,u) p E M, E M..
Immersion
A smooth map f M + R" whose differential df is injective every E M.
Embedding
An immersion that is also homeomorphism onto its image ((M).
Vector Field (X)
A smooth map (section) .M + TM such that X(p) E T.M
for all p M.
X(M)
The space of all vector fields on the manifold M.
Riemannian Metric (g)
smooth assignment of an inner product 9p T.M x TnM R to each tangent space T..M Also denoted (-, )p
Riemannian Manifold ((M,g)
A smooth manifold M equipped with Riemannian metric
n--uclideen Metric
A metric induced by a smooth embedding & .M + R²" via the pullback o, (v, u) (dolp (v), p (u)).
Levi-Civita Connection
The unique affine connection on x(M) that is torsion-free and metric-compatible.
Geodesic (?)
A smooth curve 1 + M whose velocity vector (t( satisfies the geodesic equation (t =
Exponential Map (exp,,)
A map exp, Tp.M M defined by exp,, == (1), where y the unique geodesic with (()) p and 7((0)
Retraction (Ret)
smooth map Ret TM + M satisfying Ret,(0) p and dRetplo = idy,M. It approximates the exponential map.
Gradient (Vf)
The vector field fp(p) := (df), p,),,) where is the musical isomor- phism M > TpM induced by the metric g.
Riemannian Stochastic Op- timization Problem
minppe f(p) = [f== {)], where f(: M + R a smooth function relying drawn from the data distribution E.
Symmetric zeroth-order es- timator
v, where μ is the perturbation stepsize and is uniformly sampled from the unit ball in T Ty.M.
(p)
2μ
SGD update rule
Pt+1 Ret,, (n f((: Er)), where 7 is the learning rate.
39


Reviews:
[{"review_id": "B3HmeZjSe4", "reviewer": "Reviewer_Jc8x", "summary": "This paper studies Riemannian zeroth-order optimization under geodesically incomplete metrics—a setting often overlooked in existing analyses that assume global geodesic completeness. The authors introduce a novel concept of structure-preserving metrics and develops an intrinsic two-point zeroth-order gradient estimator to resolve this problem. Convergence guarantees for Riemannian SGD with this estimator are provided, achieving rates comparable to those in geodesically complete settings.Empirical results on synthetic and mesh optimization tasks validate the theoretical findings.", "questions": "See \"weakness\" part.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper investigated a novel problem setting—— geodesically incomplete.\n2. This paper introduces some relative novel conceptions,especially structure-preserving metric and has elegant theoretical contribution.\n3. This paper has a great presentation with well writing.", "weaknesses": "It seems that Theorem 2.2 only provide the existence of structure-preserving metric. Can the authors give a guidance to construct such structure-preserving metrics?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wH9duPrlIN", "reviewer": "Reviewer_iA9C", "summary": "This paper considers zeroth-order optimization over Riemannian manifold that are embedded in Euclidean space. Existing notions of zeroth-order derivatives on Riemannian manifold critically rely on the exponential map. Nevertheless, this might be ill-defined when the inherited Euclidean metric is not complete. To address this issue, the authors show that there always exists a different non-Euclidean and structure-preserving metric. Using the new metric, they construct a sampling-based zeroth-order gradient estimator and establish theoretical convergence guarantee for this algorithm.", "questions": "No further questions.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "This paper has a clear structure is quite well-written. From a contribition perspcyive, the idea of using a structure-preserving metric is novel and mathematically elegant when the inherited Euclidean metric is not complete. The following theoretical analysis appears sound and follows standard arguments.", "weaknesses": "However, my major concern lies in the contribution aspect. Although the paper provides an interesting solution, the issue of metric incompleteness is quite specialized and rarely arises in practical Riemannian optimization scenarios. Researchers in this area are typically more interested in improved convergence rates or more efficient algorithms, which this work does not provide. Therefore, I tend to reject the paper due to the limited significance of its contribution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BwmgzGumWl", "reviewer": "Reviewer_sZeP", "summary": "This paper studies Riemannian zeroth-order optimization on geodesically incomplete manifolds. As the main contribution, It constructs structure-preserving metrics that are geodesically complete, conformally equivalent to the original one, while ensuring that every stationary point under the new metric remains stationary under the original metric.  A symmetric two-point zeroth-order estimator was developed with MSE analyzed. The paper then establishes convergence guarantees for SGD with this intrinsic estimator.  The proposed theory and methods are validated via synthetic experiments, and a practical mesh optimization task as well.", "questions": "1. The paper states that “the rescaling method even leads to divergence for the logistic loss objective.” Is this divergence inevitable (i.e., an inherent failure mode of the rescaling method), or could it be due to unlucky hyperparameter choices?\n\n2. In the geodesically complete case, can the more general framework presented in this paper imply the known optimal bound?\n\n3. Does there exist a class of \"structure-preserving\" transformations more general than conformal scaling that can simultaneously guarantee geodesic completeness and be more compatible with gradient / estimator perturbations; Is it feasible to use the \"structure-preserving\" idea as a preprocessing / transform to \"straighten\" a difficult-to-optimize manifold into an (computationally) equivalent problem that is more favorable for zeroth-order methods?\n\n4. Theorem 2.2 provides an existence result and a constructive proof for the conformal coefficient $h$. How can this construction be implemented in practice?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Connecting the geometric issue of **metric incompleteness** with Riemannian zeroth-order optimization is a worthwhile direction.  \n\n2. The intrinsic-view MSE analysis and the unified sampling perspective are promising.  \n\n3. The paper is generally clearly presented adn the theoretical development is solid.", "weaknesses": "## Major comments\n\n1. The Step 1 of Algorithm 1 requires eigen-decomposition of the metric matrix $A = Q\\Lambda Q^\\top$, and every sample requires generating Gaussians and computing the acceptance probability $\\sqrt{\\dfrac{v^\\top A^2 v}{\\lambda_{\\max}}}$. Eigen-decomposition is computationally expensive in high dimensions, and I guess this is why the experiments focus on relatively simple tasks. \n\n2. Assumption C.3 requires the objective function to have bounded third- and fourth-order derivatives which is a strong condition. As acknowledged by the authors, such a condition is less common in the literature.  The authors emphasized in the appendix that this assumption has also been used by Alimisis et al. (2021, Assumption 1), I fail to identify an equivalent assumption in that reference.\n\n## Minor comments\n1. It is suggested to add a short section in the appendix that lists and explains the frequently used notation. Although most symbols are defined at first use, having a centralized symbol table would greatly ease reading when symbols reappear across the paper.  \n\n2. In the Introduction section, it is desirable to provide some examples to help readers understand in which scenarios Riemannian zeroth-order optimization is preferable to conventional (Euclidean) zeroth-order optimization.\n \n3.  Concerning experiments, some more details about the settings of hyper-parameters should be provided for the sake of reproducibility.\n\n4. In Theorem 2.3, the term $\\dfrac{6}{d} + \\dfrac{8}{d}$ looks a bit weird — is there any typo here?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "D9f56LprFZ", "reviewer": "Reviewer_VPTy", "summary": "The authors consider the minimization of a real-valued smooth function on a manifold, using only function evaluations. Contributions are as follows:\n- To point an issue with existing zero-th order optimization methods, as follows. The exponential mapping is often defined on an open subset of the tangent space that contains the zero vector. There exist algorithms that rely on the exponentail mapping, but do not account for this fact. Unless the manifold is proved to be geodesically complete, that is the exponential mapping is defined for any tangent vector, such algorithms are ill-defined (theoretical guarantees are not applicable, possible breakdown in practice).\n- A theoretical guarantee that the metric of a geodesically incomplete manifold may modified such that it becomes geodesically complete, while preserving the stationary points of any function to be minimized (Th. 2.2).\n- a zero-th order gradient estimator that relies on two function evaluations, two geodesic computation, and a sampling routine on the tangent space. A mean-squared error analysis is provided (Th.  2.3).\n- an  analysis of the complexity of SGD with zero-th order gradient estimator for finding $\\epsilon$-stationary points (Th. 2.5),  and an analysis of the complexity of SGD and mapping from one metric to a geodesically complete metric, and back to the original metric (Coro 2.6).", "questions": "- Why does the scheme of (eq. 2, $f(Exp(\\mu v) - f(Exp(-\\mu v)))$) differs from the two references (l. 46 -  $f(Exp(\\mu v) - f(x))$)?\n- Why does working with a non-euclidean metric precludes the use of previous estimators, such as in Li et al, 2023b?\n- What are situations where the manifold and its metric are not geodesically complete?\n- l. 198-202: why do existing results in Riemannian optimization may not be applied with arbitrary metrics, different from the ones of the ambiant space?\n- What are examples with relevant applications of manifold that are not geodesically complete?\n- The same notation, $\\hat{\\nabla} f(p)$, is used for the gradient estimator when sampling with exponential mapping and retraction (eq. 2 and 5). To which estimator do Theorems 2.3 and 2.5 apply.\n- As far as I understand, while Th. 2.2 guarantees existence of a structure preserving metric, this proof does not help in designing the structure preserving metric in a practical situation, as for instance the example of fig 1 or the experiments. Am I correct in this? The applicability of the method thus fully relies on the practitionner. Can this be addressed by your theory?\n- Experiments, sec 3.1 and 3.3. Can you report experiments relative to time, and with variance indicators for 3.3?\n- Appendix C (proofs):\n  - Lemma C.5, can you review the proof with more details?\n  - What is the use of lemmas C.9, C.10, C.13?\n  - l. 1072: why introduce $g$\\# but not use it?\n  - l. 1106: what is $n$?\n  - l. 1107: the Ricci tensor is nowhere defined, nor discussed except in this lemma. This lemma is not used is the appendix as far as I can see. Why is that?", "limitations": "", "rating": 2, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- the paper identifies a reasonable issue in existing works,\n- the result Th. 2.3 account for manifolds curvature, a difficult topic, in a reasonably clear way,\n- one experiment (sec. 3.3) stems from a real-world application.", "weaknesses": "The following aspects motivate my current assessment of the paper.\n- The scheme of (eq. 2, $f(Exp(\\mu v) - f(Exp(-\\mu v)))$) differs from the two references (l. 46 -  $f(Exp(\\mu v) - f(x))$). Yet, this difference is not acknowledged, motivated, nor discussed.\n- I fail to see why does working with a non-euclidean metric precludes the use of previous estimators, such as in Li et al, 2023b.\n- The work focuses on situations where the manifold and its metric are not geodesically complete. However, the introduction does not provide any practical situations where this situation occurs. The only example in the paper is optimization on the simplex (fig. 1, and experiments), or a union thereof. I am not convinced this example alone is sufficient to generate interest from the iclr optimization community.\n- l. 198-202: I fail to see why existing results in Riemannian optimization may not be applied with arbitrary metrics, different from the ones of the ambiant space.\n- The only example with relevant application of a manifold that is not geodesically complete is the unit simplex (Fig 1), or a collection thereof (numerical experiments).\n- It is not clear to me that the proposed procedure to sample uniformly on an ellipse is novel (a contribution of this work) or classical. I am not an expert of distribution sampling, but I suspect the rejection sampling on an ellipse is classical.\n- The same notation, $\\hat{\\nabla} f(p)$, is used for the gradient estimator when sampling with exponential mapping and retraction (eq. 2 and 5). It is thus not clear to which estimator Theorems 2.3 and 2.5 apply.\n- As far as I understand, while Th. 2.2 guarantees existence of a structure preserving metric, this proof does not help in designing the structure preserving metric in a practical situation, as for instance the example of fig 1 or the experiments. The applicability of the method thus fully relies on the practitionner.\n- Experiments, sec 3.1: authors compare SGD with rescale sampling and rejection sampling in terms of SGD iterations. This does not account for the additional complexity of the rejection sampling, relative to the rescale sampling. As such, the plots does not inform on the performance of the complete method (SGD + gradient sampling scheme).\n- Experiments, sec 3.3: authors compare variants of SGD relative to the number of iterations, with four different strategies for the update. Besides, there is one trajectory for each method, which does not provide any information on the variance of the methods.  Again, I believe this does not inform on the performance of the overall methods, which I believe is ulitmately the metric of interest.\n- Appendix C (proofs) shows several serious issues:\n  - Lemma C.5: in view of $f*$ definition, there holds, by definition, $B = 0$.\n  - Lemma C.5, l. 1005-1008: the proof is not detailled enough to be checked.\n  - the lemmas C.9, C.10 are not referenced anywhere. Lemmas C.9, C.10, and C.13 are not used anywhere as far as I can see.\n  - l. 1072: introduces $g$\\# but does not uses it.\n  - l. 1106: $n$ is used but not defined there, nor in the referenced Lemma C.8.\n  - l. 1107: the Ricci tensor is nowhere defined, nor discussed except in this lemma. This lemma is not used is the appendix as far as I can see.\n\nMinor points\n- l. 42: what is nondifferentiable modules?\n- l. 90: writing issue \"while maintain\"\n- l. 127: syntax issue \"present main\"\n- l. 134: syntax issue \"we establishes\"\n- l.155-160 : ($\\epsilon$)stationary point not defined\n- l. 255: syntax issue in \"metric $g$ the choice\"\n- l. 303: I would appreciate that the theorems assumptions are stated in the main body rather than supplementary\n- l. 305: the SGD dynamics do not solve (1) (i.e. converge to the global minimizers of $f$ on M)\n- l. 372: $g_A$ not defined\n- l. 408: syntax issue \"aligns our\"\n- l. 939, last sentence of section: syntax issue.\n- l. 1015, lemma C.6: syntex \"there exists a function [...] is proper\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)