Landing with the Score: Riemannian Optimization through Denoising
Authors: Andrey Kharitenko
Zebang Shen
Riccardo De Santi
Niao He
Florian Dorfler
Venue: ICLR 2026 Poster (2026)
Keywords: Riemannian Optimization
Manifold Learning
Diffusion Models
Primary area: optimization

Abstract: Under the \emph{data manifold hypothesis}, high-dimensional data concentrate near a low-dimensional manifold. We study Riemannian optimization when this manifold is only given implicitly through the data distribution, and standard geometric operations are unavailable. 
This formulation captures a broad class of data-driven design problems that are central to modern generative AI. Our key idea is a \emph{link function} that ties the data distribution to the geometric quantities needed for optimization: its gradient and Hessian recover the projection onto the manifold and its tangent space in the small-noise regime. This construction is directly connected to the score function in diffusion models, allowing us to leverage well-studied parameterizations, efficient training procedures, and even
pretrained score networks from the diffusion model literature to perform optimization. On top of this foundation, we develop two {efficient} inference-time algorithms for optimization over data manifolds: \emph{Denoising Landing Flow} (DLF) and \emph{Denoising Riemannian Gradient Descent} (DRGD). We provide theoretical guarantees for approximate feasibility (manifold adherence) and optimality (small Riemannian gradient norm). We demonstrate the effectiveness of our approach on finite-horizon reference tracking tasks in data-driven control, illustrating their potential for practical generative and design applications.

Published as a conference paper at ICLR 2026
LANDING WITH THE SCORE:
RIEMANNIAN OPTIMIZATION THROUGH DENOISING
Andrey Kharitenko
Zebang Shen
Riccardo De Santi ETH Zurich, ETH AI Center x redesantiii meant@@ehhz.cc
ETH Zurich,
ETH Zurich, Institute for Machine Learning zebang- shen i nf. wtw.. .ch
Institute for Machine Learning
Niao He
Florian Dörler
ETH Zurich,
ETH Zurich,
Institute for Machine Learning niao he@ inff..... inf.ethe wettt.tcc ch
Automatic Control Laboratory
doer fl er@centrol ee wee.eth..cc
ABSTRACT
Under the data manifold hypothesis, high-dimensional data concentrate near a low-dimensional manifold. We study Riemannian optimization when this mani- fold is only given implicitly through the data distribution, and standard geometric operations are unavailable. This formulation captures a broad class of data-driven design problems that are central to modern generative AI. Our key idea is link function that ties the data distribution to the geometric quantities needed for opti- mization: its gradient and Hessian recover the projection onto the manifold and its tangent space in the small-noise regime. This construction is directly connected to the function in diffusion models, allowing us leverage well-studied pa- rameterizations, efficient training procedures, and even pretrained score networks from the diffusion model literature perform optimization. On top of this foun- dation, we develop two efficient inference-tine algorithms for optimization over data manifolds: Denoising Landing Flow (DLF) and Denoising Riemannian Gra- dient Descent (DRGD). We provide theoretical guarantees for approximate feasi- bility (manifold adherence) and optimality (small Riemannian gradient norm). We demonstrate the effectiveness of our approach on hoite-horizoo reference tracking tasks in data-driven control, illustrating their potential for practical generative and design applications.
INTRODUCTION
Riemannian optimization Boumal 2023) Absil et al. (2008); Hu al. (2020) considers minimizing an objective function f Rd over an explicitly known embedded submanifold M C Rdd,
min f(x).
#EM
Problem (T) is ubiquitous in fields of machine learning and control and encompasses problems such as independent component analysis Nishimori (1999), low-rank matrix completion Vandereycken (2013), training of orthogonally normalized neural networks Bansal et al.) (2018), the control of rigid bodies Duong et al. (2024). as well as sensor network localization Patwari & Hero (2004) and many others. Compared to general constrained optimization, Riemannian optimization promises the advantage of exploiting the natural geometry of the problem, producing feasible iterates and increased numerical robustness [ouuma] (2023). thus allowing for early termination and making more suitable for real-time implementation.
In contrast to the above classical setup, in this work we focus on the setting where the manifold M is given implicitly through finite set of samples from underlying population data distribution Adata supported on M. This perspective is especially relevant in view of the data manifold hypoth- esis Goaiza-Gannem al. 2024), which posits that many real-world data sets lie (approximately)

---

Published as a conference paper at ICLR 2026
on manifold with dimension much smaller than that of the ambient space Fefferman et al. (016)) Importantly, such data manifolds are not just cow-dimensional geometric structures - they also cap- ture rich semantic meaning. For instance, the image manifold corresponds to photo-realistic images Pope et al., the system behavior manifold to dynamically feasible input-output trajectories |Willems) & Polderman (1997). while the manifold of airfoils represents aerodynamically viable shapes (Zheng et alll]l (2025). The optimization problem (I) in this implicit setting thus encompasses a broad class of modern tasks, including airfoil- and ship hull design (Chen et al. 2025 Bagazinski & Ahmed, 2023), additive manufacturing Peng et al. (2023), reinforcement learning (Lee & Choi 2022). and Bayesian inverse problems (Chung et 2002b). In this data-driven regime, methods from classical smooth Riemannian optimization cannol be applied directly, since they rely on explicit manifold operations such as tangent-space projection, retraction, or exponential maps (Boumal et al.] 2019 Boumal 2023). Graph-based Riemannian optimization has recently been proposed in (Wang et al., 2025). but the setup is quite different, as to be claborated in Section 1.2 Furthermore, while there
has been extensive research on manifold learning Meilä & Zhang (2024 Lin & Zha (2008); Cay- ton et al. ) Belkin & Niyogi (2005) in the past, none of these works addressed (I) from an optimization point of view and focused instead on learning the manifold geometry, rather than rep- resentation models that are suitable to be included as a constraint in an optimization problem.
Instead, in this work, we propose novel data-driven approach recover the fundamental operations needed for optimization on manifolds: Starting first from the population data distribution (daaaa we smooth it with a Gaussian kernel to obtain
Po -N o21) Pdatas
and define the associated finkfunction
(x) 1 ||-||2 o² (o) x x EER
(3)
We show that, as the smoothing parameter a decreases, the gradient Vl, recovers the projection back to the manifold, while the Hessian V recovers the projection onto its tangent space. These results reveal that core ingredients of Riemannian optimization such as retraction and gradient computation can be implemented directly from the derivative information of la, which itself can be constructed from Mdata-
Given the above novel theoretical findings, a key practical challenge is how to access (even approximately) the gradient VL, and Hessian 2r when, in real-world applications, only samples from (ıdara are available. A crucial observation is that pa coincides with the marginal distribution of the Sariance-Exploding SDE (VE-SDE Note that the gradient of logPo. i.e. the score function, is central to diffusion models. It is parameterized by neural network and learned directly from samples of jdata. typically via denoising score matching (Vincent) 20IT). We can hence exploit this well-developed toolkit for practical implementation of our idea: Let s(r, be pretrained score network. We use v(x) = x + ( σ) to approximately represent Vfg. while the Hessian V'E can be recovered by computing its Jacobian. This approach offers two key advantages: (i) strong inductive biases can be incorporated through the neural network parameterization, and (ii) efficient, well-established training techniques from the diffusion model literature can be directly leveraged. Taken together, the theoretical link between Vl, and manifold geometry, and the practical machinery for learning the score, form the foundation of the paradigm shift: from classical smooth Riemannian optimization with explicit manifold knowledge to a data-driven framework where geometry is recovered from samples-thur enabling principled manifold optimization in generative and design-driven applications.
Building on these insights, we propose the first score-based framework for optimization over data manifolds. We derive two algorithms: denoising landing flow (DLF) and denoising Riemannian gradient descent (DRGD). Both rely on an approximate score network a ( from samples of Mdata using standard diffusion-mode training. DRGD uses the learned gradient and its Jacobian (for a small, fixed a) to mimic Riemannian gradient descent with retraction, while DLF performs gradient flow on a penalized objective, relaxing feasibility at intermediate iterates. We establish non-asymptotic convergence to approximate stationary points as - O. Importantly, our methods require only inference of the neural network and gradients with respect to its inputs- not
VlggPo Po learned
Analogous extensions to VP-SDE and DDPM can also be similarly derived.
2

---

Published as a conference paper at ICLR 2026
Car position
0050
75
-100
-175
trjjectory Vinmam data
-200
Reforence Urajectory
0.50
0.75
1.50
3.75 2,00
Figure 1: Optimized trajectory (orange) on the system trajectory manifold for the unicycle car model that is desired to track the reference trajectory (black, thin). The closest tracking trajectory in the set of available manifold samples is given in (green, dotted). See Section for more details.
with respect to the network parameters. Thus, pretrained score network is already available for given task, no additional training is required to enable Riemannian optimization on the corre- sponding data manifold. Viewed from this perspective, our approach can also be interpreted as an inference-time algorithm, aligning with growing trend in modern machine learning research.
1.1 OUR CONTRIBUTIONS
We next summarize our main contributions and outline the structure of the paper.
Link function and data-driven manifold operations. Building on the smoothed data distribution Po and the associated link function & introduced above, show in Section B that its gradient and Hessian recover, in the small-o regime, the projection onto the manifold and its tangent space. Further, with the help of score learning in diffusion model training, requirements of classical Riemannian optimization (which assumes explicit manifold knowledge) and the emerging need optimize over data manifolds that are only implicitly available.
bridge the gap between the
First score-based algorithms for optimization over data manifolds. We propose two algo- rithms denoising landing flow (DLF) in Section 4 and denoising Riemannian gradient descent (DRGD) in Section 5 the best of our knowledge, the first in the literature that exploit opera- tions enabled by pretrained score function. Our methods require only inexpensive inference- queries and pack-propagatton with respect to the input of the neural network, making them com- putationally efficient and readily applicable when pretrained scores are available.
Non-asymptotic guaranties. We establish non-asymptotic convergence guarantees for both al- gorithms, showing approximate feasibility (outputs close the data manifold) and approximate optimality (small Riemannian gradient norm) as + 0 (Theorems 3 and 5). A key technical ingredient is a uniform control of how Vl, and 2 approximate the ideal manifold operations at finite (Theorem
Finally, in Section we validate approach classical Riemannian optimization benchmarks and a data-driven optimal control problem for reference trajectory tracking. We demonstrate that methods can generate feasible points on the manifold with objective values significantly lower than those observed in the training data (Figure[]) illustrating how strong inductive biases of modern deep networks can be harnessed for constrained optimization.
1.2
RELATED WORK
In this section we briefly give an overview over the related work. More connections our approach previous work can be found in Appendix
3

---

Published as a conference paper at ICLR 2026
Smooth Riemannian Optimization. Riemannian optimization (RO) was originally developed un- der the assumption that the constraint manifold is explicitly known, either through closed -oomm de- scriptions such matrix manifolds via nonlinear equality constraints (Sato [021 [Bommal| (023).
A prototypical algorithm in this setting is Riemannian gradient descent, which updates by taking a gradient step along the tangent space and then retracting back to the manifold (Boumal al. 2019).
Such algorithms guarantee that the iterates remain on the manifold at every iteration and are thus re- ferred to as feasible methods. In contrast, infeasible methods, where the iterates are not constrained to stay on MM, have also been studied for example, augmented Lagrangian approaches (Xie & Wright 202T). However, these typically involve solving complicated optimization subproblems in each inner loop, making them computationally expensive in practice. More recently, a new class of landing-type algorithms has emerged: instead of enforcing feasibility at every step, they regularize the objective function with the distance to the manifold and then perform gradient flow or descent on this regularized objective. Such methods have demonstrated strong empirical performance, offering promising alternative to classical approaches (Ablin & Peyré 2022 Scheechmaan et al. 2023
Optimization over graph Recently, Wang et al. (2225) propose a graph-based strategy for man- ifold denoising, i.e. given observations sampled from the data manifold, for a new query point, find the closest point on the data manifold. The mixed-order optimization algorithm in this work consists of two types of updates: (1) By extending the concept of tangent space in classical smooth manifolds to the graph setting, first-order steps can be taken over the tangent bundle graph; (2) By enumerating the points connected by the "zero-order edge", zero-th order steps can be taken avoid local minima and achieve a strong result of convergence to global minima. We highlight that their setup is quite different from ours: (i) Their method is purely non-parametric and hence cannot ex- ploit the recent development of deep learning; (ii) The proposed mix-order algorithm is more like a discrete search method. As an example, their strong global convergence result is built on the exis- tence of zero-order edges, which "teleport" the iterate out of the local minima automatically. This is not the typical behavior of continuous optimization method, either deterministic or stochastie; (iii) Their method is purely designed for the denoising task while our focus is for minimizing general smooth objectives over data manifold. We highlight that due to the requirement of their analysis, for different objectives, the zero-order edges need to be reconstructed to avoid local minima.
Manifold Learning. High-dimensional data in modern machine learning often exhibit an intrin- sic lower-dimensional structure. Such structure is of central importance: it enables more efficient representation and compression of data, facilitates interpretability by revealing meaningful semantic organization, and provides foundation for designing algorithms that exploit geometry rather than ambient dimensionality. The task of uncovering this structure is commonly referred to as manifold estimation or manifold learning, with large body of work devoted to this goal.
Non-parametric approaches include Isomap (Tenenbaum al. 2000 Laplacian Eigenmaps
(Belkin & Niyogil (2003), and Locally Linear Embedding (LLE) (Roweis & Saul 2000), Dif- fusion Map (Coifiman & Lafon 2006). among many subsequent developments, e.g. (Zhou al.,
2020
Yao et al. 2023).
Parametric manifold learning techniques such generative adversarial networks (GAN) and various types of autoencoders (AE) are used to estimate data manifolds M by learning map ei Z - M that parametrizes M through some lower-dimensional "latent space" Z and which can be interpreted as a learned coordinate chart on the former Goodfellow et al. 2014); [Berahmand et al.) (2024). In the case of GAN this map corresponds to the generator, while for AE it corresponds to the decoder. For both GAN and AE, it has been noted that vb is not always a local chart in the sense Lee et al. (2022) and an additional line of work considered enforcing the latter property for the decoder through regularization or architectural constraints Sorrenson ett al. (2023); Kumar et al. (2020). Manifold learning flows (M-flows) Brehmer & Cranmer (2020) address this problem by parametrizing vo in such way that & is injective and sa b can be efficiently computed and have been used for manifold learning and density estimation. Having obtained a chart b. one can locally trivialize the problem equation T to min o 1h, variants of which have also been known under the name of latent space optimization Tripp et al. (2020); Notin et al. (002)) We provide a detailed comparison of this paragraph to our work in Appendix
differential-geomees
4

---

Published as a conference paper at ICLR 2026
Comparison with Pull-back Gradient Flow. When the map e is learned, it induces pull-back metric on the latent space Z. Smooth Riemannian optimization algorithms on M can then be equiv- alently formulated in Z; for instance. Riemannian gradient flow M corresponds to pull-back gradient flow on 2. While we are not aware of prior work that explicitly exploits this idea to address optimization over data manifolds, we provide comparison with this approach in Appendix an2 highlight several nontrivial challenges that should not be overlooked.
Diffusion models. Diffusion models have achieved remarkable success in generative modeling, where the goal to generate new samples consistent with an underlying data distribution (datt- A central ingredient of these methods is the learning of the score function - the gradient of the log- density of the diffused distribution Pa Tang & Zhao (2255) [ong et al. (2020). In practice, the score is parameterized by neural network whose architectural design has been extensively studied, and its training is carried out using eell-established techniques such as denoising score matching (Song et al. 2020). This combination of principled theory and mature practice has made diffusion models one of the most effective tools for data-driven generative modeling.
When the data distribution resides on manifold, a recent observation in the diffusion model com- munity is that the score is asymptotically orthogonal to the manifold surface Stanczuk et al. (2024). This observation has been exploited for the estimation of the manifold dimension. See (Kamkari ett al. 2024) for the same task. Furthermore, Ventura et al. (2024) has shown that in the case of linear manifolds (i.e. affine subspaces), the Jacobian of the score scaled by the diffusion tempera-
ture asymptotically approximates the projection of the manifold onto the normal space and used it to study the geometric phases of diffusion models.
Recent work on the statistical complexity of diffusion models under the manifold hypothesis pro- vides indirect evidence that these models capture geometric information about the underlying data manifold (Oko et al., 2023] Tang & Yang. 2024). In particular, the sample complexity required to learn the data distribution Mdata depends only on the intrinsic dimension of the manifold, rather than on the ambient dimension.
Two Formulations: Optimization and Posterior Sampling, To avoid possible confusion, we stress the distinction between our optimization formulation and the posterior sampling literature. In short, our optimization formulation in eq. ( enforces the manifold constraint directly. This en- sures (approximate) feasibility at the final step and guarantees that the optimization process remains semantically meaningful, which is not guaranteed by the sampling formulation, as discussed below.
Across the literatures Classifier (-ereefee(-FFree Guidance in Diffusion Models (Dhariwal Nichol 2021 Ho & Salimans| 2022) and the Plug-and-Play Framework in Bayesian Inverse Problems ( (Venkatakrish- nnn et al. Laumont et al. 2022 Pesme et al. 2025 Graikos et al. 2022 Chung et al.) 2022a).
unifying perspective is the task of sampling from posterior distribution of the form
Ppost Pere exp
where r denotes a cost function to be minimized, corresponding to our objective f, and a > 0 is a lemperature-like parameter (Domingo-Enrich) et al. 2024). In Classifier Guidance, encodes a classifier signal (c.g., specified by a prompt) and in Bayesian inverse problems, r is the negative log-likelihood of observations. Meanwhile, Ppre serves prior distribution, typically derived from a large-scale pretrained generative model, which is expected to capture the semantic structure of the data manifold. Sampling from Ppost thus aims to balance semantic plausibility with low cost.
If Ppre were exactly supported on the data manifold M and + 0, this formulation would reduce to the constrained optimization problem eq. (T) In practice, however, the situation is very different: the distribution induced by pretrained diffusion model is not concentrated on a low-dimensional manifold but rather has support of non-zero Lebesgue measure in the ambient space-indeed, in many cases, essentially the full space (due to the noisy generation process of Ppre). As result, when a is set too small, the posterior Ppost becomes dominated by the exponential tilt exp(-r/a). pushing samples into regions far from the true data manifold M and thereby losing semantie meaning. Consequently, these sampling-bars based frameworks must carefully tune a to trade off semantic fidelity (staying close M) against optimization quality (achieving low r).
5

---

Published as a conference paper at ICLR 2026
2 PRELIMINARIES
Here we introduce some preliminary notation and concepts to state our results. We refer the reader the Appendices and D for more details.
2.1 MANIFOLDS AND DISTANCE FUNCTIONS
Let M S _ be a k-dimensional embedded compact Co-submanifold without boundary. For any point E .M we denote by T M S Rd and Np M _ R the tangent and normal spaces of M at p. respectively, and their orthogonal projections by PT, M and PN, M. For C²-submanifold, then there exists radius TM such that every point in ((TM), where
T(r) = {x Rd
dist(x,, M) < <}}, dist(z, inf || pll: PEM
has unique projection (2x) E .M. Sets of the form T(T) with E (0, TM) are called tubular neighborhoods of M. The squared distance function is defined by d(x) = 2 dist((IM)) and and are both differentiable on T(TM) with 2 |11 -=))]1 = d(xr) and x (x) E Nx(a) M and d' (x) and r' (1: given explicitly in Appendix B. Finally, for function f M - R
we denote the Riemannian gradient by grad fp). which in the case of f Rd - R is given by
grad M (() PT, M (p()..
2.2 THE STEIN SCORE FUNCTION AND SCORE-BASEE DIFFUSION MODELS
For Borel probability measure denote its Gaussian blurring by Pa = NO 2 μ for 0 with Pu = μ.. The score function logp, of Pa and its Jacobian allow for the following interpreta- tion (see Jaffer & Gupta (1972), also Appendix Ch:
x + 0 Pog (x) = (.) = Euz,o + 0² log
σ² Cov(Vz,)), (4)
where l, is the link function 3 and Vs.o is the posterior distribution observing under the noise model Po and prior μ.. The representation for Evg.o has also been known under the name of Tweedie's formula Robbins (1992); Efron (201I). The score function has gained recent attention due its use in score-based diffusion models in the field of generative modelling. Specifically in the so-called variance exploding (VE) diffusion scheme one seeks to learn V log Po for differ- ent noise scales a via neural network s(,,a) rr) by minimizing the conditional score matching loss (ssM(s(, r)), which attains its unique minimum in s(., V Po. Then s(, is used for sampling from y by following a particular reverse-time SDE or ODE flow the noise scale o (see Appendix D for more details).
SCORE AS RETRACTION AND TANGENT SPACE PROJECTION
In this section give another interpretation to the score function V log p,(x) lo 42 log Por (x). Namely, it has been already observed in Stanczuk et al. (2024) that the score function for a + asymptotically orthogonal to the tangent space of the data manifold. Our first contri- bution is showing that when the support of the distribution is manifold M and μ is absolutely continuous w.r.t. its volume measure, then both quantities (4) approximate the projection operator (x) and its Jacobian nt (x) uniformly on tubular neighborhoods of M. Formally we establish the following
and its Jacobian
Theorem 1 (Main). Let M S gd be a compact, embedded C CB-submanifold and μ P(RR)) dp Borel probability measure with supp M and μ & Volm such that a E (3M). Then for any T E (0, TM) there exist some constants K = K(T,M, > and T = M(T,M, > depending t, M and p such that
1 σ
Ev x(x)] {o]log(o)) and
Cov(Va.o) x(x)) S {Ko(lo
(5)
for all a E (0, and x E T(T).

---

Published as a conference paper at ICLR 2026
The proof is deferred to JppendixxEannd is based careful non- con-ayympooiic estimate of the Laplace integral method. As a consequence and together with fact that T r²) coincides with Pr, M for E M (see Appendix(B) we obtain the following result.
Corollary 2. Let M and p be as Theorem 1 and suppose that E M. Then lim 1 + o log (x() = Pr, M R-0
In view of Theorem I let us abbreviate 4 into
(a) - log Po
To =a( = r o Vlgo
P=(=)=I+o log Po (x)
(6)
with the limiting cases dp = d, mo = π and Po(z)= (x) We stress that the expressions in 6 are
defined for all x E Rd, whereas d. श and n' in are only sensible in a tubular neighborhood T of .M.
Thus, in theory, a well-trained diffusion model score s(:oo) and its Jacobian s' s'(., allow us to
approximate the closest-point projection (xx) and the tangent space projection PT. . as a -
arbitrarily well via the operator v(z) = x + o²s s o) and its Jacobian, respectively.
4 DENOISING RIEMANNIAN GRADIENT FLOW WITH LANDING
In this section we show how for some smooth f E C (R"). Assuming that we have access to sufficiently accurate estimate of the score in form of a vector function
the score function from Section 3 for Riemannian optimization of
C. (IR: IE*) such that vzx) x + = o) and [[v - To (z)]l S e and [[v((x) - (x))l S E for x E T(T)
(7)
for some E TM), we propose for 1 2 0 the demoising landing flow (DLF)
3 v'
+ n[v(x) x).
(8)
Po (x) (i.e. = 0 in flow (8 the gradient flow + (T (2) - x) with Fx [T + (z) (9) and the dynamics in (9) consists of two parts: An approximate projection Po (z) (To (x)) of the gradient ( (z)) and an approximate landing term (7 (2)) - x) corresponding to the penalty function dy (z). In the further case of 0 and ==)) MM the flow (9) reduces to the ordinary Riemannian gradient flow, which has been extensively studied Helmke & Moore (001)):(mbrrsii
In the exact case v(x) = Th (z) and v(xx) (
== - = -
al. 2005). Interestingly, when 0, but only <x)) E TTr), then (9) reduces
= - -H=² grad (T(())) + [)((()))- -)),
(10)
with 11,1g linear operator on Tzx M given in Appendix In particular the two terms in (10
belong to T((x)) M and N N=(a) M., respectively, and are orthogonal to each other, implying that the distance between x and M is non-increasing, which allows for perfect landing on the manifold
via similar arguments as in Ablin & Peyré (202);; Schechtman et al. (2023), see Theorem 22 in
Appendix F3 For a > 0 or e O, the two summands in (9) and (8) in general not perpendicular to each other and the landing is not exact. However, using Theorem we show the following.
Theorem 3. Consider the flow (8) and E (0, TM) Ser
|| fll and L Lip( ff). Sup-
pose that for some e 0nd
(T, p)) with e+ K(, M, (oo(lo)(o)) < min (T, 1+0/m
the function v satisfies (7). Then for any x0) E T(T) the solution sr() to o exists for all 2 0 and is contained in T (T). Moreover, every accumulation point x. of this flow satisfies
dist distMM (z,) S To :0 1 ( (c
+1
(e+ K(T, ((.K(.,))oolllg))))
for the projection p. T(x.) it holds that
llgrad M f(p))ll < 2( [2 + C 2m) + (1 + C/n)//MM 1-T/TM
c) (c+ (11)
Thus, Theorem 3 shows that can still use the flow (8) for fixed 0 to converge approxi- mate critical points of the objective f. at which the approximation error and norm of the Riemannian gradient are both Ō(a) plus the score error
Remark 4. We can evaluate the right hand side of the flow { in single forward-backward pass of the network v. Namely, given an input x, we compute and store p v(z) by a forward pass of v, while keeping the computational graph of v(x). Then evaluate g f and build the computational graph of y = (v(x), g), while detaching g. Finally we backpropagate on x in y to obtain v (x)g v (())V(((())
7

---

Published as a conference paper at ICLR 2026
5 DENOISING RIEMANNIAN GRADIENT DESCENT
In a practical implementation one has to consider discretized version of the flow (8 natural alternative is to study the following approximate version of the Riemannian gradient descent Absil ct al. 2008);; Boumaa (2023)
Fk++ = v( - 7vv' (2t (xx))
(12)
which we term the denoising Riemannian gradient descent (DRGD), Here acts approximate retraction and v' as an approximate projection onto the tangent space, We obtain the following convergence result for this algorithm.
Theorem 5. Let T and
and C [V and Lip(V)). D -IflTT))
Lu Lo= = SC (2(
22-- + TM)
TM -1) + +L.
Suppose that for some (0, (T, M, p)) and 0 with = K K(T, M, ())o((og(()) S T/2 the function v satisfies 7 and that the step-size is constrained by E [Tmin Ymax] with
C(4+T) Then for any (TT/2) the iterates of 12 satisfy {= S TT()) T(T/22 and for the projection pt (xxx) the following average-on of-rradient-nom condition holds:
Tuiin
mmax
min
Lo
N Igrad f(pa)1 4D/N + (8C26/TJJ + 2(20 + Toxe (C +))))
Ymin La %max)
In particular there exists at least accumulation point H. T(c) of (DA) such that its projection = xxx)) satisfies
Ilgrady f(p.)11 (CC2)/TM+2200 + Love L)))e'
Vmin (1 Lo 7mmax)
and O. recover, up constants, the classical result iterates of Remark 6. we descent with known manifold M and non-convex objective f with Lipschitz the Riemannian gradient
gradient (Boumal) 2023 Corollary 4.9) (see Appendix A.3
Remark Note that both Theorem and Theorem require the rather strong T e--appoxximation assumption (7) on and its Jacobian. The analysis under weaker LL2bound is out of score for this paper and left for future work.
NUMERICAL EXPERIMENTS AND APPLICATIONS
In this section we provide some numerical results for our proposed algorithms, namely the denoising landing flow 8 (more precisely, the discretized version (5I) and the denoising Riemannian gradient
descent 12
6.1 OPTIMIZATION ON ORTHOGONAL GROUP (nn)
In this section we evaluate flow on a synthetic example. In order to compare the error of our method classical Riemannian optimization techniques, consider distributions supported manifolds that are the focus of study in Absil et al. (2008); [[]]/Buu]] (2023). Specifically we consider
the orthogonal group manifold M O C 1xx" with Vol,M being the uniform volume
measure and Brockett's 's cost function Helmke & Moore (2012) defined by
min f(X) tr(AXQx XEO(n)
where A, Q E Saxn are given. We consider the cases {10,20) and assume that we are given set Dirain C M of Naata = 20000 data points from μ and train the score function s with denoising score matching (see Appendix D for diffusion models and Appendix H.I.I for implementation de- tails). In Figure 2 (left) we compare the evolution of the objective value our approximation to the exact landing flow (1Op for different noise levels 0 and dimension n. We observe that can obtain objective values with cost lower than the best possible point in the training set and that the accuracy improves - O.
(8)
8

---

Published as a conference paper at ICLR 2026
Dust kinding
Figure 2: Objective value vs. flow time & for the orthogonal manifold for n = 10 (left, different > 0) and for n 20 (right, 0.05). Here "exact landing" refers equation g with exact operations = π and v' = T'.
6.2 REFERENCE TRACKING VIA DATA-DRE DRIVEN CONTROL
Problem definition: In this example we consider applying our method to the control of discrete- time dynamical systems on finite horizon. Specifically we assume that we are given a discrete-time state-space system
Tk+1 = f ux): yk: = g( ux),
N₁ O...., 1,,
(13)
on a finite time horizon NA with state 2i R== input un RR=", output yk RR and fixed initial state zn = 0 E R The task is to find inputs u (wo,.... WN,-1) such that the corresponding
outputs y (yo. (....... yn.) closely track a prespecified reference trajectory (ro...., ,N,)) by
solving the optimal control problem
min (u.yleMMo
f(u,y)
(14)
with tracking objective
N
f(u, y)
Σ
ug Rug + (yk Tk) T Tx) + (yN, T TNN))
(15)
k==0
for positive-definite weight matrices R gnuxna and Q SuyXns and feasible input-output set
{ (u, E (R x (IR")) Nut1 exists (xo,..... NN..) R
N
Mio =
with (u, x, y) satisfying 13
Under smoothness assumptions on the dynamics f and g. the set Mio. being a graph of a smooth map, is (non-compact) embedded smooth submanifold of (R"")) (IR")) No The problem
14 is ubiquitous receding horizon control applications such as model predictive control (MPC)
and used for e.g. autonomous driving Vu al. (2021). motion planning Cohen et al. (2020), op-
timizing HVAC system energy efficiency (Serale et al. (2018) and inventory control Kostič (2009). In many of these applications the dynamics (13) governing the system are not known explicitly. In-
stead, in data-driven control Dörfler (2023a by: Markovsky et al. (2023) one assumes that (13) is
given implicitly by a finite number of measured input-output trajectories
Derain = {(u,, y:) i 1..., Naata) C M1o
where the input u is persistently exciting Willems et al. (2005), e.g. given by (white) noise Ljung (1999). In particular the so-called system behavior manifold Mio Willems & Polderman (197)) is given by samples from a distribution j on its in- and outputs and fits precisely into framework of data-driven Riemannian optimization I A similar setup with observable state = x has been considered in the domain of reinforcement learning JJanerr et al. (2022). We test proposed denoising Riemannian gradient descent two classical systems from the

---

Published as a conference paper at ICLR 2026
First pendulum angle
Car position
6:60
1.20
-0.25
-0.00
PondJun stgle
0.50
-4.7
0.25
-1.00
-1.2-
-0.25
Optinieed Ergecdory Tue tigedaary
Opcinseed trpecdrry Tue Bia
-5.00
Minmul toin data
Knnmum train data Reference trjjectory
Reforence Injectory
Figure 3: Denoising Riemannian gradient descent: Angle of the first pendulum (left) and unicycle position (right) with the optimized output trajectory y" (blue, dashed), the true system trajectory ytuue (orange), the initial trajectory yo (green, dotted) and the reference trajectory r (red)
control domain: The discretized double pendulum system and the unicycle car model La Valle 2006) (see Appendix H.2.T detailed information the systems and the particular choice of r, R and Q in (15), each on a horizon of N₁. = 100. To apply our proposed methods, we train diffusion model (see Appendix H.2.2 for implementation details)
the measured trajectories Dtrain and solve [14p via the denoising Riemannian gradient descent to obtain a solution (u² y²). As initial values we take the trajectories from the training set that minimize the objective cost f. ie. (No, yo) = arg min (u.)) f(u, y). Note that, as seen in Section in general the final iterate will not exactly lie on the input-output manifold, i.e. (u",y") $ .Mio To account for this deviation we back-test our generated input trajectory u* by implementing it on the true underlying system 13 to obtain the real output y tre
depict the evolution of the ob- jective value w.r.t. iteration count and in Figure 3 we depict the final optimizing trajectories y" and gruu We can observe that the error yyy y tru is small, which shows that (u", yy" is close the true system behavior Mio- Moreover, we can see that the trajectory y " tracks the corresponding reference much better than the train set minimum yo, which shows generalization capability of our diffusiin model. We from Figure (left) that the current objective can depart from the true objective significantly. This is due to the iterates deviating from Mio- In this example the algorithm (DRGD) recovers and we have found it to be robust w.r.t. moderate deviations from the manifold. Note that we have set iteration budget of Niter = 3000 and Niter 2500, respectively, while the objective is still decreasing. Accelerating the convergence of DRGD is left for future work.
Results and discussion: In Figure (in Appendix H.2.3
CONCLUSION AND FUTURE WORK
In this paper, show that the denoising score and its Jacobian allow to perform manifold opera- tions such as the closest-point and tangent space projection without the explicit knowledge of the manifold. We then propose landing flow for the corresponding manifold-constraim nian optimization problem and show that its limit points approximate critical points of the original problem. Moreover, we investigate the approximate version of the Riemannian gradient descent and provide bound on the average-gradee graegee-rrdieennoooom of its iterates, which converges to zero as the manifold operations become more exact. We apply both algorithms on known manifolds and to finite-horizon reference tracking in the domain of data-driven control. Future work will consist of deriving error bounds for this flow when the denoising score is trained with a non-zero errerror as well as the study of more sophisticated classical Riemannian optimization algorithms such as New- ton and trust region methods when using the approximate manifold operations with the trained score accelerate convergence.
Rieman-
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS AND DISCLOSURE OF FUNDING
We thank the reviewers for their constructive suggestions. The work is supported by ETH research grant, Swiss National Science Foundation (SNSF) Project Funding No. 200021-207343, and SNSF Starting Grant. Andrey Kharitenko is supported as a part of NCCR Automation, a National Cen- tre of Competence in Research, funded by the Swiss National Science Foundation (grant number 5INF40 225155). Riccardo De Santi is supported by the ETH AI Center through an ETH AI Center doctoral fellowship.
REFERENCES
Theagenis Abatzoglou. The minimum norm projection on C2-manifolds in Rn. Trunsuctions of the American Mathematical Society, 243:115-122, 1978.
Pierre Ablin and Gabriel Peyré. Fast and accurate optimization on the orthogonal manifold without retraction. In International Conference on Artificial Intelligence and Statistics, pp. 5636-5657. PMLR, 2022.
P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithmss on matrix manifolds. Princeton University Press, 2008.
Guillaume Alain and Yoshua Bengio. What regularized auto-encoders learn from the data mata-gennrating distribution. The Journal of Machine Learning Research, 15(1):3563-3593, 2014.
Felipe Alvarez, Jerome Bolte, and Olivier Brahic. Hessian riemannian gradient flows in convex programming. SIAM journal control and optimization, 43(2):477-501, 2004.
Shun-Ichi Amari. Natural gradient works efficiently in learning. Neural computation, 10(2):251- 276, 1998.
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer, 2005.
Noah Bagazinski and Faez Ahmed. Shipgen: diffusion model for parametric ship hull generation with multiple objectives and constraints. Journal of Marine Science and Engineering, 11(22): 2215, 2023.
Nitin Bansal, Xiaohan Chen, and Zhangyang Wang. Can we gain more from orthogonality regu- larizations in training deep networks? Advances in Neural Information Processing Systems, 31, 2018.
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373-1398) 2003.
Mikhail Belkin and Partha Niyogi. Towards theoretical foundation for Laplacian- based manifold methods. In International conference on computational learning theory, pp. 486-500. Springer, 2005.
Kamal Berahmand, Fatemeh Daneshfar, Elaheh Sadat Salehi, Yuefeng Li, and Yue Xu. Autoen- coders and their applications in machine learning: survey. Artificial intelligence review, 57(2): 28, 2024.
Nicolas Boumal, An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.
Nicolas Boumal, Pierre-Antoine Absil, and Coralia Cartis. Global rates of convergence for noncon- vex optimization on manifolds. IMA Journal of Numerical Analysis, 39(1): 1011)11-3,, 2019.
Johann Brehmer and Kyle Cranmer. Flows for simultancous manifold learning and density estima- tion. Advances in neural information processing systems, 33:442-453, 2020.
Paul Breiding and Nick Vannieuwenhoven. The condition number of Riemannian approximation problems. SIAM Journal on Optimization, 3191):1049-1077, 2021.

---

Published as a conference paper at ICLR 2026
Lawrence Cayton al. Algorithms for manifold learning. Univ. of California at San Diego Tech. Rep, 221-17))1,,005, 2005.
Long Chen, Emre Oezkaya, Jan Rottmayer, Nicolas R Gauger, Zebang Shen, and Yinyu Ye. Adjoint- based aerodynamic shape optimization with manifold constraint learned by diffusion models. arXiv preprini arXiv:2507.25445 2025.
Hyungjin Chung, Jeongsol Kim, Michael T' Mccann, Marc Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. arXiv preprint (tti::2209.14688. 2022a.
Hyungjin Chung. Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye. Improving diffusion models for inverse problems using manifold constraints. Advances in Neural Information Processing Systems, 35:25683-25696, 2022b.
Mitchell Cohen, Khairi Abdulrahim, and James Richard Forbes. Finite-horizon LQR control of quadrotors on SE2(3). IEEE Robotics and Automation Letters, 5(4):5748-57) 57544 2020.
Ronald R Coifman and Stéphane Lafon. Diffusion maps. Applied and computational harmonic analysis, 21(1):5-30, 2006.
Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. Advances in neural information processing systems, 34:8780-8777 00-894, 2021.
Vincent Divol. Measure estimatiin manifolds: an optimal transport approach. Probubility Theory and Related Fields, 183(1):581-647,22 2022.
Carles Domingo-Enrich, Michal Drozdzal, Brian Karrer, and Ricky TQ Chen. Adjoint matching: Fine-tuning flow and diffusion generative models with memoryless stochastic optimal control. arXiv preprint htth::400900886.. 2024.
Florian Dörfler. Data-driven control: Part one of two: A special issue sampling from a vast and dynamic landscape. IEEE Control Systems Magazine, 43(5):24 27, 2023a.
Florian Dörfler. Data-driven control: Part two of two: Hot take: Why not go with models? IEEE Control Systems Magazine, 43(6) 2023b.
Thai Duong, Abdullah Altawaitan, Jason Stanley, and Nikolay Atanasov. Port-Hamiltonian neural ODE networks Lie groups for robot dynamics learning and control. IEEE Transactions on Robotics, 2024.
Bradley Efron. tion, 106(496):16021166) 14. 2011.
Tweedie's formula and selection bias. Journal of the American Statistical Associa-
Bálint Farkas and Sven-Ake Wegner. Variations Barbalat's lemma. The American Mathematical Monthly, (2):8)8 830, 2016.
Charles Fefferman, Sanjoy Mitter, and Hariharan Narayanan. Testing the manifold hypothesis. Journal of the American Mathematical Society, 29(4)9933-1049 2016.
Ian J Goodfellow, Jean Pouget-Abadee Mehdi Mirza, Bing Xu, David Warde-Farley. Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in neural information processing systems, 27, 2014.
Alexandros Graikos, Nikolay Malkin, Nebojsa Jojic, and Dimitris Samaras. Diffusion models as plug-and-play priors. Advances in Neural Information Processing Systems, 35:14715-14728. 2022.
Enkelejd Hashorva, Dmitry Korshunov, and Vladimir Piterbarg. On Laplace asymptotic method with application random chaos. 2015.
Uwe Helmke and John B Moore. Optimization and dynamical systems. Springer Science & Business Media, 2012.
Jonathan Ho and Tim Salimans. arXiv:2207. /2,988 2022.
Classifier-free
diffusion guidance.
arXiv preprint
12

---

Published as a conference paper at ICLR 2026
Jiang Hu, Xin Liu, Zai- Wen Wen, and Ya-Xiang Yuan. A brief introduction manifold optimiza- tion. Journal of the Operations Research Society of China, 8(2):199-2 248, 2020.
Chii-Ruey Hwang. Laplace's method revisited: weak convergence of probability measures. The Annals of Probability, pp- 1177-1182, 1980.
Tadeusz Inglot and Piotr Majerski. Simple upper and lower bounds for the multivariate Laplace approximation. Journal of Approximation Theory, 186:1-11, 2014.
Amin G Jaffer and Someshwar C Gupta. On relations between detection and estimation of discrete time processes. Information and Control, 20(1))46-54, 1972.
Michael Janner, Yilun Du, Joshua B Tenenbaum, and Sergey Levine. Planning with diffusion for ficxible behavior synthesis. arXiv preprint S2Kiv::205009991,, 2022.
Hamid Kamkari, Brendan Ross, Rasa Hosseinzadeh, Jesse Cresswell, and Gabriel Loaiza-Ganen A geometric view of data complexity: Efficient local intrinsic dimension estimation with diffusion models. Advances in Neural Information Processing Systems, 37:38307-38354, 2024.
Hassan K Khalil and Jessy W Grizzle. Nonlinear systems, volume 3. Prentice hall Upper Saddle River, NJ, 2002.
Konstantin Kostić. Inventory control as a discrete system control for the fixed-order quantity system. Applied Mathematical Modelling, 33(11):4201-4214 2009.
Abhishek Kumar, Ben Poole, and Kevin Murphy. Regularized autoencoders via relaxed injective probability flow. In International conference artificial intelligence and statistics, pp. 4292- 4301. PMLR, 2020.
Tomasz M Lapinski. Multivariate Laplace approximation with estimated error and application to limit theorems. Journal of Approximation Theory, 248: 105305, 2019.
Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. Bayesian imaging using plug & play priors: when langevin meets tweedie. SIAM Journal on Imaging Sciences, 15(2):701-737, 2022.
Steven M Vallee. Planning algorithms. Cambridge University Press, 2006.
Kyowoon Lee and Jaesik Choi. Local manifold approximation and projection for manifold aware diffusion planning. arXiv preprint arXiv:2506. 008677 2025.
Yonghyeon Lee, Sangwoong Yoon, Minjun Son, and Frank C Park. Regularized autoencoders for isometric representation learning. In International Conference on Learning Representations, 2022.
Tong Lin and Hongbin Zha. Riemannian manifold learning. IEEE transactions on pattern analysis and machine intelligence, 30(5):796-809, 2008
Lennart Ljung. System identification (2nd ed.): theory for the user. Prentice Hall PTR, USA, 1999. ISBN 0136566952.
Gabriel Daaiz--Ganem, Brendan Leigh Ross, Rasa Hosseinzadeh, Anthony Caterini, and Jesse Cresswell. Deep generative models through the lens of the manifold hypothesis: A survey and new connections. Transactions on Machine Learning Research, 2024.
Piotr Majerski. Simple error bounds for the multivariate Laplace approximation under weak local assumptions. arXiv preprint 151iv:/5/1..000 2015.
Ivan Markovsky, Linbin Huang, and Florian Dörfler. Data-driven control based on the behavioral approach: From theory to applications in power systems. IEEE Control Systems Magazine, 43 (5):28-68) 68, 2023.
Marina Meilä and Hanyu Zhang. Manifold learning: What, how, and why. Annual Review of Statistics and Its Application, 11(1):393-417, 2024.
13

---

Published as a conference paper at ICLR 2026
Yasunori Nishimori. Learning algorithm for independent component analysis by geodesic flows orthogonal group. In IJCNN' 99. International Joint Conference on Neural Networks. Proceedings (Cat. No. 99CH36339), volume 2. pp. 333-938. IEEE, 1999.
Pascal Notin, José Miguel Hernández- Lobato,, and Yarin Gal. Improving black-box optimization in vae latent space using decoder uncertainty. Advances in Neural Information Processing Systems, 34:802-814, 2021.
Kazusato Oko, Shunta Akiyama, and Taiji Suzuki. Diffusion models minimax optimal distri- bution estimators. In International Conference on Machine Learning, pp. 26517-26582. PMLR, 2023.
Neal Patwari and Alfred O Hero. Manifold learning algorithms for localization wireless sensor networks. In 2004 IEEE international conference acoustics, speech, and signal processing, volume 3, pp. i1--877. IEEE, 2004.
Paraskevas Pegios, Aasa Feragen, Andreas Abildtrup Hansen, and Georgios Arvanitidis. Counter- factual explanations via riemannian latent space traversal. In NeurIPS 2024 Workshop on Sym- metry and Geometry in Neural Representations.
Bo Peng, Ye Wei, Yu Qin, Jiabao Dai, Yue Li, Aobo Liu, Yun Tian, Liuliu Han, Yufeng Zheng, and Peng Wen. Machine leaming-enabled constrained multi-objective design of architected materials. Nature Communications, 14(1)66630, 2023.
Scott Pesme, Giacomo Meanti, Michael Arbel, and Julien Mairal. Map estimation with denoisers: Convergence rates and guarantees. arXiv preprint arXiv:2507. 15397, 2025.
Phil Pope, Chen Zhu, Ahmed Abdelkader, Micah Goldblum, and Tom Goldstein. The intrinsic dimension of images and its impact on learning. In International Conference on Learning Repre- sentations.
Herbert E Robbins. An empirical bayes approach to statistics. In Breakthroughs in Statistics: Foundations and basic theory, pp. 388-394, Springer, 1992.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi- cal image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pp. 234-241. Springer, 2015.
Sam T Roweis and Lawrence Saul. Nonlinear dimensionality reduction by locally linear embed- ding. science, 29065500) 2000.
Hiroyuki Sato. Riemannian optimization and its applications, volume 670. Springer, 2021.
Sholom Schechtman, Daniil Tiapkin, Michael Muehlebach, and Eric Moulines. Orthogonal direc- tions constrained gradient method: from non-linear equality constraints to stiefel manifold. In The Thirty Sixth Annual Conference Learning Theory, pp. 1228-1258. PMLR, 2023.
Gianluca Serale, Massimo Fiorentini, Alfonso Capozzoli, Daniele Bernardini, and Alberto Bempo- rad. Model predictive control (mpc) for enhancing building and hvac system energy efficiency: Problem formulation, applications and opportunities. Energies, 11(3):631, 2018
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. 2020.
Peter Sorrenson, Felix Draxler, Armand Rousselot, Sander Hummerich, Lea Zimmermann, and UlI- rich Köthe. Lifting architectural constraints of injective flows. arXiv preprint arXiv:2306.014 018433 2023.
Jan Pawel Stanczuk, Georgios Batzolis, Teo Deveney, and Carola-Bibian Schönlieb. Diffusion models encode the intrinsic dimension data manifolds. In Forty-first International Conference Machine Learning, 2024.
14

---

Published as a conference paper at ICLR 2026
Rong Tang and Yun Yang. Adaptivity of diffusion models manifold structures. In International Conference on Artificial Intelligence and Statistics, pp. 1648-1656. PMLR, 2024.
Wenpin Tang and Hanyang Zhao. Score-based diffusion models via stochastic differential equations. Statistic Surveys, 19:28-66 64, 2025.
Joshua B Tenenbaum, Vin de Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290655000;2331-22333, 2000.
Austin Tripp, Erik Daxberger, and José Miguel Hernández-Lobato. Sample-efficient optimization in the latent space of deep generative models via weighted retraining. Advances in Neural Infor- mation Processing Systems, 33:1 1259-1 2020.
Bart Vandereycken. Low-rank matrix completion by Riemannian optimization. SIAM Journal on Optimization, 23(2):1214-2236, 2013.
Singanallur V Venkatakrishhan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In 2013 IEEE global conference on signal and information process- ing, pp. 945-948. IEEE, 2013.
Enrico Ventura, Beatrice Achilli, Gianluigi Silvestri, Carlo Lucibello, and Luca Ambrogioni. Man- ifolds, random matrices and spectral gaps: The geometric phases of generative diffusion. arXiv preprint arXiv:2410.05888 2024.
Pascal Vincent. connection between score matching and denoising autoencoders. Neural compu- tation, 23(7):1661-1674 2011.
Trieu Minh Vu, Reza Moezzi, Jindrich Cyrus, and Jaroslav Hlava. Model predictive control for autonomous driving vehicles. Electronics, 10(21):2593, 2021.
Shiyu Wang, Mariam Avagyan, Yihan Shen, Arnaud Lamy, Tingran Wang, Szaboles Marka,
Zsuzsanna Marka, and John Wright. Fast, accurate manifold denoising by tunnlling riemannian optimization. In Forty- International Conference on Machine Learning, 2025.
Jan C Willems and Jan W Polderman. Introduction to mathematical systems theory: a behavioral approach, volume 26. Springer Science & Business Media, 1997.
Jan Willems, Paolo Rapisarda, Ivan Markovsky, and Bart LM De Moor. A note on persistency of excitation. Systems & Control Letters, 54(4):325 2005.
Yue Xie and Stephen J Wright. Complexity of proximal augmented lagrangian for nonconvex opti- mization with nonlinear equality constraints. Journal of Scientific Computing, 86(3):38, 2021.
Zhigang Yao, Jiaji Su, Bingjie Li, and Shing-Tung Yau. arXiv:2304.07680 2023.
Manifold fitting.
arXiv preprint
Boda Zheng, Abhijith Moni, Weigang Yao, and Min Xu. Manifold learning for aerodynamic shape design optimization. Aerospuce, 12(3):258, 2025.
Yufan Zhou, Changyou Chen, and Jinhui Xu. Learning manifold implicitly via explicit heat-kernel learning. Advances in Neural Information Processing Systems, 33:477-487, 2020.
A CONNECTIONS AND COMPARISONS TO PRIOR WORK
A.1 CONNECTION TO DENOISING AUTOENCODERS
In the following we show how our results connect to learning denoising auto-encoders. The follow- ing theorem is shown in Alain & Bengio 2014):
15

---

Published as a conference paper at ICLR 2026
Theorem 8 (Theorem in Alain & Bengio 2014 Let μ be a probability measure R and consider a > O the reconstruction problem
min TDAE(F) for LDAE(F) Buwn e> r(x: + ae)ll
If is absolutely continuous with (Lebesgue) density p. the optimal reconstruction function r² sat- isfies
(x) Ee- (oep)pp((eee((-re- σe) o+t x V + olo²)
Ee-Nio -oo) oe)
at any I E Rd with p(x) > O.
To see how this result relates Theorem T. note that function r here precisely our "approximate projection" no, since the first above expression corresponds to the posterior of observing æ under the convolution the prior p with (00,, I), i.c in our notation
r (x) EV3,o (x))
Therefore, if j has the Lehssgue-density p (i.e. the data manifold would have full dimension == d), then Alain & Bengio (2014) show that
o² lo = V logp(æ) - o
In partieular, p(x) 0, then the last equation states that
lim 2 Por = a
Theorem D on the other hand deals with the case k d, i.c. when the manifold is truly lower- dimensional. In this case, has no density p w.r.t the Lebesgue measure on Rd and it holds that
lim + o V log T(x),
where π denotes the projection operation onto the data manifold MM. While Alain & Bengio (0014) does consider manifold support in their discussion (see (Alain Bengio 2014, Section 3.4, Figure 5)), note that the population data measure in this case is not absolute continuous w.r.f. the Lebesgue measure on Rd. and hence log p(z) has no meaning. Consequently, results of Alain Bengio 2014), c.g. the expansion above regarding r not applicable. Theorem | can therefore be seen a (uniform) generalization of Alain Bengio (2014) an order of Ö(a), establishing the asymptotic convergence of r and its Jacobian for + 0 as
(z: (x) + Ō()), (+t)())= <') Ō(6).
A.2 TRIVIALIZATIONS AND NATURAL GRADIENTS USING PARAMETRIC MANIFOLD LEARNING
In this section we compare our approach to latent space optimization. For this purpose, we first describe the latter problem in terms of trivializations of (I) using local parametrizations: Given finite atlas A (1,, V of M with open V; C R and U. C M and charts di; Ut - V, clearly holds
min f(x) min min (()x.
TEM
i= zEV,
In the following we will focus inner subproblem by fixing single chart (æ.U, V) and consider
min f(x) min zEU
min f()(2))
(16)
sev
While for the above equality we only need té to be surjective (i.e. (VV) = 44, for equiva- lence of first-order stationary points of both problems this, note that stationarity x. E 24 for the left hand side in equation [6] equivalent to grad v f(x.) PT., M Vf(x.) and stationarity of any 2. V 1-102 for the right
need s to be submersion: To
16

---

Published as a conference paper at ICLR 2026
hand side is equivalent to (z.) V((())) = Vf(x.) = 0. In particular, for these two notions to be equivalent a-priori for any f and 2., the columns of ez(),) must span Tz, M, which implies that d is submersion. In particular, if db is chart, then the stationary points of both sides in equation are in one-to-one correspondence and only then both problems are equivalent from an first-order optimization point of view.
Two possible ways to solve the right hand side in (6) are the latent space gradient flow (LSGF) and the pull-back Riemannian gradient flow (PBRGF) given by
غ = - V ()) = - (=)) ff(())))
(17)
2 = (s)
(18)
respectively, where (2) = x is the metric tensor of và at E V. Here PBRGF is obtained from LSGF via change in inner product on the tangent space Tw(=) M and is equivalent to the Riemannian gradient flow in the sense that for (z( it holds that
= - -rad (z)
While (17 and 18
both possible approaches solving ①I and have been applied in the past Pegios et al. Tripp et al. (2020), there are few challenges present these methods:
Parameterizing the full atlas of manifold is hard: The atlas A of a manifold typically consists of multiple charts (1,,, Us Vo). The number of charts as well as their span across the manifold is typically unknown and hard to estimate a-priori.
Intrinsic dimension of the manifold unknown: The dimension of the latent space Z must correspond to the manifold dimension, the estimation of which is not a trivial task.
Generative model training: The implementation of (88) or (17) require & to be invertible and the inverse coordinate map b to be well approximated by a neural network. It is not straightforward to align this requirement with the standard training objectives commonly used in generative models, such as the likelihood-based loss in VAEs, the denoising score- matching loss in diffusion models, or the adversarial loss in GANs. Further approaches such M-flows Breehmer & Cranmer (2020) require & to be essentially normalizing flow, which is fundamentally incompatible with common architectures such as U-Nets or transformers.
Expensive update rules for PBRGF:
The update rule of PBRGF requires (i) computing the matrix (z),), which in turn involves the explicit evaluation of the Jacobian (f)x). and (ii) inverting Gg)( Both operations become computationally expensive when vi is implemented by a neural network even for a moderate intrinsic dimension k.
Limited latent space validity region: Since Z is representing the local coordinates of single chart, the validity of the model vb is only given in a bounded domain of 2. Opti- mizing outside of this domain requires changing charts. since otherwise the obtained latent points will have no meaning in the data-space, Estimating the validity domain of t is difficult a- a-priori Notin et (202))
In the above light, the advantages of the algorithms DLF and DRGD described in Sections and are as follows:
Easy to parameterize with neural network: The central object in approach is the projection operator onto the data manifold M, which is provably C2 in the tubular neigh- borhood of any O manifold. Assuming M is compact, standard universal approximation results ensure that the target map π can be well-approximmted by a single neural network
Dimension-agnostte algorithms: Both DLF and DRGD require no knowledge about the manifold dimension k.
A Alignment with denoising score matching: DLF and DRGD use an approximation of π via the score function of Pa (up to scaling and a residual term of 2)), so the denoising score- matching loss used to train diffusion models is directly aligned with our target quantity.
The Fuler discretization of the continuous dynamics (18) is also known natural gradient descent Amari
1998
Theorem 1) in the classical Riemannian optimization literature
17

---

Published as a conference paper at ICLR 2026
Our approach imposes architectural constraints: in our experiments, for example, successfully use both MLPs and U-Nets.
Simple and efficient update rule The update rules (8) and 12 are highly compatible with modern hardware and software stacks and requires only a single forward-backwarc pass of the network per iteration (see Remark
Convergence guarantees in the embedding space: Theorem 3 and provide explicit convergence analysis with guarantees on both optimality and feasibility, stated directly for the original data manifold in the data space, where convergence is semantically meaningful.
A.3 RECOVERY OF ITERATE COMPLEXITY OF RIEMANNIAN GRADIENT DESCENT
In this section we show that the case of for the Riemannian gradient descent. To see this, let us recall the standard complexity guarantee for Riemannian gradient descent under general smooth non-convex objective f:: In (Boumal 2023 Corollary 4.9) it is stated that (taking (2) = πτ + v) as the retraction defined on some ball can be picked uniformly over x M) if f т satisfies the
Theorem 5 recovers the corresponding result
(x) M, where for us pullback L-Lippshitz condition
L ||ull for all (x(x) Z
f(T(x + v)) f(x) S (
VS(x)+
then for step-size 7h t the RGD generates iterates (IR, Uk (1x that satisfy 22L(f(ro) f.)
min || grad f(x))l
k=
VN
provided that Up E (T) for all h: = , 1.
Now we set the following Theorem 5 in our paper: $ n, o = 0. The pullback Ly-Lipschitz condition is satisfied under our assumption L = Lip( V f) for different constant Lo = Lip( V( T)), which can be estimated in terms of L, C, t and M (see the proof of Theorem 5). The rate obtained in Theorem 5 reads then
4D
min || grad f(x))11 N Nomin (1 7max/
where 7min (1 - Lo Tmax/ / can be lower-bounded by constant constructed from and Lo- More- over, 2D is simply an estimation of f(zo) - f. via the triangle incquality (see the proof of Theorem 5). Thus in this case both theorems guarantee an ((1/VN) wwr---ooooooooaaooeeee convergence rate with similar constants.
B MORE ON MANIFOLDS AND DISTANCE FUNCTIONS
In addition to the notation introduced in Section 2.1 we note that TM is also the largest T 2 0 such that the map [p. v) E N M| ||a||| < r} - Rd (p, v) + + v is diffeomorphism. By tubular neighborhood of radius т E (0, TM] we mean a set of the form T(r) = {p + pe ||: pll denote the distance function that d(x:) 2 (x)) The second fundamental form of .M at a point E M will be denoted by I, and a symmetric bilinear map II, T x T,M + N, M intrinsic the manifold M. Fixing some u N₂ .MM also define the directed second fundamental form np Mx T M + R: (v, + (U, w), ") .M' The Weingarten map Sp at point E M in the direction E M defined the unique self-adjoint linear operator SW Tp M - T M
E M, e
M,M,llu
7). Moreover, for
Rd let (
such that (w, (1)) T M
IA (v, w) for all v. E Tp M. A useful operator that has been studied
in Abatzoglou 1978): Breiding & Vannieuwenhoven (001); Alvarez et al. 2004 is
H2- IT. M + (x())= T (())) M + m(z)) M.
(19)
In Breiding & Vannieuwenhoven/ (202T) it has been shown (see Lemma 10) that H. is invertible on T. Now we have the following useful identities that hold for x E T
n (x) = TT H-1 PT
V d((x) = x - T ((x)
2 d(x) = I - <'(x) = -- I M H PTic M
18

---

Published as a conference paper at ICLR 2026
For any p MM the map prp M - MpM: g Pr, p) is a local diffeomorphism p with inverse up defined on B TM/1 (0) =- Brmm4(0) n TT M.. Following Divol 2022 we define Mx(T, M) as the set all C* -anmfiofodds M as above such that TM > sud pllce I M. For the class Mx(T, M). manifold M E Mx(T, and E M, we denote by t always the inverse of the orthogonal projection prp (also called Monge graph chart) restricted to the particular neighborhood B
which will be (isometrically) identified with the ball M) (0) _ IR. We make frequent use of the following useful result from Divol 2022). Lemma 9 (Lemma A.1 in Divol (2022). Suppose M E [M,T,, and p e M. Then up B TTM
M)/4 (0) M is well-defined, -m-smooth and the following holds:
(i) For all r S min(TM, M3/4 /4 it holds that B, (r(p) (p (0)) (.)
For z E TT
M/44 it holds that || S (1)p - pll S 8||=11/7.
(ii) There exists a map W, : BTM M)/4
N, with W₂ = and such that vn(z) = - (z) and [x)())] & M|| for all E B Biin M))
(iii) For Ga, B BT.M
min (TM (0) + 2 v det det (=) (x)T((())))))))))))))))))))))()) it holds that Gv, (0)
and VG (0) O.
Note that for the graph chart vi we always have
(0) II, M. -= Lp(,)),
(20)
and hence (1,) S and [ || 1/TM-
B.1 PROPERTIES OF
We study the invertibility and boundedness of the operator 19 For this purpose, let us recall first the definition of the (normalized) curvature radius of M at p the direction of u E N,
TH I|oll (u. v) max(eig(S" (Sy//lll U (0)).
p(p u)
max HET, I (2,)20
||1|2
If 5% has only non-positive eigenvalues, then M is curved away from the unit vector u and thus the curvature radius is infinite. Moreover, we define the (normalized) curvature of M to be
xpp (u) = max eig eig(S%)/lll) for u E N,
and the maximal curvature by
RM = max M (M.
(M))
(21)
We have the following Lemma 10. The operator 19 is invertible on T(TM). Moreover,
19 satisfie
[[P(())] ||- (1
[<<x) - x11
p( (2)))
1 - - [x)
||
(1)]l/TT
and if sb is local parametrization of M with (0) = p. then
Po(x)= = =0) (0(( (0) +
Σ
(p - x),, v2, (T)
i=1
In particular, for any T E (0, TM), it holds
sup (T)(r)
H
< x..
Theorem 4.3) the quantity ||==1|| || has heen shown to equal the condition number of certain critical point problem associated with M.
In (Breiding& & Vannieuwenhoven 2021
19

---

Published as a conference paper at ICLR 2026
Proof. In (Breiding & Vannieuwenhoven 2021 Lemma A.2) the following condition has been tablished: Let = {(a. E R. x XM a - p E M). Then S is diffeomorphic to the normal bundle N . via the diffcomorphism d NM - S (v. + (p - IN, (M), Con- sider the operator IT : S- + TR" : (a.p) + a and the domain where its differential is invertible W = W= p). E S TT' (a.p): ((aa)) S RR" invertible). Then H2 is invertible iff (xr, (x)) E W. But II + N R (v.p) + p IN, (N)) being a diffeomorphism (and thus having an inveritble differential) is precisely the condition in the definition of the tubular neighborhood TT The formulas for (a()r) in local coordinates well as (x)]l are given in (Abatzoglou, 1978 Theorem 4.1, Corollary 4.1). To see that (x is bounded on T(T) for any T É (M, it sufficies to note that in the tubular neighborhood T (TMM) we always have [[= -x] (x)) (a)) and
that TT) is compact subset thereof. The second inequality follows from 1/p(p, u) S K.M-
Now let us derive bounds for the quantity [P)(((())) Po(x)]] when T(TM).
Lemma 11. If E T(TM), then
( S NM (z) [x)) (1
p(T(2)), ( T Txx)) )
-
[x(x) - x ==11
||0 - -<))]l
x 1-- [x)x(x)) 5
|| - - +(xx(())) ||TMM
|| [==))=m
1- 1 xxx))//
Proof. We clearly have for x ET
Po( - Po(x) II, ( M
Moreover, (I - H-, )) Tx))) M + +Tx(s) M is symmetric with eigenvalues
eig(T-H,71 { S T+ 15 eig( eig(S { -I| æ01s
[I(()) x1I 1 eig(S%(a))
where x-π(x) |--)( Thus
5
[P((T())) P
max ÇEeig(S
1 - || π (<)]1<< ||0- x())]
< <||$ (1
| x1| -
(()(T),; - -
| ||- (x)];
which shows the first inequality. The second and third inequalities follow from 1/p(p, u) S KM S 1/TM.
0
The next lemma establishes a bound on the Lipschitz-constannt of Po of some M E &M&(T, M) in terms of M and TM-
Lemma 12. If M E Mc(T,M). then
[P)(())]| s ( (1- =ET(F)
2 3 (1++ + +M) T-T/TM TM TM
1
2
sup
Proof. Let x E T(r) where r E (0, TM). First let relate the quantity T(x)) and fixed chart
& : V + U with ((x) E U for all x E W for some (small enough) open set W C Rdd The map
hz(2) = 2|1 (22 attains its minimum in some ze(2) E V and if (2,2x)
12(2)
(=))((())) - z). then (z)()),,) = O. Moreover, we have T(x) = (2 (x)) and hence for
w E Rdd
pf(x)[w, = ="> "" (z
w|],
and hence
20

---

Published as a conference paper at ICLR 2026
By the implicit function theorem bound the derivatives of 2 in terms of derivatives of FV. Indeed, we have (here again E Rd and E R are plader vectors to express the differentials)
= Fg (x) x) (v)e + F (2s (I)) w),
0 1 (z 2) (x) wIl F. (z) x) (T)(w)), (=)[[]]
(22)
+ 2F (z)) (x) zn (x) wi + (z (T), (,,
with
(2, = (())e) (0(z) z, v" (=)(e,,
(2,,)) (=)
(1,,,,) = 3 (x + (0(2) - x,
(x, u, = (=)
(2, u. =0 0.
Now we set 3 vp for the graph chart defined on BT
MJ/ (0), where p == x(x). this
case (x) and || (0. x) IIl
[))(())] (1 (MM -1 by Lemma 10 Solving for
z's (x) and 2 (x) in 22 using 20
with |_ S 1/ TM and taking norms yields then
11 1- T/T/MM
3
_ -=(()))) TM
[z)(()) T-7/TM
(TM TM
+ +
Plugging this back into the upper bound of [x)]( and using once again finishee the proof
The next lemma establishes a lower bound on the distance between point in tubular neighbor- hood and any other point that is sufficiently away from 7 (x) uniformly in r.
Lemma 13. Let T E (0, TM) and let > O. Then
inf dist ((s) distm(x)² *ET(r) 2
2 TT - TT -=1 > 0.
2
TM +
Proof. Let 7' = I+TM (0. TM). Take x E T(r) and set p x(x) and y = p + 7'n with z-p
2
7 Tz-pl Then E T(TM) with (y) = 7' and thus By) n M {p}. Then, since
M (x) C Rdd (B) U B,(2)), sufficient to show
inf ||x gl22 |r pl
TM - T
TM + T
To see this, note that E Rdd (v) U (p)) implies 11y- 9/l 2 r' = || pll and ||g pll 2 n
as well as (p q.y - q) > lp - qll². Then, abbreviating L
| 1-
E (0, 1), we have x
ty + (1 - t)p p + tly p) and
|| - 211 - ||I - 2| = [ty + (1 -11 p|
2 [t(y - q) + (1 t) -) ly - g
= (1 - - + 28(1 - t) (p 4. g)
2 (1 - t) all²
TM ,
TM + T
21

---

Published as a conference paper at ICLR 2026
B.2 DENSITIES ON MANIFOLDS
:- -u S M, the volume measure VolM is Given manifold M as in Section B and any chart : uniquely defined by
f 1
Volm(E)
Go(z) for E Cu Borel measurable.
Here (z) = Vee (z) and VolM is independent of the chart 1h. If P(R)) is absolutely continuous w.r.t. Vol then for the density u(y) avolm (v) on M we have the local representations
de
Yp
m(b(z)) d(v Volm) (=),
X(=) = G=(=)u(v(z))
# with À = - the pullback under the chart vb with density A(z) = d
(2) w.r.t. the Lehesgue
amy
measure on V. In particular when .M E (M M) with graph chart v/p, we have
A(0) = u(p), VA(0)= = grad E p M
(23)
C
MORE ON THE STEIN SCORE FUNCTION
In this section derive the representations (4). For 2 let V(0,02Ta denote the Gaussian distribution with mean zero and variance r² and for a > the Gaussian (heat) kernel in the ambient space Rd by
1 /20
Po (x)
Zo
Zo
(2
Then Pe = NN, (0) * has for o > 0 a fully supported C density, denoted by Pa as well, given by
1 Po v) du(y), e Rd.
Po (z(
Let Vr be the posterior of observing under Pa with prior μr, i..
Po (x) / Po(x - d)d))y) S Rd Borel.
V2 (E)
We have the following representation
Lemma 14. For each representations (4 hold.
Proof. Clearly
(po ( / (I (e-y)po(x du()),, ZER
σ2 J
1 σ [t((t)] 0 M (x (1 w)(x- v) po - 1) du(s))
V px
E Rd
and hence
V log (og) Vpo(x)
(x(x) σ² (I - -EVzo),
which yields the first representation in 4. Further
V'p((x)
1 2
(x v)(x- -v) ())
Po
JM
22

---

Published as a conference paper at ICLR 2026
and hence
V Po (x) = P2po(t) (x)
(x
Po (2)
1
/ (x -1((z v) duz
1
0²
M
(x
(EV,,)) Rvs,a)
09
1
(Jer yu dva.o (EVz, (BVe.)T)
1+ 0²
σ²
1
1
I Cov(us,a)
which yields the second representation in
D VARIANCE-EXPLODNS DIFFUSION MODELS
In the variance-exppa variance-exploding scheme of score- score-based diffusion models Song et al. (2020); Tang & Zhao 2025). one considers the following stochastic differential equation on the finite interval [0, 7]:
X, We, LE [0, T) Xo p.
(24)
The solution X, of this SDE is distributed according to X, Po()), where o²(1) = =22 To sample from j one exploits the fact that the reverse SDE
dx, = 2(T V)V -)) 22(T -jdW,, [E(0,]) Xo Po(T) (25)
satisfies X Po(T-t) and in particular Xo Po = μ. Here the initial distribution is approxi-
mated by Po(T) = (00 IT) and unknown score V VloPoor-- is learned by minimizing the conditional score matching loss defined
LCSM (8) Fu--Uoi F Be (1)2((s V log (2: mx))?
(26)
and which admits the unique minimizer (x(x)= Vlogp=(r). The loss LeSM be evaluated because V Po Pa (x (p) = (2 mo) is known explicitly.
E PROOF OF THEOREM
E. NONASYMPTOTIC LAPLACE METHOD
In this section we derive some non-asympiotic еπоо estimates for the Laplace method Hashorval al. (005)) Hwang 1980), which concerns itself with the asymptotic of integrals of the form
[ f(e)e M2=) dz (o (a-00)
for some functions f.h: V + R from open set V C R where h is non non-negative and attains V. While there have been results providing such an error
a unique minimum in, say, = = estimation to the first order expansion Inglot & Majerski (2014): Majerski (2015); Lapinski (201)) we need an error estimate up to the second order expansion, as provided Theorem T6] For the sake of completeness
also include the statement and proof for the first order expansion in Theorem 15
Before we state the results, we need to introduce some notation and conventions. For any f e C and h C² (V), respectively we write
Sr(2) f(z) f(0) == o(1)
(x)=)= M(=) - h(0) Vh(0)]
1
(²h(0)= =20 o([[=113
2
23

---

Published as a conference paper at ICLR 2026
for its first and second order remainder terms and abbreviate
(w) [y(ow) and X (w) Xh for > 0.
Let us note the following: If additionally / CV)) and g E (V))), then (y(2) O[I and
(2(2) O(1==19) for 2 + 0, i.. there exist some 77 and C > 0 such that (0) C V and
[s{(2) - C||=|| and (x) - C||2|| for all E (0).
(27)
In partieular for 0 and B [0, 0o) [0, oc) it holds
sup ||e]] [min(n/a,8(())
[s7 (w) - Co,3(a) and
sup (wv///aminn///,33(
[g(ww] S CoB(o)²
(28)
For a non-negative function h with a global minimum at = - 0 we define
(6)
inf h(z) - h(0).
= 126
We say that h satisfies a local quadratic growth condition in (0) S V if there exists a c > 0 with
h(2) h(0)= 2 c||=|² for all E B,,())
(29)
and that admits minimum separation outside of (0) in V if there exists a A > 0 with
h(2) h(0) 2 A > 0 for all E VA (0)
(30)
Clearly, if h satisfies 27 b. then h also satisfies (2) (for potentially smaller 7)) On the other hand,
the global assumption (30) cannot be inferred from properties of around = = 0 alone. Together,
29
and (30) imply that 7h be bounded below by
(6) 2 min(c6",AN)
(31)
which is crucial to obtain quantitative bound on the convergence of the Laplace method. While we state the following results in this full generality, our application we have (0) V and thus 30 is not needed, with 27 and (29) holding globally on V and implying that
(d) 2 co².
(32)
Theorem 15 (First Order Laplace method). Let h E C²(V) with Vh(0) =0 and f E (V) for some open B, (0) C V C Rt. Let E r0.r) and C, A be such that 27). 29 and (30) hold.
c,
Then for all a (0, with a = min(y? (log(2)(22C))
with it holds
M0) del
/ f(=)e Me) de s(0) s F.
(33)
Zo
with E (Ph(0) and Ea(o;f,h,r) h)r) O d for 0 given in 34
Proof.
Without loss of generality we can assume h(0) = 0. For brevity we also write ZFF
E for the normalizing constant. Let a : 10, + [0, oc) be any function and denote B (O(o)(O).
Za
We can split
2 /, f(z)e -2)=) dz=
1
1 1
f(z)e- -2 dz +
S f(z)e (2hh()) dz.
2 JVB(O)
2 B(o)
We will craft a 10, x0) + 10, oc) in such a way that the first integral vanishes and the second converges f(0) for + O. For the first integral have
f(z)e 12=()) dz - z
- 2 flllay)
2 .(B())
For the second integral let us write f(z) = f(0) + (2) and split
1 (
f(z)e 12)(2) dz 1(0) e M(2) dz + 3 B(o)
1
(z()e m2)) dz
E B(o)
23 B(a)
I:](o) J(o).
24

---

Published as a conference paper at ICLR 2026
First let estimate the difference between I(o) and f(0). We have, using the substitution and abbreviation (y(a) = (V B(o))/o,
(T(a) -(0)
S(0)
1
h(rm)) dw -1
zF By(o)
(1
f(0)
(c - dw
cbw Ew dw
zF
RM Bulo)
f(0)
+1(0)))(1(())))
Estimating 11(0 we obtain
[T(0)] - z
sup (x -1 1/. wEBv(o)
(R² (n) S R B ( Then, via the substitution To estimate [2( let us note that w we obtain
SeNmoaiio e - -- Ew dw
( e 3
S L
|2((0]] S
Ew dw
JB (ain
S 2Tp Go (Amin() min(r,
< _ k2T)) 221/2(R (a
e b lull du
with the Gaussian tail
Go(R) (2T))//2 (221BB((0)
e Ellll du.
Next we estimate J(a via
/ (rlowwe im Em e (x]()) dw |(())s 2 By(o)
1
(C((oo)))
le xx[((w) 1/ -
sup (WBvva)
sup wEBv(a)
In all, obtain the following error estimate
C M())
25 / f(z)e f(z)e- Mx(=) d f(0)
1
+ f(0) zF
([T(()) +
75
e xx((w) -1 -1/ Go Aminn(2) min((;;
f(0)
sup wEBu(o)
Ky (r(ww)))
sup le xx 1| 11+11 (y()
sup Bv
Thus we observe that we need to select in such way that for +
i (z
- exp(==((a((o))o3)
(ii) Go (Amin (E) minfra((o)) - 0
25

---

Published as a conference paper at ICLR 2026
(iii) suPweBr(oi (w)) +
(iv) By + 0
Let us make the Ansatz a = = a for some 3 (0, oc) -> oo) with 8(a) - x for o + 0.
Taking B(o) [og(o)), noting that le
s 2/2l for [ar| log 2 as well as using (31). (28) and
the fact that [loog(o)), (logg())) v 0/12 when (0, 1). yields the following estimates
1 exp(
0² ((a(o))) S
( 2))
minfr,a(o))
Go
min
Go (Amia (E) min(ro [log(r))))
sup w (y(c)
xf. (())
1/52
S2 we Bvbr(a)
sup [x((w) 2Co 2Co
sup we (v()
(s((ow) S CaB(a) = Collog(o)),
when 0 <o <<ss== min (log (12))) Plugging all these estimates back yields 33 with
1 E(a)
exp( spp(-minflloge Ao -0-
z
((0) (2Col Go Go (E) minfro (log(o)l)))
(34)
Collog(o)[(1 (Collog(o))))
= 0(o log(o))) Ö(o).
The next result provides an nonasymptotic estimate of the Laplace method for the second order expansion. We will need version that allows for a non-zero gradient of f. as long as it is contained in the subspace ((g-99) 1 q E R"] S Rd x R
Theorem 16 (Second Order). Let f.h E (V x V) be such that f(0) O and Vh(0) 0 for some open B, (0) C V 1. for some r > Further, suppose that h is additively separable as M(2) T(2)) + T((2)) and V
9 for some q E R Additionally, let rj E (0, and
C, A > 0 be such that 27 29 and (30 hold, where this time (7 applies also for h = f. Then
for E (0,7) with = mi, (log(2)((2C)))) it holds
e (()) Vdet Zo
f(z)e- M2=) b dz o² A(f.h)
S EE(oo,
(35)
Jxv
where (2)(0) with E2(o; 0(0'(log(o))5) for 0. Here the second order coefficient h) is given by
A2(f,h) tr (r(2-~~2)(0))
(36)
and Ea(a; J,,,r) is given in 37
Proof. The proof is similar to the proof of Theorem 15 with the modification for the product space
R T TR.. Again we can assume h()) = Let a : [0, oo) ~> [0, o0) be any function and
denote B(o)= = Ba(a)(0) (0) x Ba(m) C R" x R" and D = x V. Write
-2 75 | ((z)e 1 dz
1
1
[(=)e M+)) dz + 0² 25 /
1
[(2)e M2)) M)))))))))))) dz.
02 75 (iB((o)
JYnBIO)
26

---

Published as a conference paper at ICLR 2026
We will carefully craft such that lim,-0 a(o)) and the first and second integral converge to and 2 (rEE-V''(00) respectively, with explicit error bounds. For the first one we have
(a(o)))
1 1 ta Zg ((B(())
/ f(z)e het) d2
0z
Now consider the second integral. By substitution B(o))/o Bv(o) By(a) with Bv(a) = (V (V (a(o)())) as well as x2 (w) = 2 X for E 4f, we have
and abbreviating By(o) =(V (Vn
1 1
I
f(z)e f())==the -()))=+((s)) dz
a² Z
(nā(a)
1
1 f(0) T (25(() xf(w))
h(ow) dw
2 Bola)
2
=: (o) + I(o) + J(o).
Note that (o) vanishes due to separability of h and the particular form of Vf()):
o 02; = (Bplo) ff(0) w e- - Mow dw
BButoo I JBvio) (qJux W1 gJw,y (Twyee" h(ow)) = M(rww)) dwi dw2
Next, let consider the difference between (() and 1 V22
f(0)) given by T T2 f(0) - Jt Iw xx(w) dw tr (r(
w
By(o)
2
v2 f2((0) - Iu' xū(w) 1) dw
<
2 w
7 (a)
22 [(0) e - Tow dw
+
7F R2 R2mi,BBl0 (p())
2
[((( + |2(0)))
Then
(E)
sup
le (i
-11]
zF
Ren
2/1""
2
1
w V22
-e- - Ew du
1/ 1/31|
le xxx(1) -11
sup WE Bx(a)
Furthermore, by a variable substitution == / w we obtain
[ <
[
1/2 1/2
llll du
Z1 R3H (PPA///BB500)
< 2 [O))|E [0))|E-1/233 1/2V2 (0))=111))))
where
M(a) = Z,
21la1 du. | b1lall
By(a)
27

---

Published as a conference paper at ICLR 2026
Noting that for Σ= = (2)(0) it holds B)0))(()=>> =\l(BBuoo)
Bu(o) x
Bu(o) and B Xumia (22) min minfraa(e)))000
Bu(o) // (V (V Ba(o)
and hence
L(o) 22G (Amin(>) (Σ) minlm;
with the Gaussian second moment tail
G2(R) = (2T)n/2 /B(BBR(0)
| - =1llll du.
Now consider the final expression
/ x(())-- dw
[J(o]]
7 Bolo)
bx((w\))
/ Zi oo(o)
(w du
sup weBp(o)
Let us further estimate
e (w)) dw = 2 Bz(
/
e - Tim xime (h) dw
z BB(o)
le (e=x(((x) (w0) -11) 1.
sup wEBp(o)
Hence we see that we need to pick a in such a way that for o - 0
() (25)-10 exp((>x + 0.
(ii G2 Ammin (2)) min( +
a
(iil) sup B (w)l (w) +
(iv) sup suP b (w)) +
Again take a(o) = ()) for B(a) log[o)]. This yields
exp( गe 7h (a(()))) S
Do 2));
min(r,
G2 Amin
< G2 (Amin (E) min(iro- [log(a)))
sup (c)
(ww) 1/ 2 sup wEBv(o)
[xx(w) < 2Co = 2C)log(o))
[x] < Co C|log(a))
sup weBv(o)
when 0 < σ m7n (log(2)/(2())
Combining all these estimates one obtains, in
the proof of Theorem 15 that (35) holds with
det Σ
F(()
Zao2 exp(- exp(=min(clog(o) Ao ())))))]][|
+ 120 (oCollog))))) G2 (Amin (E:) min(ra [log(o))]
+ (1+
= (() Ö())
(37)
28

---

Published as a conference paper at ICLR 2026
E.2 LOCAL VERSION OF THEOREM |
First we prove local version of Theorem
Theorem 17. Assume that there exists an subset U S M that C d--iffeemorphic an open subser ofV S I& with (0) C V via vb V + U for some 0 S k S d and that μ < Volu, where Volu is du
the volume measure U. Moreover, assume that p(-) Vooli (2) Let z Ra be any point with (0) =: p - arg min pim || - -pll and infyemyull y11² -1 - pll > 0, and such that
(p - x). (1(0) 0. Σ = / (0) (0) + Σ
(38)
2=1
Then there exists some σ > 0 such that for all a E (0,7)
(i) it holds that
2 (o h x|| + h) + (llp =l| + (I; E))
[Vu,o - pll -
H(P)
(ii) it holds for Po = =0)= that
1 Cov(Va,)) Po
S 4d
E2 h,r))
n(p)2
x²
[[[luul2+ + Iulll) T);
Fa h,r, + (o;; Σ)
o² + 12 m(p)
|Poll
p(p)²
with h, fo and fi given in 39p.
J2 and he in 40 Mell
Bymally x||', T(o;2)
Ve to expl 2).
P1 h), max I=1....d
En Ea Ji; h,r)
E2 h,r)
max 1.j=l....d
Ez(o; f2 h,r),
with E1 (o; h,, given in 34 and F2(a;F2',h,, h,r) in 37 respectively. Remark 18. This theorem only requires M to be locally manifold, namely at U.
Proof. We denote by == -b1 be corresponding pullback measure on V and denote its positive density again by E (VV). Note that X(0) p(p) by (23). Moreover, let v and Po (0)= (0) for the rest of the proof. First we need investigate the non-asymptotic convergence of following two integrals:
/
Po (z)
Po (z - du(y)
MM
X(2)
c- 2 [+=1==((=))1
dz / Polo - v) du(x))
Zo
So
Ro
and
Po (z) - 2z) / (3) xx -)) (dddd
(EVs,o
M
2 ln ([=)] dz + I
{(=)(1(=)) -x)
(u x)polo 1) du(v)
Ze
Mie
R1
29

---

Published as a conference paper at ICLR 2026
We can write each S and R₂
= = Za / fi(z)e = dz, R4 = aa=Jmna Jm 9:(0)oal y) du(z),
where
h(2) || +z))]?? fo(2) = 1(2), fi(z) (z))(((2) - x)
(39)
as well as Bo(y) = 1 and (v) y Then fo, fi E (V) and
h(0) 2 ||011? Vh(0) 0, <h((0) = E,
where the second equality is due to p (0) being the closest point from U æ. Also obviously fo(0) = X(0) and (0) (p x). Applying Theorem 15 component-wise, we obtain for (0. with o. (n C) given as in Theorem 15 with n). 0 dependent on fi and h, the
estimate
ez Ilvll det ES, f((0) s DD(o),
with
Do(o) = Ei Jo, Di(o) = ((:): = Vd max
Ei,o; I=.....
and where E₁ is given in Theorem 15 and fi is the I-th component of fi. Furthermore, estimate for i = 0,
Rall = Mun alu)po(x - v) saanoo)
Z Inne (ilu)eatell 1 an
S z [1ul 25.
J JMMzd (19(
Zo
where we have used that || y/1 2 8 + ||ua12 for all y E EM \.. Thus we get the estimate
2025
lea
|vll R. I
Zo
where the quantity ||a]|= = fuully - x du(y) denotes the i-th centered moment of j. Let us denote
1lll /det Σ and express
5, = f((0/J + F/J
Now we can estimate the distance to p by
S1 + R
R(0)J+R/J+ Rı
Evso - -p So + Ro
(p x) (o/) Ro
(p x)
[1(0) + Fi + JR
(p-x) (p-x)+11 (p x))
T To (p- x)
(p (p-x) x))
fo(0) + Fo + JRo
1+
1 To
where T; = -F; +
JJR.)A(O) (J&)
Let us bound
|T
1
(D(()) + Vdct
|| - 700 0.
A(0)
Zo
30

---

Published as a conference paper at ICLR 2026
Then Ty < 1/2 o h(f minfou, ō) with ō depending Do. E and 6, For this E (0,7) we obtain
Iu - pll
2 æ|I + 2/
Ei fo. h r)ppp x|| + (aEi(a; Si, 4,) + (llv x1| + Ilalli) del E 2070
2
X(0)
This proves (i). The proof of (ii) is analogous. For (z, 2) E V x V define
72(2, 2)
z) (0(2) 2)((
R(2, 2) 1 + -()
(40)
Then E (V > x V). () = 0 and
E
f(,000 0) ||+||
Vh(0, 0) = 00
v²T(0,0) = - = diagdaa E)
J1 Let l2 be the (j, 1) entry of 2- Then it is elementary, but very tedious exercise show that and h satisfy the conditions of Theorem 16 and that moreover the second order coefficient is given by
Ax h) = 2000(000000000000000000000 = (Po))
ie.. Ag = (A2212 M)) Po. Now split
Po (x)² Cov(Vz,))
Po (x: (Cov(()))) - Rvz,o (Bu/,,)) ))
= - (1y 1 x1 - x) (y - x)) () - x) ))pa s) D @ m)
= -rav Jz(2, 2) 1 2 (+=) ( ||= (()))()))1 d(2,3)
z3"
82
IMUMMAn (1y=))(2 x²) - (y- 2) ()) d(H @ ()
R₂
Applying Theorem 16 component-wise to fa obtain for (0,02), with 02 o2(17, C) and 77 C > 0 dependent Fz and h, that
(erz
|ull?
det 2 o² A2
[d d max Ls=t.....
E2(o; (2'
||aal
with Ez Fa(o:TTh,TTT r) given in Theorem 37
The term R2 the other hand can be estimated again
by
222
.28
|Ral S 311lppe
Jncentra (w-=x)) =)T - (v (v-x)(0 x) [ddddpp ().
72
Se 22 - -22
(1lm)12 + [|all?)
22
For S, it holds
52 o J2
31

---

Published as a conference paper at ICLR 2026
and thus
x² Cov(Ve.) - -Po 82/0² R2/0² Po
(So + Ro)2
1(0) Po/ + G2/(02 + R2/0² - Po
(A(0)/J + Fo/J -Ro)²
Po(x)++T Po
+ (1+Tb)²
T. (2 + Po
(1+ To)
with T- = + JRa) Again we have the estimate
d E2(a) |T < X(O)² o²
(det () l2 + ||all7) 2
X(0)2
Zz022
and for a E (0, (with a as before, guaranteeing (Tol < 1/2)
92 Cov(Vz,)= - Po
S 2| + P1
s 4d Eσ(0) X(0)2 0²
[A(det (1llll + [1]l1l]) e -
Do(a) + det Ee
202
+
- -22
Zo
Poll
X(0)2
Zz02
X(0)
E.3 BOUNDS IN TERMS OF DISTRIBUTION AND MANIFOLD PARAMETERS
In this section bound the constants appearing in Theorem 17 in terms of parameters of the distribution and the manifold M MT, when the chart is given by the graph chart
Moi(fM..M// M144 (0), E T(r), cording the proof of Theorem 17 need to consider the following maps
u V M with V TTp M
(0, TM) and p T(x) Ac-
h V R:
2 1 - (
(41)
as well as
fo V - : X(2),
(42)
f1 V Rdd :d:224 + ((=)(1)<)) 2),
(43)
Jz : V x V- - Rdddd (2, + {(=)A - (1))(())) - xx) - (1)((2) - x) 1))) (44)
where A = So # = prp #p.
E.3.1 CONDITIONS
AND FOR
First we show how to express conditions 27 and
29 in terms of the manifold parameters for
given in (41). To 27. note that for E V have by the chain and product rule well as
Lemma 9 (1) that
h(())l S 31
+ || <+>)llll"""
3MM² + ([æ - pl + #
3M2 (TM 8 I1=111MM
< (3M + FTMMM.
32

---

Published as a conference paper at ICLR 2026
Thus 27 holds with, say TM.M (M TMM and min M/4. Next, for (29) and
c(l)
? min [TM, M//4 compute for E B, (0)
1
<(=) h(0) 2
inf
h(3) h(0) -
inf
x
1 -
pll
TE BEV (0)
TP
TD M (0))
2
2
1
inf
||
all
12
pl2
B (p) 2
2
1 TM - 7
TM + 1|=112
In the first inequality have used 1/ (V)) B M (0)) S Bijll| (p), which follows from Lemma
||=||
9 (i), whereas the last incquality we have used Lemma 13 and æ e T(r). Thus (2 holds on
with FM-T and implies (32)
FM
E.3.2 LOWER BOUND ON Amin (2) AND UPPER BOUND ON det Σ
We first consider Amin (E) for E = (() with h given in 41 and when E T(T). By Lemma
10 and (00) we directly obtain
Amin (Σ) ||
11 pl|
TM
Due to 20 and the definition of the shape operator Sp. an upper bound on / Σ is given by
det E Vdel(IT,.M+SS))
p(p,p x) 1
|p x1
k/2
(1+
(++m)
k/2
< 2//
TM
E.3.3 CONDITION (27 FOR Fz
An upper bound for || T [| in terms of M and the third derivatives of du V where
(2,)), (()())) (v) =) (v) (3) x)),
can be obtained by direct differentiation. We don't pursue the complete derivation here and instead say that (7) holds with, say, c22) M,
E.3.4 UPPER BOUND ON (l) (M And 122ll AND [| (0)
We clearly have ( A(V) < 1. Next, due to l1y æ|1 < diam(M) llp xI1 s diam(M) for M, we have
l & [[fall = / [Iv -1- =1(dp())= = |ell s diam(M) +T,
and
Th FFlll S (1
-
Juxe [1(y x)(U - - -y-- [| ())
< ||al|2 + ||1
< (ddiam(MM + 7)2
Furthermore for v == (%2 E R** x & we have
T" 0),
X(0)2 (201e) 201 + (p- 2) (1,((1 02) 1p u2))
+ +2 V\(0) T (0) + 02) (p - z) - -2)
Hence, using 23 and the fact that ||p x [ w2)l - 5 for any unit vectors 01, U2 T MM we obtain
TM
T ( (0)1 (0. < 6u(p)2 87p(p)]grrdd m(p))
33

---

Published as a conference paper at ICLR 2026
E.4 PROOF OF THEOREM
We apply Theorem 17 to every point & E (T) and - T(x) with the graph chart &b = vp- remains to provide universal constants for the bounds in Theorem 17 (i) and (ii). First bound
(1; So: h, with Ej given in (34) By Section E.3 we can take TTM.M = 7 min(TM. M3/4,
CTM.M 2 TM+
TM-T and CTM.M max(CCl)
TM.M M..dpt and hence, due to [l S
24/2
exp (-CTM. log(a)))
Z
+ u(p) 20 (C Go((1 min (llog(()))))))
TM
+ +Crm.M.no + ",,,. o
Similarly, due to 11 & 11lll1. we have the bound
2
Ei(o;f1)
exp (- log(o))) (diam(MM) +
Z
+ u(p) M2 Go((1 - TM
_
min min( (log(r))))
+ C MG 20 0,00
Further, due to 1172 (( & |11/|2 + 1/l] and the bound on V (0) we have
2²/2
Ez(a:
exp(-CrM.M (diamm))))) + +)²
Zo
+ u(pp) (p(p) + 8TlIgrad p(D)
kC (Crm.Mmolloo(o))
T/TM
G2 (1- (1-2) - _ min (min(rru.mo-! [log((]])
TM
+ M log(o))))) + 2C, +2T (og(o))))
Now by Lemma 13 we can pick 8 IMA 2c 2 TM.
and thus
2k/2
CMM.MTMM TM.
T(o;E) <
Zo
exp
x²
Finally, |Poll < (1-7/TM) Note that since u(.)) (9(M) and supp μ M,, there exist positive lower and upper bounds on u((-) and its derivatives on M. This shows that there and o depending only M and diam(M) (as k) and μ such that ( hold and finishes the proof.
E.5 A USEFUL COROLLARY TO THEOREM
The following elementary consequence of Theorem 1 will be useful in the proofs of the results from Section
Corollary 19. Suppose that for some v Rd Rdd and T E (0, TM) and S T have
|| - m (x))l < for all T(T).
(45)
Then v(((T)) C T(e). Moreover, for x E T(E) we have [|v(x) xll & 2e for E T(r) we have [v(x) x S 2r.
Proof. From 45 and the fact that 7 (xx) M we clearly have v(x) T(e) for any T(T).
ie. vT)(r) T(c). The last claims follows then from the triangle inequality [[(()) - x
[[v[z) - [<)] [[())) - xl|.
34

---

Published as a conference paper at ICLR 2026
F APPROXIMATE RIEMANNIAN GRADIENT FLOW WITH LANDING
In this section C | (lTmml)l and L. = Lip (T)T()))) denote the supremum of V and the Lipschitz constant of Vf on ((T)) for some manifold M.. We will often assume the following approximation condition function v v CCCC (R(:
||v()) (+))] <.,
[|v((x) Po(x))l < for all x TTT).
(46)
We will also abbreviate the (negative) right hand side of the dynamics ( by
G²(x)) = v' (z) + 7( v(z))
(47)
F.1 STATIONARY POINTS OF AND OPTIMALITY CRITERIA
We need analyze the meaning of approximate stationary points of { for the optimization problem
Lemma 20. Suppose that E (0, TM) and a > 0 are such that (46) is satisfied. Moreover, suppose that vTT(r) T(r). Suppose that 7 (0, 5 >> 0 and that v. (T() is 8-approximate stationary point of i.e.
|G)((,))] 88
(48)
Then for pa Tz.)) it holds that
Igrad f(p.)l n) +2 F/TM C+8.
(49)
1- --//TM
Proof. We have the estimates
- x.)Il S ne.
and
< [|B86fcc))///(ttt))--4/////1/11
||Po(z.) [ftl::////fvvoo.--VVVVVVVVVvvvv (Ro.))llVVVTT(()
where have used that (o(I(x))) is orthogonal projection and hence has unit norm. By
Lemma 25 applied to 21 = v ((v(z. and 22 = (v((z.) I=), Lemma II and the fact
that Pat F. (T)(())) grad f (Ir))) the inequality 49 follows.
F.2 LIE DERIVATIVE OF MANIFOLD DISTANCE
Lemma 21. Suppose that T E (0, and a > 0 is such that (46 is satisfied. Then it holds
- -(Vd(x), G*(()) S d(x) + cCC + 11) / 2d(x) for all TT(T).
Proof. We have for T(T)
(x ---a0 (x), -G)(z))
2n) + (r(x) -
+n (())
< -
(n]lm(x) +
< 2 (|V////1ll//
35

---

Published as a conference paper at ICLR 2026
F.3 EXACT LANDING WITH
Here we establish the following noiseless version of Theorem 3
Theorem 22. Consider the flow (8) for a and di Z with r00) E [TT) for some (0, TM). Then the solution æ(() exists for all times 2 0 and is contained in (()). Moreover every ac- cumulation point v. of this flow satisfies grad M f(x.)
O and is at most T away from the point p. (xx)) .M with llgrad f(p.)ll I Lt. Moreover if n >> O. then I, Pa E -M and grad f(p.) = 0, i.e. every accumulation point is critical.
Proof. Note first that
dt d(x) = (or (x)) = - 7//0 T(x)]| = 22))()),
i.c d(x(t)) = e- 2nd (0)) and hence x(t) E T for all t 0. Moreover, if > 0 then d(x(t)) + for t - cc. Let us now consider ff o Lyapunov function for (9) with o We observe that
d
dt f(((z) = ff((())) T (o(x)V
(Praa)) M Vf(x), H 2 Pr, M Vf(z))
= Ilradd (x)
with ||ol|3,, = =111 u1 T M This shows that f(((t) is non-increasing for - oc. Since T bounded, lim f(x(t)) is finite. Moreover, we have
/ Igrad ((xx(())) di s(x(0)) <0.
(50)
Notice that (=(() | 2 0} S M being compact. Since T + R: : Igrad f()) Hr
continuous, it is uniformly continuous ((dist, This, together with 50
implies, by
Barbalat's lemma Farkas & Wegner (2016), that lim Iggrad (x((t)) Aoco
U. Clearly any
accumulation point æ. of (x(t) 2 0j satisfies H grad f(x.) O, i.c. grad fxx Then (I..) E M is an accumulation point of (π(x(t)) 1t2 0} S M and
[[rrad (T(x.)) = |grad f(T(x.)) grad f(x.)) s L[[ - xx11==11 LT.
F.4 PROOF OF THEOREM 3
Note first that, by Theorem if hold with e + c', then are satisfied with ='+ (olooo))())) for a (0, (T, M, n.). In the following will write analyze the manifold distance of (8) We have by Lemma 21 whenever xr() E T(T) that
for the latter quantity. First let us
d ã7 d(>)
(dd(),,)(0 < 22)) v(C+))))))
where the right hand side is non-positive iff
(x) > 2 11
+1 =: To
Note that < < e < 2T (1 (1+1 C/n). In particular if ( (a(()) [To, r), then [(dist (xrist (x()))
invariant w.r.t. the flow (8) and (t) exists for all 2 0. Moreover, by similar argument in the standard proof of Lyapumov's 's direct method Khalil & Grizzle (2002). for each > there exists some > such that for all t > T it holds that x(t) E T + 8).
is gradient field, i..
Vg for some function g C¹ (RR), then V(x) f(v(x)) (|[==]|//2 g(x)) satisfies
(v(r) = (z). Similar as in the proof Theorem 22 we can take G, as a Lyapunov function to
obtain
d -
V(x) [c)())))
di
36

---

Published as a conference paper at ICLR 2026
along the dynamics 8). By similar argument (xx) is non-increasing and every accumulation
point of (x(t)
t 0) 0} satisfies (x.)) 0 and belongs to T(7)). Note that the condition
v(T(r)) - T() in Lemma 20 is satisfied when e S т. Applying Lemma 20 with = 0 and == To yields 11 and proves Theorem 3 for the case when is gradient field. Now consider
the more general case when v is not necessarily a gradient field. In this case consider V(x) = f(v(z)) + (x) obtain
d
dt V(x) (V=((())V(((()) + n/ -(===(=))),))))))))) (())).
S
s <-1G 2 + + G.) (z)]l
Thus, if [C](z)] 2 TEE have di Vxx) < 0. By barrier function argument Khalil & Grizzle 2002) this implies that every accumulation point 2, of (=(t) t 2 0} satisfies [G](z+)] S 7/ and belongs to ((r)). Again applying Lemma 20 with d = THE and F To yields and finishes the proof.
G DISCRETIZED RIEMANNIAN GRADIENT FLOW AND DESCENT
In this section we provide proof of Theorem 5 as well as analysis of the approximate Rieman- nian gradient descent and the discretized landing flow. As before L Lip(V V (lT(TM)) denote the supremum of V and the Lipschitz constant of Vj on ((TM) for some manifold M. Moreover, the following bounds will be useful:
| [VLl((r and
Lemma 23. IfT E (0, and (46) holds for some e 0, then
[I S 1 TT/TM =ET(r)
sup
Additionally, ife E (0, then [G][()]] s C]]v(+)]] (1)]= +)ll((z) x|| S C(e + 1
+ 2n/ for x E T(T)
T/T/M
Proof. Via triangle inequality, Lemma 10 and the definition of C.
G.1 DISCRETIZED APPROXIMATE RIEMANNIAN GRADIENT FLOW
In this section we analyze the discretized version of 9 specifically the corresponding gradient
descent
Øh+1 24 7 (Z))
(51)
for some sequence of step sizes J The selection of 7h can be inferred from any standard analysis of gradient descent for F7 to guarantee that all accumulation points of the resulting se- quence {ar k= stationary points of FV, Since we can only interpret stationary points of F7 in terms of our original problem (T) when they are contained in tubular neighborhood T(T) for some the step-size ensure {xx ) T(r).
(0, TM) (see Lemma 20).
drive conditions
Theorem 24. Let T E (0, T22) and o > be such that 46) holds for some e (0, .(+ッ) Then if [0, Ttubular] with
+ - ( n)e
Atubular (c, (E, T, = min
2(C + 2) 2
4(C(E + +) +) 317r)²
and xo E (()), the iterates of the discretized
51 belong to T(T).
Proof. First we show that It E T(T) implies 06+1 E T(r). If ss T(7/3), then, since
nullV F2 (xx))l I r/2 by Lemma 23 it follows J&+1 E T((). Hence assume ay E T(() (T(/33), ie. 22d(zz) 7/2. Let us write
1
d(xx41) = d(xx) 7h (Vd(zz), (xw)) + Tlll( Po) Ir
37

---

Published as a conference paper at ICLR 2026
where we have used V d I Po and 24+1 T(3r/2), since again 2 VV (zx) (xx)ll s T/2. Then by Lemma 21
(V d(r),, V=) < -2))(()) + <(C + (C++))2d(())
Now, Lemma
23 and the fact that (1 - (1-37 v 4 imply
8 = 2 11(I - Po)1 VVF < 4(C( + 4) + 377)²
Therefore we have
d(xx+1) d(xx)) S (- + (CC+))+))))
+ (2 + 7*.6),
where right hand side is non- non-positive for all d 2 2/2
1 2 +
2(0 CC n)TTE + n6 $o.
or, equivalently,
1 / 1
- - 1
-2(
7
8
4
2
If these satisfied, then d(æ)41) - dxxx) S /2/8 and therefore xx+1 T(T). In all, f0
T(T) provided that Zo E T(T).
G.2 PROOF OF THEOREM 5
Let us write d' = c K(T, M, (lggo) 1/2.. First we note that {ax S T(c') soon as To E (T//2), because v(((T) C T(')) and æx - <wv' E (T),, since
since
Y S 20 + (1-11170 T/TM)-11 S 2z (xx)]
where first inequality is due to (1 T/TMM) -1 < and the second inequality due to Lemma 23 For brevity let us denote yk Fy - -/ v (xg V (xxx) and 2 3* Po(x)) ( Also let
Lo = Lip(V( T)) ). By Lemma 12 and (1 T/TM)- -1 s 2 it follows that
2 _ ++M + -1
3
Lo < (T) (1- (1-7//m)-12 ' s 8C
+ 2L.
TM
TM
Then, since 11a x (IHi)| S for k: 2 1, we have
f(zx+1) - fxxx = f(v(yx)) - f(xx) f((((z)) - -((r( + +f((()))) +((( f(
f(*x)
+ vi Lo (((())) 1
S n (H)
+ - - o L| - (xx)]])
Lo
2 ( 1
(2C' Lo L))e'
Now let 7 E ('min 7maxx S (0. 2.) Then summing over 1,...,, N yields
1 I'mnn 2 7maax Lo)
1 IP(((>
f(x1) f(EN+1) + (2C + Love(C + L))e' L))e
N
N
Now note that by Lemma and the fact that (1 (1-6/TM))+ < it holds
[lgrad
- (i
which implies
1 N
4D <
Tuun 2 Ymax Lo)
Igrad v f(T(xx)))
N -C +22C L07 + (C+ L))e'
This finishes the proof.
38

---

Published as a conference paper at ICLR 2026
G.3 AUXILIARY RESULTS
G.3.1 A LEMMA ON NORMS OF ORTHOGONAL VECTORS
Lemma 25. Let L E R be orthogonal and > O. If there exists some E R with | 2 S
and |/u =|1 S then ||a||| S 2e, |yl| S 2e and ||=|| S V2.. on the other hand for some à > 0
there exist some 21, 22 E IR" with || - 21 S e, yyy - 22 S 6 and ||| - 22ll S 8, then || - +
Ilyll 2x and ||=||| & v 22 + 8/v2.
Proof. The first claim follows by geometric considerations. The second follows from the first by noting that the midpoint = (21 + /2) satisfies || zı S 6/2 and || - 22ll - 6/2 and hence || I e+6/2 and | |y I e+6/2, i.e. | & 2e+6, |llll s 2e+6 and || S v2c+8/v2 O
G.3.2 GAUSSIAN TAIL BOUNDS
The following elementary lemmas give simplifications for some of the constants appearing in proof of Theorem and and Theorem 16
Lemma 26. holds for R > 0 that
Go(R) (2T)a/22 J (R(BB(())
e =1llll du S 2c R2
and for R 2 2n that
22
| du S l(m/2) c R22² -
G2(R) (2T) 3z Bx())
In particular Go(llog(o)) ((o') and G2 ([log(o) (aa)) for any | as 0
Proof. See Majerski (2015).
H FURTHER NUMERICAL EXPERIMENTS AND IMPLEMENTATION DETAILS
H.1 OPTIMIZATION OVER (nn)
H.1.1 IMPLEMENTATION DETAILS
In our experiments, we discretize the flow equation g using the Euler scheme with a step size of lstpp = . 10-1 and set the landing gain 1) 3 103.
Data generation: In experiments, we take Q symmetric A E gnxa with (0, -)-erriss.
and randomly sampled
Score architecture: We use the following score architecture:
So (X) (X).
with (a1 x.) Raxn and (y1 Vn) = 5, (X), where
y/ MLP (1; a)) for -,...., n,
and MLP fully connected multi-layer perceptron with ReL.U activation function, / layers of width w. The features r = r(X)= e Rm are
(X) = tr(Q,XK,XT), j. m
where Qj K2 E Bxn are learnable weight matrices (shared for all i = 1,...,n). For n 10 we take 4, w 512, m 128 and for n 20 we take 4, w = 2048, m = 512.
Diffusion and training parameters: We train minimizing equation 26 with the Adam opti- Tni with T 3 and 10 and cosine learning rate scheduling from 1r = 10 3 to 1r = 5 10-5, We use Nepochs = 10000 and Nopochs = 50000 epochs for 10 and 20, respectively.
mizer, early stopping (i.e.
39

---

Published as a conference paper at ICLR 2026
H.2 DATA-DRIVEN REFERENCE TRACKING
H.2.1 BENCHMARK SYSTEMS AND TRACKING GOALS
Benchmark systems: We consider two classical benchmark systems: The double pendulum and the unicycle car model Valle (2006). For the double pendulum the state is x (01, 41, W2) with W, = 6, gravity g. masses m1. m2, lengths 41, 12, dampings d1, d2, control torque u applied at joint 1 (first pendulum), and AO == 02 01.
0 Wl
02 W2,
181 02
DØ =, 0 (()))
T = (0)
with
(m1 + m2) (2ll122 cos A
m2l1122 cos A m2l3
dh 0
M(0))
D
d2)
m2l112 sin A w.? m2l1/2 sin A w
(m1 + m2 m2)gl1 sin 01
C(O,0) 0)
G(0)
m2gl2 sin 02
We pick the output y (01,02) and set my = 1 =g 1, m2, unicycle car model the dynamics is given by x = (x. 8) with
0.5 and d1, d2 0.1. For the
x = e cos( f)
s Ô = w
and the input u = (v, w) and output y x (x, y, 0). Here (x, y), v, is the car's position, velocity, angle and angular velocity, respectively.
Tracking goals: For the double pendulum system the goal is to track reference trajectory r via the first joint angle 0 and we pick the optimal control objective f to be 15 with
Q
10 0
100=
= 0.01.
while for the unicycle model the goal is track positional reference (xx.Ty),ie i.e. we pick
10 0
(00
00.01 0
Q -
0
10
0
R =
0.01
0
Here R is some small penalty on the input u to keep it bounded during the optimization.
H.2.2 IMPLEMENTATION DETAILS
Discretization: We discretize both continuous-time dynamics à: and discretization step At to obtain the discrete-time dynamics (13)
fcont (I. u) via the RK4-method
ka f Jennt (I, u)
At (k1 + 2kg + 2kg + ka) where
h/2 Scont (x + Atki, u)
f(z, u) = x
ka Scont (x + Alkiz, w.
Aia Scont (x + Atk3, u)
We use At = 0.1 for the double pendulum and At = 0.05 for the unicycle model.
Data generation: To generate trajectories, we use i.i.d. random inputs Mi Unif[- for the double pendulum and पı = (Uk, wN) ~ Unif(0, 1] Ø (0, and horizon of N₁, 100 for both system. We use Naata = 50000 trajectories for the double pendulum and Naata 20000 trajectories for the unicycle model.
Score architecture: The score architecture a 1-dimensional version of the standard UNet architecture [Ronneberger al. (2015) with sin-cos time-embedding [Song al. (2020) and
40

---

Published as a conference paper at ICLR 2026
residual connections, where the different input-, state- and output dimensions concatenated and treated as additional channels. The down- and upsampling convolutions are done w.r.t. the temporal dimension and channels.
Diffusion and training parameters: Same as in Appendix H.I.I with this time Nepochs = 50000 training epochs. For the DRGD step-size we pick fixed step-size of 0.001.
H.2.3 EXPERIMENTS
In Figure 4 we present the objective value evolution for the experiment from Section 6.2
coot
cost Mninum train data
Mininum train data
Objective Value
Objective Value
10--
16'
500
1000
1500
2000
2500
3000
500
1000
1500
2000
2500
lterations
lterations
Figure 4: Denoising Riemannian gradient descent: Objective value f vs the iteration count j for double pendulum (left) and unicycle model (right). Note the logarithmic scale on the y-axis. The current cost (blue, dashed) the objective value f(u;, yj) at the current (in general infeasible (u,, y,) & Mioo iterate, while the true cost (orange) is the value suuy y grue). with guue obtained by simulating 13) with input u
In Figure we show optimized trajectories for two other reference trajectories. Note that we have set our iteration budget at N 4000, while the objective is still decreasing. How to accelerate the denoising Riemannian gradient descent without losing feasibility is core question for future work.
position
position
2.50
Reference
1.5
1.50
Mirium data Belenence
0.25
1.00
Figure 5: Denoising Riemannian gradient descent: Unicycle car position (right) with the optimized output trajectory y" (blue, dashed), the true system trajectory y irue (orange), the initial trajectory yo (green, dotted) and the reference trajectory r (red)
41


Reviews:
[{"review_id": "o5Wz4CjOcO", "reviewer": "Reviewer_NnVy", "summary": "The goal of the paper is creating a novel framework that ties theoretical Riemannian optimization tools (retraction, computing the Riemannian gradient, projection onto a tangent space, etc) to practical methods that rely solely on the data distribution (e.g., diffusion models). The proposed work is impactful and impressive in numerous ways. I highly recommend this paper to be accepted.", "questions": "1. It makes sense that authors need a $C^2$-submanifold. Do you have any intuition on what this means in terms of limitations, both practical and theoretical?\n\n2. How does the proposed approach compare to well-known baselines, such as VAEs and GANs for modeling manifolds? This would help position the work in the literature.\n\n3. Could the proposed approach converge slowly (in terms of convergence rate) in certain scenarios?\n\n4. Is the approach easily extendable to incorporating methods like heavyball (momentum) into the optimization process? If so, what are some potential challenges that you foresee?\n\n5. Are there certain manifolds where the proposed framework could struggle?", "limitations": "", "rating": 10, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "Strong points: the paper is well-written, with clearly made arguments, well-structured problem, meaningful theoretical as well as experimental results. The paper provides a comprehensive literature review and acknowledges its own limitations.\n\nThe paper supports results with rigorous mathematical proofs, as well as intuition-based explanations, which I appreciated very much. This work has many, many strengths. To name a few: this work provides a rigorous theoretical analysis of how Riemannian optimization tools can in fact approximate practical operations such as projection (for denoising purposes), and how this can be done by learning from data. This works provides a massive step forward in this literature by connecting two fields that have been struggling to converge -- theory of high-dimensional geometry and diffusion models.\n\nThe significance of this work is huge. It makes a rigorous connection between theoretical and practical fields (Riemannian optimization and diffusion models), which is very high-impact and rare to come across. Creating such impactful and meaningful work requires a vast and deep knowledge of the field, and one can dream to be at that level as a researcher!\n\nAuthors make several innovations. Impressively, authors prove that the quantities in equation 4 converge to manifold operations as $\\sigma \\rightarrow 0$. For example, the proposed link function provides a way to approximately project points onto the unknown manifold wihtout explicitly knowing its equation, using the score function learned from data. Authors also clearly demonstrate how to use these operations for optimization by proposing two concrete algorithms that tell us how exactly their work can be implemented. Lastly, one non-obvious strength is in experiments, where the authors show that even though theoretical results are strong when $\\sigma \\rightarrow 0$, experiments work even when $\\sigma$ is nonzero.", "weaknesses": "The authors heavily rely on the claim that classical algorithms require knowledge of the manifold to be able to use Riemannian optimization tools, however, this claim is not necessarily true. Plenty of works, especially [1], develop provable, feasible, accurate, learnable, and inference-time efficient algorithms that create machinery to use Riemannian optimization tools without knowledge about the manifold, while having access to noisy samples only. In fact, the noise level can also be relatively high. Authors are also missing another key work [2] which provides a provable and accurate denoising method without knowing the underlying manifold. While I find the results of the paper impressive and useful, if one were to make the paper even stronger, I would recommend re-positioning / re-phrasing the argument that \"classical Riemannian optimization cannot be applied directly, since they rely on explicit knowledge of the manifold.\"\n\nLastly, this work heavily relies on the fact that a diffusion model will learn the score function well. If this is not the case, the proposed framework will not be useful.\n\n[1] Wang, Shiyu, Mariam Avagyan, Yihan Shen, Arnaud Lamy, Tingran Wang, Szabolcs Márka, Zsuzsa Márka, and John Wright. \"Fast, Accurate Manifold Denoising by Tunneling Riemannian Optimization.\" Proceedings of the 42 nd International Conference on Machine\nLearning (ICML), Vancouver, Canada. PMLR 267, 2025.\n\n[2] Yao, Zhigang, Jiaji Su, Bingjie Li, and Shing-Tung Yau. \"Manifold fitting.\" arXiv preprint arXiv:2304.07680 (2023).", "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": "waLT33uME3", "reviewer": "Reviewer_eZg1", "summary": "The paper studies constrained optimization on data manifolds when the manifold is only available implicitly via samples. It smooths the data distribution with a Gaussian of scale σ and defines a “link function” $\\ell_\\sigma$ whose gradient and Hessian can be related to manifold operations. Theorem 1 shows that as $\\sigma \\to 0$, $x + \\sigma^2\\nabla log p_\\sigma(x)$ uniformly approximates the closest-point projection $\\pi(x)$ to the manifold, and $I + \\sigma^2\\nabla^2 log p_\\sigma(x)$ approximates the Jacobian of that projection with error $O(\\sigma|log \\sigma|^3)$ on a tubular neighborhood. Corollary 2 then states that, on the manifold, $I + \\sigma^2\\nabla^2 log p_\\sigma(x) \\to P_{TM}$. By this process, the denoising score approximates manifold projection/retraction.\n\nBuilding on this, the paper proposes two algorithms that only require inference-time access to a trained score network sσθ and its Jacobian–vector products: (i) Denoising Landing Flow (DLF), a continuous-time flow combining a Riemannian-like step with a “landing” term toward the manifold (eq. 8), and (ii) Denoising Riemannian Gradient Descent (DRGD), a discrete method that applies an approximate tangent step followed by an approximate retraction given by the Tweedie mean (eq. 12). The paper provides non‑asymptotic guarantees: for DLF, accumulation points are near the manifold and have small Riemannian gradient norm with bias proportional to the geometry-approximation error (Theorem 3); for DRGD, an average Riemannian gradient-norm bound with similar bias terms (Theorem 5). Theorem 20 shows that with exact manifold operators ($\\sigma=0$) the flow reduces to classical exact Riemannian gradient flow. Experiments on illustrate the approach. Figures 1 and 2 show optimized trajectories tracking references more closely than the best sample from the training set.", "questions": "1. The theorems assume a fixed small $\\sigma$ to control the geometry error. In practice, do you anneal $\\sigma$ during optimization? Or do you pick a $\\sigma$ that trades off stability and bias? Some concrete empirical guidance would be helpful.\n2. DRGD performs “tangent step + retraction” using s and s′. If we hypothetically replace $s$ with the exact projection $\\pi$ and $s'$ with $P_{T_x M}$, do your bounds recover the standard complexity guarantees for Riemannian GD (up to constants), or is there an intrinsic loss because you work in ambient coordinates with the landing penalty? It would be helpful to position Theorem 5 against known rates. \n3. For real networks, $s$ may not be conservative. Apart from the barrier‑type argument in Appendix E.4, are there practical regularizers or architectural choices (e.g., gradient‑field parameterizations) that make s approximately integrable?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The link function viewsystematically connect $\\nabla \\log p_\\sigma$ and $\\nabla^2 \\log p_\\sigma$ to retraction and tangent projection, respectively, with explicit uniform error bounds (eq. 5). This is a useful and reusable idea. \n\n* Theorem 1 leverages Laplace-method estimates to obtain uniform $O(\\sigma|log \\sigma|^3)$ control on a tubular neighborhood; the constants are traced to curvature/tubular radius and density smoothness (Appendix D). This is technically solid and extends beyond heuristic “score is normal to the manifold” statements. \n\n* Both DLF and DRGD require s(x) and Jacobian–vector products s′(x)v (Remark 4), which makes them practical with modern autodiff. \n\n* The analysis is structured so that (a) convergence of optimization with approximate manifold operators is established assuming an $L_{\\infty}$ approximation error, and (b) Theorem 1 supplies such an error by letting $\\sigma$ be small. This cleanly decomposes the problem. \n\n* The paper clearly distinguishes constrained optimization from posterior-guided sampling and explains why small “temperatures” in the latter can drift off the true manifold—motivating the explicit feasibility machinery here (Section 1.2).", "weaknesses": "* The dependence on $\\sigma \\to 0$ is strong in practice. Theorem 1 and Corollary 2 fundamentally rely on $\\sigma$ being small, and all optimization guarantees carry a bias term $\\sigma|log \\sigma|^2$. Training accurate scores precisely in the small‑$\\sigma$ regime is hard, so the theoretically clean limit may be brittle in practice. The paper acknowledges this by proposing a landing flow, but still the guarantees rely on $\\sigma$ shrinking (Theorem 3). \n\n* Uniform $L_\\infty$ approximation is a strong assumption. Convergence theorems (7) require s to approximate both $\\pi_\\sigma$ and its Jacobian uniformly over a neighborhood. Score models are typically trained in L2 (denoising score matching). The authors do note this (Remark 6) and defer L2‑type analyses to future work.\n\n* Comparison to classical Riemannian GD is limited. With $\\sigma = 0$ and exact operators, DLF reduces to Riemannian gradient flow with exact landing (Theorem 20), which is consistent with classical results; however, for $\\sigma > 0$ the paper provides average gradient‑norm bounds with additive bias rather than full complexity guarantees comparable to known rates for feasible Riemannian methods. A more direct side‑by‑side empirical and theoretical comparison to “Riemannian GD with a known retraction” would be  useful.\n\n* The experiments demonstrate improvement over the best training sample, but comparisons to explicit‑manifold baselines (where available) are limited to an “exact landing flow”. There are no comparisons to standard Riemannian optimizers using the true retraction on the same objective, nor ablations for the effect of $\\sigma$ and Jacobian accuracy in higher‑dimensional settings.", "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": "Hb6X23OW88", "reviewer": "Reviewer_xQd1", "summary": "The paper targets the problem of optimizing objective function on data manifolds, i.e., manifolds that are given only implicitly through the data distribution. The paper proposes to do so via Riemannian optimization, by obtaining the required quantities (Riemannian gradient and Hessian, projection onto the tangent space) through a link function associated with the data distribution which is associated with a diffusion score function.", "questions": "- It is unclear to me why it is necessary to learn a network $s$. Could the gradient $\\nabla \\ell_{\\sigma}$ and Hessian be computed solely from equation 3 instead?\n- How is the network $s$ trained? Which training data are needed?\n- How does the proposed approach compare to Riemannian optimization on learned data manifold endowed with a Riemannian metric?\n- How is the probability measure obtained in general? How is it obtained in the problems considered in the experiments?\n- How does the approach perform for more complex and higher dimensional optimization problems?\n- How should $\\sigma$ be chosen in practice?\n- What does “exact landing” refers to in Figure 3?\n- How does the proposed approach compare to gradient descent on the O(n) manifold? I know that the paper considers the manifold as unknown, but providing such a comparison would still be valuable to understand the influence of the manifold learning component.\n- It would be better if figures corresponding to main results, e.g., Figure 3, would appear in the main text. If space is an issue, I would suggest to reduce the introduction (which takes 4 pages out of 9).", "limitations": "", "rating": 2, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper targets an interesting problem, which is relevant to the field of manifold learning. The approach nicely combines manifold learning with elements from recent works in diffusion models and score matching. The paper is generally clearly written.", "weaknesses": "The paper does not mention an important category of related works in manifold learning, where the data manifold is characterized via the pullback metric of the decoder of a (deep) generative model (Tosi et al, Arvanitidis et al, Syrota et al, Diepeveen et al). In such cases, the data manifold is endowed with a Riemannian pullback metric, which can directly be used for computing explicit manifold operations for Riemannian optimization. In this case, the motivation of the paper (“In this data-driven regime, methods from classical Riemannian optimization cannot be applied directly, since they rely on explicit manifold operations such as tangent-space projection, retraction, or exponential maps”) does not hold anymore, as such quantities can be computed based on the Riemannian metric. I think that these approaches must be discussed, and some of the claims of the paper revised according to this perspective. \n\nIn addition, it is unclear to me if the proposed approach really qualifies as Riemannian optimization. In Riemannian optimization (Absil et al, Boumal, as cited in the paper), the objective function is defined on the manifold of interest and the optimization strictly takes place on the manifold, i.e., all iterate x lie on M. In contrast, in equation 10, the second term belongs to the orthogonal space to the manifold, and pushes the iterate towards the manifold, i.e., the iterates do not necessarily belong to the manifold of interest. \n\nThe experiments provided in the paper are simple and conducted only on low-dimensional spaces. I would like to see how the approach performs on more complex and higher dimensional optimization problems. The fact that score functions and diffusion models scale well to high dimensions could allow the approach to scale well, and could be presented as a strength if shown experimentally.\n\nAll in all, the performance of the approach seems very dependent on the value of $\\sigma$. In Fig. 3, only $\\sigma=1e^{-5}$ leads to a good optimization, with significant differences compared to $\\sigma=1e^{-4}$. It is unclear how the approach would perform for different values of $\\sigma$ in more complex optimization problems.\n\nGiven that the link function seems to be available based on the data distribution, which is assumed to be given, it is unclear to me why it is necessary to additionally learn a network $s$ to represent the gradient of the link function. \n\n\nReferences: \n- Tosi, A., Hauberg, S., Vellido, A. & Lawrence, N. D. Metrics for Probabilistic Geometries. in UAI (2014).\n- Arvanitidis, G., Hansen, L. K. & Hauberg, S. Latent Space Oddity: on the Curvature of Deep Generative Models. in International Conference on Learning Representations (ICLR) (2018).\n- Syrota, S., Zainchkovskyy, Y., Xi, J., Bloem-Reddy, B. & Hauberg, S. Identifying metric structures of deep latent variable models. in ICML (2025). \n- Diepeveen, W., Batzolis, G., Shumaylov, Z. & Schönlieb, C.-B. Score-based pullback Riemannian geometry: Extracting the Data Manifold Geometry using Anisotropic Flows. in ICML (2025).", "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": "VG4HuTxBR0", "reviewer": "Reviewer_rhkJ", "summary": "In this paper, the author's introduce a new way to harness data manifolds described using score-based diffusion models. By introducing a so-called link function, whose gradient is related to the score function and thus to the projection onto the data manifold, they derive Riemannian optimization algorithms that work directly on the score-based description of the manifold. The authors analyze the expected performance of their algorithms under accuracy assumptions on the underlying score functions. Finally, they demonstrate the applicability of their method in two synthetic examples.", "questions": "- the authors should please check again their usage of \\citep and \\citet. In some places, it seems to be off.\n- I would maybe switch (8) and (9), so that the method is first introduced with the more accurate \\pi before going to the score-based approximation.\n- One reference I missed is \"What regularized auto-encoders learn from the data-generating distribution\"  by Guillaume Alain and Yoshua Bengio. It also analyses denoising autoencoders (although not considering an underlying manifold) and it would be nice to set it into relation with the present work.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "To me, the main strength of the paper is that it uses a representation of data manifold that is widely used in current machine learning (i.e. score-based models) to introduce practically usable algorithms on said manifolds. There is, of course, a lot of papers currently that try to achieve something to this end. Here, I enjoyed the clear combination of geometry and probabilistic components in the construction of the method. The resulting algorithms seem natural and actually usable in practice.\n\nFurthermore, I appreciate the provided theoretical results, which give an insight into what to expect of the method. At the same time, the authors do not oversell the impact of their results \n\nFinally, I found the paper to be well written. I could read the paper in one go without any road blocks even though I am not an expert on optimization. After this, I believe to have understood the method to a point where I could now implement it myself, which is huge plus in my book. I also enjoyed the discussion of related work and clear positioning of the paper.", "weaknesses": "The primary weakness of the presented, from my perspective, is the experimental validation. It is limited to two synthetic examples, one on O(n) and one discrete data-driven control. Since the authors themselves make the point that score-based models have become ubiquitous in modern machine learning, I would have hope that one of these would also be considered here. For example, for image generation models one could try to modify a given image to mislead a given classifier. There are certainly other creative examples one could come up with. However, I understand that this requires a different expertise and computational resources, and this is not a deal-breaker for me.", "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)