Geometric regularization of autoencoders via observed stochastic dynamics
arXiv:2604.16282v1 [cs.LG] 17 Apr 2026
Sean Hill∗
Felix X.-F. Ye†
Abstract Stochastic dynamical systems with slow or metastable behavior evolve, on long time scales, on an unknown low-dimensional manifold in high-dimensional ambient space. Building a reduced simulator from short-burst ambient ensembles is a long-standing problem: local-chart methods like ATLAS suffer from exponential landmark scaling and per-step reprojection, while autoencoder alternatives leave tangent-bundle geometry poorly constrained, and the errors propagate into the learned drift and diffusion. We observe that the ambient covariance Λ already encodes coordinate-invariant tangent-space information, its range spanning the tangent bundle. Using this, we construct a tangent-bundle penalty and an inverse-consistency penalty for a three-stage pipeline (chart learning, latent drift, latent diffusion) that learns a single nonlinear chart and the latent SDE. The penalties induce a function-space metric, the ρ-metric, strictly weaker than the Sobolev H 1 norm yet achieving the same chart-quality generalization rate up to logarithmic factors. For the drift, we derive an encoder-pullback target via Itô’s formula on the learned encoder and prove a bias decomposition showing the standard decoder-side formula carries systematic error for any imperfect chart. Under a W 2,∞ chart-convergence assumption, chart-level error propagates controllably to weak convergence of the ambient dynamics and to convergence of radial mean first-passage times. Experiments on four surfaces embedded in up to 201 ambient dimensions reduce radial MFPT error by 50–70% under rotation dynamics and achieve the lowest inter-well MFPT error on most surface–transition pairs under metastable Müller–Brown Langevin dynamics, while reducing end-to-end ambient coefficient errors by up to an order of magnitude relative to an unregularized autoencoder.
Keywords: autoencoders, manifold learning, stochastic differential equations, geometric regularization, Itô calculus, generalization bounds MSC codes: 60H10, 58J65, 68T07
1
Introduction
Many stochastic systems of scientific interest evolve, on long time scales, on an unknown lowdimensional invariant manifold M ⊂ RD [9, 22, 38, 5]. The effective dynamics on M arise, for example, after averaging or homogenization of a multiscale (fast–slow) SDE [30, 38], and the practical goal is to learn a reduced simulator that faithfully reproduces the slow dynamics, including drift, diffusion, and passage-time statistics, directly from data. Because the full system is stiff, long direct simulations are prohibitively expensive. A data-driven alternative is to run short-burst ensembles from strategically chosen initial conditions: each burst is cheap, embarrassingly parallel, and reveals local effective dynamics and local geometry at the initial point [12, 48]. The ATLAS ∗ †
Department of Mathematics and Statistics, University at Albany, SUNY ([email protected]). Department of Mathematics and Statistics, University at Albany, SUNY ([email protected]).
1
framework [12, 48] organizes this idea into a complete simulator. A well-distributed network of landmarks {xi }m i=1 is placed on M , and short-burst ensembles at each landmark yield local estimates of the ambient drift and covariance {b(xi ), Λ(xi )}m i=1 with statistical guarantees [12]. A Gaussian-kernel smoother then interpolates these estimates into continuous fields, and the simulator steps in ambient coordinates using linear tangent-plane models at the landmarks. However, three additional computational bottlenecks arise in the ATLAS simulator: 1. Curse of dimensionality. The number of landmarks m and the number of neighbors used in kernel smoothing both scale exponentially with the manifold dimension d. The root cause is the reliance on linear local models, namely tangent ellipsoids glued together to cover M . 2. Re-projection onto M . Because ATLAS simulates in ambient RD , every time step requires re-projection onto M via weighted oblique affine projections at neighboring landmarks; without this, sample paths quickly drift off the manifold due to accumulation of linearization error. 3. Per-step simulation cost. Each step blends D × D covariance matrices from neighboring landmarks at O(mD2 ) cost; in our experiments ATLAS is already computationally infeasible at D=201 with only m=200 landmarks and d=2. Autoencoders [29, 37, 39, 16, 8, 24, 42] address all three ATLAS bottlenecks: an encoder π : RD → Rd and decoder ϕ : Rd → RD replace exponentially many linear patches with a single nonlinear learned chart, and simulating intrinsically in d dimensions eliminates re-projection entirely. Several lines of work build on this idea. Autoencoder-based methods have been combined with deterministic latent dynamics [8, 34, 32, 16, 6] and, more recently, with latent SDEs [20, 31], but none use observed dynamics to regularize the chart geometry; the contractive autoencoder [40] and the pullback metric of [3] regularize geometry but without a dynamics-derived tangent target. Outside the autoencoder framework, the equation-free approach [25] learns SDE coefficients in externally supplied reduction coordinates, and diffusion-map methods [10, 15, 14] estimate generator eigenfunctions from dense trajectory data in a spectral regime. Our pipeline instead targets the sparse-landmark setting (N =50) and requires a decoder for ambient-space reconstruction; neither the equation-free nor the spectral setting is directly comparable. When the manifold is known, Riemannian diffusion models [13, 23] and normalizing flows [35] define dynamics intrinsically. No existing method jointly learns a geometrically regularized chart and SDE coefficients from sparse ambient observations. Training with reconstruction loss alone, however, controls only the L2 error of the chart map, leaving the first-order geometry of tangent spaces poorly constrained [16], especially in the sparsedata regime where only tens of observations are available. The resulting geometric errors amplify when the chart is used to simulate dynamics: the Itô correction depends on both the decoder Jacobian and Hessian [1, 2], so tangent-bundle and curvature inaccuracies propagate directly into the recovered drift and diffusion, and the chart cannot extrapolate reliably beyond its training set. Existing autoencoder theory [42, 33] achieves a squared generalization error of order m−2/(d+2) log4 m from point-cloud data but does not address this geometric gap. Sobolev training yields provably better generalization [46, 47], but requires chart-derivative labels; semi-supervised methods [43] improve chart quality when function values are observed, but treating drift and covariance as regression targets in latent coordinates yields chart-dependent quantities. Neither is directly applicable. Our framework addresses this geometric gap by exploiting a key observation: the ambient covariance Λ already encodes coordinate-invariant geometric information. Its range spans the tangent space [11], and the resulting tangent projector P = Dϕ g −1 Dϕ⊤ is invariant under reparameterization. Penalizing the discrepancy between P and the data-derived projector regularizes first-order geometry without chart-derivative labels and defines a function-space metric ρ strictly between L2 and the Sobolev H 1 norm. For the latent drift, Itô’s formula applied to the learned encoder yields an exact pullback target, avoiding the systematic bias of the decoder-side formula. This paper develops the single-chart theory for this program; extending to a multi-chart atlas is a natural next 2
step discussed in Section 6. Contributions. 1. A three-stage pipeline whose geometric losses (tangent-bundle, drift, diffusion) are coordinateinvariant: it extracts geometric penalties from Λ, fits the latent drift via an encoder-pullback target, and learns the diffusion under metric-weighted losses. 2. A generalization guarantee showing that an idealized ρ-ERM achieves the same rate as Sobolev H 1 training for chart quality, with controlled propagation to SDE coefficients and weak convergence (§4). Outline. Section 2 derives geometric conditions and the encoder-pullback drift target from SDE coefficients. Section 3 introduces the ρ-metric and geometric penalties. Section 4 establishes generalization bounds, the decoder-side bias analysis, and error propagation to dynamics. Section 5 reports experiments; Section 6 discusses extensions.
2
Geometry from Stochastic Dynamics
We work in the single-chart setting throughout: the manifold M is assumed to be covered by one coordinate chart. This restricts the theory to manifolds diffeomorphic to a bounded open subset of Rd , but the penalties and pipeline extend chart-wise to an atlas (Section 6). Let M be a smooth, connected, d-dimensional Riemannian submanifold of RD (d < D), parameterized by a smooth embedding ϕ : Ω → RD (the decoder ) from a fixed open, bounded domain Ω ⊂ Rd , with encoder π = ϕ−1 : M → Ω. We write Dϕ for the Jacobian, g = Dϕ⊤ Dϕ for the induced metric, and P = Dϕ g −1 Dϕ⊤ for the orthogonal projector onto Tϕ(z) M .
2.1
Itô transformation rules
Expressed in Itô form in ambient coordinates, the effective dynamics on M take the form dXt = b(Xt ) dt + η(Xt ) dWt ,
Xt ∈ M ⊂ RD ,
where b is the effective drift, η is the ambient diffusion coefficient, and Λ = ηη ⊤ is the rank-d covariance. Within a chart (Ω, ϕ), the local-coordinate SDE takes the form dZt = µ(Zt ) dt + σ(Zt ) dBt ,
Zt ∈ Ω ⊂ Rd .
We write Σ = σσ ⊤ for the local covariance and define the Itô correction vector q i (Σ) := ⟨Σ , ∇2 ϕi ⟩F for i = 1, . . . , D. Lemma 2.1 (Local to ambient). The ambient drift and covariance are Λ = Dϕ Σ Dϕ⊤ .
b = Dϕ µ + 21 q(Σ),
Proof. Apply Itô’s formula to each component X i = ϕi (Z). The first-order term gives Dϕi µ, and the second-order quadratic-variation term gives 21 ⟨Σ , ∇2 ϕi ⟩F . Stacking over i = 1, . . . , D yields the ambient drift b. The diffusion coefficient transforms by the chain rule as Dϕ σ, so the covariance is Λ = Dϕ Σ Dϕ⊤ . To invert this transformation and recover the latent drift µ from the ambient data (b, Λ), we need an identity that relates the Itô correction in the two charts. 3
Lemma 2.2 (Quadratic covariation under chart change). With ϕ, π, Λ, Σ as above, assume rng(Λ(x)) ⊂ Tx M for all x ∈ M . Then ⟨Λ , ∇2 π⟩F = −Dπ ⟨Σ , ∇2 ϕ⟩F . Proof. Differentiate the identity π ◦ ϕ = idΩ . At first order, Dπ Dϕ = Id . Differentiating once 2 i ∂ϕl ∂ϕk ∂π i ∂ 2 ϕk more gives ∂x∂ kπ∂xl ∂z m ∂z j + ∂xk ∂z j ∂z m = 0. Contracting both sides with a d × d PSD matrix Smj and recognizing the Frobenius inner products gives ⟨Dϕ S Dϕ⊤ , ∇2 π i ⟩F = − Dπ ⟨S , ∇2 ϕ⟩F i . Stacking over i gives the identity for any PSD S. Tangentiality (rng Λ ⊂ Tx M ) ensures Λ = d (PSD since Σ = σσ ⊤ ); S = Σ is therefore an admissible Dϕ Σ Dϕ⊤ with Σ = Dπ Λ Dπ ⊤ ∈ S+ choice and recovers the stated identity. Lemma 2.3 (Ambient to local). With ϕ, π, b, Λ, q as above, assume rng(Λ(x)) ⊂ Tx M for all x ∈ M . Then Z = π(X) solves a local SDE with covariance and drift h i Σ = Dπ Λ Dπ ⊤ , µ = Dπ b − 21 q(Σ) . Proof. Apply Itô’s formula to Z = π(X). The drift of Z is µ = Dπ b + 12 ⟨Λ , ∇2 π⟩F and the local covariance is Σ = Dπ Λ Dπ ⊤ . Applying Lemma correction yields µ = 2.21 to the2 Hessian 1 1 2 2 Dπ b + 2 ⟨Λ , ∇ π⟩F = Dπ b − 2 Dπ ⟨Σ , ∇ ϕ⟩F = Dπ b − 2 ⟨Σ , ∇ ϕ⟩F , as claimed. Computing the latent drift target. The local-coordinate coefficients of Lemma 2.3 coincide with the optimal mean-square approximation in the Itô-jet framework [1, 2]. The lemma gives two equivalent representations of the latent drift µ: a decoder-side form µ = Dπ[b− 12 q(Σ)] that requires the decoder Hessian via q(Σ), and an encoder-side form from the intermediate step of the proof that uses only the encoder derivatives: µ = Dπ b + 12 ⟨Λ , ∇2 π⟩F .
(1)
For a learned encoder πθ , the Jacobian Dπθ and Hessian ∇2 πθ are available via automatic differentiation, so the latent drift target µ(zi ) = Dπθ (xi ) b(xi ) + 21 ⟨Λ(xi ) , ∇2 πθ (xi )⟩F
(2)
is computable without touching the decoder. This encoder-pullback drift trains the latent drift network in Stage 2 (Algorithm 1). For an imperfect autoencoder the two formulas diverge: the decoder-side target carries a deterministic bias for any fixed learned chart (Proposition 4.6), making the choice between them consequential for trajectory fidelity.
2.2
Geometric consistency and coordinate invariance
The Itô rules above show how to convert between ambient and latent SDE coefficients given a chart. We now turn to the question of what the ambient covariance Λ reveals about the geometry of M itself. The relation Λ = Dϕ Σ Dϕ⊤ (Lemma 2.1) implies that the range of Λ must span the tangent space defined by ϕ; we call this the geometric consistency condition (GC). d Since Σ ∈ S++ and Dϕ has rank d, Λ has rank d with rng(Λ) = rng(Dϕ) = Tx M . Writing Λ = Q diag(λ1 , . . . , λD ) Q⊤ with λ1 ≥ · · · ≥ λd > 0 = λd+1 = · · · = λD , the top-d eigenvectors Ud := Q1:d ∈ RD×d satisfy Ud⊤ Ud = Id and P = Ud Ud⊤ is the orthogonal projection onto Tx M . Thus the tangent projector P (x) can be recovered from the ambient covariance Λ(x) by spectral truncation, connecting observable dynamics data to the tangent bundle of M . For (GC) to serve as a training penalty, it must not depend on the parameterization of Ω.
4
Lemma 2.4 (Coordinate invariance). Let α : Ω̃ → Ω be a diffeomorphism and ϕ̃ = ϕ ◦ α. (i) The orthogonal projector P = Dϕ g −1 Dϕ⊤ is invariant: P̃ = P . (ii) The Itô drift µ is not a geometric vector, but the difference of any two Itô drifts sharing a common local covariance Σ transforms as a contravariant vector: µ̃1 − µ̃2 = Dα−1 (µ1 − µ2 ). In particular, ∥µ1 − µ2 ∥g̃2 = ∥µ1 − µ2 ∥2g . (iii) For any symmetric (2,0)-tensor difference ∆Σ = Σ1 − Σ2 , the metric-weighted norm Tr (g ∆Σ)2 is invariant. Proof. (i) The chain rule gives D[ϕ̃] = Dϕ Dα, so g̃ = Dα⊤ g Dα and P̃ = D[ϕ̃] g̃ −1 D[ϕ̃]⊤ = Dϕ g −1 Dϕ⊤ = P . (ii) Under α, the Itô drift acquires a Hessian correction µ̃ = Dα−1 µ + 12 ⟨Σ , ∇2 α−1 ⟩F . Both µ1 and µ2 receive the same correction, so µ̃1 − µ̃2 = Dα−1 (µ1 − µ2 ). Then ∥µ̃1 − µ̃2 ∥g̃2 = (µ1 −µ2 )⊤ (Dα)−⊤ Dα⊤ g Dα Dα−1 (µ1 −µ2 ) = ∥µ1 − µ2 ∥2g . g = Dα−1 ∆Σ (Dα)−⊤ and the metric as g̃ = (iii) The covariance difference transforms as ∆Σ g = Dα⊤ g ∆Σ (Dα)−⊤ . Squaring and taking the trace, the similarity Dα⊤ g Dα. Multiplying: g̃ ∆Σ 2 g factors cancel: Tr (g̃ ∆Σ) = Tr Dα⊤ (g ∆Σ)2 (Dα)−⊤ = Tr (g ∆Σ)2 . Consequently, both (GC) and the tangent-bundle penalty derived from it in Section 3 are chartinvariant, and the metric-weighted losses for drift (Stage 2) and diffusion (Stage 3) in Algorithm 1 are coordinate-invariant.
3
The ρ-Metric and Geometric Penalties
In this section we introduce the function-space metric that underlies our geometric losses. Section 3.1 motivates a well-conditioned Jacobian assumption and defines the ρ-metric. Section 3.2 specifies the empirical ρ-loss, details the data requirements, and shows that the full pipeline is efficiently computable from dynamics observations.
3.1
The ρ-metric on E d (s)
Our metric will penalize the tangent-projector mismatch ∥Pϕ −Pψ ∥F . A natural question is whether this quantity is controlled by the Jacobian error ∥Dϕ − Dψ∥F . In general the answer is no: when the smallest singular value of Dϕ degenerates, an arbitrarily small Jacobian perturbation can rotate the tangent space by a large angle. This motivates restricting attention to charts whose Jacobians have a uniform lower singular-value bound. Let Rd (s) := {X ∈ RD×d : σmin (X) ≥ s} denote the set of full-rank matrices with controlled smallest singular value. We first record two identities about the orthogonal projection map f (X) = X(X ⊤ X)−1 X ⊤ that will be used throughout. Lemma 3.1 (Frobenius distance of orthogonal projections). Let P1 and P2 be D × D, rank-d orthogonal projections. Write Pi = Hi Hi⊤ with Hi⊤ Hi = Id , and let Ni be D × (D − d) with orthonormal columns spanning (rng Pi )⊥ . Then 2 ⊤ 2 ⊤ 2 ⊤ 2 1 2 ∥P1 − P2 ∥F = d − ∥H1 H2 ∥F = ∥N1 H2 ∥F = ∥N2 H1 ∥F .
Proof. Expanding gives ∥P1 − P2 ∥2F = 2(d − ⟨P1 , P2 ⟩F ) since ∥Pi ∥2F = d. The cyclic property of the trace yields ⟨P1 , P2 ⟩F = ∥H1⊤ H2 ∥2F . For the normal form, use I = Hi Hi⊤ + Ni Ni⊤ to decompose H2 and note ∥H2 ∥2F = d.
5
The lower singular-value bound is necessary: without it, the √ projection map is not√Lipschitz. Setting ϕ(t) = ϵ t e1 and ψ(t) = ϵ t e2 gives ∥Dϕ − Dψ∥F = 2 ϵ but ∥Pϕ − Pψ ∥F = 2, so the ratio blows up as ϵ → 0. On Rd (s), however, the projection map is Lipschitz: Lemma 3.2 (Lipschitz property of the projection map). Let f (X) = X(X ⊤ X)−1 X ⊤ be the orthogonal projection onto rng X. Then for all X, Y ∈ Rd (s), √ 2 ∥f (X) − f (Y )∥F ≤ ∥X − Y ∥F . s Proof. Write gY = Y ⊤ Y . Since f (X)X = X, we have (I − f (X))X = 0, so by Lemma 3.1 applied −1/2 −1/2 −1/2 with H2 = Y gY , 12 ∥f (X) − f (Y )∥2F = ∥(I − f (X)) Y gY ∥2F = ∥(I − f (X)) (Y − X) gY ∥2F . −1/2
Since I − f (X) is a contraction and ∥gY s−2 ∥X − Y ∥2F .
∥2 = σmin (Y )−1 ≤ s−1 , the right-hand side is at most
With this Lipschitz bound in hand, we define the chart class and the metric. Let E d (s) := {ϕ ∈ H 1 (Ω; RD ) : Dϕ(z) ∈ Rd (s) a.e. on Ω} denote the set of well-conditioned charts, i.e. those whose Jacobians have σmin ≥ s at almost every point. The requirement is mild: any C 1 immersion [28] ϕ⋆ : Ω̄ → RD satisfies s0 := minΩ̄ σmin (Dϕ⋆ ) > 0 by the extreme value theorem, so ϕ⋆ ∈ E d (s0 ) automatically. Every Monge patch ϕ(z) = (z ⊤ , f (z))⊤ satisfies σmin (Dϕ) ≥ 1 because Dϕ = [Id ; Df ]. A smaller s0 enlarges the Lipschitz constant in Lemma 3.2 and hence the generalization bounds. In practice the assumption is not explicitly enforced during training; we verify it post-hoc in Section 5. To use E d (s) as a hypothesis class for ERM, we first verify that it is closed under H 1 limits. Lemma 3.3. E d (s) is closed in H 1 . Proof. Let ϕk ∈ E d (s) with ϕk → ϕ in H 1 . By definition of the H 1 norm, ∥Dϕk − Dϕ∥L2 ≤ ∥ϕk − ϕ∥H 1 → 0, so Dϕk → Dϕ in L2 (Ω; RD×d ). By the Lp subsequence theorem, there exists a subsequence Dϕkj → Dϕ pointwise almost everywhere. Weyl’s singular-value perturbation inequality gives |σmin (A) − σmin (B)| ≤ ∥A − B∥F for any matrices A, B of the same size, so σmin is 1-Lipschitz on RD×d . Therefore σmin (Dϕkj (z)) → σmin (Dϕ(z)) for almost every z ∈ Ω. Since σmin (Dϕkj (z)) ≥ s for all j, the limit satisfies σmin (Dϕ(z)) ≥ s almost everywhere, hence ϕ ∈ E d (s). With the chart class in place, we combine the pointwise L2 reconstruction error and the tangentprojector discrepancy from Lemma 3.2 into a single metric. Lemma 3.4 (The ρ-metric). For ϕ, ψ ∈ E d (s), define ρ(ϕ, ψ)2 := ∥ϕ − ψ∥2L2 + 21 ∥Pϕ − Pψ ∥2F,L2 , where Pϕ = f (Dϕ) is the orthogonal projection onto rng(Dϕ) as in Lemma 3.2. Then the following hold: (i) ρ is a metric on E d (s). (ii) ∥ϕ − ψ∥L2 ≤ ρ(ϕ, ψ) ≤ Cs ∥ϕ − ψ∥H 1 with Cs = max(1, s−1 ). (iii) ρ is not equivalent to the H 1 norm. Proof. (i) Non-negativity and symmetry are immediate. If ρ(ϕ, ψ) = 0 then both terms vanish; in particular ∥ϕ − ψ∥L2 = 0, so ϕ = ψ a.e. on Ω and hence Dϕ = Dψ a.e., giving ϕ = ψ in H 1 ∩ E d (s). For the triangle inequality, let ϕ, ψ, χ ∈ E d (s) and write a1 = ∥ϕ − ψ∥L2 , b1 = √12 ∥Pϕ − Pψ ∥F,L2 ,
6
and define (a2 , b2 ), (a3 , b3 ) analogously for (ψ, χ) and (ϕ, χ). The L2 triangle inequality gives a3 ≤ a1 + a2 and b3 ≤ b1 + b2 , so q p ρ(ϕ, χ) = a23 + b23 ≤ (a1 +a2 )2 + (b1 +b2 )2 ≤ ρ(ϕ, ψ) + ρ(ψ, χ), where the last inequality is the triangle inequality for the Euclidean norm on R2 . (ii) The lower bound ∥·∥L2 √≤ ρ is immediate from the definition. For the upper bound, Lemma 3.2 gives ∥Pϕ (x) − Pψ (x)∥F ≤ s2 ∥Dϕ(x) − Dψ(x)∥F pointwise a.e. Squaring, integrating over Ω, and combining with the L2 term yields ρ(ϕ, ψ)2 ≤ ∥ϕ − ψ∥2L2 +
1 ∥Dϕ − Dψ∥2L2 ≤ Cs2 ∥ϕ − ψ∥2H 1 . s2
Pd a (iii) Let z = (z1 , . . . , zd ) ∈ Ω and define ϕ(z) = j=1 zj ej and ϕk (z) = ϕ(z) − k cos(kz1 ) e1 for 0 < a < 1 − s. Then ∂1 ϕk (z) = (1 + a sin(kz1 )) e1 and ∂j ϕk (z) = ej for j = 2, . . . , d, so σmin (Dϕk ) ≥ 1 − a > s and ∥Dϕk ∥F is uniformly bounded. Hence ϕk ∈ E d (s). Both Pϕk and Pϕ equal the constant orthogonal projection onto span{e1 , . . . , ed }, so ∥Pϕk − Pϕ ∥F,L2 = 0. Moreover ∥ϕk − ϕ∥L2 → 0 as k → ∞ (the oscillation has amplitude a/k). Hence ρ(ϕk , ϕ) → 0. However, Dϕk − Dϕ has first column a sin(kz1 ) e1 and all other columns zero, so Z Z a2 a2 a2 2 2 2 sin (kz1 ) dz = |Ω| − ∥Dϕk − Dϕ∥L2 (Ω) = a cos(2kz1 ) dz → |Ω| > 0, 2 2 Ω 2 Ω by the Riemann–Lebesgue lemma. Thus ϕk ̸→ ϕ in H 1 , so ρ is strictly weaker than H 1 . In summary, ρ sits strictly between L2 and H 1 on E d (s): it controls tangent-space alignment without requiring full derivative matching. This is precisely the gap exploited by the tangent-bundle penalty.
3.2
Efficient computation from dynamics data
The pointwise ρ-loss at a latent coordinate z is ℓρ (z; θ) := ∥ϕθ (z) − ϕ⋆ (z)∥2 + 21 ∥Pϕθ (z) − Pϕ⋆ (z)∥2F , where ϕ⋆ is the target chart and ϕθ is the decoder parametrized by θ. For a probability measure R ν on Ω, the population risk is Rν,ρ (θ) := Ω ℓρ (z; θ) dν(z), and the empirical risk minimizer is i.i.d. 1 Pm θ̂S,ρ ∈ arg minθ∈Θ m i=1 ℓρ (zi ; θ) with zi ∼ ν. Since ℓρ is coordinate-invariant (Lemma 2.4), Rν,ρ is reparametrization-invariant whenever ν is√the pullback of a fixed measure on M , in particular for the normalized Riemannian volume dν⋆ ∝ det g⋆ dz. For the generalization analysis in Section 4, we take ν to be Lebesgue on [0, 1]d ; the two are equivalent up to bounded constants since g⋆ has bounded singular values on compact Ω. The target chart ϕ⋆ is unknown in practice. The training data {xi , b(xi ), Λ(xi )}m i=1 from the ATLAS exploration phase provide two ingredients for realizing the ρ-loss; the drift b(xi ) enters in Stage 2 (below). Writing zi = πθ (xi ) and P̂i = Ud Ud⊤ for the rank-d spectral projector of Λ(xi ) (Section 2.2): (i)
• Reconstruction LR = ∥ϕθ (zi ) − xi ∥2 : a surrogate for the ρ-loss term ∥ϕθ (z) − ϕ⋆ (z)∥2 , since xi = ϕ⋆ (zi⋆ ) for the true latent zi⋆ = π⋆ (xi ). The two coincide when πθ = π⋆ . 7
Algorithm 1 Three-stage pipeline: geometric chart learning and latent SDE estimation D with drift b(x ) and covariance Λ(x ) {e.g. from Require: Ambient observations {xi }m i i i=1 ⊂ R short-burst ensembles [12, 48]} Require: Penalty weights λT , λF ; epochs E Ensure: Chart (πθ , ϕθ ); latent SDE coefficients (µ̂ω , σ̂ψ ) 1: Stage 1: Chart learning {Geometric autoencoder} 2: Train (πθ , ϕθ ) in two phases (warmup E1 , fine-tune E2 epochs), both minimizing L = LR + λT LT + λF LF via (3) 3: Stage 2: Encoder-pullback drift fitting {Frozen chart} 4: Encode zi ← πθ (xi ); compute encoder-pullback target µ(zi ) = Dπθ (xi ) b(xi ) + 1 2 π (x )⟩ via (2) ⟨Λ(x ) , ∇ i θ i F P 2 2 5: Train µ̂ω to minimize i ∥µ̂ω (zi ) − µ(zi )∥g(zi ) 6: Stage 3: Latent diffusion fitting {Frozen chart} 7: Compute target Σ(zi ) = Dπθ (xi ) Λ(xi ) Dπθ (xi )⊤ P ⊤ 2 8: Train σ̂ψ to minimize i Tr (g(zi ) [σ̂ψ (zi )σ̂ψ (zi ) − Σ(zi )])
(i)
• Tangent-bundle LT = 12 ∥Pϕθ (zi )− P̂i ∥2F : replaces the unknown Pϕ⋆ with the data-derived P̂i . (i)
We augment the ρ-loss with an inverse-consistency penalty [17], LF = ∥Dπθ (xi ) Dϕθ (zi ) − Id ∥2F , which enforces the first-order inverse identity Dπ Dϕ = Id , ensuring that the encoder and decoder are consistent as approximate chart maps. The condition Dπ Dϕ = Id is coordinate-invariant (it holds in all parameterizations or none), though the Frobenius penalty itself is not; LF is a practical regularizer, not part of the ρ-metric or the generalization theory. The full training objective for Stage 1 (chart learning) is m
1 X (i) (i) (i) L(θ) = LR + λT LT + λF LF . m
(3)
i=1
Stages 2 and 3 freeze the chart and fit latent SDE coefficients using coordinate-invariant losses (Lemma 2.4 (ii)–(iii)). Stage 2 regresses a drift network µ̂ω (z) onto the encoder-pullback target (2) under ∥·∥2g ; Stage 3 regresses a diffusion network σ̂ψ (z) onto the pulled-back covariance Σ(zi ) := Dπθ (xi ) Λ(xi ) Dπθ (xi )⊤ under Tr (g · )2 . The geometric losses are coordinate-invariant: LR is computed in ambient RD ; LT uses the ambient projector (Lemma 2.4 (i)); Stage 2 and 3 use metric-weighted norms (Lemma 2.4 (ii)–(iii)). The law of the learned ambient process is therefore independent of the latent parameterization. The auxiliary penalty LF enforces a coordinate-invariant condition (Dπ Dϕ = Id ) but its Frobenius value is chart-dependent. The full pipeline is summarized in Algorithm 1. Naı̈vely, the tangent-bundle penalty forms D × D projection matrices (O(D2 d) time, O(D2 ) memory). The following proposition shows that all stages can be evaluated with O(Dd2 ) flops and O(Dd) memory. Proposition 3.5 (Computational cost). Every loss and target in the three-stage pipeline can be evaluated in O(Dd2 ) flops per sample, with O(Dd) for the drift loss, without forming any D × D matrix. Specifically, with g = Dϕ⊤ Dϕ ∈ Rd×d : (R) reconstruction is O(D); (T) tangent-bundle is O(Dd2 ); (F) inverse-consistency is O(Dd2 ); (Stage 2) the encoder-pullback target is O(Dd2 ) via O(d2 ) Hessian–vector products at O(D) each; the drift loss itself is O(Dd) given the precomputed Jacobian Dϕ ∈ RD×d ; (Stage 3) pulled-back covariance target and diffusion loss are each O(Dd2 ). 8
Proof. Write Ud ∈ RD×d for the top-d eigenvectors of Λ, P = Ud Ud⊤ , P̂ = Dϕ g −1 Dϕ⊤ , and C = Dϕ⊤ Ud ∈ Rd×d . (R) LR = ∥ϕθ (z) − x∥2 is a single D-vector difference and dot product: O(D). (T) Since P̂ and P are both orthogonal projections of rank d, idempotency gives P̂ 2 = P̂ , Expanding: ∥P̂ − P ∥2F = Tr(P̂ 2 ) − 2 Tr(P̂ P ) + Tr(P 2 ) = 2d − 2 Tr(P̂ P ). Substituting P̂ = Dϕ g −1 Dϕ⊤ and P = Ud Ud⊤ and applying the cyclic property of the trace: Tr(P̂ P ) = Tr Dϕ g −1 Dϕ⊤ Ud Ud⊤ = Tr g −1 (Dϕ⊤ Ud )(Dϕ⊤ Ud )⊤ = Tr(g −1 CC ⊤ ), so
P 2 = P , and Tr(P̂ ) = Tr(P ) = d.
2 −1 ⊤ 1 . 2 ∥P̂ − P ∥F = d − Tr g C C
(4)
Forming C = Dϕ⊤ Ud costs O(Dd2 ); g −1 CC ⊤ is d × d, so the trace is O(d3 ) ⊂ O(Dd2 ). When the inverse-consistency penalty drives Dπ Dϕ ≈ Id , we have Dπ ≈ g −1 Dϕ⊤ and the loss simplifies to d − Tr(Dπ Ud Ud⊤ Dϕ), which is the form used in Algorithm 1. (F) LF = ∥Dπ Dϕ − Id ∥2F : Dπ ∈ Rd×D times Dϕ ∈ RD×d is a d × d product costing O(Dd2 ). (Stage 2) Target. The encoder-pullback drift (2) contains the ⟨Λ , ∇2 π j ⟩F . Since Pd Itô correction 2 j Λ has rank d on the manifold, its spectral decomposition Λ = m=1 λm um u⊤ m gives ⟨Λ , ∇ π ⟩F = Pd ⊤ 2 j ⊤ 2 j m=1 λm um (∇ π ) um , where each um (∇ π ) um is a single forward-over-reverse Hessian–vector product (HVP), avoiding the full d × D × D encoder Hessian. Over j = 1, . . . , d encoder outputs this is d2 HVPs at O(D) each, totalling O(Dd2 ). Loss. ∥r∥2g = r⊤ g r = ∥Dϕ r∥2 rewrites the metric-weighted norm as an ambient Euclidean norm. Computing Dϕ r is a matrix–vector product in O(Dd); the squared norm costs O(D). (Stage 3) The pulled-back covariance target Σ(zi ) = Dπ Λ Dπ ⊤ is naively O(dD2 ) because Λ ∈ RD×D . Since rnk(Λ) = d, write Λ = Ud diag(λ1 , . . . , λd ) Ud⊤ and set B = Dπ Ud ∈ Rd×d at cost O(Dd2 ); then Σ = B diag(λ) B ⊤ in O(d3 ). For the loss, g = Dϕ⊤ Dϕ is precomputed in O(Dd2 ); the product g ∆Σ is a d × d matrix multiplication costing O(d3 ). Squaring (g ∆Σ)2 and taking the trace are likewise O(d3 ) ⊂ O(Dd2 ).
4
Generalization Theory
The previous section constructed a pipeline whose geometric losses (LR , LT , Stage 2, Stage 3) are coordinate-invariant and depend only on point-cloud positions, local covariance eigenvectors, ambient drift, and the learned chart; no ground-truth Jacobians or curvature labels are required. Yang et al. [46] establish optimal generalization rates for deep networks trained in the full Sobolev H k norm, where the training loss requires labels for all derivatives up to order k. Our ρ-loss is strictly less supervised than H 1 training, which would require full Jacobian labels Dϕ⋆ (xi ). A natural question is whether this weaker supervision comes at a statistical cost. Perhaps surprisingly, the answer is that no rate is lost. We show that ρ-ERM achieves the same rate (up to logarithmic factors) as full H 1 training. The key new ingredient is a bridge lemma (Lemma 4.4) that converts a W 1,∞ approximation guarantee, available from [46], into a ρ-approximation guarantee. Once this bridge is in place, the estimation error is controlled by a covering-number argument over the first-order feature class, which reduces to the entropy bounds already established in [46]. This improves on the chart autoencoder rate of Liu et al. [33], which achieves m−2/(d+2) log4 m from reconstruction loss alone: geometric supervision via the tangent projector replaces the fixed smoothness n = 2 in their bound with the true regularity n of the target chart. Section 4.2 then analyzes the systematic bias that arises when the decoder-side drift formula is used in place of the encoder-pullback target. We begin with three standing assumptions.
9
Assumption 4.1 (Target chart regularity). Let Ω ⊂ Rd be a bounded Lipschitz domain and let ϕ⋆ : Ω → RD be the ground-truth chart. Assume ϕ⋆ ∈ E d (s0 ) ∩ W n,∞ (Ω; RD ) for some s0 > 0 and n > 1, and after rescaling, ∥ϕ⋆ ∥W n,∞ (Ω) ≤ 1. The well-conditioned chart class E d (s0 ) from Section 3.1 ensures that σmin (Dϕ⋆ ) ≥ s0 a.e. on Ω, so the ρ-metric is well-defined on a neighborhood of ϕ⋆ . Assumption 4.2 (Hypothesis class). Fix constants B ≥ 1 and s ∈ (0, s0 ]. For each architecture (N, L) (width N , depth L), let H(N, L) = {ϕθ : Ω → RD : θ ∈ ΘN,L } be a DeNN hypothesis set. The constrained class used for ERM is n o HB,s (N, L) := ϕθ ∈ H(N, L) : ∥ϕθ ∥W 1,∞ (Ω) ≤ B and ϕθ ∈ E d (s) . Since HB,s (N, L) ⊂ E d (s), the ρ-metric and the projector Pϕθ are well-defined for every member. Assumption 4.3 (W 1,∞ approximation). Assume the DeNN architecture is chosen as in [46] (width N = O(log L), depth L, parameter count W = O(L(log L)3 )). By the W 1,∞ approximation result (Proposition A.2 and Theorem 3.2 of [46]), there exists a constant Capp > 0 depending only on (d, D, n) and fixed architectural constants such that, for each (N, L), one can find θ♯ ∈ ΘN,L with −2(n−1)/d ϕθ♯ ∈ HB,s (N, L) and ∥ϕθ♯ − ϕ⋆ ∥W 1,∞ (Ω) ≤ Capp (logWW )2 . Once the right-hand side is at most s0 /2, Weyl’s perturbation bound [21] gives ϕθ♯ ∈ E d (s0 /2), so membership in the constrained class is automatic. The singular-value bound defining E d (s) depends on the choice of coordinates on Ω; the latent domain with its Euclidean structure is fixed as part of the statistical setup. As noted in Section 3, the ρ-loss and Stage 2/3 losses are coordinate-invariant (Lemma 2.4), and under bounded-distortion reparametrisations the hypothesis class is stable: ϕ ∈ E d (s) ⇒ ϕ ◦ ψ ∈ E d (as) when σmin (Dψ) ≥ a.
4.1
Main generalization results
The first ingredient is a bridge lemma showing that W 1,∞ -closeness to a well-conditioned target implies ρ-closeness. Lemma 4.4. Let Ω ⊂ Rd be measurable with 0 < |Ω| < ∞ and let ϕ⋆ : Ω → RD be a C 1 chart such that for all x ∈ Ω, the Jacobian J⋆ (x) := Dϕ⋆ (x) ∈ RD×d has full rank d and s0 ≤ σmin (J⋆ (x)) ≤ σmax (J⋆ (x)) ≤ r0 for some constants 0 < s0 ≤ r0 < ∞. For any C 1 map ϕ : Ω → RD with Dϕ(x) −1 full rank, define Pϕ (x) := Dϕ(x) Dϕ(x)⊤ Dϕ(x) Dϕ(x)⊤ and ρ(ϕ, ϕ⋆ )2 := ∥ϕ − ϕ⋆ ∥2L2 (Ω) + 1 2 2 ∥Pϕ − Pϕ⋆ ∥F,L2 (Ω) . Assume there exists a neural-network realizable ϕ = ϕN,L satisfying ∥ϕ − ϕ⋆ ∥W 1,∞ (Ω) ≤ δ,
(5)
and suppose δ ≤ s0 /2. Then (i) For all x ∈ Ω, s20 ≤ σmin (Dϕ(x)) ≤ σmax (Dϕ(x)) ≤ r0 + δ, so Dϕ(x) is full rank and Pϕ is well-defined. (ii) There is a constant Cρ (s0 , Ω) := |Ω| 1 + s42 such 0 that ρ(ϕ, ϕ⋆ )2 ≤ Cρ (s0 , Ω) ∥ϕ − ϕ⋆ ∥2W 1,∞ (Ω) ≤ Cρ (s0 , Ω) δ 2 . (6) Consequently, if an approximation theorem yields δ ≲ N −2(n−1)/d L−2(n−1)/d for ϕ⋆ ∈ W n,∞ , then ρ(ϕ, ϕ⋆ )2 ≲ N −4(n−1)/d L−4(n−1)/d . In the DeNN regime where N = O(log L) and the parameter count satisfies W = O(N 2 L log L), the above can be rewritten (up to logarithmic factors) − 4(n−1) d as ρ(ϕ, ϕ⋆ )2 ≲ (logWW )2 . 10
Proof. Write J(x) := Dϕ(x) and E(x) := J(x) − J⋆ (x). (i)) For each fixed x ∈ Ω, Weyl’s perturbation bound for singular values gives |σi (J(x)) − σi (J⋆ (x))| ≤ ∥E(x)∥2 , for each singular value index i. Since ∥E(x)∥2 ≤ ∥E(x)∥F ≤ ∥E∥F,∞ ≤ δ s0 by (5), we obtain σmin (J(x)) ≥ σmin (J⋆ (x)) − δ ≥ s0 − δ ≥ 2 , and similarly σmax (J(x)) ≤ σmax (J⋆ (x)) + δ ≤ r0 + δ. This proves (i). (ii)) Let f (X) := X(X ⊤ X)−1 X ⊤ be the projector map. By the Lipschitz√bound for the projector on √ 2 2 2 matrices with σmin (·) ≥ s0 /2 (Lemma 3.2), we have ∥f (X)−f (Y )∥F ≤ s0 /2 ∥X−Y ∥F = s0 ∥X− Y ∥F . Applying this with X =√J(x) and Y = J⋆ (x) yields,√ for all x ∈ Ω, ∥Pϕ (x) − Pϕ⋆ (x)∥F = ∥f (J(x)) − f (J⋆ (x))∥F ≤ 2s02 ∥J(x) − J⋆ (x)∥F ≤ 2s02 δ. Therefore ∥Pϕ − Pϕ⋆ ∥2F,L2 (Ω) = R 8 2 2 Ω ∥Pϕ (x) − Pϕ⋆ (x)∥F dx ≤ |Ω| · s2 δ . 0
Also by (5), ∥ϕ − ϕ⋆ ∥2L2 (Ω) ≤ |Ω| · ∥ϕ − ϕ⋆ ∥2L∞ (Ω) ≤ |Ω| δ 2 . Combining the two bounds, ρ(ϕ, ϕ⋆ )2 = ∥ϕ − ϕ⋆ ∥2L2 (Ω) + 12 ∥Pϕ − Pϕ⋆ ∥2F,L2 (Ω) ≤ |Ω|δ 2 + 21 |Ω| s82 δ 2 = |Ω| 1 + s42 δ 2 , which is (6) 0 0 and completes the proof. With this bridge in hand, we state the main result for an oracle ρ-ERM: ℓρ (z; θ) evaluates ϕθ against ϕ⋆ at latent points z ∼ ν known to the oracle, a standard device in covering-number oracle inequalities; quasi-uniform ATLAS designs [18] achieve comparable rates. The implemented loss differs in two ways. First, the learned encoder replaces the true latent coordinate. On the wellconditioned class, both the reconstruction and projector terms in ℓρ are Lipschitz in the evaluation point, so the practical ρ-risk exceeds the oracle risk by a term controlled by E∥πθ (ϕ⋆ (z)) − z∥2 ; when Dπ Dϕ ≈ Id (enforced by LF ), this gap is small because the encoder is an approximate local inverse. Second, if P̂i is formed from an empirical covariance Λ̂i , the Davis–Kahan sin Θ theorem yields ∥P̂i − Pϕ⋆ (zi )∥F ≤ C∥Λ̂i − Λi ∥op /(λd − λd+1 ), where λd −λd+1 is the eigengap separating tangent from normal eigenvalues; thus the projector error vanishes when the covariance estimator is accurate and the eigengap stays bounded away from zero. Theorem 4.5 (Generalization in ρ has the same order as H 1 ). Under Assumptions 4.1–4.3, let θ̂ be the ρ-ERM over HB,s0 /2 (N, L) with {zi }m i=1 drawn i.i.d. from a probability measure ν on Ω. −2(n−1)/d Assume N, L are large enough that Capp N L−2(n−1)/d ≤ s0 /2. Then there exists a constant C > 0 depending only on (d, D, n, B, s0 ) such that " # 2 L2 log L log N 4(n−1) 4(n−1) N 2 2 E Rν,ρ (θ̂) ≤ C N − d L− d + log m , (7) m i.i.d.
where the expectation is over (zi )m i=1 ∼ ν. In particular, in the DeNN regime N = O(log L) and parameter count W = O(N 2 L log L), the bound (7) can be rewritten (up to logarithmic factors) as " # − 4(n−1) W2 d W E Rν,ρ (θ̂) ≤ C (log W )2 + log m , (8) m and optimizing over W yields the same sample-only rate (up to logs) as H 1 -Sobolev training: 2(n−1)
− 2(n−1)+d
E Rν,ρ (θ̂) ≲ m
.
The ρ-metric discards tangential stretching information and is therefore strictly weaker than H 1 (Lemma 3.4 (iii)). That it achieves the same rate reflects two complementary facts. On the approximation side, the well-conditioning assumption σmin (Dϕ) ≥ s makes the map Dϕ 7→ Pϕ 11
Lipschitz (Lemma 3.2), so any W 1,∞ approximant at scale δ automatically gives ρ-error at scale δ 2 with no loss in exponent. On the estimation side, the ρ-loss is a Lipschitz function of the same firstorder features (ϕ, Dϕ) that determine the H 1 covering number, so the entropy bound is inherited unchanged. In short, ρ uses less supervision than H 1 because the covariance data supplies the tangent-space information that would otherwise require Jacobian labels, yet on the well-conditioned class both exponents remain first-order. ⋆ Proof. Let θ̂ be the ρ-ERM and write Rν,ρ := inf{Rν,ρ (θ) : ϕθ ∈ HB,s0 /2 (N, L)}. Let Lρ := {ℓρ (·; θ) : ϕθ ∈ HB,s0 /2 (N, L)}. Since the squared loss is uniformly bounded, 0 ≤ ℓρ ≤ Bρ2 by eq. (11) below, we may apply a fast-rate covering-number oracle inequality for ERM with bounded square loss [18, 4, Theorem 11.3]. Concretely, there exist constants C, c0 > 0 such that 1 ⋆ E Rν,ρ (θ̂) ≤ Rν,ρ + C log N c0 /m, Lρ , m + 1 . (9) m ⋆ ≤ R (θ̄) for any fixed comparator θ̄ in the class, it remains to (i) choose θ̄ to control Since Rν,ρ ν,ρ the approximation error and (ii) bound log N (·, Lρ , m). By Assumption 4.3, there exists θ♯ such that ϕθ♯ ∈ HB,s0 /2 (N, L) and ∥ϕθ♯ − ϕ⋆ ∥W 1,∞ (Ω) ≤ 2(n−1)
2(n−1)
Capp N − d L− d . Denote this bound by δ; by hypothesis δ ≤ s0 /2. The proof of Lemma 4.4 gives the pointwise bound ℓρ (z; θ♯ ) ≤ (1 + 4/s20 ) δ 2 ; integrating against the probability measure ν yields 4(n−1) 4(n−1) Rν,ρ (θ♯ ) ≤ (1 + 4/s20 ) δ 2 ≤ C N − d L− d . Let θ̄ := θ♯ . Then Rν,ρ (θ̄) ≤ C N −
4(n−1) d
L−
4(n−1) d
.
(10)
First note ℓρ (z; θ) is uniformly bounded: since ∥ϕθ ∥L∞ ≤ B and ∥ϕ⋆ ∥L∞√≤ 1, we have ∥ϕθ − ϕ⋆ ∥ ≤ B+1, and since Pϕθ , Pϕ⋆ are orthogonal projections, ∥Pϕθ − Pϕ⋆ ∥F ≤ 2 d, hence 0 ≤ ℓρ (z; θ) ≤ Bρ2
for all z ∈ Ω, θ,
(11)
for a constant Bρ = Bρ (d, B). Next, by the Lipschitz bound for the projector map on the set {J : σmin (J) ≥ s0 /2}, ∥Pϕθ (z) − Pϕθ′ (z)∥F ≤ sC0 ∥Dϕθ (z) − Dϕθ′ (z)∥F for a.e. z. Using the elementary inequality ∥a∥2 − ∥b∥2 ≤ (∥a∥ + ∥b∥) ∥a − b∥ and the uniform bounds in W 1,∞ , there exists C = C(d, D, B, s0 ) such that for a.e. z and all θ, θ′ , |ℓρ (z; θ) − ℓρ (z; θ′ )| ≤ C ∥ϕθ (z) − ϕθ′ (z)∥ + ∥Dϕθ (z) − Dϕθ′ (z)∥F . (12) n Define the first-order feature map class F := z 7→ (ϕθ (z), ∂1 ϕθ (z), . . . , ∂d ϕθ (z)) : ϕθ ∈ o HB,s0 /2 (N, L) . Then (12) implies the uniform covering numbers satisfy N ε, Lρ , m ≤ N c ε, F, m for some constant c = c(d, D, B, s0 ) > 0, where N (·, ·, m) denotes the uniform covering number on m sample points. Bounding N (ε, F, m) is exactly the “first-derivative” (k=1) instance of the Sobolev-loss generalization machinery in [46, Theorem 3.2]. That result bounds the covering number of the full DeNN derivative class; since HB,s0 /2 (N, L) is a subclass, the same entropy bound applies a fortiori. Moreover, [46, Remark 3.1] explains that the same entropy bound holds for general i.i.d. sampling measures. Specializing their result to the first-order feature class F yields, for ε of order 1/m, log N (ε, Lρ , m) ≤ C N 2 L2 log2 L log2 N · log m, 12
(13)
where C = C(d, D, B, s0 ). Inserting (13) into (9) gives ⋆ E Rν,ρ (θ̂) ≤ Rν,ρ + C
N 2 L2 log2 L log2 N log m. m
(14)
⋆ ≤ R (θ̄), combining (10) with (14) yields (7). In the DeNN regime N = O(log L) Since Rν,ρ ν,ρ and W = O(N 2 L log L), (8) follows up to logarithmic factors. Balancing the two terms in (8) gives 2(n−1)
− 2(n−1)+d
the sample-only rate m
4.2
(up to logs), matching the Sobolev H 1 training rate.
Bias analysis of decoder-side drift fitting
Theorem 4.5 guarantees that the oracle ρ-ERM chart converges at the optimal rate. The next question is how to extract latent SDE coefficients from this chart. The latent drift admits two equivalent representations as derived in Section 2.1: the encoder-pullback form (1) using (Dπ, ∇2 π), and the decoder form of Lemma 2.3 using the metric pseudo-inverse g −1 Dϕ⊤ and the decoder Hessian. For an ideal autoencoder (π = ϕ−1 ) these coincide; for a learned autoencoder they diverge, and the decoder-side target carries a systematic bias. Proposition 4.6 (Decoder-side bias decomposition). Let (ϕθ , πθ ) be a C 2 encoder–decoder pair −1 for the decoder-side latent covariance, with Dϕθ of full column rank. Write Σ̂ = g −1 Dϕ⊤ θ Λ Dϕθ g 1 q(Σ̂) for the corresponding Itô correction, µdec := g −1 Dϕ⊤ θ [b − 2 q(Σ̂)] for the decoder-side drift target, and µenc for the encoder-pullback target (2). Then µdec − µenc = (g −1 Dϕ⊤ − Dπ)(b − 21 q) − 21 ⟨Σ̂ , D2 (π◦ϕ)⟩F − 21 ⟨Λ − DϕΣ̂Dϕ⊤ , ∇2π⟩F , | {z } | {z } | {z } (II)
(I)
(15)
(III)
where (I) is the pseudo-inverse ̸= encoder gap, (II) the cycle-Hessian bias, and (III) the covariance mismatch. Proof. Write ϕ, π, g, Σ̂, q for ϕθ , πθ , gθ , Σ̂θ , q(Σ̂), suppressing θ. Subtracting µenc from µdec and splitting the b-dependent and Itô-correction terms gives µdec − µenc = (g −1 Dϕ⊤ − Dπ) b − 12 g −1 Dϕ⊤ q + ⟨Λ , ∇2 π⟩F . For the bracket, write g −1 Dϕ⊤ q = Dπ q + (g −1 Dϕ⊤ −Dπ) q and evaluate Dπ q via the secondderivative chain rule applied to h := π ◦ ϕ: D2 hjab = (∇2 π j )kl (Dϕ)ka (Dϕ)lb + (Dπ)jk (∇2 ϕk )ab . Contracting with Σ̂ab and stacking over j: ⟨Σ̂ , D2 (π ◦ ϕ)⟩F = ⟨Dϕ Σ̂ Dϕ⊤ , ∇2 π⟩F + Dπ q. Substituting back: g −1 Dϕ⊤ q + ⟨Λ , ∇2 π⟩F = (g −1 Dϕ⊤ −Dπ) q + ⟨Σ̂ , D2 (π ◦ ϕ)⟩F + ⟨Λ − Dϕ Σ̂ Dϕ⊤ , ∇2 π⟩F . Collecting the (g −1 Dϕ⊤ −Dπ) terms yields the three-term decomposition. Interpretation. All three terms are biases in the decoder-side formula for a fixed learned chart. Terms (I) and (III) are first-order: LF reduces Term (I) by enforcing Dπ ≈ g −1 Dϕ⊤ on rng(Dϕ). By the Itô decomposition, b − 12 q(Σ) = Dϕ µ is exactly tangential; when (GC) also holds, Σ̂ ≈ Σ so q(Σ̂) ≈ q(Σ) and b − 12 q(Σ̂) is approximately tangential, making LF effective on the relevant subspace. The diagnostic E = ∥(I − Pθ )(b − 21 q)∥2 in Section 5 measures this normal-space leakage, 13
while (GC) reduces Term (III) by aligning rng(Λ) with rng(Dϕ). Term (II), the cycle-Hessian bias governed by D2 (π◦ϕ), is genuinely second-order: it requires π◦ϕ to be locally affine, which the firstorder penalties T and F do not enforce. Empirically, ∥D2 (π ◦ ϕ)∥F ≈ 0.9–1.3 across trained T+F models, so this term need not be small. Directly penalizing ∥D2 (π ◦ ϕ)∥2F during training degrades reconstruction and tangent alignment, because the penalty competes with the chart-fitting objective when the autoencoder is far from an exact inverse. By contrast, the encoder-pullback target (2) is exact by Itô’s formula applied to the learned encoder, bypasses the decoder-side decomposition and the amplification of Term (II), and in Section 5 consistently improves trajectory error across all N ×D configurations.
4.3
From coefficient convergence to weak convergence
The convergence chain has three links, each proved conditionally on the preceding one. 1. Chart-level generalization (Theorem 4.5): the ρ-ERM controls reconstruction and tangentbundle error. 2. Coefficient propagation (Theorems 4.9 and 4.10 below): assuming ϕn → ϕ in W 2,∞ with σmin (Dϕn ) ≥ s > 0, the induced local coefficients (Σ∗n , µ∗n ) converge uniformly on compacta, and the reconstructed ambient coefficients recover (b, Λ). 3. Weak convergence (Theorem 4.12 below): given uniform coefficient convergence, Stroock– Varadhan yields weak convergence of the learned dynamics. The gap in this chain is between links (1) and (2): the ρ-metric is weaker than H 1 and does not control second derivatives, so W 2,∞ convergence of the training sequence is assumed as a sufficient condition, not derived from the training objective. Uniform Hessian boundedness is automatic for networks with bounded depth and weight norms (the hypothesis class of Assumption 4.3), but convergence of Hessians requires additional structure not imposed by the current penalties. We first state a technical lemma on Frobenius inner products used in both directions. Lemma 4.7. Let K ⊂ Rd be compact, ϕ ∈ C 2 (K; RD ), and ϕn ∈ C 2 (K; RD ) with d ) and A : Rd → S d supn supK ∥∇2 ϕin ∥F ≤ CK for each i. Let A ∈ C 1 (Rd , S++ n ++ satisfy 2 i (a) An → A uniformly on K and (b) ⟨A , ∇ ϕn ⟩F → ⟨A , ∇2 ϕi ⟩F uniformly on K. Then ⟨An , ∇2 ϕin ⟩F → ⟨A , ∇2 ϕi ⟩F uniformly on K. Proof. By bilinearity and an add-zero trick, ⟨An , ∇2 ϕin ⟩F − ⟨A , ∇2 ϕi ⟩F = ⟨An − A , ∇2 ϕin ⟩F + ⟨A , ∇2 ϕin − ∇2 ϕi ⟩F . The second term vanishes uniformly by (b). By Cauchy–Schwarz the first is bounded by CK supK ∥An − A∥F → 0 by (a). Ambient-to-local convergence. The next lemma controls the local covariance in terms of the decoder and ambient covariance. Lemma 4.8. Let K ⊂ Ω be compact. Let Σ, Σ̂ denote local-coordinate d × d covariance fields, and let Λ, Λ̂ denote the corresponding ambient D × D covariance fields related via Σ = (Dϕ)† Λ((Dϕ)† )⊤ and Σ̂ = (Dϕ̂)† Λ̂((Dϕ̂)† )⊤ . Suppose ϕ, ϕ̂ satisfy σmin (Dϕ), σmin (Dϕ̂) ≥ s > 0 and σmax (Dϕ), σmax (Dϕ̂) ≤ R on√ K, with ∥Σ∥F , ∥Σ̂∥F ≤ c̃0 . Then ∥Σ − Σ̂∥F ≤ C0 ∥Dϕ − Dϕ̂∥F + √ C1 ∥Λ − Λ̂∥F , where C0 = 2 c̃0s 2 [s−2 + 2 R2 s−4 ] and C1 = d/s2 . Proof. Using Σ = (Dϕ)† Λ((Dϕ)† )⊤ , expand Σ − Σ̂ into a Λ-difference term and a pseudo-inversedifference term. The is locally Lipschitz on √ well-conditioned matrices: ∥(Dϕ)† − √ pseudo-inverse † −2 2 −4 † (Dϕ̂) ∥F ≤ [s + 2R s ]∥Dϕ − Dϕ̂∥F , and ∥(Dϕ) ∥ ≤ d/s. Substituting yields the stated constants.
14
For the drift, µ − µ̂ = Dϕ† b − Dϕ̂† b̂ + 12 (Dϕ̂† q̂ − Dϕ† q) gives, by the same well-conditioning, ∥µ − µ̂∥ ≤ C̃0 ∥Dϕ − Dϕ̂∥F + C̃1 ∥b − b̂∥2 + C̃2 ∥q̂ − q∥2 .
(16)
Theorem 4.9. Let K ⊂ Ω be compact, ϕ ∈ C 3 (Ω; RD ) with σmin (Dϕ) ≥ s > 0 on K, b ∈ D ). Denote by (µ, Σ) the local coefficients induced by the chart ϕ and the C 1 (Ω; RD ), Λ ∈ C 1 (Ω; S++ ambient data (b, Λ) via µ = (Dϕ)† (b− 21 q) and Σ = (Dϕ)† Λ((Dϕ)† )⊤ , and analogously (µ̂n , Σ̂n ) from (ϕ̂n , b̂n , Λ̂n ). Let ϕ̂n satisfy σmin (Dϕ̂n ) ≥ s > 0 on K with supn supK ∥∇2 ϕ̂in ∥F ≤ CK , and assume d ), ϕ̂n → ϕ in W 1,∞ (K), b̂n → b, Λ̂n → Λ uniformly on K, and for all bounded A ∈ Cb1 (Ω, S++ 2 i 2 i ⟨A , ∇ ϕ̂n ⟩F → ⟨A , ∇ ϕ ⟩F uniformly on K. Then µ̂n → µ and Σ̂n → Σ uniformly on K. Proof. By Lemma 4.8, Dϕ̂n → Dϕ and Λ̂n → Λ in L∞ (K) imply Σ̂n → Σ. For the drift, Lemma 4.7 applied with An = Σ̂n , A = Σ and the uniform Hessian bound gives q̂n → q uniformly. Inserting into (16) together with the W 1,∞ and b-convergence assumptions yields µ̂n → µ. Local-to-ambient convergence. The reverse direction closes the loop: once Stage 2/3 surrogates fit the converging targets, the reconstructed ambient coefficients recover (b, Λ). Theorem 4.10. Let K ⊂ Ω be compact, ϕ ∈ C 3 (Ω, RD ) with Dϕ of rank d, µ ∈ C 1 (Ω, Rd ), d ), and define b = Dϕ µ + 1 q and Λ = Dϕ ΣDϕ⊤ . Suppose (ϕ , µ , Σ ) satisfy, Σ ∈ C 1 (Ω, S++ n n n 2 d ), uniformly on K: ϕn → ϕ, Dϕn → Dϕ, µn → µ, Σn → Σ, and for all bounded A ∈ Cb1 (Ω, S++ 2 i 2 i 2 i ⟨A , ∇ ϕn ⟩F → ⟨A , ∇ ϕ ⟩F . Assume supn supK ∥∇ ϕn ∥F < ∞ and that µn , Σn are uniformly bounded on K. Setting bn = Dϕn µn + 21 qn and Λn = Dϕn Σn Dϕ⊤ n , we have bn → b and Λn → Λ uniformly on K. Proof. Write Dϕn µn − Dϕµ = (Dϕn − Dϕ)µn + Dϕ(µn − µ); both terms vanish uniformly by hypothesis. Lemma 4.7 with An = Σn gives qn → q, hence bn → b. For the covariance, the identity ⊤ ⊤ Λn − Λ = Dϕn [(Σn − Σ)Dϕ⊤ n + Σ(Dϕn − Dϕ) ] + (Dϕn − Dϕ)ΣDϕ with uniform bounds on Dϕn , Σ yields Λn → Λ. The following proposition bridges the two theorems from a single joint hypothesis. Proposition 4.11. Fix a compact K ⊂ Ω. Suppose ϕn ∈ C 2 , ϕ ∈ C 3 , and µn , µ, Σn , Σ are C 1 with ∥ϕn − ϕ∥W 2,∞ (K) + ∥µn − µ∥L∞ (K) + ∥Σn − Σ∥L∞ (K) → 0 and inf n inf K σmin (Dϕn ) ≥ s > 0. Then the hypotheses of both Theorems 4.9 and 4.10 hold on K. In particular, bn → b and Λn → Λ uniformly on K. Proof. W 2,∞ (K) convergence gives ∇2 ϕin → ∇2 ϕi uniformly, implying uniform Hessian boundedness and Cauchy–Schwarz: |⟨A , ∇2 ϕin − ∇2 ϕi ⟩F | ≤ ∥A∥F,∞ ∥∇2 ϕin − ∇2 ϕi ∥F,∞ → 0. Uniform convergence on compact K gives uniform boundedness of µn , Σn . Direct comparison with the hypotheses of the two theorems completes the argument. Given link (2) above, the following theorem completes the convergence chain. Theorem 4.12 (Weak convergence). Let Z (n) be processes solving the martingale problems for (µn , Σn ) on Ω, and set X (n) := ϕn (Z (n) ). Let Z solve the martingale problem for (µ, Σ), and set X := ϕ(Z). Assume there exists a compact set K ⊂ Ω such that (n)
P Zt
∈ K for all t ∈ [0, T ] = 1,
P Zt ∈ K for all t ∈ [0, T ] = 1
for every n. Assume further that (µ, Σ) are continuous on K with Σ uniformly positive definite, and that µn → µ, Σn → Σ, and ϕn → ϕ uniformly on K. Finally, assume that for each starting point 15
z ∈ K, the martingale problem for (µn , Σn ) admits a solution Pzn , and the target martingale problem (n) for (µ, Σ) admits a unique solution Pz . Fix a common initial point z0 ∈ K and set Z0 = Z0 = z0 . Then L(Z (n) ) ⇒ L(Z) on C([0, T ]; Rd ), L(X (n) ) ⇒ L(X) on C([0, T ]; RD ). Proof. Since Z (n) and Z remain in K almost surely up to time T , we may extend (µn , Σn ) and (µ, Σ) from K to bounded continuous coefficients on Rd (preserving positive semidefiniteness) without changing the laws of the processes on [0, T ]; see Appendix A.1 for details. The uniform convergence on K ensured by the coefficient-level convergence assumptions transfers to these extensions, so the Stroock–Varadhan theorem [44] yields Pzn ⇒ Pz on C([0, T ]; Rd ). For the ambient processes, uniform convergence ϕn → ϕ on K implies (n)
(n)
sup∥ϕn (Zt ) − ϕ(Zt )∥2 ≤ sup ∥ϕn (x) − ϕ(x)∥2 −→ 0 t≤T
x∈K
deterministically. Applying the continuous mapping theorem to the fixed continuous map γ 7→ ϕ◦γ and combining with the display above gives ϕn (Z (n) ) ⇒ ϕ(Z) on C([0, T ]; RD ). (n)
Corollary 4.13 (MFPT convergence). Under the hypotheses of Theorem 4.12, let τr = inf{t ≥ (n) 0 : ∥Xt − X0 ∥2 ≥ r} ∧ T and τr the corresponding hitting time for X, both capped at T . Assume the limiting process crosses the sphere {∥x − X0 ∥ = r} transversally almost surely: writing τr+ = inf{t : ∥Xt − X0 ∥ > r} ∧ T , P τr < T and τr = τr+ = 1. (n)
Then τr
(n)
⇒ τr in distribution and E[τr ] → E[τr ].
Proof. The hitting-time functional γ 7→ inf{t : ∥γ(t) − γ(0)∥ ≥ r} ∧ T is continuous at every path that crosses the sphere immediately (i.e. τr = τr+ < T ); see [45, Thm. 13.6.1]. By the transversality assumption the limit X lies a.s. in the continuity set, so the continuous mapping theorem gives (n) (n) τr ⇒ τr . Since τr ∈ [0, T ], the family is uniformly bounded and hence uniformly integrable; convergence of means follows.
5
Experiments
The convergence chain of Sections 3–4 has three links: (i) ρ-metric chart quality → (ii) coefficient convergence → (iii) weak/MFPT convergence. This paper controls link (i) via Theorem 4.5; link (ii) is established in Theorems 4.9 and 4.10 under a W 2,∞ chart-convergence assumption (Proposition 4.11); link (iii) is Theorem 4.12 and Corollary 4.13. We validate each link experimentally. The ablation study in Section 5.2 measures chart quality via reconstruction error, tangent-space error, and tangent-space fidelity E, as well as end-to-end coefficient accuracy via Eb and EΛ . Section 5.3 tests the full chain from chart to dynamics via MFPT under both smooth and metastable regimes. Post-training diagnostics verify the well-conditioning assumption σmin (Dϕ) ≥ s of Assumption 4.1.
5.1
Setup
Surfaces and embedding. We evaluate on four Monge-patch surfaces ϕ(u, v) = (u, v, f (u, v))⊤ with (u, v) ∈ [−1, 1]2 : paraboloid (f = u2 + v 2 , positive Gaussian curvature), hyperbolic paraboloid (f = u2 − v 2 , negative curvature), quartic dome (f = (u2 +v 2 ) − (u4 +v 4 )/2, sign-changing curvature), sinusoidal (f = sin(u + v), intrinsically flat). Each surface is embedded in RD by appending 16
KF Fourier coordinate pairs, giving D = 3 + 2KF ; the Monge-patch bound σmin (Dϕ) ≥ 1 is preserved. We test D ∈ {11, 201}. Two dynamics. We consider two latent SDEs with complementary challenges. (i) Overdamped Langevin (MB). The Müller–Brown potential [36] is rescaled to (u, v) ∈ [−1, 1]2 via x = 2.25 u−0.25 and y = 2.25 v +1.0, divided by V0 =200, giving drift µ = −∇VMB and isotropic √ diffusion σ = 2kB T I2 with kB T =0.10. The rescaled potential has three metastable wells and two saddle-point transition channels; this tests whether the learned SDE captures inter-well transition rates. (ii) Rotation drift with state-dependent diffusion (Rot). The drift and diffusion are µ(u, v) = (−v, u)⊤ and 1 + u2 /4 u+v σ(u, v) = . 0 1 + v 2 /4 This exercises all three pipeline stages with anisotropic, position-dependent noise. Full SDE parameters and MB coefficients are in Appendix B.1. Training and penalties. All experiments use 10 random seeds and follow Algorithm 1. Rotation uses N =50 training points, the sparse-data regime where geometric penalties matter most; MB uses N =200, enough to populate the three metastable wells. The ablation study in Section 5.2 compares six conditions: baseline (reconstruction loss only), T (tangent-bundle penalty LT from Section 2.2, λT =1), F (inverse-consistency penalty LF =∥Dπ Dϕ−Id ∥2F , λF =1), C (contractive penalty [40], LC =∥Dπ∥2F , λC =0.01), T+F (both LT and LF ), and ATLAS [48] (the local-chart baseline of Section 1, using Gaussian-kernel blending of drift and diffusion at N landmarks with oracle coefficients). Architecture details are in Appendix B.3. Oracle coefficients. All methods receive exact ambient drift b(xi ) and covariance Λ(xi ) computed from the known SDE, isolating the effect of chart-quality regularization from coefficient-estimation noise. In practice, (b, Λ) would be estimated from short-burst trajectory data [12, 48]; quantifying the resulting perturbation to the pipeline losses is an open problem. Landmarks and baselines. Landmarks are placed via a greedy δ-net in the induced Riemannian metric: candidates are accepted greedily so that consecutive landmarks are separated by at least δ, with δ chosen by binary search to yield approximately N points. This quasi-uniform design follows the ATLAS construction [12] and provides good geometric coverage at a given landmark budget. ATLAS is included in the coefficient ablation but excluded from MFPT tables because its ambientspace simulation becomes infeasible at D=201 and its local-chart interpolation degrades under sparse landmarks. Statistical tests. Unless otherwise noted, significance is assessed by one-sided paired Wilcoxon signed-rank tests for directional AE-vs-baseline comparisons. Table 2 (rotation MFPT) uses paired t-tests with common Brownian noise and reports means. No multiple-comparison correction is applied; marginal p-values near 0.05 in Table 3 should be interpreted as directional evidence. Metrics. We evaluate at five levels: (i) chart quality: reconstruction error ∥x − ϕ(π(x))∥2 and tangent-space error ∥Pθ − P ∥2F on 500 held-out test points, reported as per-seed medians; (ii) tangent-space fidelity: E := ∥(I − Pθ )(b − 21 q)∥2 , the normal-space residual of the Itô-corrected drift; (iii) end-to-end coefficient quality: Eb := ∥Dϕ µ̂ + 12 q(σ̂σ̂ ⊤ ) − b∥2 (ambient drift error) and EΛ := ∥Dϕ σ̂σ̂ ⊤ Dϕ⊤ − Λ∥2F (ambient covariance error), reported as per-seed medians over 200 heldout evaluation points; (iv) trajectory quality: we use two MFPT observables. Under rotation drift, radial MFPT at radius r measures the mean first-passage time to ambient distance r from the start, averaged over 500 trajectories with T =2 and ∆t=0.01; under the transversality condition of Corollary 4.13, convergence of means follows from weak convergence. Under MB Langevin, we report per-pair inter-well MFPT for the W0→W1 and W0→W2 transitions, using 2000 trajectories with
17
T =50 and ∆t=0.005. The learned and ground-truth processes are driven by the same Brownian increments (common random numbers) so that the reported difference isolates model error from Monte Carlo variance. By the Kramers rate [27, 19] k ∼ exp(−∆V /kB T ), inter-well MFPT is exponentially sensitive to drift errors in the barrier region, making it a stringent empirical stress test beyond the fixed-horizon guarantee of Theorem 4.12. Both MFPT metrics are reported as relative error vs. ground-truth simulation on the true manifold; further algorithmic details (core-set well assignment, dwell confirmation, censoring) are in Appendix B.4. (v) Extrapolation: reconstruction error at increasing distance δ beyond the [−1, 1]2 training domain, evaluated on 10 seeds at D=11.
5.2
Chart and coefficient ablation
Table 1 reports chart quality, measured by reconstruction error, tangent error, and E on the paraboloid at D=11 and D=201 under both dynamics with 10 seeds each; full results on all four surfaces are in Appendix B.5. Table 1: Ablation study on the paraboloid, 10 seeds, medians reported. Bold = best per column. ∗∗ p<0.01 vs. baseline, paired Wilcoxon test. E: tangent-space fidelity; E , E : true end-to-end Λ b ambient drift and covariance errors from the full pipeline. F and C coincide at MB D=201: without tangent alignment, both collapse across all 10 seeds to a degenerate chart with near-zero decoder-Jacobian singular values, underscoring the need for the tangent-bundle penalty. Rotation drift (N =50) D=11
ATLAS baseline T∗∗ F C T+F∗∗
D=201
Rec.
Tang.
E
Eb
EΛ
.008 .001 <.001 .017 .001 .001
.185 .029 .002 .490 .055 .005
.104 .536 .032 5.60 .736 .063
.10 3.67 3.51 5.46 2.62 .41
2.95 4.54 1.48 22.7 6.27 1.19
Rec.
Tang.
E
Eb
EΛ
.001
.139 .093 .004 2.00 1.64 .035
.140 1.67 .320 11.7 10.2 .540
.14 7.54 9.35 11.7 12.2 1.86
3.15 5.79 2.34 53.5 48.6 3.83
<.001 <.001
.006 .006 .001
MB Langevin (N =200) D=11
ATLAS baseline T∗∗ F C T+F∗∗
D=201
Rec.
Tang.
E
Eb
EΛ
Rec.
Tang.
E
Eb
EΛ
.001
.012 .003 <.001 .015 .430 <.001
.773 .006 <.001 .039 1.64 <.001
.465 .029 .017 .088 3.95 .011
<.001
<.001 <.001
.014 .007 <.001 2.00 2.00 <.001
.779 .014 <.001 11.3 11.3 <.001
.511 .040 .024 10.8 10.8 .022
<.001
.001 <.001 .006 .062 <.001
<.001 <.001 <.001
.002 <.001
<.001
.002 .002 <.001
<.001 <.001
.426 .426 <.001
The tangent-bundle penalty T directly minimizes the projector term of the ρ-metric (Section 2.2). At D=11, T alone reduces tangent error 15× (from .029 to .002) and E by 17×. Under MB Langevin (N =200), T and T+F both achieve E < 10−3 , while baseline sits at .006 and ATLAS at .77, two to three orders of magnitude worse. Adding F yields the best reconstruction, extrapolation (Figure 1), and end-to-end drift accuracy Eb (0.41 vs T’s 3.51 at D=11), because the encoder-pullback target requires the inverse consistency that F enforces. F alone and C penalize Dπ without a dynamics-informed tangent target, so they do not reduce the projector term that drives the ρ-metric. Under MB at D=201, both collapse to a degenerate AE with E > 10 and tangent error saturating at 2.0, indicating that the learned tangent space is es18
Reconstruction error
sentially unrelated to the true one. ATLAS [48] interpolates oracle coefficients via Gaussian kernels and achieves the lowest Eb under rotation (0.10–0.14), confirming that its coefficient approximation is locally accurate. However, its tangent-space fidelity is limited by kernel smoothing: E=0.15 under rotation (vs 0.053 for T) and E=0.77 under MB (vs < 10−3 for T+F). The three computational bottlenecks from Section 1 persist: exponential landmark scaling, per-step re-projection onto M , and O(mD2 ) simulation cost, making ATLAS infeasible at D=201. The columns Eb and EΛ report the true ambient drift and covariance errors of the full three-stage pipeline. Under rotation, T+F reduces Eb by 9× at D=11 and 4× at D=201 relative to baseline. T alone achieves good EΛ but worse Eb than baseline, because the encoder-pullback drift target assumes Dπ Dϕ ≈ Id ; without F enforcing this, the targets are biased and the drift network learns biased coefficients. This validates the complementary roles of T and F: T aligns the tangent space, F ensures the encoder–decoder consistency that the pullback formula requires. Under MB, all AE conditions achieve very small Eb and EΛ , with T+F yielding the best drift accuracy (Eb =.011 at D=11). Post-training diagnostics confirm that σmin (Dϕθ ) ≥ s > 0 holds across all T+F models, validating the hypothesis of Assumption 4.1. Figure 1 shows that T+F consistently extrapolates best, with error growing 2–3× slower than baseline at δ=0.3. F alone produces the worst extrapolation: without tangent alignment, the decoder diverges off-manifold.
Paraboloid
1.5
2
Quartic dome
1.0
1 0
Hyp. paraboloid
0.5 0.1
0.2
Distance
0.3
0.2
0.2
0.1
0.1
0.0
0.0 0.0
0.3
0.3
0.0
0.1
0.2
Distance
0.3
0.0
0.1
0.2
Distance
0.3
0.0
Sinusoidal Baseline T F T+F
0.0
0.1
0.2
Distance
0.3
Figure 1: Reconstruction error vs extrapolation distance δ beyond the [−1, 1]2 training domain (D=11, 10 seeds, mean ± std). T+F extrapolates best on all four surfaces.
5.3
From coefficients to dynamics
We now test whether improved chart and coefficient quality translates to better dynamics. The ablation above shows that F alone and C fail to learn a meaningful tangent space, and ATLAS cannot simulate at D=201, so the dynamics comparison focuses on baseline, T, and T+F. Rotation drift Table 2 reports radial MFPT relative error at r=2 across all four surfaces. T+F reduces error by 50–70% relative to baseline and is significant (p<0.05) on all surfaces at both dimensions, with T alone showing intermediate gains at D=201. Radial MFPT is a smooth functional with no barrier sensitivity, so the Eb improvement from the ablation translates directly to dynamics, consistent with Corollary 4.13. Metastable dynamics Inter-well MFPT goes beyond the fixed-horizon weak convergence of Theorem 4.12: the Kramers rate [27, 19] k ∼ exp(−∆V /kB T ) makes it exponentially sensitive to drift errors in the barrier region. We simulate 2000 trajectories starting near well W0 with paired 19
Table 2: Radial MFPT relative error (%) at r=2 under rotation drift (N =50, δ-net landmarks, 10 seeds, mean). ∗∗ p<0.01, ∗ p<0.05 vs. baseline, paired t-test with common noise. D=11
D=201
Surface
base
T
T+F
base
T
T+F
Paraboloid Hyp. parab. Quartic dome Sinusoidal
39.8 49.9 25.7 16.9
37.6 32.2∗ 27.1 16.1
16.0∗ 18.7∗∗ 11.0∗∗ 6.0∗∗
38.3 51.9 29.3 25.9
38.9 32.3∗ 14.8∗∗ 19.3∗∗
19.2∗∗ 20.6∗∗ 12.1∗∗ 7.7∗∗
Brownian noise to isolate model error from Monte Carlo variance and report per-pair MFPT for the W0→W1 and W0→W2 transitions; see the Appendix B.4 for details. Table 3 shows that T+F achieves the lowest MFPT error on most surface–transition pairs, though several p-values are near 0.05 and no multiple-comparison correction is applied, so these results should be interpreted as directional evidence rather than definitive. On the paraboloid, W0→W1 error drops from 5.1% to 2.2% at D=11 (p<0.05) and W0→W2 from 4.3% to 1.7% at D=201 (p<0.05). On the hyperbolic paraboloid at D=201, T+F reduces W0→W1 from 6.6% to 4.1% (p<0.05) and T alone reaches 2.9% on W0→W2. Quartic dome and sinusoidal results are directionally consistent, with T+F achieving the lowest error on most transitions and all conditions in the 1–5% range. Preliminary oracle-drift experiments (replacing the learned Stage 2 with exact latent drift) suggest that, at D=201, the remaining error is limited by drift learning near the saddle rather than by autoencoder quality. Table 3: Per-pair MFPT error (%) under MB Langevin (N =200, δ-net landmarks, 10 seeds, median). Only W0→W1 and W0→W2 reported, since all trajectories start at W0. Bold = lowest error per row. ∗ p<0.05 vs. baseline, one-sided paired Wilcoxon with common noise. D=11 Surface
6
Pair
base
T ∗
D=201 T+F
base
T
T+F
∗
Paraboloid
W0→W1 W0→W2
5.1 2.7
4.4 1.3
2.2 2.2
4.1 4.3
4.9 2.6
3.2 1.7∗
Hyp. parab.
W0→W1 W0→W2
5.3 3.3
4.8 2.1
4.3 2.7
6.6 6.8
4.8 2.9
4.1∗ 4.0∗
Quartic dome
W0→W1 W0→W2
4.3 2.1
4.5 2.3
3.7 1.4∗
5.5 4.4
4.9 4.3
4.3 2.9
Sinusoidal
W0→W1 W0→W2
2.5 3.1
4.3 1.8
2.1 1.2∗
2.8 2.5
2.9 2.1
2.3 1.2∗
Discussion and Future Work
We have developed a geometric regularization framework for autoencoders trained on data from non-singular diffusions on Riemannian submanifolds. The ρ-metric controls first-order geometry through the tangent-bundle penalty, achieves the same generalization rate as Sobolev H 1 training for chart quality in an oracle ρ-ERM (Theorem 4.5), and under stronger W 2,∞ regularity, chartlevel errors propagate controllably to the ambient SDE coefficients as shown in Section 4.3. The 20
encoder-pullback drift target (2) provides exact latent drift via Itô’s formula applied to the learned encoder, avoiding the systematic bias of the decoder-side formula (Proposition 4.6). The two penalties T and F serve complementary roles. T is chart-invariant by construction, built from Λ (Lemma 2.4), and controls tangent-space orientation. F enforces the coordinateinvariant condition Dπ Dϕ = Id ; the Frobenius penalty ∥Dπ Dϕ − Id ∥2F itself is not invariant under reparametrization, but this is moot in practice because the latent coordinates are learned jointly and only one chart is ever trained. Together, T+F produces charts with well-aligned tangent bundles and accurate encoder–decoder inverses, which determine the quality of the induced latent SDE. To assess the downstream effect on dynamics, we use MFPT, which depends only on the law of the learned process and avoids the coupling dependence of pathwise metrics. Under rotation drift, T+F translates directly to better dynamics on all four surfaces. Under metastable Langevin dynamics, T+F achieves the lowest MFPT error on most surface–transition pairs. T+F also extrapolates 2–3× better than baseline beyond the training domain (Figure 1), a prerequisite for multi-chart extensions. Future work. The single-chart theory extends to a finite atlas via chart-wise penalties weighted by a partition of unity, requiring control of transition-map Lipschitz constants. The encoder Hessian gap not targeted by T+F could be addressed by second-order drift objectives such as penalizing pullback-drift roughness [7] or operator-aligned regression preserving dominant timescales [26]. Finally, replacing oracle (b, Λ) with estimates bootstrapped from short trajectory bursts [12, 48] and quantifying the resulting perturbation to the ρ-loss is an important open problem.
A
Deferred proofs
A.1
Error propagation: additional details
The coefficient-propagation theorems and bridging proposition are stated and proved in the main text (Section 4.3). We record here only the discussion of the W 2,∞ assumption and the pathwise Lp error bound that complements the weak-convergence result. Discussion of the W 2,∞ assumption. The bridging proposition requires ϕn → ϕ in W 2,∞ . Theorem 4.5 controls the ρ-risk (reconstruction + tangent-bundle error), which lies strictly between L2 and H 1 ; it does not imply W 1,∞ or W 2,∞ convergence. The additional second-order control is needed because the Stage 2 drift target involves ⟨Σ, ∇2 ϕ⟩, so uniform convergence of the decoder Hessian is required to pass to the limit. Three remarks are in order. 1. Uniform Hessian boundedness. For feed-forward networks with bounded weights and bounded activation derivatives (true of tanh, sigmoid, and softplus), all second derivatives of ϕθ are automatically bounded by a constant depending only on the weight norms and depth. Thus the uniform Hessian bound supn supK ∥∇2 ϕ̂in ∥F ≤ CK holds whenever the depth and weight-norm bounds are uniform in n (as is the case for the hypothesis class of Assumption 4.2). 2. W 2,∞ convergence. Uniform Hessian boundedness does not imply that Hessians converge. In general, the ρ-loss minimizers ϕn need not converge in W 2,∞ even if they converge in W 1,∞ . We therefore present W 2,∞ convergence as a sufficient condition for the propagation chain, not as a consequence of the training objective. 3. Modular design. The convergence pipeline is deliberately modular: Theorem 4.5 controls chart quality (link (i)), the bridging proposition controls coefficient propagation under W 2,∞ (link (ii)), and Stroock–Varadhan gives weak convergence (link (iii)). Closing the gap between links (i)
21
and (ii)—showing that ρ-minimizers enjoy W 2,∞ convergence, or finding a weaker sufficient condition—is an open problem noted in Section 4.3. Weak convergence via Stroock–Varadhan. The classical Stroock–Varadhan theorem [44] provides weak convergence of solutions to the martingale problem under uniform convergence of coefficients. D and b : RD → RD be locally bounded Theorem A.1 (Stroock–Varadhan). Let Λ : RD → S+ measurable functions which are continuous in x, and assume that for each x ∈ RD the martingale problem for (b, Λ) starting from x has exactly one solution Px . Suppose that for each n ≥ 1, D and b : RD → RD are measurable functions satisfying, for all R > 0: Λn : RD → S+ n sup sup ∥bn (x)∥2 + ∥Λn (x)∥F < ∞ n≥1 ∥x∥≤R
and lim sup
n→∞ ∥x∥≤R
∥Λn (x) − Λ(x)∥F + ∥bn (x) − b(x)∥2 = 0.
Let Pxn be a solution to the martingale problem for (bn , Λn ) starting from x. Then Pxn → Px weakly as measures on path space. Coefficient extension. To apply Theorem A.1, it suffices to extend (µn , Σn ) and (µ, Σ) from K to bounded continuous coefficients on Rd ; Stroock–Varadhan requires only continuity and local boundedness. A bounded linear extension operator E : C(K) → Cb (Rd ) extends each field, and composing d preserves positive semidefiniteness the matrix extensions with the nearest-point projection onto S+ and compact-uniform convergence. Since the processes remain in K a.s., their laws are unchanged. The proof of the weak convergence theorem appears in Section 4.3. The pathwise error bound below complements the weak convergence result; its proof uses standard BDG and Gronwall tools. Theorem A.2 (Lp pathwise error bound). Let Xt and X̂t be strong solutions of dX = µ(X) dt + σ(X) dW and dX̂ = µ̂(X̂) dt + σ̂(X̂) dW driven by the same Brownian motion W , with X0 = X̂0 and σ = Σ1/2 , σ̂ = Σ̂1/2 , both principal square roots. Assume there exists a compact set K ⊂ Rd such that P Xs , X̂s ∈ K for all s ∈ [0, T ] = 1. Suppose ∥µ̂ − µ∥L∞ (K) ≤ ϵµ ,
∥Σ̂ − Σ∥L∞ (K) ≤ ϵΣ .
Assume µ, µ̂ are Cµ -Lipschitz and σ, σ̂ are Cσ -Lipschitz on K. Assume moreover that there exists λ0 > 0 such that λmin (Σ(x)), λmin (Σ̂(x)) ≥ λ0 for all x ∈ K. Then for any p ≥ 2 and 0 ≤ t ≤ T , there exists C2 = C2 (p, d, T, Cµ , Cσ ) such that 1 p E (∆∗t )p ≤ C2 ϵpµ + ϵ , p/2 Σ 2p λ0
(17)
where ∆∗t = sups≤t ∥Xs − X̂s ∥. Proof. Write ∆s = Xs − X̂s . The BDG inequality [41] and Hölder’s inequality give, for C = C(p, d, T ), Z t h i E (∆∗t )p ≤ C E ∥µ̂(X̂s ) − µ(Xs )∥p + ∥σ̂(X̂s ) − σ(Xs )∥pF ds. (18) 0
22
The drift and diffusion coefficient errors decompose via triangle inequality into approximation error plus Lipschitz feedback: ϵΣ + Cσ ∥x̂ − x∥, ∥σ̂(x̂) − σ(x)∥F ≤ 2√ λ
∥µ̂(x̂) − µ(x)∥ ≤ ϵµ + Cµ ∥x̂ − x∥,
0
where Cµ , Cσ are Lipschitz constants on K and the diffusion bound uses the perturbation inequality for positive-definite square roots with eigenvalues ≥ λ0 . Substituting into (18) and applying (a + b)p ≤ 2p−1 (ap + bp ) yields Z t ∗ p ∗ p ϵpΣ ′′ p ′ E (∆s ) ds + C T ϵµ + E (∆t ) ≤ C . p/2 0 2p λ0 Gronwall’s lemma [41, Lemma 11.11] gives (17). Practical limitations of the pathwise bound. Theorem A.2 is mathematically correct but provides an upper bound that is not tight in practice. We conducted a controlled diagnostic experiment to test whether improving the coefficient errors ϵµ , ϵΣ translates to pathwise trajectory improvement; the results reveal three structural reasons why the bound is loose in the learned-chart setting. (i) Gronwall amplification. The constant C2 in (17) arises from Gronwall’s lemma and grows as ′ eC T , where C ′ depends on the Lipschitz constants Cµ , Cσ of both the true and learned coefficients. In the end-to-end pipeline, the learned coefficients are obtained by inverting the ambient SDE through the learned chart via the Jacobian, Hessian, and metric inverse, each of which has its own Lipschitz constant. For a neural network with Tanh activations on a compact domain, these constants are moderate but their product enters the exponential, making C2 substantially larger than 1 even at T = 1. As a result, a 30–60% reduction in ϵµ can be multiplied by a Gronwall factor large enough that the bound does not meaningfully tighten. (ii) Chart reconstruction error dominates. The bound (17) assumes that the two processes X, X̂ evolve in the same coordinate system and differ only in their coefficients. In the learnedchart pipeline, however, the ground-truth process is simulated in the true local coordinates (u, v), while the learned process evolves in the latent coordinates z and is mapped to ambient space via the learned decoder ϕθ . Even if the latent-space coefficients were exact, the trajectory error would include a chart reconstruction component: at each time step the decoded position ϕθ (πθ (x)) deviates from the true position x by the autoencoder’s reconstruction error. This reconstruction error accumulates over T /∆t Euler steps and, in our experiments, accounts for the majority of the observed trajectory error. (iii) Magnitude mismatch. The coefficient errors ϵµ , ϵΣ enter the bound at the same power p, but in practice the two terms contribute very differently: the absolute magnitudes can differ by an −p/2 order of magnitude, so the bound treats them symmetrically up to the factor λ0 , and even a large relative improvement in the smaller term has negligible effect on the bound’s value. The pathwise bound establishes the consistency of the learned dynamics: if ϵµ , ϵΣ → 0, then ∆∗t → 0. However, the multiplicative constants and the chart-reconstruction floor prevent it from being an actionable predictor of trajectory quality at finite sample size. The distributional convergence of the weak convergence theorem (Theorem 4.12), which avoids these pathwise difficulties, provides a more practically relevant guarantee.
23
B
Additional experiments
B.1
Müller–Brown potential
The standard Müller–Brown potential [36] is VMB (x, y) =
4 X
Ai exp ai (x − x0i )2 + bi (x − x0i )(y − yi0 ) + ci (y − yi0 )2 ,
i=1
with the coefficients in Table 4. We rescale to (u, v) ∈ [−1, 1]2 via the affine map x = 2.25 u − 0.25,
y = 2.25 v + 1.0,
and divide by V0 = 200 to obtain the rescaled potential Ṽ (u, v) = VMB (x(u), y(v))/V0 . The overdamped Langevin SDE on the Monge-patch surface is p u kB T = 0.10. d = −∇Ṽ dt + 2kB T I2 dWt , v Table 4: Müller–Brown potential coefficients.
B.2
i
Ai
ai
bi
ci
x0i
yi0
1 2 3 4
−200 −100 −170 15
−1 −1 −6.5 0.7
0 0 11 0.6
−10 −10 −6.5 0.7
1 0 −0.5 −1
0 0.5 1.5 1
Landmark sampling
Training landmarks are sampled via greedy δ-net subsampling in local coordinates, following the ATLAS construction [12]. A large uniform candidate pool is generated in the training domain ([−1, 1]2pfor rotation, [−0.55, 0.55]2 for MB), and each point is accepted if a local metric approximation ∆u⊤ g(ū) ∆u, where g is the induced metric evaluated at the midpoint ū, to all previously accepted points exceeds δ. The separation δ is chosen by binary search to yield approximately N landmarks. The candidate pool size is max(10 000, 100N ).
B.3
Architecture and training
The encoder πθ : RD → Rd and decoder ϕθ : Rd → RD are feedforward networks with two hidden layers and tanh activations; width is 64 at D=11 and 256 at D=201. The Stage 2 drift network µ̂ω : Rd → Rd has hidden layers [64, 64] under rotation and [256, 256, 256] under MB Langevin. The Stage 3 diffusion network σ̂ψ : Rd → Rd×d has hidden layers [64, 64] in all conditions. All hidden layers use tanh activations. Hyperparameters differ by dynamics. Under rotation drift: Stage 1 trains for 500 epochs (Adam, lr = 0.005, batch size 20); Stages 2 and 3 each train for 300 epochs (Adam, lr = 0.001) with the AE frozen. Under MB Langevin: Stage 1 trains for 4000 epochs (Adam, lr = 0.005, batch size 20); Stages 2 and 3 each train for 3000 epochs (Adam, lr = 0.001) with the AE frozen.
24
B.4
MFPT computation
We describe the trajectory simulation and MFPT extraction procedure in detail, as these methodological choices affect the reported results. Ground-truth simulation. The reference SDE is integrated in the true local coordinates (u, v) using Euler–Maruyama with step size ∆t = 0.005 (MB) or ∆t = 0.01 (rotation), driven by i.i.d. Gaussian increments ∆Wk ∼ N (0, Id ). Learned simulation. The learned latent SDE (µ̂ω , σ̂ψ ) is integrated from encoded initial conditions z0 = πθ (x0 ) using the same Brownian increments ∆Wk as the ground truth via common random numbers. This paired-noise design reduces Monte Carlo variance in the MFPT error estimate, so that the reported difference isolates model error from sampling noise. The learned latent trajectory is decoded to ambient coordinates via ϕθ , and the first two components (u, v) are used for well assignment. MB initial conditions. All 2000 trajectories start near well W0: x0 ∼ N (zW0 , 0.01 I2 ) clipped to [−0.55, 0.55]2 , where zW0 ≈ (−0.137, 0.196) in rescaled MB coordinates. The three well centers in rescaled (u, v) are W0 = (−0.137, 0.196), W1 = (0.388, −0.432), W2 = (0.089, −0.237). Core-set well assignment (MB). Each time step is assigned to well i ∈ {0, 1, 2} if the decoded (u, v) lies within Euclidean radius rcore = 0.08 of the i-th well center. Steps outside all cores are labeled −1 for unassigned. The dwell confirmation (ndwell = 10 consecutive steps in a core) is applied during passage extraction, not during per-step assignment. Pairwise MFPT extraction. For each trajectory, we scan the well-assignment sequence and record confirmed passages: a passage from well i to well j begins when the trajectory is in core i and ends at the first confirmed arrival in core j, i.e. ndwell consecutive steps in j. The passage time is backdated to the start of the dwell window. Multiple passages per trajectory are collected since the scan resumes from the arrival time. The pairwise MFPT τij is the mean over all recorded i → j passages across 2000 trajectories. Error metric. Since all trajectories start near W0, only the first-passage transitions W0→W1 and W0→W2 have model-independent initial conditions. Follow-on transitions (e.g. W1→W0) start from wherever the learned dynamics places the trajectory after the first arrival, introducing a model-dependent selection bias that contaminates the comparison: the entrance distribution in well j differs between GT and learned dynamics, so subsequent MFPTs out of j are not directly comparable. We therefore report the per-pair relative error MFPT(i→j) = err
|τ̂ij − τij | , τij
separately for W0→W1 and W0→W2 in Section 5.3. Pairs with non-finite values, indicating no observed transitions, are excluded and flagged. Censoring. Trajectories whose decoded (u, v) ever leaves [−1, 1]2 or contains non-finite values are censored: their well assignments are set to −1 for all time steps, excluding them from all MFPT metrics. The censoring is applied identically to GT and learned trajectories as a matched treatment. The exit fraction is reported per condition. Radial MFPT (rotation drift). Under rotation dynamics, we instead measure the mean time for a trajectory to first reach ambient distance r from its own starting point. The same paired-noise and censoring procedures apply.
B.5
Full ablation across surfaces
Tables 5 and 6 extend the ablation in Section 5.2 (paraboloid only) to all four surfaces (10 seeds, medians), including end-to-end coefficient errors Eb and EΛ . 25
Table 5: Full ablation under rotation drift (N =50, 10 seeds, medians). Bold = best per column. D=11
D=201
Surface
Cond.
Tang.
E
EΣ
Eb
EΛ
Tang.
E
EΣ
Eb
EΛ
Paraboloid
ATLAS baseline T F C T+F
.19 .029 .002 .49 .055 .005
.10 .54 .03 5.60 .74 .06
.85 .63 .04 20.4 1.53 .13
.10 3.67 3.51 5.46 2.62 .41
2.95 4.54 1.48 22.7 6.27 1.19
.14 .093 .004 2.00 1.64 .035
.14 1.67 .32 11.7 10.2 .54
1.07 3.01 .15 53.5 41.8 1.02
.14 7.54 9.35 11.7 12.2 1.86
3.15 5.79 2.34 53.5 48.6 3.83
Hyp. parab.
ATLAS baseline T F C T+F
.16 .018 .059 .34 .039 .149
.11 .39 8.27 1.98 .57 1.29
.82 .54 1.07 19.7 .83 3.08
.11 4.82 48.1 2.10 2.81 5.89
2.65 3.27 9.14 22.7 4.47 7.65
.13 .086 .158 1.99 1.34 .023
.14 1.85 5.65 5.05 4.21 .27
.93 2.20 3.66 53.1 47.6 .69
.14 10.9 31.5 5.07 5.05 .96
2.49 6.70 18.5 53.1 49.2 2.15
Quartic dome
ATLAS baseline T F C T+F
.075 .077 .010 .37 .152 .009
.30 1.07 .18 2.64 1.61 .17
.31 .66 .14 6.24 1.96 .12
.30 2.08 2.17 2.88 2.80 .45
.54 1.06 .77 6.43 3.74 .42
.060 .162 .024 2.00 1.91 .036
.38 2.16 1.37 5.37 4.97 .71
.44 2.08 .49 20.8 19.9 .52
.38 3.72 2.39 5.41 5.41 1.22
.71 3.12 .79 20.8 19.9 1.19
Sinusoidal
ATLAS baseline T F C T+F
.092 .038 .002 .20 .093 .001
.08 .32 .01 .93 .49 .02
.23 .47 .04 3.21 1.33 .02
.08 1.29 1.25 .95 .67 .07
.45 .61 .15 3.40 1.59 .22
.071 .144 .006 2.00 1.16 .023
.11 1.57 .17 5.17 3.04 .34
.28 2.19 .11 32.0 14.3 .51
.11 3.19 2.10 5.23 4.99 .67
.48 2.28 .20 32.0 16.3 .69
References [1] John Armstrong and Damiano Brigo. Intrinsic stochastic differential equations as jets. Proceedings of the Royal Society A, 474(2210):20170559, 2018. [2] John Armstrong, Damiano Brigo, and Emilio Ferrucci. Projections of SDEs onto submanifolds. Information Geometry, 7(Suppl 1):397–427, 2024. [3] Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg. Latent space oddity: on the curvature of deep generative models. In International Conference on Learning Representations (ICLR), 2018. [4] Peter L. Bartlett, Olivier Bousquet, and Shahar Mendelson. Local Rademacher complexities. The Annals of Statistics, 33(4):1497–1537, 2005. [5] Nils Berglund and Barbara Gentz. Noise-Induced Phenomena in Slow-Fast Dynamical Systems: A Sample-Paths Approach. Probability and its Applications. Springer, 2006. [6] Tom Bertalan, Felix Dietrich, Igor Mezić, and Ioannis G. Kevrekidis. On learning Hamiltonian systems from data. Chaos, 29(12):121107, 2019. [7] Andreas Bittracher, Péter Koltai, Stefan Klus, Ralf Banisch, Michael Dellnitz, and Christof Schütte. Transition manifolds of complex metastable systems: Theory and data-driven computation of effective dynamics. J. Nonlinear Sci., 28(2):471–512, 2018. 26
Table 6: Full ablation under MB Langevin (N =200, 10 seeds, medians). Bold = best per column. D=11
D=201
Surface
Cond.
Tang.
E
EΣ
Eb
EΛ
Tang.
E
EΣ
Eb
EΛ
Paraboloid
ATLAS baseline T F C T+F
.001 .002 < .001 .01 .44 < .001
.43 .006 < .001 .04 1.59 < .001
< .001 < .001 < .001 .003 .18 < .001
.43 .03 .02 .08 5.98 .01
< .001 < .001 < .001 .004 .18 < .001
.001 .003 < .001 2.00 2.00 < .001
.43 .007 < .001 11.3 11.3 < .001
< .001 < .001 < .001 .42 .42 < .001
.43 .04 .02 11.3 11.3 .01
< .001 < .001 < .001 .42 .42 < .001
Hyp. parab.
ATLAS baseline T F C T+F
.001 .001 < .001 .01 .44 < .001
.38 .003 < .001 .03 1.22 < .001
< .001 < .001 < .001 .003 .18 < .001
.38 .03 .01 .07 5.52 .01
< .001 < .001 < .001 .003 .19 < .001
.001 .007 < .001 2.00 1.99 < .001
.33 .02 < .001 12.0 11.9 < .001
< .001 < .001 < .001 .42 .42 < .001
.33 .09 .02 12.0 12.0 .01
< .001 .002 < .001 .42 .42 < .001
Quartic dome
ATLAS baseline T F C T+F
.001 .002 < .001 .01 .39 < .001
.36 .003 < .001 .03 1.40 < .001
< .001 < .001 < .001 .002 .13 < .001
.36 .03 .02 .05 5.11 .01
< .001 < .001 < .001 .002 .13 < .001
.001 .004 < .001 2.00 2.00 < .001
.32 .008 < .001 11.0 11.0 < .001
< .001 < .001 < .001 .37 .37 < .001
.32 .04 .03 11.0 11.0 .01
< .001 < .001 < .001 .37 .37 < .001
Sinusoidal
ATLAS baseline T F C T+F
.001 .001 < .001 .01 .05 < .001
.38 .002 < .001 .02 .10 < .001
< .001 < .001 < .001 .003 .02 < .001
.38 .03 .02 .05 4.07 .02
< .001 < .001 < .001 .003 .03 < .001
.001 .004 < .001 2.00 2.00 < .001
.34 .008 < .001 13.3 13.8 .001
< .001 .001 < .001 .68 .68 < .001
.34 .06 .02 13.3 13.9 .02
< .001 .001 < .001 .68 .68 < .001
[8] Kathleen Champion, Bethany Lusch, J Nathan Kutz, and Steven L Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45):22445–22451, 2019. [9] G.S. Chirikjian. Stochastic Models, Information Theory, and Lie Groups, Volume 1: Classical Results and Geometric Methods. Applied and Numerical Harmonic Analysis. Birkhäuser Boston, 2009. [10] Ronald R. Coifman, Ioannis G. Kevrekidis, Stéphane Lafon, Mauro Maggioni, and Boaz Nadler. Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems. Multiscale Modeling & Simulation, 7(2):842–864, 2008. [11] Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006. [12] Miles Crosskey and Mauro Maggioni. ATLAS: a geometric approach to learning highdimensional stochastic systems near manifolds. Multiscale Model. Simul., 15(1):110–156, 2017. [13] Valentin De Bortoli, Emile Mathieu, Michael Hutchinson, James Thornton, Yee Whye Teh, and Arnaud Doucet. Riemannian score-based generative modelling. In Advances in Neural Information Processing Systems, volume 35, pages 2406–2422, 2022. [14] Felix Dietrich, Alexei Makeev, George Kevrekidis, Nikolaos Evangelou, Tom Bertalan, Sebastian Reich, and Ioannis G. Kevrekidis. Learning effective stochastic differential equations from 27
microscopic simulations: Linking stochastic numerics to deep learning. Chaos, 33(2):023121, 2023. [15] Nikolaos Evangelou, Felix Dietrich, Eliodoro Chiavazzo, Daniel Lehmberg, Marina Meila, and Ioannis G. Kevrekidis. Double diffusion maps and their latent harmonics for scientific computations in latent space. Journal of Computational Physics, 485:112072, 2023. [16] Daniel Floryan and Michael D Graham. Data-driven discovery of intrinsic dynamics. Nature Machine Intelligence, 4(12):1113–1120, 2022. [17] Hastings Greer, Roland Kwitt, François-Xavier Vialard, and Marc Niethammer. ICON: Learning regular maps through inverse consistency. In Proc. IEEE/CVF Intl. Conf. Computer Vision (ICCV), pages 3396–3405, 2021. [18] László Györfi, Michael Kohler, Adam Krzyżak, and Harro Walk. A Distribution-Free Theory of Nonparametric Regression. Springer Series in Statistics. Springer, 2002. [19] Peter Hänggi, Peter Talkner, and Michal Borkovec. Reaction-rate theory: fifty years after Kramers. Rev. Mod. Phys., 62(2):251–341, 1990. [20] Ali Hasan, João M. Pereira, Sina Farsiu, and Vahid Tarokh. Identifying latent stochastic differential equations. IEEE Transactions on Signal Processing, 70:89–104, 2022. [21] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 2 edition, 2013. [22] Elton P. Hsu. Stochastic Analysis on Manifolds, volume 38 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 2002. [23] Chin-Wei Huang, Milad Aghajohari, Avishek Joey Bose, Prakash Panangaden, and Aaron Courville. Riemannian diffusion models. In Advances in Neural Information Processing Systems, volume 35, 2022. [24] George A Kevrekidis, Mauro Maggioni, Soledad Villar, and Yannis G Kevrekidis. Thinner latent spaces: Detecting dimension and imposing invariance through autoencoder gradient constraints. arXiv preprint arXiv:2408.16138, 2024. [25] Ioannis G. Kevrekidis, C. William Gear, James M. Hyman, Panagiotis G. Kevrekidis, Olof Runborg, and Constantinos Theodoropoulos. Equation-free, coarse-grained multiscale computation: Enabling microscopic simulators to perform system-level analysis. Communications in Mathematical Sciences, 1(4):715–762, 2003. [26] Stefan Klus, Feliks Nüske, Péter Koltai, Hao Wu, Ioannis Kevrekidis, Christof Schütte, and Frank Noé. Data-driven model reduction and transfer operator approximation. J. Nonlinear Sci., 28:985–1010, 2018. [27] H. A. Kramers. Brownian motion in a field of force and the diffusion model of chemical reactions. Physica, 7(4):284–304, 1940. [28] John M. Lee. Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics. Springer, 2nd edition, 2012.
28
[29] Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 404:108973, 2020. [30] Ben Leimkuhler and Charles Matthews. Molecular Dynamics: With Deterministic and Stochastic Numerical Methods, volume 39 of Interdisciplinary Applied Mathematics. Springer, 2015. [31] Xuechen Li, Ting-Kam Leonard Wong, Ricky T. Q. Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), volume 108 of PMLR, pages 3870–3882, 2020. [32] Alec J. Linot and Michael D. Graham. Deep learning to discover and predict dynamics on an inertial manifold. Physical Review E, 101(6):062209, 2020. [33] Hao Liu, Alex Havrilla, Rongjie Lai, and Wenjing Liao. Deep nonparametric estimation of intrinsic data structures by chart autoencoders: Generalization error and robustness. Applied and Computational Harmonic Analysis, 68:101602, 2024. [34] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 9(1):4950, 2018. [35] Emile Mathieu and Maximilian Nickel. Riemannian continuous normalizing flows. In Advances in Neural Information Processing Systems, volume 33, pages 2503–2515, 2020. [36] K. Müller and L. D. Brown. Location of saddle points and minimum energy paths by a constrained simplex optimization procedure. Theoret. Chim. Acta, 53:75–93, 1979. [37] Samuel E Otto and Clarence W Rowley. Linearly recurrent autoencoder networks for learning dynamics. SIAM Journal on Applied Dynamical Systems, 18(1):558–593, 2019. [38] Grigorios A. Pavliotis and Andrew M. Stuart. Multiscale Methods: Averaging and Homogenization, volume 53 of Texts in Applied Mathematics. Springer, 2008. [39] Erez Peterfreund, Ofir Lindenbaum, Felix Dietrich, Tom Bertalan, Matan Gavish, Ioannis G Kevrekidis, and Ronald R Coifman. Local conformal autoencoder for standardized data coordinates. Proceedings of the National Academy of Sciences, 117(49):30918–30927, 2020. [40] Salah Rifai, Pascal Vincent, Xavier Muller, Xavier Glorot, and Yoshua Bengio. Contractive auto-encoders: explicit invariance during feature extraction. In Proceedings of the 28th International Conference on Machine Learning (ICML), ICML’11, page 833–840, Madison, WI, USA, 2011. Omnipress. [41] L. C. G. Rogers and David Williams. Diffusions, Markov Processes and Martingales: Volume 2, Itô Calculus. Cambridge Mathematical Library. Cambridge University Press, 2 edition, 2000. [42] Stefan Schonsheck, Jie Chen, and Rongjie Lai. Chart auto-encoders for manifold structured data, 2019. [43] Stefan C. Schonsheck, Scott Mahan, Timo Klock, Alexander Cloninger, and Rongjie Lai. Semi-supervised manifold learning with complexity decoupled chart autoencoders, 2022. 29
[44] Daniel W. Stroock and S. R. Srinivasa Varadhan. Multidimensional Diffusion Processes, volume 233 of Grundlehren der mathematischen Wissenschaften. Springer-Verlag, Berlin, 1979. [45] Ward Whitt. Stochastic-Process Limits: An Introduction to Stochastic-Process Limits and Their Application to Queues. Springer, New York, 2002. [46] Yahong Yang and Juncai He. Deeper or wider: A perspective from optimal generalization error with sobolev loss. In Proceedings of the 41st International Conference on Machine Learning (ICML), volume 235 of PMLR, pages 56109–56138, 2024. [47] Yahong Yang and Juncai He. Deep neural networks with general activations: Superconvergence in sobolev norms, 2025. [48] Felix X.-F. Ye, Sichen Yang, and Mauro Maggioni. Nonlinear model reduction for slow–fast stochastic systems near unknown invariant manifolds. Journal of Nonlinear Science, 34(1):22, 2024.
30