Paper deep dive
Generalization Bounds for Flow-matching Generative Models for Intrinsically Low-dimensional Data
Saptarshi Chakraborty, Quentin Berthet, Peter L. Bartlett
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 10/6/2026, 10:19:03 AM
Summary
This paper establishes finite-sample generalization bounds for flow-matching generative models, demonstrating that they adapt to the intrinsic low-dimensional structure of data distributions. The authors derive high-probability error bounds in Wasserstein-p distance, showing that convergence rates depend on the intrinsic dimension rather than the ambient dimension, thereby mitigating the curse of dimensionality. This provides a theoretical explanation for the empirical success of flow matching on structured high-dimensional data like images and molecular geometries.
Entities (10)
Relation Signals (7)
Saptarshi Chakraborty → affiliatedwith → University of Michigan
confidence 99% · Affiliation: Department of Statistics, University of Michigan
Quentin Berthet → affiliatedwith → Google DeepMind
confidence 99% · Affiliation: Google DeepMind
Peter L. Bartlett → affiliatedwith → UC Berkeley
confidence 99% · Affiliation: Department of Statistics, University of California, Berkeley
Peter L. Bartlett → affiliatedwith → Google DeepMind
confidence 99% · Affiliation: Google DeepMind
Flow Matching → uses → Wasserstein-p Distance
confidence 95% · We derive finite-sample error bounds on the learned generative distribution, measured in the Wasserstein-p distance
Flow Matching → adaptsto → Intrinsic Dimension
confidence 92% · the convergence exponent depends on the intrinsic rather than ambient dimension.
Flow Matching → mitigates → Curse of Dimensionality
confidence 90% · Our results demonstrate that flow matching naturally adapts to the intrinsic geometry of data and mitigates the curse of dimensionality
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Despite the remarkable empirical success of flow-matching models, their statistical generalization guarantees remain underdeveloped. Existing analyses often impose restrictive assumptions on the estimated velocity field and yield convergence rates that fail to reflect the intrinsic low-dimensional structure common in real data, such as natural images and molecular geometries. In this work, we study the statistical generalization of flow-matching models for learning an unknown distribution $P_{\mathrm{data}}$ from finitely many samples. We derive finite-sample error bounds on the learned generative distribution, measured in the Wasserstein-$p$ distance, for all $p\geq 1$. Specifically, given $n$ i.i.d. samples from $P_{\mathrm{data}}$, we show that, for every $d>d_p^\ast(P_{\mathrm{data}})$ and appropriately chosen network architectures and hyperparameters, the learned distribution $\widehat{P}^{\mathrm{FM}}$ satisfies $ \mathbb{W}_p(\widehat{P}^{\mathrm{FM}},P_{\mathrm{data}}) \lesssim n^{-1/d}+n^{-1/(2p)}\bigl(\log(1/\xi)\bigr)^{1/(2p)}$ with probability at least $1-\xi$, where $d_p^\ast(P_{\mathrm{data}})$ denotes the Wasserstein-$p$ dimension of the target measure. Our results demonstrate that flow matching naturally adapts to the intrinsic geometry of data and mitigates the curse of dimensionality, as the convergence exponent depends on the intrinsic rather than ambient dimension. These guarantees remain meaningful in high-dimensional regimes and provide a theoretical explanation for the empirical success of flow matching on structured data distributions under substantially more relaxed assumptions than those in existing analyses.
Tags
Links
- Source: https://arxiv.org/abs/2610.02663v2
- Canonical: https://arxiv.org/abs/2610.02663v2
Trouble viewing inline? Open PDF directly →
Full Text
165,076 characters extracted from source content.
Expand or collapse full text
Generalization Bounds for Flow-matching Generative Models for Intrinsically Low-dimensional Data Saptarshi Chakraborty †thanks: email: saptarsc@umich.edu Affiliation: Department of Statistics, University of Michigan Quentin Berthet †thanks: email: qerthet@google.com Affiliation: Google DeepMind Peter L. Bartlett †thanks: email: peter@berkeley.edu Affiliation: Department of Statistics, University of California, Berkeley Affiliation: Department of Electrical Engineering and Computer Sciences, UC Berkeley Affiliation: Google DeepMind Abstract Despite its remarkable empirical performance, the theoretical properties of flow-matching in terms of generalization accuracy remain largely unexplored. The current state-of-the-art generalization analyses not only impose stringent assumptions on the estimated velocity field, but also scale poorly with the ambient dimension—precisely the high-dimensional setting where flow matching is most commonly applied. In particular, the data distributions targeted by flow matching models, such as natural images, molecular geometries, etc., are commonly hypothesized to concentrate near a low-dimensional structures embedded in the high-dimensional ambient data space. However, existing analyses fail to reflect this geometry in their derived convergence rates. In this paper, we bridge the gap between the theory and practice of flow matching by establishing generalization bounds on the Wasserstein-p distance between the learned distribution and the target measure, for any p≥1p≥ 1. Specifically, given n independent and identically distributed samples from a target distribution PdataP_data and for every d>dp∗(Pdata)d>d _p(P_data), we show that for appropriate network and hyperparameter choices, p(P^FM,Pdata)≲n−1/d+n−12p(log(1/ξ))12pW_p\! ( P^FM,\,P_data ) n^-1/d+n^- 12p\! ( (1/ξ) ) 12p with probability at least 1−ξ1-ξ, where dp∗(Pdata)d_p (P_data) denotes the Wasserstein-p dimension of the target measure. Our results demonstrate that flow matching naturally adapts to the intrinsic geometry of the data, effectively mitigating the curse of dimensionality: the convergence exponent depends only on the intrinsic dimension, dp∗(Pdata)d_p (P_data), rather than the ambient dimension of the feature space, yielding convergence rates that remain meaningful in high-dimensional regimes and offering a principled theoretical explanation for the strong empirical performance of flow matching on structured data distributions with significantly relaxed assumptions, compared to the existing literature. 1 Introduction Learning high-dimensional probability distributions from a handful of samples is a foundational problem in machine learning, with applications spanning image synthesis (Rombach et al., 2022; Saharia et al., 2022), scientific modeling (Watson et al., 2023; Geffner et al., 2025; Corso et al., 2023), and generative artificial intelligence. The dominant paradigm for this task has shifted decisively from learning a single generator/decoder (Kingma and Welling, 2014; Goodfellow et al., 2020) toward continuous-time generative models such as diffusion models and score-based methods that learn to reverse a fixed noising process (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song and Ermon, 2019; Song et al., 2021), while optimal transport-based approaches construct more direct paths between distributions (Korotin et al., 2023; Bunne et al., 2022). Flow matching (Lipman et al., 2023; Albergo and Vanden-Eijnden, 2023; Liu et al., 2023; Albergo et al., 2025; Chen et al., 2018) has emerged as a particularly compelling framework within this landscape, offering simulation-free training, straight probability paths, and strong empirical performance. The deterministic formulation of flow matching endows the resulting models with several attractive computational advantages, including stable optimization, flexible control over discretization during sampling, and seamless integration with modern Ordinary Differential Equation (ODE) solvers. In practice, flow-matching methods have demonstrated state-of-the-art performance across a range of challenging high-dimensional applications, including text-to-image generation, video synthesis, and molecular design (Ma et al., 2024; Bose et al., 2024; Esser et al., 2024). However, a rigorous understanding of when and why flow-matching generative models succeed remains incomplete, as prior theory has focused on sampling guarantees. Early work on score-based models established such guarantees under stochastic sampling (Block et al., 2020; Bortoli, 2022; Lee et al., 2022), with later results showing that an ϵε-approximate oracle score yields KL error (poly(D,log(1/δ))ϵ−2)O(poly(D, (1/δ))ε^-2) under mild regularity (Chen et al., 2023a; Lee et al., 2023; Albergo and Vanden-Eijnden, 2023; Albergo et al., 2025; Lu et al., 2022; Chen et al., 2023b), where D is the ambient dimension and δ the early stopping time. Stochastic interpolants bridge diffusion and flow matching (Albergo et al., 2025; Albergo and Vanden-Eijnden, 2023): Albergo and Vanden-Eijnden (2023) bound the 22-Wasserstein distance between observed and estimated flow ODE endpoints via the 2L_2-error of the estimated velocity field, assuming Lipschitzness in time, while Albergo et al. (2025) analyze the KL error of flow matching under oracle approximations. For the deterministic sampler, Chen et al. (2023d) derive polynomial discretization bounds that do not extend to learned scores, and Chen et al. (2023c) obtain polynomial guarantees with learned scores via interleaved predictor-corrector updates. Most recently, Chen et al. (2026) established the first polylog(1/ϵ)polylog(1/ε)-sampling guarantee, with ~(dim¯M(Pdata)polylog(1/ϵ)) O( dim_M(P_data)polylog(1/ε)) steps, where dim¯M(Pdata) dim_M(P_data) is the Minkowski dimension (Section 2.3) of the target PdataP_data. All of these works, however, assume oracle access to an ϵε-accurate velocity or score estimate, without establishing its learnability. Despite this progress, existing guarantees for flow matching suffer from the curse of dimensionality: their rates depend exponentially on the ambient dimension D, and thus, become uninformative precisely in the high-dimensional regimes where flow matching is most used. For example, Gao et al. (2024) derived 2W_2 convergence guarantees for flow matching under the assumption that the target distribution is log -concave and admits a bounded density, obtaining 2(P^FM,Pdata)=~(n−1/(D+5))EW_2( P^FM,P_data)= O (n^-1/(D+5) ), when the estimated vector velocity field is uniformly Lipschitz. Zhou and Liu (2025) removed this Lipschitzness assumption on the learned velocity field at the cost of a slower rate, establishing 2(P^FM,Pdata)=~(n−13(D+5))EW_2( P^FM,P_data)= O (n^- 13(D+5) ). Under additional Besov-smoothness assumptions on the target density, Fukumizu et al. (2025) showed that for any δ>0δ>0, 2(P^FM,Pdata)=(n−s+(2κ)−1−δ2s+D)EW_2( P^FM,P_data)=O (n^- s+(2κ)^-1-δ2s+D ), where s is the smoothness parameter of the density and κ≥1/2κ≥ 1/2 is a tunable model parameter. This sits in sharp contrast with the widely held belief that high-dimensional data concentrate near low-dimensional structures (Pope et al., 2020), suggesting the D-dependence is an artifact of analysis rather than intrinsic to the problem—as in classical nonparametric statistics, where kernel estimators and Gaussian process regression achieve minimax rates depending only on intrinsic dimensionality (Bickel and Li, 2007; Kim et al., 2019). Dimension-adaptive, minimax-optimal guarantees have since been established for GANs (Dahal et al., 2022; Huang et al., 2022; Chakraborty and Bartlett, 2025), Wasserstein Autoencoders (Chakraborty and Bartlett, 2024), and diffusion models (Chen et al., 2023d; Oko et al., 2023; Tang and Yang, 2024; Chakraborty et al., 2026), but whether analogous results hold for flow matching remains open. While recent studies have begun exploring flow matching on manifolds (Benton et al., 2024; Pi et al., 2026) and diffusion models for multi-modal low-dimensional data (Wu and Cai, 2026), existing guarantees either rely on restrictive manifold reach and density lower bounds or analyze in-expectation 1W_1/2W_2 metrics. Kumar et al. (2026) partially address this, establishing 2W_2 rates under oracle access to the ODE sampler for data on a smooth, lower-bounded-reach manifold with Hölder-smooth density; however, their analysis relies on strong, unverifiable geometric conditions and excludes even simple distributions (Gaussian, uniform, degenerate). In contrast, our framework provides the sharpest known non-asymptotic, high-probability concentration guarantees for arbitrary p≥1p≥ 1 without requiring manifold reach or global Lipschitz regularity on the learned vector field. Summary of Contributions The main contributions of this paper are as follows. • High-probability Wasserstein guarantees for every p≥1p≥ 1. We develop a general theoretical framework for analyzing the statistical performance of flow-matching generative models. For arbitrary p≥1p≥ 1, we derive non-asymptotic, high-probability bounds on the Wasserstein-p error p(P^FM,Pdata).W_p\! ( P^FM,P_data ). Specifically, when one has access to i.i.d. samples x1,…,xn∼Pdatax_1,…,x_n P_data, we show that, for any d>dp∗(Pdata)d>d _p(P_data), with probability at least 1−ξ1-ξ, p(P^FM,Pdata)≲n−1/dp∗(Pdata)+n−1/(2p)(log(1/ξ))1/(2p),W_p\! ( P^FM,P_data ) n^-1/d_p (P_data)+n^-1/(2p) ( (1/ξ) )^1/(2p), under appropriate choices of the velocity-field network class, discretization scheme, and early-stopping rule. In contrast to existing analyses (Benton et al., 2024; Gao et al., 2024; Zhou and Liu, 2025; Kumar et al., 2026), which are largely expectation-based and often restricted to p∈1,2p∈\1,2\, our results hold for every p≥1p≥ 1 and provide guarantees at arbitrarily high confidence. We also give sufficient conditions for selecting the relevant hyperparameters as functions of the sample size and the population distribution. • Minimal distributional and model assumptions. Our guarantees require no density, smoothness, or Lipschitz assumptions on either the data-generating distribution or the velocity-field class. Despite this generality, the resulting rates adapt to the intrinsic dimension of the target distribution, providing a rigorous statistical explanation for the effectiveness of flow matching on structured high-dimensional data, including images, videos, and molecular geometries that may concentrate near low-dimensional subsets of the ambient space. • Minimax rates. For targets supported on d~ d-dimensional regular sets, such as compact d~ d-dimensional manifolds, compact convex sets of dimension d~ d, and self-similar sets of similarity dimension d~ d, our bound gives (n−1/d)O(n^-1/d) for every d>d~d> d and every 1≤p<d~/21≤ p< d/2, matching the minimax rate for distribution estimation on such sets up to an arbitrarily small slack in the exponent. 2 Background and Problem Setting 2.1 Notations and Definitions Bϱ(x,r)B_ (x,r) denotes the open ball of radius r around x w.r.t. the metric ϱ . For a probability measure γ, supp(γ)=x:γ(Bϱ(x,r))>0, for all r>0supp(γ)=\x:γ(B_ (x,r))>0, for all r>0\. (⋅) 1(·) denotes the indicator function. We write An≲BnA_n B_n (equivalently An=(Bn)A_n=O(B_n)) if there exists C>0C>0, independent of n, such that An≤CBnA_n≤ CB_n; An≾BnA_n B_n (equivalently An=~(Bn)A_n= O(B_n)) if An≤Bn⋅logC(en)A_n≤ B_n· ^C(en) for some C>0C>0, i.e. up to poly-log factors in n; and An≍BnA_n B_n if both An≲BnA_n B_n and Bn≲AnB_n A_n hold. For k∈ℕk , [k]=1,…,k[k]=\1,…,k\. For a function of t, ξt _t, we write ξ˙t=dξt/dt ξ_t=d _t/dt. For a distribution μ, M(μ)=Law(Z∥Z∥≤M) T_M(μ)=Law(Z 1\\|Z\|≤ M\) where Z∼μZ μ. divx(⋅)div_x(·) denotes the divergence operator w.r.t. x. We also recall the following standard definitions. Definition 1 (Covering Number). For a metric space (,ϱ)(S, ), the ϵε-covering number with respect to (w.r.t.) ϱ is (ϵ;,ϱ)=infn∈ℕ:∃x1,…,xn such that ∪i=1nBϱ(xi,ϵ)⊇.N(ε;S, )= \n :∃\,x_1,…,x_n such that _i=1^nB_ (x_i,ε) \. Definition 2 (Neural networks). Let L∈ℕL and Nii∈[L]⊂ℕ\N_i\_i∈[L] . An L-layer neural network f:ℝD→ℝNLf: R ^D→ R ^N_L is f(x)=AL∘σL−1∘AL−1∘⋯∘σ1∘A1(x),f(x)=A_L _L-1 A_L-1 … _1 A_1(x), (1) where Ai(y)=Wiy+biA_i(y)=W_iy+b_i, with Wi∈ℝNi×Ni−1W_i∈ R ^N_i× N_i-1, bi∈ℝNib_i∈ R ^N_i, N0=DN_0=D, and each σj _j applied component-wise. Here Wi\W_i\ are weights, bi\b_i\ biases, and σi1≤i≤L−1\ _i\_1≤ i≤ L-1 activation functions. Define: depth ℒ(f):=LL(f):=L; number of weights (f)=∑i=1LNiNi−1W(f)= _i=1^LN_iN_i-1; maximum weight ℬ(f)=max1≤j≤L(‖bj‖∞∨‖Wj‖∞)B(f)= _1≤ j≤ L(\|b_j\|_∞ \|W_j\|_∞); and maximum output (f)=supx∈ℝD∥f(x)∥∞K(f)= _x∈ R ^D\|f(x)\|_∞. Let σi1≤i≤L−1(L,W,B,K)=f of the form (1):ℒ(f)≤L,(f)≤W,ℬ(f)≤B,(f)≤K. _\ _i\_1≤ i≤ L-1(L,W,B,K)= \f of the form e1:L(f)≤ L,W(f)≤ W,B(f)≤ B,K(f)≤ K \. When σj(x)=x∨0 _j(x)=x 0 for all j∈[L−1]j∈[L-1], we write ℛ(L,W,B,K)RN(L,W,B,K) for σi1≤i≤L−1(L,W,B,K)NN_\ _i\_1≤ i≤ L-1(L,W,B,K). Definition 3 (Wasserstein p-distance). For probability measures μ,νμ,ν on ℝD R ^D, the Wasserstein-p distance (under ℓ2 _2) is p(μ,ν)=(infγ∈Γ(μ,ν)(X,Y)∼γ‖X−Y‖2p)1/p,W_p(μ,ν)= ( _γ∈ (μ,ν)E_(X,Y) γ\|X-Y\|_2^p )^1/p, where Γ(μ,ν) (μ,ν) is the set of couplings of μ and ν. 2.2 Flow-matching The fundamental problem in generative modeling is to learn an unknown distribution PdataP_data on a data space ⊆ℝDX ^D from n i.i.d. observations x1,…,xn∼Pdatax_1,…,x_n P_data, where D can be extremely large—e.g., D=49,152D=49,152 for 128×128128× 128 color ImageNet images. Flow matching (Lipman et al., 2023) constructs a continuous probability path Ptt∈[0,1]\P_t\_t∈[0,1] from an easy-to-sample reference P0P_0 (typically standard Gaussian) to the target P1≡PdataP_1≡ P_data, governed by a time-dependent vector field v(⋅,⋅):ℝD×[0,1]→ℝDv(·,·): R ^D×[0,1]→ R ^D whose induced flow transports P0P_0 to PdataP_data: starting from X0∼P0X_0 P_0, X˙t=v(Xt,t),X0∼P0. X_t=v(X_t,t), X_0 P_0. (2) Since many velocity fields transport P0P_0 to P1P_1, a canonical choice arises from the linear interpolant I(x0,x1,t)=αtx1+βtx0I(x_0,x_1,t)= _tx_1+ _tx_0, with αt,βt∈[0,1] _t, _t∈[0,1] satisfying α0=β1=0 _0= _1=0, α1=β0=1 _1= _0=1; this gives I(x0,x1,0)=x0I(x_0,x_1,0)=x_0 and I(x0,x1,1)=x1I(x_0,x_1,1)=x_1. The process Xt0≤t≤1\X_t\_0≤ t≤ 1 constructed as Xt=I(X0,X1,t)X_t=I(X_0,X_1,t) for independent X0∼P0X_0 P_0, X1∼P1≡PdataX_1 P_1≡ P_data, with marginal law Pt=Law(Xt)P_t=Law(X_t). Since PdataP_data is unknown, it is replaced in practice by Pnemp=1n∑i=1nδxiP^emp_n= 1n _i=1^n _x_i, so X1∼PnempX_1 P^emp_n. Among velocity fields consistent with Pt\P_t\, flow matching targets the conditional expectation v⋆(x,t)=(X˙t∣Xt=x),v_ (x,t)=E ( X_t X_t=x ), satisfying the continuity equation ∂tϱt+divx(v⋆ϱt)=0 _t _t+div_x(v_ _t)=0, t∈[0,1)t∈[0,1), where ϱt _t is the density of PtP_t. Since v⋆(⋅,t)v_ (·,t) may explode as t→1t→ 1, the flow is instead run until t=1−δt=1-δ. Theorem 4 shows the resulting ODE X˙t=v⋆(Xt,t),X0∼P0,t∈[0,1−δ] X_t=v_ (X_t,t), X_0 P_0, t∈[0,1-δ] (3) is well-defined: for any initialization the solution exists and is unique, and XtX_t has density ϱt _t. Theorem 4. For all t∈[0,1−δ]t∈[0,1-δ], PtP_t admits a density, ϱt _t, w.r.t. the Lebesgue measure and satisfies the continuity equation (in the weak sense), ∂tϱt+divx(ϱtv⋆(⋅,t))=0, _t _t+div_x( _tv_ (·,t))=0, where v⋆(x,t)=βt−1β˙tx+βt−1(α˙tβt−β˙tαt)(∑i=1nexp(−‖x−αtxi‖22βt2))−1∑i=1nxiexp(−‖x−αtxi‖22βt2)v_ (x,t)= _t^-1 β_tx+ _t^-1( α_t _t- β_t _t) ( _i=1^n (- \|x- _tx_i\|^22 _t^2 ) )^-1 _i=1^nx_i (- \|x- _tx_i\|^22 _t^2 ). Let yx′(t)=v⋆(yx(t),t), with yx(0)=x,y_x (t)=v_ (y_x(t),t), with y_x(0)=x, (4) be the Initial Value Problem (IVP) with initial point x. Then (4) has a unique solution for every x∈ℝDx∈ R ^D, and Pt=(Yt)♯P0P_t=(Y_t)_ P_0, where Yt:x↦yx(t)Y_t:x y_x(t). Flow matching trains a neural network v:ℝD×[0,1]→ℝDv:R^D×[0,1] ^D to approximate v⋆v_ by minimizing ℒ(v)=∫01‖v(Xt,t)−v⋆(Xt,t)‖2t.L(v)= _0^1E\|v(X_t,t)-v_ (X_t,t)\|^2dt. A key observation of Lipman et al. (2023) is that minimizing ℒ(v)L(v) is equivalent (Lemma 11) to minimizing the conditional flow matching objective ℒCFM(v)=∫01‖v(Xt,t)−Xt˙‖2t=∫01‖v(Xt,t)−(α˙tX1+β˙tX0)‖2t,L_CFM(v)= _0^1E\|v(X_t,t)- X_t\|^2dt= _0^1E \|v(X_t,t)-( α_tX_1+ β_tX_0) \|^2dt, (5) which can be estimated from data without knowledge of v⋆v_ . Substituting PnempP^emp_n for PdataP_data gives the empirical objective, ℒCNF(v)=1n∑i=1n∫01[‖v(Xi,t,t)−(α˙txi+β˙tX0)‖2]t,L_CNF(v)= 1n _i=1^n _0^1E [\|v(X_i,t,t)-( α_tx_i+ β_tX_0)\|^2 ]dt, where Xi,t=αtxi+βtX0X_i,t= _tx_i+ _tX_0, X0∼P0X_0 P_0. Since this integral is not computed exactly in practice, we study its numerical approximation ℒ^n(v)=1n∑i=1n∑k=0N−1(tk+1−tk)‖v(Xi,tk,tk)−(α˙tkxi+β˙tkX0)‖2, L_n(v)= 1n _i=1^n _k=0^N-1(t_k+1-t_k)E\|v(X_i,t_k,t_k)-( α_t_kx_i+ β_t_kX_0)\|^2, (6) where 0=t0<t1<⋯<tN=1−δ0=t_0<t_1<…<t_N=1-δ partitions [0,1−δ][0,1-δ]; we use a uniform step size, tk+1=tk+ht_k+1=t_k+h for k∈[N−1]k∈[N-1]. The expectation in (6) is typically further approximated via Monte Carlo samples X0(j)∼(0,ID),X1(j)∼Unif(X1,…Xn),Xt(j)=αtX1(j)+βtX0(j),X_0^(j) (0,I_D),\,X_1^(j) (\X_1,… X_n\ ),\,X^(j)_t= _tX_1^(j)+ _tX_0^(j), giving the Monte Carlo estimate ℒ^n,mmc(v)=1m∑j=1m∑k=0N−1(tk+1−tk)‖v(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2. L^mc_n,m(v)= 1m _j=1^m _k=0^N-1(t_k+1-t_k) \|v(X_t_k^(j),t_k)-( α_t_kX_1^(j)+ β_t_kX_0^(j)) \|^2. (7) The resulting velocity field estimates are v^n=argminv∈ℒ^n(v)andv^n,mmc=argminv∈ℒ^n,mmc(v). v_n= argmin_v L_n(v) v_n,m^mc= argmin_v L_n,m^mc(v). (8) Stochastic Sampling After the velocity field is learned, one can generate new samples by following the ODE, dXt=v^(Xt,t)dt,X0∼(0,ID).dX_t= v(X_t,t)dt,\;X_0 (0,I_D). However, analyzing the above ODE is rather difficult, as stability of the flow relies critically on the regularity of the learned velocity field. In particular, standard perturbation arguments require v^(⋅,t) v(·,t) to be uniformly globally Lipschitz, and the error between trajectories generated by the learned and optimal velocity fields typically grows according to Grönwall’s inequality, leading to an exponential dependence on the Lipschitz constant. Since this Lipschitz constant may potentially blow up as t→1t→ 1, these bounds become increasingly loose near the terminal time, making it difficult to obtain meaningful control of the Wasserstein distance between the generated and target distributions. The existing literature artificially imposes a boundedness on the Lipschitz constant of the learned velocity field (Kumar et al., 2026; Gao et al., 2024). Alternatively, one can perform the generation through an SDE. This is because the Fokker-Plank equation corresponding to (3), ∂tϱt= _t _t= −divx(ϱtv⋆)=−divx(ϱt(v⋆+12σ2∇logϱt))+12σ2Δϱt, -div_x( _tv_ )=-div_x ( _t (v_ + 12σ^2∇ _t ) )+ 12σ^2 _t, -5.0pt which is satisfied by the SDE, dXt=(v⋆(Xt,t)+σ2∇logϱt(Xt))dt+σdWt.dX_t= (v_ (X_t,t)+ σ2∇ _t(X_t) )dt+ _t. A formal statement that ϱt _t satisfies the above SDE appears in Appendix A.3. Here, σ>0σ>0. Since, v⋆(x,t)=α˙tαtx+βt2(α˙tαt−β˙tβt)∇logϱt(x),v_ (x,t)= α_t _tx+ _t^2 ( α_t _t- β_t _t )∇ _t(x), One can generate samples by following the SDE, dXt=(atv⋆(Xt,t)−btXt)dt+σdWt,0≤t≤1−δ,dX_t= (a_tv_ (X_t,t)-b_tX_t )dt+σ\ dW_t, 0≤ t≤ 1-δ, -2.0pt (9) where, let at=1+(2βt(α˙tβt−β˙tαt))−1σ2αta_t=1+ (2 _t( α_t _t- β_t _t) )^-1σ^2 _t and bt=(2βt(α˙tβt−β˙tαt))−1σ2α˙tb_t= (2 _t( α_t _t- β_t _t) )^-1σ^2 α_t. After the velocity field is learned, one generates new samples through the Euler–Maruyama discretization of the SDE as, dX^t=(atkv^(X^tk,tk)−btX^t)dt+σdWt,tk≤t<tk+1, X_t= (a_t_k v( X_t_k,t_k)-b_t X_t )dt+σ\ dW_t, t_k≤ t<t_k+1, (10) where v v is either v^n v_n or v^n,mmc v_n,m^mc. The distribution of X^1−δ X_1-δ is denoted as P^1−δ(v^):=Law(X^1−δ) P_1-δ( v):=Law ( X_1-δ ). Letting B(s,t)=∫stbrrB(s,t)= _s^tb_rdr, the solution to (10) satisfies X^tk+1=e−B(tk,tk+1)X^tk+atkv^(X^tk,tk)∫tktk+1e−B(s,tk+1)s+σξk∫tktk+1e−2B(s,tk+1)s, X_t_k+1=e^-B(t_k,t_k+1) X_t_k+a_t_k v( X_t_k,t_k) _t_k^t_k+1e^-B(s,t_k+1)ds+σ\ _k _t_k^t_k+1e^-2B(s,t_k+1)ds, where ξk∼(0,ID) _k (0,I_D). Such a data generation procedure has been used by practitioners such as Singh and Fischer (2024); Ma et al. (2024), who report empirical superiority compared to deterministic sampling in flow-matching. The goal of this paper is to analyze the distance of P^1−δ(v^)=Law(X^1−δ) P_1-δ( v)=Law ( X_1-δ ), for both the estimators defined in (8). 2.3 Intrinsic Data Dimension Several notions of the effective dimension of a probability measure formalize the intuition of intrinsic dimension, most characterizing the growth rate of covering numbers on the logarithmic scale. Let (,ϱ)(S, ) be a Polish space endowed with the ℓ∞ _∞ metric and let μ be a probability measure on S. Recall the (ϵ,τ)(ε,τ)-covering number of μ (Posner et al., 1967), ϵ(μ,τ)=inf(ϵ;S,ϱ):S⊆,μ(S)≥1−τ, N_ε(μ,τ)= \N(ε;S, ):S ,\,μ(S)≥ 1-τ\, the minimal number of ϵε-balls needed to cover a set carrying at least 1−τ1-τ of the mass of μ. A natural measure of dimensionality is the upper Minkowski dimension dim¯M(μ)=lim supϵ↓0(log(1/ϵ))−1log(ϵ,supp(μ),ℓ∞), dim_M(μ)= _ε 0 ( (1/ε) )^-1 (ε;\,supp(μ),\, _∞), which depends only on the metric entropy of the support and so requires no smoothness or manifold structure, accommodating irregular geometries such as fractals. Its statistical implications are well established: covering numbers of smooth function classes are governed by the Minkowski dimension of their domain (Kolmogorov and Tikhomirov, 1961), and deep networks, WGANs and Wasserstein Autoencoders all attain rates reflecting it rather than the ambient dimension (Nakada and Imaizumi, 2020; Huang et al., 2022; Chakraborty and Bartlett, 2024). Being a property of the support alone, however, it ignores how mass is distributed and can be as large as D even when μ concentrates near a low-dimensional region (Chakraborty and Bartlett, 2025). Weed and Bach (2019) therefore refined Dudley’s (Dudley, 1969) entropic dimension into a notion sensitive to the mass distribution of μ, and showed that it governs the convergence rate of the empirical Wasserstein distance. Definition 5 (Wasserstein Dimension (Weed and Bach, 2019)). For p>0p>0, the p-upper Wasserstein dimension of μ is dp∗(μ)=infs>2p:lim supϵ↓0(log(1/ϵ))−1logϵ(μ,ϵsps−2p)≤s,d _p(μ)= \s>2p: _ε 0 ( (1/ε) )^-1 N_ε\! (μ,\,ε sps-2p )≤ s \, and the lower Wasserstein dimension is d∗(μ)=limτ↓0lim infϵ↓0logϵ(μ,τ)log(1/ϵ).d_ (μ)= _τ 0\, _ε 0 N_ε(μ,τ) (1/ε). Weed and Bach (2019) established that these jointly control the convergence rate of the empirical Wasserstein distance, n−1/d∗(μ)≲p(μ^n,μ)≲n−1/d∗p(μ)n^-1/d_ (μ) \,W_p( μ_n,μ) n^-1/d _p(μ). The upper Wasserstein dimension has recently been used to bound error rates for GANs (Chakraborty and Bartlett, 2025) and score-matching diffusion models (Chakraborty et al., 2026). 3 A Proof of Concept Experiment (a) Results on the basenji class (b) Results on the bubble class (c) Results on the goldfish class Figure 1: Average generalization error (in terms of FID scores) for different values of n for flow-matching. The error bars denote the standard deviation out of 1010 replications. Before presenting our theoretical analysis, we give an experiment illustrating that Flow Matching (FM) error rates depend primarily on the intrinsic dimension of the data, not the ambient dimension. We establish finite-sample generalization bounds using the Fokker-Planck equivalent SDE, bypassing artificial Lipschitz assumptions on the velocity field. Empirically, probability-flow ODE simulations confirm that our predicted intrinsic-dimensional error decay holds under standard deterministic sampling. Since intrinsic dimensionality is difficult to control in real image data, we follow the synthetic construction of Pope et al. (2020) and Chakraborty and Bartlett (2024): we use a pre-trained BigGAN (Brock et al., 2019) (128-dim latent space, 128×128×3128× 128× 3 output, trained on ImageNet (Deng et al., 2009)) to generate 11,00011,000 images each for the basenji, bubble, and goldfish classes, fixing all but d∈10,100d∈\10,100\ latent coordinates to zero so that generated data lie on a d-dimensional manifold. Images are downsampled to 28×2828× 28 via bilinear interpolation. For each d, we train a standard UNet FM model with the usual objective ‖vt−vθ(xt,t⋅1000)‖22E\|v_t-v_θ(x_t,t· 1000)\|_2^2, varying the training set size over 1000,2000,…,10000\1000,2000,…,10000\ (Adam, lr 10−410^-4, batch size 6464, 5050 epochs; NVIDIA RTX 6000 Blackwell GPUs). At evaluation, we simulate the learned probability-flow ODE (100100 steps) to draw 1,0001,000 samples and compute FID (Heusel et al., 2017) against a held-out set of equal size (indices 10,00010,000–11,00011,000), averaged over 1010 repetitions. The empirical results, presented in Fig. 1, clearly validate our hypothesis. As the number of training samples increases from 10001000 to 1000010000, the FM model trained on data with intrinsic dimension d=10d=10 consistently achieves substantially lower FID scores than the model trained on data with intrinsic dimension d=100d=100. The gap persists across large sample sizes, and in fact widens slightly as more data become available, indicating that lower intrinsic-dimensional structure yields faster error decay under flow-based generative learning. This behavior aligns precisely with our theoretical predictions: the intrinsic dimension, rather than the ambient pixel dimension, governs the sample complexity and achievable error rates of flow-matching generative models. All code pertaining to the experiment is publicly available through https://github.com/saptarshic27/fmintrinsic. 4 Theoretical Analyses 4.1 Main Result and Inference To begin our theoretical analyses, we make the following standard assumption (Huang et al., 2022; Oko et al., 2023; Chakraborty and Bartlett, 2025) that the data is independent and identically distributed and supported on the unit hypercube, [0,1]D[0,1]^D. Assumption 1. We assume that x1,…,xnx_1,…,x_n are independent and identically distributed according to the distribution PdataP_data, such that Pdata([0,1]D)=1P_data ([0,1]^D )=1. Our approach relaxes several restrictive assumptions found in recent literature. Unlike Kumar et al. (2026), we do not require the existence of densities, smoothness, or a support characterized as a smooth manifold with bounded reach. Additionally, we avoid the uniform score-Lipschitzness and log-concavity of the target measure imposed by Gao et al. (2024) and Zhou and Liu (2025). Further, unlike Fukumizu et al. (2025), we do not require that PdataP_data to have a lower-bounded smooth density. Design Choices To ensure that the interpolating path is regular enough, we impose the following design choices on αt _t, and βt _t. Design choices D 1(a)–(c) are satisfied by the commonly used interpolation schedules, (αt=tCLOSE( _t=t, OPENβt=1−t) _t=1-t), (αt=tCLOSE( _t=t, OPENβt=1−t2) _t= 1-t^2). D 1. We assume that for any t,s∈[0,1−δ]t,s∈[0,1-δ], the following hold: (a) |β˙tβt|,|α˙tβt−β˙tαtβt|,at2,|bt|,βt−1 | β_t _t |,\ | α_t _t- β_t _t _t |,\ a_t^2,\ |b_t|,\ _t^-1, |β¨tβt−(β˙tβt)2||βt−4(3αtβtα˙tβ˙t−βt2α˙t2− 2αt2β˙t2)| | β_t _t- ( β_t _t )^2 |\; | _t^-4 (3 _t _t α_t β_t\;-\; _t^2 α_t^2\;-\;2 _t^2 β_t^2 ) |, |α¨t−αtβ¨tβt−1−α˙tβ˙tβt−1+αtβ˙t2βt−2| | α_t- _t β_t\, _t^-1- α_t β_t\, _t^-1+ _t β_t^2\, _t^-2 |, |αt(α˙tβt−β˙tαt)2βt−4|≲δ−κ1 | _t ( α_t _t- β_t _t )^2\, _t^-4 | δ^- _1. (b) |αt−αs|,|βt−βs|,|at−as|≲δ−κ2|t−s|| _t- _s|,\ | _t- _s|,\,|a_t-a_s| δ^- _2|t-s|. (c) 1−α1−δ,β1−δ≲δκ31- _1-δ,\ _1-δ δ _3. Note that for the linear interpolation, (αt=tCLOSE( _t=t, OPENβt=1−t) _t=1-t), it suffices the above hold with κ1=4 _1=4, κ2=2 _2=2 and κ3=1 _3=1, where as for the Variance Preserving (VP) choice (αt=tCLOSE( _t=t, OPENβt=1−t2) _t= 1-t^2), the above hold with κ1=3 _1=3, κ2=0.5 _2=0.5 and κ3=0.5 _3=0.5. Under Assumption 1, we establish high-probability upper bounds on the Wasserstein-p distance between the flow-matching estimators P^1−δ(v^n) P_1-δ( v_n) and P^1−δ(v^n,mmc) P_1-δ( v^mc_n,m) and the target distribution PdataP_data. These bounds depend on the training sample size through a rate exponent governed solely by the upper p-Wasserstein dimension of PdataP_data, rather than the ambient dimension D. Consequently, with data-driven choices of the discretization step-size h, early stopping time δ, and velocity-field function class V, minimizing the flow-matching MSE loss yields an estimator that automatically adapts to the low-dimensional structure of PdataP_data. Theorem 6 derives a high-probability bound in terms of the Wasserstein-p distance. Theorem 6 (Main Theorem). Suppose that d>dp∗(Pdata)d>d _p (P_data ) and Assumption 1 holds. Then under Design Choices 1, for any p≥1p≥ 1, there exists a constant n0∈ℕn_0 , such that if n≥n0n≥ n_0, and ≡ℛ(W,L,B,K)V (W,L,B,K ) with δ≍n−1dκ3δ n^- 1d _3, h≍n−2p(2p+1)d−10κ1+2κ2dκ3h n^- 2p(2p+1)d- 10 _1+2 _2d _3, M≍n2pdM n 2pd, then, with probability at least 1−ξ1-ξ, p(Pdata,M(P^1−δ(v^n)))≲n−1/d+n−12plog12p(1/ξ). _p (P_data, T_M ( P_1-δ( v_n) ) ) n^-1/d+n^- 12p 12p(1/ξ). (11) Here, W≍np(2p+1)(D+2)d+(20+7D)κ1+4κ22dκ3logD/2+1nW n p(2p+1)(D+2)d\,+\, (20+7D) _1+4 _22d _3 ^D/2+1n, L≍lognL n, B≍n(2p(2p+1)d+20κ1+4κ22dκ3)∨(3κ1)∨1dκ3lognB n ( 2p(2p+1)d+ 20 _1+4 _22d _3 )\, \, (3 _1) 1d _3 n, and K≍n(3κ1)∨1dκ3lognK n (3 _1) 1d _3 n. Further, if m≍np(2p+1)(D+6)d+(24+7D)κ1+4κ2+4((6κ1)∨2)2dκ3log2n(logD/2+3n+log(1/ξ)),m n p(2p+1)(D+6)d+ (24+7D) _1+4 _2+4 ((6 _1) 2 )2d _3 ^2n ( ^D/2+3n+ (1/ξ) ), with probability at least 1−ξ1-ξ, p(Pdata,M(P^1−δ(v^n,mmc)))≲n−1/d+n−12plog12p(1/ξ). _p (P_data, T_M ( P_1-δ( v^mc_n,m) ) ) n^-1/d+n^- 12p 12p(1/ξ). (12) Inference for Data Supported on a Manifold Recall that a set A is d~ d-regular w.r.t. the d~ d-dimensional Hausdorff measure ℍd~H d if ℍd~(Bϱ(x,r))≍rd~H d(B_ (x,r)) r d for all x∈x∈ A (Weed and Bach, 2019, Definition 6). By Weed and Bach (2019, Proposition 8), if supp(Pdata)supp(P_data) is d~ d-regular and Pdata≪ℍd~P_data d, then dp∗(Pdata)=d~d_p (P_data)= d for every p<d~/2p< d/2; in particular this holds for every compact d~ d-dimensional differentiable manifold, nonempty compact convex set of dimension d~ d, its relative boundary, and self-similar sets of similarity dimension d~ d. For such supports and any 1≤p<d~/21≤ p< d/2, Theorem 6 applies with every d>d~d> d, giving (n−1/d+n−1/2plog1/2p(1/ξ))O (n^-1/d+n^-1/2p ^1/2p(1/ξ) ); as d↓d~d d this matches the minimax rate for distribution estimation on d~ d-dimensional manifolds up to the slack d−d~d- d (Chakraborty, 2025), and improves on Kumar et al. (2026), whose comparable guarantee is not minimax optimal and additionally requires bounded reach and Hölder-smooth densities. Comparison with Existing Bounds Unlike Gao et al. (2024); Zhou and Liu (2025), whose (n−1/(D+5))O(n^-1/(D+5))-type rates depend on the log-concavity or Lipschitzness of the score/velocity field and scale with the ambient dimension D, and unlike Fukumizu et al. (2025), who require Besov smoothness of the density, our bound in Theorem 6 holds for every p≥1p≥ 1, is non-asymptotic and high-probability rather than in-expectation, and imposes no density, smoothness, or Lipschitzness condition on PdataP_data or the velocity class. The claim concerns the exponent of n; the prescribed architecture still depends on D, and more than polynomially (see below). This places flow matching on the same footing as the dimension-adaptive guarantees recently established for GANs (Chakraborty and Bartlett, 2025) and diffusion models (Chakraborty et al., 2026), while relaxing the compact-support and manifold-structure assumptions common to both lines of work. Choice of Early Stopping and Discretization Theorem 6 sets the early-stopping offset δ≍n−1/(dκ3)δ n^-1/(d _3), which by D 1(c) keeps the truncation error p(P1−δ,Pdata)=(δκ3)W_p(P_1-δ,P_data)=O(δ _3) on the same order as the sampling error, while preventing v⋆v_ from blowing up as t→1−t→ 1^-. Because the sampler is an SDE with a nondegenerate diffusion coefficient, its increments over a step are of mean square order h rather than h2h^2, and the cumulative Euler–Maruyama error is (δ−(10κ1+2κ2)h)O(δ^-(10 _1+2 _2)h) (Lemma 8); matching this to the sampling error gives the step size h≍n−2p(2p+1)d−10κ1+2κ2dκ3h n^- 2p(2p+1)d- 10 _1+2 _2d _3 above, and a number of partitions N≍h−1N h^-1 that still grows only polynomially in n. The proof also fixes internally the approximation accuracy ϵε and the excursion probability η, both polynomial in n; the latter sets the radius R≍Dlog(1/η)≍DlognR D (1/η) D n of the box on which v⋆v_ is approximated. Network and Monte Carlo Scaling The velocity network width W≍np(2p+1)(D+2)d+(20+7D)κ1+4κ22dκ3logD/2+1nW n p(2p+1)(D+2)d+ (20+7D) _1+4 _22d _3 ^D/2+1n is polynomial in n for each fixed D, with an exponent that grows linearly in the ambient dimension D and shrinks as the intrinsic dimension d increases, while depth L≍lognL n suffices; the caps B and K carry an extra DlognD n, the magnitude v⋆v_ attains on the truncation box. The polylogarithmic factors reflect the cost of covering the box [−R,R]D[-R,R]^D on which v⋆v_ is approximated, and do not affect the exponent of n, which is what governs the rate. The Monte Carlo sample size m≍np(2p+1)(D+6)d+(24+7D)κ1+4κ2+4((6κ1)∨2)2dκ3log2n(logD/2+3n+log(1/ξ))m n p(2p+1)(D+6)d+ (24+7D) _1+4 _2+4 ((6 _1) 2 )2d _3 ^2n ( ^D/2+3n+ (1/ξ) ) inherits this (D,d)(D,d)-dependence, inflated by the covering number of V and by K4K^4 from the sub-exponential concentration in Lemma 14; larger D/dD/d ratios thus demand proportionally more generated samples, mirroring the scaling seen for GANs (Chakraborty and Bartlett, 2025; Huang et al., 2022) and diffusion models (Chakraborty et al., 2026). 4.2 Proof Sketch Error Decomposition To derive an upper bound on the excess risk, we first decompose the KL estimation error into a statistical error due to velocity field estimation and a numerical error resulting from the forward Euler-Maruyama discretization of the SDE. This decomposotion is shown in Lemma 7 with proof appearing in Appendix B.1. Through out the remainder of the paper, all expectations are taken conditioned on the data i.e., (⋅∣x1,…,xn)E(· x_1,…,x_n). Lemma 7. Under Assumption 1 and Design Choices D 1, for any estimator v^∈ v of v⋆v_ , KL(P1−δ∥P^1−δ(v^))≤2σ2∑k=0N−1∫tktk+1at2∥v⋆(Xt,t)−v⋆(Xtk,tk)∥2dt (P_1-δ\,\|\, P_1-δ( v) )≤\; 2σ^2 _k=0^N-1 _t_k^t_k+1a_t^2\,E \|v_ (X_t,t)-v_ (X_t_k,t_k) \|^2dt +2σ2∑k=0N−1∫tktk+1(at−atk)2∥v⋆(Xtk,tk)∥2dt+1σ2∑k=0N−1hkatk2∥v⋆(Xtk,tk)−v^(Xtk,tk)∥2, + 2σ^2 _k=0^N-1 _t_k^t_k+1(a_t-a_t_k)^2\,E \|v_ (X_t_k,t_k) \|^2dt+ 1σ^2 _k=0^N-1h_k\,a_t_k^2\,E \|v_ (X_t_k,t_k)- v(X_t_k,t_k) \|^2, where every expectation is taken under the law of the process Xtt∈[0,1−δ]\X_t\_t∈[0,1-δ] solving (9). Using Lemma 7, we aim to bound the terms individually in the subsequent paragraphs. Discretization Error The first in Lemma 7, ∑k=0N−1∫tktk+1‖atv⋆(Xt,t)−atkv⋆(Xtk,tk)‖2t _k=0^N-1 _t_k^t_k+1E\|a_tv_ (X_t,t)-a_t_kv_ (X_t_k,t_k)\|^2dt captures the discretization error incurred when the continuous-time SDE generated by the velocity field v⋆v_ is approximated using the forward Euler-Maruyama scheme. This term measures the cumulative discrepancy between the velocity evaluated at the left endpoint of each interval and the true velocity along the continuous trajectory. Consequently, it quantifies the bias introduced by replacing the exact SDE dynamics with piecewise-constant updates over the partition tjj=0N\t_j\_j=0^N. The following lemma shows that this error decreases as the uniform discretization step size h tends to zero. Lemma 8. Under Assumption 1 and Design Choices D 1, with uniform step size h≤1h≤ 1, the following hold, all expectations being taken under the law of the process Xt\X_t\ solving (9): ∑k=0N−1∫tktk+1at2‖v⋆(Xt,t)−v⋆(Xtk,tk)‖2t≲D3(σ2∨1)δ−10κ1h, _k=0^N-1 _t_k^t_k+1a_t^2\,E \|v_ (X_t,t)-v_ (X_t_k,t_k) \|^2dt\; \;D^3 (σ^2 1 )\,δ^-10 _1\,h, (13) ∑k=0N−1∫tktk+1(at−atk)2‖v⋆(Xtk,tk)‖2t≲Dδ−2(κ1+κ2)h2. _k=0^N-1 _t_k^t_k+1(a_t-a_t_k)^2\,E \|v_ (X_t_k,t_k) \|^2dt\; \;D\,δ^-2( _1+ _2)\,h^2. (14) In particular, writing c1:=10κ1+2κ2c_1:=10 _1+2 _2, the sum of the two is ≲D3(σ2∨1)δ−c1h D^3(σ^2 1)\,δ^-c_1h. Approximation Error A key component in controlling the first term in Lemma 7 is the approximation error of the function class V used to approximate v⋆v_ . Neural network approximation theory has been studied extensively: Cybenko (1989) and Hornik (1991) established universal approximation for sigmoid-type single-hidden-layer networks on compact domains, and more recent work has characterized the expressive power and approximation rates of deep networks (Yarotsky, 2017; Petersen and Voigtlaender, 2018; Shen et al., 2022; Nakada and Imaizumi, 2020; Schmidt-Hieber, 2020). Lemma 9. For any ϵ∈(0,1)ε∈(0,1) and any R≥1R≥ 1, there exists a network v0∈ℛ(W,L,B,K)v_0 (W,L,B,K) with supx∈[−R,R]D‖v0(x,ti)−v⋆(x,ti)‖≤ϵ _x∈[-R,R]^D\|v_0(x,t_i)-v_ (x,t_i)\|≤ε for all i=0,…,Ni=0,…,N, and (v0)≲NRDδ−3κ1Dϵ−D(log(1/δ)+log(1/ϵ)+logR),W(v_0) NR^Dδ^-3 _1Dε^-D ( (1/δ)+ (1/ε)+ R ), ℒ(v0)≲log(1/δ)+log(1/ϵ)+logRL(v_0) (1/δ)+ (1/ε)+ R, ℬ(v0)≲δ−(3κ1)∨1RD+N,B(v_0) \;δ^-(3 _1) 1R D+N, and (v0)≲δ−(3κ1)∨1RD.K(v_0) δ^-(3 _1) 1R D. Further, v0(x,t)=0v_0(x,t)=0 if x∉[−(R+1),R+1]Dx ∈[-(R+1),R+1]^D, and Lip(v0(⋅,t))≲δ−3κ1RLip(v_0(·,t)) δ^-3 _1R for all t∈[0,1−δ]t∈[0,1-δ], where the implied constants depend only on D. The proof (detailed in Appendix C.1) then proceeds to control the first term in Lemma 7. For the estimator v^n v_n, this term can be upper bounded by a function of the approximation error, as shown in Lemma 9. For the estimator v^n,mmc v^mc_n,m, controlling this term becomes more nuanced, as this involves controlling both an approximation error and a generalization error that accounts for estimating the expectation via Monte-Carlo in Equation (7) as shown in Lemma 14. Finally, we use the bound on KL(P1−δ,P^1−δ(v^))KL (P_1-δ, P_1-δ( v) ) to bound p(P1−δ,P^1−δ(v^))W_p (P_1-δ, P_1-δ( v) ) via truncation as shown in Lemma 21. 5 Conclusion and Discussions In this work, we developed a theoretical understanding of flow matching for data distributions whose complexity is substantially lower than their ambient dimension. We established the first intrinsic-dimension-adaptive, high-probability generalization bound for flow matching in Wasserstein-p distance, for any p≥1p≥ 1: with appropriate choices of network architecture, early stopping, discretization scheme and truncation level, given n i.i.d. samples from PdataP_data, and for every d>dp∗(Pdata)d>d _p(P_data), with probability at least 1−ξ1-ξ, p(Pdata,M(P^1−δ(v^n)))≲n−1/d+n−12plog12p(1/ξ),W_p (P_data, T_M( P_1-δ( v_n)) )\; \;n^-1/d+n^- 12p 12p(1/ξ), where dp∗(⋅)d _p(·) is the upper Wasserstein-p dimension (Weed and Bach, 2019) and M T_M truncates at level M≍n2p/dM n^2p/d. Since d may be taken arbitrarily close to dp∗(Pdata)d _p(P_data), the exponent is governed by the intrinsic dimension up to an arbitrarily small slack. Our analysis dispenses with the uniform Lipschitz assumption on the learned velocity field standard in state-of-the-art analyses—the regularity of v⋆v_ itself is established rather than assumed—and prescribes network sizes polynomial in n for each fixed D. The resulting rates avoid the ambient-dimensional curse of dimensionality and, for intrinsically low-dimensional distributions—including those on compact smooth manifolds—approach the corresponding minimax rates, giving a theoretical basis for the empirically observed effectiveness of flow matching in high dimensions with low-dimensional structure. We also address several algorithmic considerations—discretization, early stopping, and truncation of the generated law—with theoretically motivated prescriptions. Our bounds imply an early stopping time δ=Θ(n−1dκ3)δ= (n^- 1d _3), ensuring the learned velocity field stays close to the ground truth while avoiding blow-up as t→1−t→ 1^-, and a truncation level M=Θ(n2p/d)M= (n^2p/d), balancing the tail of P1−δP_1-δ against the transport cost it induces. For numerical integration on [0,1−δ][0,1-δ], we recommend a uniform step size h=Θ(n−2p(2p+1)d−10κ1+2κ2dκ3)h= (n^- 2p(2p+1)d- 10 _1+2 _2d _3 ). This is dictated by the Brownian component of the sampler: over a step of length h the increments of the generating SDE have mean square of order h rather than h2h^2, so the accumulated Euler–Maruyama error is (δ−(10κ1+2κ2)h)O(δ^-(10 _1+2 _2)h), and the step size above is what places it on the same order as the velocity estimation error. For the velocity network, ReLU feed-forward architectures suffice, with depth L≍lognL n and W≍np(2p+1)(D+2)d+(20+7D)κ1+4κ22dκ3logD/2+1nW n^\, p(2p+1)(D+2)d\,+\, (20+7D) _1+4 _22d _3 ^D/2+1n weights. The exponent of n here is linear in D and decreasing in d. While our results give detailed statistical guarantees, the optimization error from training nonconvex, time-dependent velocity networks remains hard to control; our bounds are optimization-agnostic and can be combined with future optimization-theoretic advances. Two further directions remain open. First, our approximation analysis uses only the Lipschitz regularity of v⋆v_ ; exploiting its Gaussian-mixture structure, or additional smoothness, should yield smaller networks. Second, we assume a ReLU velocity network—standard in the statistical learning literature and compatible with existing approximation results for Hölder and Besov-type classes—whereas practical models use structured architectures such as U-Nets or transformers. These often exhibit strong empirical smoothness and hierarchical approximation properties, but rigorous guarantees for them remain unavailable. Acknowledgements We gratefully acknowledge the support of the NSF and the Simons Foundation for the Collaboration on the Theoretical Foundations of Deep Learning through awards DMS-2031883 and #814639, the NSF’s support of FODSI through grant DMS-2023505, and the ONR’s support through MURI award N000142112431. References Albergo et al. (2025) Michael Albergo, Nicholas M. Boffi, and Eric Vanden-Eijnden. Stochastic interpolants: A unifying framework for flows and diffusions. Journal of Machine Learning Research, 26(209):1–80, 2025. URL http://jmlr.org/papers/v26/23-1605.html. Albergo and Vanden-Eijnden (2023) Michael Samuel Albergo and Eric Vanden-Eijnden. Building normalizing flows with stochastic interpolants. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=li7qeBbCR1t. Benton et al. (2024) Joe Benton, George Deligiannidis, and Arnaud Doucet. Error bounds for flow matching methods. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL https://openreview.net/forum?id=uqQPyWFDhY. Bickel and Li (2007) Peter J Bickel and Bo Li. Local polynomial regression on unknown manifolds. Lecture Notes-Monograph Series, pages 177–186, 2007. Block et al. (2020) Adam Block, Youssef Mroueh, and Alexander Rakhlin. Generative modeling with denoising auto-encoders and langevin sampling. arXiv preprint arXiv:2002.00107, 2020. Bogachev et al. (2015) Vladimir I. Bogachev, Nikolai V. Krylov, Michael Röckner, and Stanislav V. Shaposhnikov. Fokker–Planck–Kolmogorov Equations, volume 207 of Mathematical Surveys and Monographs. American Mathematical Society, Providence, RI, 2015. Bortoli (2022) Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022. ISSN 2835-8856. URL https://openreview.net/forum?id=MhK5aXo3gB. Expert Certification. Bose et al. (2024) Joey Bose, Tara Akhound-Sadegh, Guillaume Huguet, Kilian FATRAS, Jarrid Rector-Brooks, Cheng-Hao Liu, Andrei Cristian Nica, Maksym Korablyov, Michael M. Bronstein, and Alexander Tong. SE(3)-stochastic flow matching for protein backbone generation. In The Twelfth International Conference on Learning Representations, 2024. URL https://openreview.net/forum?id=kJFIH23hXb. Brock et al. (2019) Andrew Brock, Jeff Donahue, and Karen Simonyan. Large scale GAN training for high fidelity natural image synthesis. In International Conference on Learning Representations, 2019. URL https://openreview.net/forum?id=B1xsqj09Fm. Bunne et al. (2022) Charlotte Bunne, Andreas Krause, and marco cuturi. Supervised training of conditional monge maps. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=sPNtVVUq7wi. Chakraborty (2025) Saptarshi Chakraborty. Minimax lower bounds for estimating distributions on low-dimensional spaces. Transactions on Machine Learning Research, 2025. ISSN 2835-8856. URL https://openreview.net/forum?id=wIgRV336hC. Chakraborty and Bartlett (2024) Saptarshi Chakraborty and Peter Bartlett. A statistical analysis of wasserstein autoencoders for intrinsically low-dimensional data. In The Twelfth International Conference on Learning Representations, 2024. URL https://openreview.net/forum?id=WjRPZsfeBO. Chakraborty and Bartlett (2025) Saptarshi Chakraborty and Peter L. Bartlett. On the statistical properties of generative adversarial models for low intrinsic data dimension. Journal of Machine Learning Research, 26(111):1–57, 2025. URL http://jmlr.org/papers/v26/24-0054.html. Chakraborty et al. (2026) Saptarshi Chakraborty, Quentin Berthet, and Peter L. Bartlett. Generalization properties of score-matching diffusion models for intrinsically low-dimensional data, 2026. URL https://arxiv.org/abs/2603.03700. Chen et al. (2026) Fan Chen, Sinho Chewi, Constantinos Costis Daskalakis, and Alexander Rakhlin. High-accuracy sampling for diffusion models and log-concave distributions. In Forty-third International Conference on Machine Learning, 2026. URL https://openreview.net/forum?id=GW3umRqsZZ. Chen et al. (2023a) Hongrui Chen, Holden Lee, and Jianfeng Lu. Improved analysis of score-based generative modeling: User-friendly bounds under minimal smoothness assumptions. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 4735–4763. PMLR, 23–29 Jul 2023a. URL https://proceedings.mlr.press/v202/chen23q.html. Chen et al. (2023b) Minshuo Chen, Kaixuan Huang, Tuo Zhao, and Mengdi Wang. Score approximation, estimation and distribution recovery of diffusion models on low-dimensional data. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 4672–4712. PMLR, 23–29 Jul 2023b. URL https://proceedings.mlr.press/v202/chen23o.html. Chen et al. (2018) Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL https://proceedings.neurips.c/paper_files/paper/2018/file/69386f6b1dfed68692a24c8686939b9-Paper.pdf. Chen et al. (2023c) Sitan Chen, Sinho Chewi, Holden Lee, Yuanzhi Li, Jianfeng Lu, and Adil Salim. The probability flow ODE is provably fast. In Thirty-seventh Conference on Neural Information Processing Systems, 2023c. URL https://openreview.net/forum?id=KD6MFeWSAd. Chen et al. (2023d) Sitan Chen, Giannis Daras, and Alex Dimakis. Restoration-degradation beyond linear diffusions: A non-asymptotic analysis for DDIM-type samplers. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 4462–4484. PMLR, 23–29 Jul 2023d. URL https://proceedings.mlr.press/v202/chen23e.html. Coddington and Levinson (1955) Earl A. Coddington and Norman Levinson. Theory of Ordinary Differential Equations. McGraw-Hill, New York, 1955. ISBN 978-0070115422. Corso et al. (2023) Gabriele Corso, Hannes Stärk, Bowen Jing, Regina Barzilay, and Tommi S. Jaakkola. Diffdock: Diffusion steps, twists, and turns for molecular docking. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=kKF8_K-mBbS. Cybenko (1989) George Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of control, signals and systems, 2(4):303–314, 1989. Dahal et al. (2022) Biraj Dahal, Alexander Havrilla, Minshuo Chen, Tuo Zhao, and Wenjing Liao. On deep generative models for approximation and estimation of distributions on manifolds. Advances in Neural Information Processing Systems, 35:10615–10628, 2022. Deng et al. (2009) Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large-scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255, 2009. doi: 10.1109/CVPR.2009.5206848. Dudley (1969) Richard Mansfield Dudley. The speed of mean glivenko-cantelli convergence. The Annals of Mathematical Statistics, 40(1):40–50, 1969. Esser et al. (2024) Patrick Esser, Sumith Kulal, Andreas Blattmann, Rahim Entezari, Jonas Müller, Harry Saini, Yam Levi, Dominik Lorenz, Axel Sauer, Frederic Boesel, Dustin Podell, Tim Dockhorn, Zion English, and Robin Rombach. Scaling rectified flow transformers for high-resolution image synthesis. In Forty-first International Conference on Machine Learning, 2024. URL https://openreview.net/forum?id=FPnUhsQJ5B. Fukumizu et al. (2025) Kenji Fukumizu, Taiji Suzuki, Noboru Isobe, Kazusato Oko, and Masanori Koyama. Flow matching achieves almost minimax optimal convergence. In The Thirteenth International Conference on Learning Representations, 2025. URL https://openreview.net/forum?id=2OMyAFjiJJ. Gao et al. (2024) Yuan Gao, Jian Huang, Yuling Jiao, and Shurong Zheng. Convergence of continuous normalizing flows for learning probability distributions. arXiv preprint arXiv:2404.00551, 2024. Geffner et al. (2025) Tomas Geffner, Kieran Didi, Zuobai Zhang, Danny Reidenbach, Zhonglin Cao, Jason Yim, Mario Geiger, Christian Dallago, Emine Kucukbenli, Arash Vahdat, and Karsten Kreis. Proteina: Scaling flow-based protein structure generative models. In The Thirteenth International Conference on Learning Representations, 2025. URL https://openreview.net/forum?id=TVQLu34bdw. Goodfellow et al. (2020) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial networks. Communications of the ACM, 63(11):139–144, 2020. Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. Gans trained by a two time-scale update rule converge to a local nash equilibrium. Advances in Neural Information Processing Systems, 30, 2017. Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 6840–6851. Curran Associates, Inc., 2020. URL https://proceedings.neurips.c/paper_files/paper/2020/file/4c5bcfec8584af0d967f1ab10179ca4b-Paper.pdf. Hornik (1991) Kurt Hornik. Approximation capabilities of multilayer feedforward networks. Neural networks, 4(2):251–257, 1991. Huang et al. (2022) Jian Huang, Yuling Jiao, Zhen Li, Shiao Liu, Yang Wang, and Yunfei Yang. An error analysis of generative adversarial networks for learning distributions. Journal of Machine Learning Research, 23(116):1–43, 2022. URL http://jmlr.org/papers/v23/21-0732.html. Karatzas and Shreve (1991) Ioannis Karatzas and Steven E. Shreve. Brownian Motion and Stochastic Calculus. Springer, Berlin, 2nd edition, 1991. Kim et al. (2019) Jisu Kim, Jaehyeok Shin, Alessandro Rinaldo, and Larry Wasserman. Uniform convergence rate of the kernel density estimator adaptive to intrinsic volume dimension. In International Conference on Machine Learning, pages 3398–3407. PMLR, 2019. Kingma and Welling (2014) Diederik P Kingma and Max Welling. Auto-encoding variational bayes. In International Conference on Learning Representations (ICLR), 2014. Kolmogorov and Tikhomirov (1961) Andrey N Kolmogorov and Vladimir Mikhaılovich Tikhomirov. ϵε-entropy and ϵε-capacity of sets in function spaces. Translations of the American Mathematical Society, 17:277–364, 1961. Korotin et al. (2023) Alexander Korotin, Daniil Selikhanovych, and Evgeny Burnaev. Neural optimal transport. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=d8CBRlWNkqH. Kumar et al. (2026) Shivam Kumar, Yixin Wang, and Lizhen Lin. Flow matching is adaptive to manifold structures, 2026. URL https://arxiv.org/abs/2602.22486. Lee et al. (2022) Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence for score-based generative modeling with polynomial complexity. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=dUSI4vFyMK. Lee et al. (2023) Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence of score-based generative modeling for general data distributions. In Shipra Agrawal and Francesco Orabona, editors, Proceedings of The 34th International Conference on Algorithmic Learning Theory, volume 201 of Proceedings of Machine Learning Research, pages 946–985. PMLR, 20 Feb–23 Feb 2023. URL https://proceedings.mlr.press/v201/lee23a.html. Lipman et al. (2023) Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. In 11th International Conference on Learning Representations, ICLR 2023, 2023. Liu et al. (2023) Xingchao Liu, Chengyue Gong, and qiang liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=XVjTT1nw5z. Lu et al. (2022) Cheng Lu, Kaiwen Zheng, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Maximum likelihood training for score-based diffusion ODEs by high order denoising score matching. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 14429–14460. PMLR, 17–23 Jul 2022. URL https://proceedings.mlr.press/v162/lu22f.html. Ma et al. (2024) Nanye Ma, Mark Goldstein, Michael S Albergo, Nicholas M Boffi, Eric Vanden-Eijnden, and Saining Xie. Sit: Exploring flow and diffusion-based generative models with scalable interpolant transformers. In European Conference on Computer Vision, pages 23–40. Springer, 2024. Nakada and Imaizumi (2020) Ryumei Nakada and Masaaki Imaizumi. Adaptive approximation and generalization of deep neural network with intrinsic dimensionality. Journal of Machine Learning Research, 21(174):1–38, 2020. URL http://jmlr.org/papers/v21/20-002.html. Oko et al. (2023) Kazusato Oko, Shunta Akiyama, and Taiji Suzuki. Diffusion models are minimax optimal distribution estimators. In International Conference on Machine Learning, pages 26517–26582. PMLR, 2023. Øksendal (2003) Bernt Øksendal. Stochastic differential equations. In Stochastic differential equations: an introduction with applications, pages 38–50. Springer, 2003. Petersen and Voigtlaender (2018) Philipp Petersen and Felix Voigtlaender. Optimal approximation of piecewise smooth functions using deep relu neural networks. Neural Networks, 108:296–330, 2018. Pi et al. (2026) Sophia Pi, Mingcheng Lu, Maojiang Su, Weimin Wu, Jerry Yao-Chieh Hu, and Han Liu. Learning manifold data with flow matching. In Forty-third International Conference on Machine Learning, 2026. URL https://openreview.net/forum?id=9mA5d705Gd. Pope et al. (2020) 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 Representations, 2020. Posner et al. (1967) Edward C Posner, Eugene R Rodemich, and Howard Rumsey Jr. Epsilon entropy of stochastic processes. The Annals of Mathematical Statistics, pages 1000–1020, 1967. Rombach et al. (2022) Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Bjorn Ommer. High-Resolution Image Synthesis with Latent Diffusion Models . In 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 10674–10685, Los Alamitos, CA, USA, June 2022. IEEE Computer Society. doi: 10.1109/CVPR52688.2022.01042. URL https://doi.ieeecomputersociety.org/10.1109/CVPR52688.2022.01042. Saharia et al. (2022) Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Raphael Gontijo-Lopes, Burcu Karagol Ayan, Tim Salimans, Jonathan Ho, David J. Fleet, and Mohammad Norouzi. Photorealistic text-to-image diffusion models with deep language understanding. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=08Yk-n5l2Al. Santambrogio (2015) Filippo Santambrogio. Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling, volume 87 of Progress in Nonlinear Differential Equations and Their Applications. Springer International Publishing, 2015. ISBN 978-3-319-20827-5. doi: 10.1007/978-3-319-20828-2. URL https://doi.org/10.1007/978-3-319-20828-2. Schmidt-Hieber (2020) Johannes Schmidt-Hieber. Nonparametric regression using deep neural networks with ReLU activation function. The Annals of Statistics, 48(4):1875 – 1897, 2020. doi: 10.1214/19-AOS1875. URL https://doi.org/10.1214/19-AOS1875. Shen et al. (2022) Zuowei Shen, Haizhao Yang, and Shijun Zhang. Optimal approximation rate of relu networks in terms of width and depth. Journal de Mathématiques Pures et Appliquées, 157:101–135, 2022. Singh and Fischer (2024) Saurabh Singh and Ian Fischer. Stochastic sampling from deterministic flow models. arXiv preprint arXiv:2410.02217, 2024. Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International conference on machine learning, pages 2256–2265. pmlr, 2015. Song and Ermon (2019) Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019. Song et al. (2021) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=PxTIG12RRHS. Tang and Yang (2024) Rong Tang and Yun Yang. Adaptivity of diffusion models to manifold structures. In Sanjoy Dasgupta, Stephan Mandt, and Yingzhen Li, editors, Proceedings of The 27th International Conference on Artificial Intelligence and Statistics, volume 238 of Proceedings of Machine Learning Research, pages 1648–1656. PMLR, 02–04 May 2024. URL https://proceedings.mlr.press/v238/tang24a.html. Tao (2011) Terence Tao. An Introduction to Measure Theory, volume 126 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 2011. ISBN 978-0-8218-6919-2. URL https://bookstore.ams.org/gsm-126. Vershynin (2026) Roman Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2 edition, 2026. Watson et al. (2023) Joseph L Watson, David Juergens, Nathaniel R Bennett, Brian L Trippe, Jason Yim, Helen E Eisenach, Woody Ahern, Andrew J Borst, Robert J Ragotte, Lukas F Milles, et al. De novo design of protein structure and function with rfdiffusion. Nature, 620(7976):1089–1100, 2023. Weed and Bach (2019) Jonathan Weed and Francis Bach. Sharp asymptotic and finite-sample rates of convergence of empirical measures in wasserstein distance. Bernoulli, 25(4A):2620 – 2648, 2019. doi: 10.3150/18-BEJ1065. URL https://doi.org/10.3150/18-BEJ1065. Wu and Cai (2026) Jingda Wu and Changxiao Cai. Diffusion models are statistically optimal for learning low-dimensional multi-modal distributions. In Forty-third International Conference on Machine Learning, 2026. URL https://openreview.net/forum?id=L5JTAPUdbQ. Yarotsky (2017) Dmitry Yarotsky. Error bounds for approximations with deep relu networks. Neural Networks, 94:103–114, 2017. Zhou and Liu (2025) Zhengyu Zhou and Weiwei Liu. An error analysis of flow matching for deep generative modeling. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=vES22INUKm. Appendices Appendix A Omitted Proofs from Section 2 A.1 Existence and Uniqueness of the True Vector Velocity Field (Proof of Theorem 4) See 4 Proof. The characteristic function of XtX_t is given by, φt(κ)=exp(ικ⊤Xt)=∫exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1). _t(κ)=E ( κ X_t)= ( κ ( _t x_1 + _t x_0 ))dP_0(x_0)dP^emp_n(x_1). For any t∈[0,1−δ]t∈[0,1-δ], ∫|φt(κ)|κ≲δ−2κ1(2π)D/2<∞ | _t(κ)|dκ δ^-2 _1(2π)^D/2<∞, by Lemma 10. Hence, applying the Fourier inversion theorem, XtX_t has a density and is given by, ϱt(x)=(2π)−D∫exp(−ικ⊤x)φt(κ)κ. _t(x)=(2π)^-D (- κ x) _t(κ)dκ. By Lemma 10, since, ∫‖∫(x1−x0)exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)‖κ \| (x_1-x_0) ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1) \|dκ ≲ δκ12Γ(D+12)Γ(D2)×δκ1(2π)D/2δκ1(D+1)+D(2π)D/2δκ1D<∞, δ _1 2 ( D+12 ) ( D2 )× δ _1(2π)^D/2δ _1(D+1)+ D(2π)^D/2δ _1D<∞, (15) Applying the Fubini-Tonelli Theorem (Tao, 2011, Corollary 1.7.23), ∂tϱt(x)= _t _t(x)= (2π)−D∫exp(−ικ⊤x)φt′(κ)κ (2π)^-D (- κ x) _t (κ)dκ = = (2π)−D∫ικ⊤(α˙tx1+β˙tx0)exp(−ικ⊤x)exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)κ (2π)^-D κ ( α_tx_1+ β_tx_0) (- κ x) ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1)dκ = = −divx((2π)−D∫(α˙tx1+β˙tx0)exp(−ικ⊤x)exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)κ) -div_x ((2π)^-D ( α_tx_1+ β_tx_0) (- κ x) ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1)dκ ) = = −divx(∫(α˙tx1+β˙tx0)((2π)−D∫exp(−ικ⊤x)exp(ικ⊤(αtx1+βtx0))κ)CLOSE -div_x ( ( α_tx_1+ β_tx_0) ((2π)^-D (- κ x) ( κ ( _tx_1+ _tx_0) )dκ ) OPENdP0(x0)dPnemp(x1)). 56.9055ptdP_0(x_0)dP^emp_n(x_1) ). (16) We note that, ∫exp(ικ⊤x)v⋆(x,t)ϱt(x)x ( κ x)v_ (x,t) _t(x)dx = = (exp(ικ⊤Xt)v⋆(Xt,t)) ( ( κ X_t)v_ (X_t,t) ) = = (exp(ικ⊤Xt)(X˙t∣Xt)) ( ( κ X_t)E ( X_t X_t ) ) = = (exp(ικ⊤Xt)X˙t∣Xt) ( ( κ X_t) X_t X_t ) = = (exp(ικ⊤Xt)X˙t) ( ( κ X_t) X_t ) = = (exp(ικ⊤(αtX1+βtX0))(α˙tX1+β˙tX0)) ( ( κ ( _t X_1 + _t X_0))( α_tX_1+ β_tX_0) ) = = ∫exp(ικ⊤(αtx1+βtx0))(α˙tx1+β˙tx0)dP0(x0)dPnemp(x1) ( κ ( _t x_1 + _t x_0))( α_tx_1+ β_tx_0)dP_0(x_0)dP^emp_n(x_1) (17) Taking an inverse Fourier transform on both sides, we get, v⋆(x,t)ϱt(x) v_ (x,t) _t(x) = = (2π)−D∫exp(−ικ⊤x)exp(ικ⊤(αtx1+βtx0))(α˙tx1+β˙tx0)dP0(x0)dPnemp(x1)κ. (2π)^-D (- κ x) ( κ ( _tx_1 + _t x_0))( α_tx_1+ β_tx_0)dP_0(x_0)dP^emp_n(x_1)dκ. (18) Note that we can apply the inverse Fourier transform since (15) is satisfied. Hence, from (16) and (18), we get, ∂tϱt(x)+divx(v⋆(x,t)ϱt(x))=0. _t _t(x)+div_x(v_ (x,t) _t(x))=0. (19) Note that since v⋆(⋅,t)v_ (·,t) is uniformly (in t∈[0,1−δ]t∈[0,1-δ]) Dδ−3κ1Dδ^-3 _1-Lipschitz, by Corollary 16, The Picard–Lindelöf theorem (Coddington and Levinson, 1955, Theorem 1.3.1) implies that (4) has a unique solution for all initialization x∈ℝDx∈ R ^D. From Santambrogio (2015, Theorem 4.4), (Yt)♯P0(Y_t)_ P_0 satisfies the continuity equation (19), which admits a unique solution. ∎ Lemma 10. For any 0≤t<10≤ t<1, the following holds, ∫|∫exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)|κ≤(2π)D/2βtD, | ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1) |dκ≤ (2π)^D/2 _t^D, (20) ∫‖∫(x1−z)exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)‖κ \| (x_1-z) ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1) \|dκ ≤ ≤ β˙t2Γ(D+12)Γ(D2)×(2π)D/2βtD+1+α˙tD(2π)D/2βtD. β_t 2 ( D+12 ) ( D2 )× (2π)^D/2 _t^D+1+ α_t D(2π)^D/2 _t^D. (21) Proof. ∫|∫exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)|κ | ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1) |dκ ≤ ≤ 1n∑i=1n∫|∫exp(ικ⊤(αtxi+βtx0))dP0(x0)|κ 1n _i=1^n | ( κ ( _tx_i+ _tx_0) )dP_0(x_0) |dκ = = 1n∑i=1n∫|exp(ιαtκ⊤xi)exp(−12βt2‖κ‖2)|κ 1n _i=1^n | ( _t κ x_i) (- 12 _t^2\|κ\|^2 ) |dκ ≤ ≤ ∫exp(−12βt2‖κ‖2)κ (- 12 _t^2\|κ\|^2 )dκ = = (2πβt2)D/2. ( 2π _t^2 )^D/2. ∫‖∫(α˙tx1+β˙tx0)exp(ικ⊤(αtx1+βtx0))dP0(x0)dPnemp(x1)‖κ \| ( α_tx_1+ β_tx_0) ( κ ( _tx_1+ _tx_0) )dP_0(x_0)dP^emp_n(x_1) \|dκ ≤ ≤ 1n∑i=1n∫‖∫(α˙txi+β˙tx0)exp(ικ⊤(αtxi+βtx0))dP0(x0)‖κ 1n _i=1^n \| ( α_tx_i+ β_tx_0) ( κ ( _tx_i+ _tx_0) )dP_0(x_0) \|dκ ≤ ≤ β˙tn∑i=1n∫‖∫x0exp(ικ⊤(αtxi+βtx0))dP0(x0)‖κ β_tn _i=1^n \| x_0 ( κ ( _tx_i+ _tx_0) )dP_0(x_0) \|dκ +α˙tn∑i=1n∫‖∫xiexp(ικ⊤(αtxi+βtx0))dP0(x0)‖dκ + α_tn _i=1^n \| x_i ( κ ( _tx_i+ _tx_0) )dP_0(x_0) \|dκ ≤ ≤ β˙tn∑i=1n∫‖∫x0exp(ικ⊤(αtxi+βtx0))dP0(x0)‖κ β_tn _i=1^n \| x_0 ( κ ( _tx_i+ _tx_0) )dP_0(x_0) \|dκ +α˙tn∑i=1n∫‖xiexp(ιαtκ⊤xi)∫exp(ικ⊤βtx0)dP0(x0)‖dκ + α_tn _i=1^n \|x_i ( _tκ x_i ) ( κ _tx_0 )dP_0(x_0) \|dκ = = β˙tn∑i=1n∫‖exp(ιαtκ⊤xi)⋅ιβtκexp(−βt2‖κ‖22)‖κ β_tn _i=1^n \| ( _t κ x_i)· _tκ (- _t^2\|κ\|^22 ) \|dκ +α˙tn∑i=1n∫‖xiexp(ιαtκ⊤xi)exp(−βt2‖κ‖22)‖dκ + α_tn _i=1^n \|x_i ( _t κ x_i) (- _t^2\|κ\|^22 ) \|dκ ≤ ≤ β˙t∫‖κ‖exp(−βt2‖κ‖22)κ+α˙tn∑i=1n‖xi‖∫exp(−βt2‖κ‖22)κ β_t \|κ\| (- _t^2\|κ\|^22 )dκ+ α_tn _i=1^n\|x_i\| (- _t^2\|κ\|^22 )dκ = = 2Γ(D+12)Γ(D2)×(2π)D/2βtD+1+α˙tn∑i=1n‖xi‖×(2π)D/2βtD 2 ( D+12 ) ( D2 )× (2π)^D/2 _t^D+1+ α_tn _i=1^n\|x_i\|× (2π)^D/2 _t^D ≤ ≤ β˙t2Γ(D+12)Γ(D2)×(2π)D/2βtD+1+α˙tD(2π)D/2βtD. β_t 2 ( D+12 ) ( D2 )× (2π)^D/2 _t^D+1+ α_t D(2π)^D/2 _t^D. ∎ A.2 Equivalence to the Conditional Flow-matching Objective Lemma 11. ℒCFM(v)=ℒ(v)+∫01‖v⋆(Xt,t)−Xt˙‖2tL_CFM(v)=L(v)+ _0^1E\|v_ (X_t,t)- X_t\|^2dt. Proof. Note that, ℒCFM(v) _CFM(v) = = ∫01‖v(Xt,t)−Xt˙‖2t _0^1E\|v(X_t,t)- X_t\|^2dt = = ∫01(‖v(Xt,t)−v⋆(Xt,t)‖2+‖v⋆(Xt,t)−Xt˙‖2+2⟨v(Xt,t)−v⋆(Xt,t),v⋆(Xt,t)−Xt˙⟩)t _0^1E (\|v(X_t,t)-v_ (X_t,t)\|^2+\|v_ (X_t,t)- X_t\|^2+2 v(X_t,t)-v_ (X_t,t),v_ (X_t,t)- X_t )dt = = ∫01‖v(Xt,t)−v⋆(Xt,t)‖2t+∫01‖v⋆(Xt,t)−Xt˙‖2t _0^1E\|v(X_t,t)-v_ (X_t,t)\|^2dt+ _0^1E\|v_ (X_t,t)- X_t\|^2dt +2∫01(⟨v(Xt,t)−v⋆(Xt,t),v⋆(Xt,t)−Xt˙⟩∣Xt)dt +2 _0^1EE ( v(X_t,t)-v_ (X_t,t),v_ (X_t,t)- X_t X_t )dt = = ℒ(v)+∫01‖v⋆(Xt,t)−Xt˙‖2t. (v)+ _0^1E\|v_ (X_t,t)- X_t\|^2dt. ∎ A.3 Equivalence of the probability-flow ODE and the SDE family Let ϱtt∈[0,1]\ _t\_t∈[0,1] be the marginal law of the stochastic interpolant Xt=αtX1+βtX0X_t= _tX_1+ _tX_0, absolutely continuous with density ϱt _t for t∈(0,1)t∈(0,1), and let v⋆(x,t)v_ (x,t) denote the interpolant velocity field. For a constant σ>0σ>0 define t(x):=v⋆(x,t)+σ22∇logϱt(x). b_t(x)\;:=\;v_ (x,t)+ σ^22∇ _t(x). (22) By Theorem 4, v⋆v_ solves, in the sense of distributions on (0,1)×ℝD(0,1)×R^D, ∂tϱt+∇⋅(v⋆ϱt)=0,ϱt|t→0+=ϱ0=(0,ID) weakly. _t _t+∇·(v_ _t)=0, _t |_t→ 0^+= _0=N(0,I_D) weakly. (23) We first record the regularity of the score, which is what makes (24) below well posed. Lemma 12 (Lipschitzness of the score). Under Assumption 1 and Design Choices D 1, for every t∈[0,1−δ]t∈[0,1-δ] the map x↦∇logϱt(x)x ∇ _t(x) is differentiable with ∇2logϱt(x)=αt2βt4Var(X1∣Xt=x)−1βt2ID,∥∇2logϱt(x)∥≲Dδ−4κ1,∇^2 _t(x)\;=\; _t^2 _t^4\,Var (X_1 X_t=x )\;-\; 1 _t^2\,I_D, \|∇^2 _t(x) \|\; \;D\,δ^-4 _1, uniformly in x∈ℝDx ^D and t∈[0,1−δ]t∈[0,1-δ]. In particular ∇logϱt(⋅)∇ _t(·) is globally Lipschitz with constant (Dδ−4κ1)O(Dδ^-4 _1), uniformly in t∈[0,1−δ]t∈[0,1-δ]. Proof. By Lemma 15, ∇logϱt(x)=−βt−2(x−αtm(x,t))∇ _t(x)=- _t^-2 (x- _tm(x,t) ) with m(x,t)=[X1∣Xt=x]m(x,t)=E[X_1 X_t=x], and ∇xm(x,t)=αtβt−2Var(X1∣Xt=x) _xm(x,t)= _t _t^-2\,Var(X_1 X_t=x). Differentiating the first identity in x and substituting the second gives the display. For the bound, Var(X1∣Xt=x)⪰0Var(X_1 X_t=x) 0 and, since X1∼PnempX_1 P^emp_n is supported in [0,1]D[0,1]^D by Assumption 1, each conditional coordinate variance is at most 1/41/4, so ‖Var(X1∣Xt=x)‖≤trace(Var(X1∣Xt=x))≤D4. \|Var(X_1 X_t=x) \|\;≤\;trace (Var(X_1 X_t=x) )\;≤\; D4. With αt≤1 _t≤ 1 and βt−1≲δ−κ1 _t^-1 δ^- _1 by D 1(a), ∥∇2logϱt(x)∥≤D4αt2βt−4+βt−2≲Dδ−4κ1\|∇^2 _t(x)\|≤ D4 _t^2 _t^-4+ _t^-2 Dδ^-4 _1. ∎ Theorem 13 (Marginal equivalence of the ODE and the SDE family). Fix δ∈(0,1/2)δ∈(0,1/2). Under Assumption 1 and Design Choices D 1, there is a pathwise unique strong solution Xtt∈[0,1−δ]\X_t\_t∈[0,1-δ] of dXt=t(Xt)dt+σdWt,X0∼(0,ID),dX_t= b_t(X_t)\,dt+σ\,dW_t, X_0 (0,I_D), (24) with Law(Xt)=ϱtLaw(X_t)= _t for every t∈[0,1−δ]t∈[0,1-δ]. Proof. Let ℒtσf:=t⋅∇f+σ22ΔfL_t^σf:= b_t·∇ f+ σ^22 f be the generator of (24), with formal adjoint (ℒtσ)∗(L_t^σ)^*. Since ϱt _t is a Gaussian smoothing of a compactly supported measure it is smooth and strictly positive on ℝDR^D for every t∈(0,1)t∈(0,1), so the following manipulations are classical. Using the score identity ϱt∇logϱt=∇ϱt _t∇ _t=∇ _t, (ℒtσ)∗ϱt= (L_t^σ)^* _t= −∇⋅(tϱt)+σ22Δϱt -∇·( b_t _t)+ σ^22 _t = = −∇⋅(v⋆ϱt+σ22∇ϱt)+σ22Δϱt -∇·\! (v_ _t+ σ^22∇ _t )+ σ^22 _t = = −∇⋅(v⋆ϱt) -∇·(v_ _t) = = ∂tϱt. \; _t _t. (25) In the above calculations, (25) follows from (23). This holds for every σ≥0σ≥ 0 simultaneously: the diffusive term σ22Δϱt σ^22 _t contributed by the noise is exactly cancelled by the score correction in (22), not merely asymptotically, and (25) does not depend on σ at all. Hence t↦ϱtdxt _t\,dx is a probability solution of ∂tμt=(ℒtσ)∗μt,t∈[0,1−δ],μ0=ϱ0=(0,ID), _t _t=(L_t^σ)^* _t, t∈[0,1-δ], _0= _0=N(0,I_D), (26) for every admissible σ, the initial condition being the weak limit in (23). By Corollary 16, v⋆(⋅,t)v_ (·,t) is (Dδ−3κ1)O(Dδ^-3 _1)-Lipschitz uniformly over t∈[0,1−δ]t∈[0,1-δ], and by Lemma 12, ∇logϱt(⋅)∇ _t(·) is (Dδ−4κ1)O(Dδ^-4 _1)-Lipschitz uniformly over the same range. Hence t(x)=v⋆(x,t)+σ22∇logϱt(x) b_t(x)=v_ (x,t)+ σ^22∇ _t(x) is globally Lipschitz in x, uniformly over t∈[0,1−δ]t∈[0,1-δ], with Lipschitz constant Λδ:=(D(1∨σ2)δ−4κ1) _δ:=O (D(1 σ^2)δ^-4 _1 ), and is of at most linear growth, since ‖t(0)‖≲δ−κ1D\| b_t(0)\| δ^- _1 D by Corollary 16 and ‖∇logϱt(0)‖≤βt−2αt‖X1‖≲δ−2κ1D\|∇ _t(0)\|≤ _t^-2 _tE\|X_1\| δ^-2 _1 D. Since σ is a fixed positive constant, the classical Itô existence–uniqueness theorem for SDEs with globally Lipschitz coefficients of linear growth (Karatzas and Shreve, 1991, Theorem 2.9, Chapter 5) gives a pathwise unique strong solution Xtt∈[0,1−δ]\X_t\_t∈[0,1-δ] of (24) for any initial law with finite second moment, in particular for Law(X0)=(0,ID)Law(X_0)=N(0,I_D). The drift b is globally Lipschitz in x uniformly in t∈[0,1−δ]t∈[0,1-δ] and of linear growth, and the diffusion coefficient σIDσ I_D is constant and nondegenerate, so the Cauchy problem (26) admits at most one probability solution on [0,1−δ][0,1-δ] with initial datum ϱ0 _0 (Bogachev et al., 2015, Ch. 9). Equivalently, and without appeal to that theory: t↦Law(Xt)t (X_t) for the pathwise-unique solution of (24) is itself a probability solution of (26) with the same initial datum, and two probability solutions of a well-posed SDE with the same initial law have equal marginals. From the above arguments, the map, t↦ϱt _t is also a probability solution of (26) with ϱt|t=0=ϱ0=(0,ID) _t|_t=0= _0=N(0,I_D). Uniqueness therefore forces Law(Xt)=ϱtLaw(X_t)= _t for every t∈[0,1−δ]t∈[0,1-δ]. ∎ Appendix B Omitted Proofs from Section 4.2 B.1 Error Decomposition See 7 Proof. Write Ω=C([0,1−δ];ℝD) =C([0,1-δ]; R ^D) for path space and let P denote the law on Ω of the solution Xt\X_t\ of dXt=(atv⋆(Xt,t)−btXt)dt+σdWt,X0∼(0,ID),dX_t= (a_tv_ (X_t,t)-b_tX_t )dt+ _t, X_0 (0,I_D), which exists and is pathwise unique by Theorem 13. Let Q denote the law on Ω of the Euler–Maruyama process (10), namely the solution of dX^t=(atkv^(X^tk,tk)−btX^t)dt+σdWt,t∈[tk,tk+1),X^0∼(0,ID).d X_t= (a_t_k v( X_t_k,t_k)-b_t X_t )dt+ _t, t∈[t_k,t_k+1), X_0 (0,I_D). Both are driven by the same diffusion coefficient σ>0σ>0 and started from the same law, so P and Q differ only through their drifts. Writing these as functionals of the path, μt(ω)=atv⋆(ωt,t)−btωt,μt(ω)=atkv^(ωtk,tk)−btωtfor t∈[tk,tk+1),μ P_t(ω)=a_tv_ ( _t,t)-b_t _t, μ Q_t(ω)=a_t_k v( _t_k,t_k)-b_t _t t∈[t_k,t_k+1), both are adapted, and the linear terms −btωt-b_t _t coincide, so that μt(ω)−μt(ω)=atv⋆(ωt,t)−atkv^(ωtk,tk),t∈[tk,tk+1).μ P_t(ω)-μ Q_t(ω)=a_tv_ ( _t,t)-a_t_k v( _t_k,t_k), t∈[t_k,t_k+1). (27) By Corollary 16, ‖v⋆(x,t)‖≲δ−κ1(‖x‖+D)\|v_ (x,t)\| δ^- _1(\|x\|+ D) uniformly on [0,1−δ][0,1-δ]; by D 1(a), supt≤1−δ|at|≲δ−κ1/2 _t≤ 1-δ|a_t| δ^- _1/2; and ‖v^(⋅,⋅)‖≤D(v^)≤DK\| v(·,·)\|≤ D\,K( v)≤ DK. Hence the drift difference (27) satisfies ∥μt(ω)−μt(ω)∥≤Cδ(1+sups≤t∥ωs∥),t∈[0,1−δ],Cδ≲δ−κ1/2D(δ−κ1+K), \|μ P_t(ω)-μ Q_t(ω) \|\;≤\;C_δ (1+ _s≤ t\| _s\| ), t∈[0,1-δ], C_δ\; \;δ^- _1/2 D\, (δ^- _1+K ), which is finite for each fixed δ>0δ>0. Since Beneš’s condition requires only linear growth of the drift difference in the sup-norm of the path, with no restriction on the constant, it holds here (Karatzas and Shreve, 1991, Corollary 3.5.16), and the stochastic exponential below is a true martingale. Girsanov’s theorem (Øksendal, 2003, Theorem 8.6.5) therefore gives ∼ Q P with dd(ω)=exp(1σ∫01−δ⟨μt−μt,dWt⟩−12σ2∫01−δ‖μt−μt‖2t), d Pd Q(ω)= ( 1σ _0^1-δ μ P_t-μ Q_t,dW_t - 12σ^2 _0^1-δ \|μ P_t-μ Q_t \|^2dt ), and taking log(d/d)E_ P ( / ), the stochastic integral being a true martingale under P, KL(∥)=12σ2∑k=0N−1∫tktk+1∥atv⋆(Xt,t)−atkv^(Xtk,tk)∥2dt.KL( P\,\|\, Q)= 12σ^2 _k=0^N-1 _t_k^t_k+1E \|a_tv_ (X_t,t)-a_t_k v(X_t_k,t_k) \|^2dt. (28) The marginal laws at the terminal time are P1−δ=∘π1−δ−1P_1-δ= P _1-δ^-1 and P^1−δ(v^)=∘π1−δ−1 P_1-δ( v)= Q _1-δ^-1, where π1−δ(ω)=ω1−δ _1-δ(ω)= _1-δ, the first identification being Theorem 13. Since the Kullback–Leibler divergence contracts under the pushforward by π1−δ _1-δ, KL(P1−δ∥P^1−δ(v^))≤KL(∥).KL (P_1-δ\,\|\, P_1-δ( v) ) ( P\,\|\, Q). (29) It remains to split (28). Adding and subtracting atkv⋆(Xtk,tk)a_t_kv_ (X_t_k,t_k) and using ‖u+v‖2≤2‖u‖2+2‖v‖2\|u+v\|^2≤ 2\|u\|^2+2\|v\|^2, KL(∥)≤ ( P\,\|\, Q)≤ 1σ2∑k=0N−1∫tktk+1‖atv⋆(Xt,t)−atkv⋆(Xtk,tk)‖2t \; 1σ^2 _k=0^N-1 _t_k^t_k+1E \|a_tv_ (X_t,t)-a_t_kv_ (X_t_k,t_k) \|^2dt +1σ2∑k=0N−1∫tktk+1atk2∥v⋆(Xtk,tk)−v^(Xtk,tk)∥2dt. + 1σ^2 _k=0^N-1 _t_k^t_k+1a_t_k^2\,E \|v_ (X_t_k,t_k)- v(X_t_k,t_k) \|^2dt. The second sum has an integrand constant in t on each [tk,tk+1)[t_k,t_k+1) and so equals σ−2∑khkatk2‖v⋆(Xtk,tk)−v^(Xtk,tk)‖2σ^-2 _kh_ka_t_k^2\,E\|v_ (X_t_k,t_k)- v(X_t_k,t_k)\|^2. For the first, adding and subtracting atv⋆(Xtk,tk)a_tv_ (X_t_k,t_k) and splitting once more, ‖atv⋆(Xt,t)−atkv⋆(Xtk,tk)‖2≤2at2‖v⋆(Xt,t)−v⋆(Xtk,tk)‖2+2(at−atk)2‖v⋆(Xtk,tk)‖2. \|a_tv_ (X_t,t)-a_t_kv_ (X_t_k,t_k) \|^2≤ 2a_t^2 \|v_ (X_t,t)-v_ (X_t_k,t_k) \|^2+2(a_t-a_t_k)^2 \|v_ (X_t_k,t_k) \|^2. Combining the three displays with (29) gives the claim. ∎ B.2 Discretization Error See 8 Proof. Throughout, t∈[tk,tk+1)t∈[t_k,t_k+1) for the relevant k, and all expectations are under the law of Xtt∈[0,1−δ]\X_t\_t∈[0,1-δ] solving (9). By Theorem 13, Law(Xs)=ϱs=Law(αsX1+βsX0)Law(X_s)= _s=Law( _sX_1+ _sX_0) for each fixed s. Since X1∼PnempX_1 P^emp_n is supported in [0,1]D[0,1]^D by Assumption 1 and X0∼(0,ID)X_0 (0,I_D), ‖Xs‖2≤2αs2‖X1‖2+2βs2‖X0‖2≲D,s∈[0,1−δ].E\|X_s\|^2≤ 2 _s^2\,E\|X_1\|^2+2 _s^2\,E\|X_0\|^2 D, s∈[0,1-δ]. (30) Combining (30) with Corollary 16, ‖v⋆(Xs,s)‖2≲δ−2κ1(‖Xs‖2+D)≲Dδ−2κ1.E \|v_ (X_s,s) \|^2 δ^-2 _1 (E\|X_s\|^2+D ) D\,δ^-2 _1. (31) We write μs(x):=asv⋆(x,s)−bsx _s(x):=a_sv_ (x,s)-b_sx for the drift of (9). By (31), D 1(a) (which gives as2≲δ−κ1a_s^2 δ^- _1 and |bs|≲δ−κ1|b_s| δ^- _1) and (30), ‖μs(Xs)‖2≤2as2‖v⋆(Xs,s)‖2+2bs2‖Xs‖2≲Dδ−3κ1.E\| _s(X_s)\|^2≤ 2a_s^2\,E\|v_ (X_s,s)\|^2+2b_s^2\,E\|X_s\|^2 D\,δ^-3 _1. (32) Integrating (9) from tkt_k to t, Xt−Xtk=∫tktμs(Xs)s+σ(Wt−Wtk).X_t-X_t_k= _t_k^t _s(X_s)\,ds+σ (W_t-W_t_k ). By ‖u+v‖2≤2‖u‖2+2‖v‖2\|u+v\|^2≤ 2\|u\|^2+2\|v\|^2, the Cauchy–Schwarz inequality applied to the drift integral, and ‖Wt−Wtk‖2=D(t−tk)E\|W_t-W_t_k\|^2=D(t-t_k), ‖Xt−Xtk‖2 \|X_t-X_t_k \|^2 ≤2(t−tk)∫tkt‖μs(Xs)‖2s+2σ2D(t−tk) ≤ 2(t-t_k) _t_k^tE\| _s(X_s)\|^2\,ds+2σ^2D\,(t-t_k) ≲Dδ−3κ1(t−tk)2+σ2D(t−tk) D\,δ^-3 _1(t-t_k)^2+σ^2D\,(t-t_k) ≲D(σ2+δ−3κ1h)(t−tk). \;D (σ^2+δ^-3 _1h )(t-t_k). (33) Note that the σ2σ^2 appearing here is exactly cancelled by the factor σ−2σ^-2 carried by Lemma 7, so the resulting contribution to the KL bound is independent of σ; see Remark 1 below. Splitting in space and in time, ‖v⋆(Xt,t)−v⋆(Xtk,tk)‖2 \|v_ (X_t,t)-v_ (X_t_k,t_k) \|^2 ≤ ≤\; 2‖v⋆(Xt,t)−v⋆(Xtk,t)‖2+2‖v⋆(Xtk,t)−v⋆(Xtk,tk)‖2. 2\,E \|v_ (X_t,t)-v_ (X_t_k,t) \|^2+2\,E \|v_ (X_t_k,t)-v_ (X_t_k,t_k) \|^2. (34) For the first summand, Corollary 16 gives Lip(v⋆(⋅,t))≲Dδ−3κ1Lip(v_ (·,t)) Dδ^-3 _1 uniformly in t∈[0,1−δ]t∈[0,1-δ]. Thus, from (33), ‖v⋆(Xt,t)−v⋆(Xtk,t)‖2≲D2δ−6κ1‖Xt−Xtk‖2≲D3δ−6κ1(σ2+δ−3κ1h)(t−tk).E \|v_ (X_t,t)-v_ (X_t_k,t) \|^2 D^2δ^-6 _1\,E\|X_t-X_t_k\|^2 D^3δ^-6 _1 (σ^2+δ^-3 _1h )(t-t_k). (35) For the second summand, the two evaluations share the spatial argument XtkX_t_k, so only the one-time marginal of XtkX_t_k is involved. By Lemma 18, ‖∂tv⋆(x,t)‖≲δ−κ1(1+‖x‖)\| _tv_ (x,t)\| δ^- _1(1+\|x\|) uniformly on [0,1−δ][0,1-δ]. Hence, by the mean value theorem and (30), ‖v⋆(Xtk,t)−v⋆(Xtk,tk)‖2≲δ−2κ1(1+‖Xtk‖2)(t−tk)2≲Dδ−2κ1(t−tk)2.E \|v_ (X_t_k,t)-v_ (X_t_k,t_k) \|^2 δ^-2 _1\,E (1+\|X_t_k\|^2 )(t-t_k)^2 D\,δ^-2 _1(t-t_k)^2. (36) Using at2≲δ−κ1a_t^2 δ^- _1, then (34), (35) and (36), and finally ∫tktk+1(t−tk)t=h2/2 _t_k^t_k+1(t-t_k)dt=h^2/2 and ∫tktk+1(t−tk)2t=h3/3 _t_k^t_k+1(t-t_k)^2dt=h^3/3, ∑k=0N−1∫tktk+1at2‖v⋆(Xt,t)−v⋆(Xtk,tk)‖2t≲ _k=0^N-1 _t_k^t_k+1a_t^2\,E \|v_ (X_t,t)-v_ (X_t_k,t_k) \|^2dt \; δ−κ1[D3δ−6κ1(σ2+δ−3κ1h)∑k=0N−1h22+Dδ−2κ1∑k=0N−1h33] δ^- _1 [D^3δ^-6 _1 (σ^2+δ^-3 _1h ) _k=0^N-1 h^22+D\,δ^-2 _1 _k=0^N-1 h^33 ] ≲ \; D3δ−7κ1(σ2+δ−3κ1h)h+Dδ−3κ1h2 D^3δ^-7 _1 (σ^2+δ^-3 _1h )\,h+D\,δ^-3 _1h^2 ≲ \; D3(σ2∨1)δ−10κ1h, D^3 (σ^2 1 )\,δ^-10 _1\,h, where the second line uses Nh=1−δ≤1Nh=1-δ≤ 1, so that ∑kh2≤h _kh^2≤ h and ∑kh3≤h2 _kh^3≤ h^2, and the last uses h≤1h≤ 1 to absorb h2h^2 into h. This is (13). To bound (14), we note that this term involves v⋆v_ at the single time tkt_k only, so (31) applies directly. By D 1(b), |at−atk|≲δ−κ2(t−tk)|a_t-a_t_k| δ^- _2(t-t_k), whence ∑k=0N−1∫tktk+1(at−atk)2‖v⋆(Xtk,tk)‖2t≲Dδ−2κ1−2κ2∑k=0N−1∫tktk+1(t−tk)2t≲Dδ−2(κ1+κ2)h2, _k=0^N-1 _t_k^t_k+1(a_t-a_t_k)^2\,E \|v_ (X_t_k,t_k) \|^2dt D\,δ^-2 _1-2 _2 _k=0^N-1 _t_k^t_k+1(t-t_k)^2dt D\,δ^-2( _1+ _2)h^2, using ∑kh3≤h2 _kh^3≤ h^2. This is (14). Since h≤1h≤ 1 and 2(κ1+κ2)≤10κ1+2κ22( _1+ _2)≤ 10 _1+2 _2, the two bounds combine as stated. ∎ Remark 1. The diffusion coefficient cannot be tuned to remove the O(h)O(h) term. Lemma 7 carries a factor σ−2σ^-2, while the leading term of (13) carries σ2σ^2 through the Brownian increment in (33); the two cancel, and the discretization contribution to KL(P1−δ∥P^1−δ(v^))KL(P_1-δ\| P_1-δ( v)) is of order D3δ−c1hD^3δ^-c_1h, uniformly in σ. B.3 Approximation Result See 9 Proof. Fix ϵ∈(0,1)ε∈(0,1) and R≥1R≥ 1. By Corollary 16, for every t∈[0,1−δ]t∈[0,1-δ], Lip(v⋆(⋅,t))=:Lt≲Lδ≍δ−3κ1,supx∈[−R,R]D‖v⋆(x,t)‖≲δ−κ1(RD+D)≲δ−κ1RD.Lip (v_ (·,t) )=:L_t\; \;L_δ δ^-3 _1,\; _x∈[-R,R]^D\|v_ (x,t)\|\; \;δ^- _1 (R D+ D )\; \;δ^- _1R D. (37) Let, K0:=cδ−(3κ1)∨1RDK_0\;:=\;c\,δ^-(3 _1) 1R D (38) with c a large enough absolute constant that K0≥1+supt∈[0,1−δ]supx∈[−R,R]D‖v⋆(x,t)‖∞K_0≥ 1+ _t∈[0,1-δ] _x∈[-R,R]^D\|v_ (x,t)\|_∞, which is possible by (37). Fix i∈0,…,Ni∈\0,…,N\ and ℓ∈[D] ∈[D]. Define gi(ℓ):[0,1]D→ℝg_i^( ):[0,1]^D→ R by gi(ℓ)(y)=12RLti(v⋆(ℓ)(2Ry−R,ti)−v⋆(ℓ)(−R,ti)).g_i^( )(y)\;=\; 12RL_t_i (v_ ^( ) (2Ry-R1,\,t_i )-v_ ^( ) (-R1,\,t_i ) ). Since x↦v⋆(ℓ)(x,ti)x v_ ^( )(x,t_i) is LtiL_t_i-Lipschitz, gi(ℓ)g_i^( ) is 11-Lipschitz on [0,1]D[0,1]^D and gi(ℓ)(0)=0g_i^( )(0)=0. Let hi(ℓ)h_i^( ) be the ε -approximating ReLU network for gi(ℓ)g_i^( ) furnished by Lemma 22, and set v~ti(ℓ)(x)= 2RLtihi(ℓ)(x+R2R)+v⋆(ℓ)(−R,ti). v^( )_t_i(x)\;=\;2RL_t_i\,h^( )_i\! ( x+R12R )+v_ ^( ) (-R1,t_i ). By construction, supx∈[−R,R]D|v~ti(ℓ)(x)−v⋆(ℓ)(x,ti)|= _x∈[-R,R]^D | v^( )_t_i(x)-v_ ^( )(x,t_i) |= 2RLtisupy∈[0,1]D|hi(ℓ)(y)−gi(ℓ)(y)| \;2RL_t_i _y∈[0,1]^D |h^( )_i(y)-g^( )_i(y) | ≤ ≤ 2RLtiε \;2RL_t_i ≲ RLδε. \;R\,L_δ\, . (39) Let v~ti v_t_i be obtained by stacking v~ti(ℓ)ℓ=1D\ v^( )_t_i\_ =1^D. By Lemma 22, ‖hi(ℓ)‖∞≤D\|h^( )_i\|_∞≤ D, so on all of ℝD R ^D (v~ti)≲D(ε−Dlog(1/ε)+D),ℒ(v~ti)≲log(1/ε)+D,ℬ(v~ti)≲RLδ+δ−κ1RD≲K0,W( v_t_i) D ( ^-D (1/ )+D ),\;L( v_t_i) (1/ )+D,\;B( v_t_i) RL_δ+δ^- _1R D K_0, and, since hi(ℓ)h^( )_i is 10D10D-Lipschitz in ℓ1 _1 and the inner affine map has gain 1/(2R)1/(2R), Lip(v~ti)≤ 2RLti⋅10D⋅12R⋅D≲D3/2Lδ.Lip( v_t_i)\;≤\;2RL_t_i· 10D· 12R· D\; \;D^3/2L_δ. (40) We define the clipping function, clK0(u)=ReLU(u+K0)−ReLU(u−K0)−K0cl_K_0(u)=ReLU(u+K_0)-ReLU(u-K_0)-K_0, so that clK0(u)=ucl_K_0(u)=u for |u|≤K0|u|≤ K_0, |clK0|≤K0|cl_K_0|≤ K_0 everywhere, and clK0cl_K_0 is 11-Lipschitz; write clK0(v~ti)cl_K_0( v_t_i) for its coordinatewise application. By (37), (38) and (39), clipping does not alter v~ti v_t_i on [−R,R]D[-R,R]^D once ε≤1 ≤ 1. Next set χ(x)=ReLU(1−∑ℓ=1DReLU(ReLU(xℓ)+ReLU(−xℓ)−R)),χ(x)\;=\;ReLU\! (1- _ =1^DReLU (ReLU(x_ )+ReLU(-x_ )-R ) ), a ReLU network of width O(D)O(D) and depth 33 with χ≡1χ≡ 1 on [−R,R]D[-R,R]^D, χ≡0χ≡ 0 off [−(R+1),R+1]D[-(R+1),R+1]^D, 0≤χ≤10≤χ≤ 1, and Lip(χ)≲DLip(χ) D. We define t−1=−1t_-1=-1 and tN+1=1t_N+1=1, let δi=12((ti−ti−1)∧(ti+1−ti)) _i= 12 ((t_i-t_i-1) (t_i+1-t_i) ) and ξa,b(x)=ReLU(x+a−b)−ReLU(x+ba−b)−ReLU(x−ba−b)+ReLU(x−a−b), _a,b(x)=ReLU ( x+aa-b )-ReLU ( x+ba-b )-ReLU ( x-ba-b )+ReLU ( x-aa-b ), so that ξi:=ξδi/2,δi/4(⋅−ti) _i:= _ _i/2,\, _i/4(\,·\,-t_i) satisfies 0≤ξi≤10≤ _i≤ 1, ξi≡1 _i≡ 1 on |t−ti|≤δi/4|t-t_i|≤ _i/4, and suppξi⊆|t−ti|≤δi/2supp _i \|t-t_i|≤ _i/2\; its weights are of order 1/δi≍N1/ _i N. Since δi≤12|ti−tj| _i≤ 12|t_i-t_j| for j≠ij≠ i, the supports of ξi\ _i\ are pairwise disjoint and ξj(ti)=0 _j(t_i)=0 for j≠ij≠ i. Let ψi(x,t)=min(χ(x),ξi(t))=χ(x)−ReLU(χ(x)−ξi(t)), _i(x,t)\;=\; (χ(x), _i(t) )\;=\;χ(x)-ReLU (χ(x)- _i(t) ), which is again a ReLU network, takes values in [0,1][0,1], equals 11 on [−R,R]D×|t−ti|≤δi/4[-R,R]^D×\|t-t_i|≤ _i/4\, and vanishes whenever χ(x)=0χ(x)=0 or ξi(t)=0 _i(t)=0. Finally, let Π denote the ReLU approximate-multiplication network on [−K0,K0]×[0,1][-K_0,K_0]×[0,1] of accuracy ϵ/(4D)ε/(4 D), of width and depth (log(K0D/ϵ))O ( (K_0 D/ε) ), and satisfying Π(u,0)=0 (u,0)=0 identically. We define v0(ℓ)(x,t)=∑i=0NΠ(clK0(v~ti(ℓ)(x)),ψi(x,t)),ℓ∈[D].v_0^( )(x,t)\;=\; _i=0^N (cl_K_0 ( v^( )_t_i(x) ),\; _i(x,t) ), ∈[D]. Fix i∈0,…,Ni∈\0,…,N\ and x∈[−R,R]Dx∈[-R,R]^D. Then ψj(x,ti)=0 _j(x,t_i)=0 for j≠ij≠ i, so those summands vanish exactly, while ψi(x,ti)=1 _i(x,t_i)=1. Hence, by (39), |v0(ℓ)(x,ti)−v⋆(ℓ)(x,ti)|≤ϵ4D+2RLtiε. |v_0^( )(x,t_i)-v^( )_ (x,t_i) |\;≤\; ε4 D+2RL_t_i . Choosing ε=ϵ/(8DRLδ) =ε / (8 D\,RL_δ ) makes the second term at most ϵ/(4D)ε/(4 D), so that the coordinatewise error is at most ϵ/(2D)ε/(2 D) and supx∈[−R,R]D‖v0(x,ti)−v⋆(x,ti)‖≤D⋅ϵ2D≤ϵfor every i=0,…,N. _x∈[-R,R]^D \|v_0(x,t_i)-v_ (x,t_i) \|\;≤\; D· ε2 D\;≤\;ε every i=0,…,N. With this ε we have ε−D≲RDδ−3κ1Dϵ−D ^-D R^Dδ^-3 _1Dε^-D and log(1/ε)≍log(K0D/ϵ)≍log(1/δ)+log(1/ϵ)+logR (1/ ) (K_0 D/ε) (1/δ)+ (1/ε)+ R, whence, summing the N+1N+1 parallel branches and recalling that D is fixed, (v0)≲ (v_0) N(ε−Dlog(1/ε)+log(K0D/ϵ)+D)≲NRDδ−3κ1Dϵ−D(log(1/δ)+log(1/ϵ)+logR), \;N ( ^-D (1/ )+ (K_0 D/ε)+D )\; \;NR^Dδ^-3 _1Dε^-D ( (1/δ)+ (1/ε)+ R ), ℒ(v0)≲ (v_0) log(1/ε)+log(K0D/ϵ)+D≲log(1/δ)+log(1/ϵ)+logR, \; (1/ )+ (K_0 D/ε)+D\; \; (1/δ)+ (1/ε)+ R, ℬ(v0)≲ (v_0) maxK0, 1/δi, 1≲δ−(3κ1)∨1RD+N, \; \K_0,\,1/ _i,\,1 \\; \;δ^-(3 _1) 1R D+N, (v0)≲ (v_0) K0≲δ−(3κ1)∨1RD. \;K_0\; \;δ^-(3 _1) 1R D. Moreover v0(x,t)=0v_0(x,t)=0 for x∉[−(R+1),R+1]Dx∉[-(R+1),R+1]^D, since then χ(x)=0χ(x)=0, hence ψi(x,t)=0 _i(x,t)=0 for every i and Π(⋅,0)=0 (·,0)=0. Finally, for fixed t the map x↦Π(u(x),ψ(x))x (u(x),ψ(x)) has |∂uΠ|≤1| _u |≤ 1 and |∂ψΠ|≲K0| _ψ | K_0, so by (40) and Lip(χ)≲DLip(χ) D, Lip(v0(⋅,t))≲D3/2Lδ+DK0≲D5/2δ−3κ1R,t∈[0,1−δ].Lip (v_0(·,t) )\; \;D^3/2L_δ+DK_0\; \;D^5/2\,δ^-3 _1R, t∈[0,1-δ]. (41) ∎ Appendix C Proof of the Main Result C.1 Proof of Theorem 6 See 6 Proof. Throughout, R≍Dlog(1/η)R D (1/η) denotes the truncation radius of Lemma 9, and K denotes the output bound of the approximating network v0v_0 of that lemma, namely K≲δ−(3κ1)∨1RD≍δ−(3κ1)∨1Dlog(1/η).K\; \;δ^-(3 _1) 1R D\; \;δ^-(3 _1) 1D (1/η). We abbreviate c1:= 10κ1+2κ2,c_1\;:=\;10 _1+2 _2, the exponent appearing in Lemma 8. The ambient dimension D and the diffusion coefficient σ are held fixed, so factors of D and σ are absorbed by ≲ in the final rate; we nonetheless track those of D explicitly below. By Remark 1, the factor σ2σ^2 carried by Lemma 8 cancels against the σ−2σ^-2 of Lemma 7, so the discretization contribution below is in fact uniform in σ. Proof for estimator v^n v_n: Since conditioning on XtkX_t_k eliminates the cross term in ‖v(Xtk,tk)−X˙tk‖2E\|v(X_t_k,t_k)- X_t_k\|^2 exactly as in the proof of the conditional flow-matching identity (Lemma 11), so ℒ^n(v)=∑k=0N−1hk‖v(Xtk,tk)−v⋆(Xtk,tk)‖2+ℒ^n(v⋆) L_n(v)= _k=0^N-1h_kE\|v(X_t_k,t_k)-v_ (X_t_k,t_k)\|^2+ L_n(v_ ), the second term being independent of v. Hence v^n v_n, the minimizer of ℒ^n L_n over V, also minimizes the left-hand side below, and ∑k=0N−1hk‖v^n(Xtk,tk)−v⋆(Xtk,tk)‖2= _k=0^N-1h_kE \| v_n(X_t_k,t_k)-v_ (X_t_k,t_k) \|^2= infv∈∑k=0N−1hk‖v(Xtk,tk)−v⋆(Xtk,tk)‖2 _v _k=0^N-1h_kE \|v(X_t_k,t_k)-v_ (X_t_k,t_k) \|^2 ≤ ≤ ∑k=0N−1hk‖v0(Xtk,tk)−v⋆(Xtk,tk)‖2 _k=0^N-1h_kE \|v_0(X_t_k,t_k)-v_ (X_t_k,t_k) \|^2 = = ∑k=0N−1hk‖v0(Xtk,tk)−v⋆(Xtk,tk)‖2(‖Xtk‖≤R) _k=0^N-1h_kE \|v_0(X_t_k,t_k)-v_ (X_t_k,t_k) \|^2 1 (\|X_t_k\|≤ R ) +∑k=0N−1hk‖v0(Xtk,tk)−v⋆(Xtk,tk)‖2(∥Xtk∥>R) + _k=0^N-1h_kE \|v_0(X_t_k,t_k)-v_ (X_t_k,t_k) \|^2 1 (\|X_t_k\|>R ) ≲ ∑k=0N−1hkϵ2+∑k=0N−1hk(DK2η+δ−2κ1D2η1/2) _k=0^N-1h_kε^2+ _k=0^N-1h_k (DK^2η+δ^-2 _1D^2η^1/2 ) (42) ≲ ϵ2+D2δ−2((3κ1)∨1)η1/2. \ ε^2+D^2δ^-2 ((3 _1) 1 )η^1/2. (43) Here (42) follows from Lemmata 9 and 17, bounding the second summand by 2‖v0‖2(⋅)+2‖v⋆‖2(⋅)2E\|v_0\|^2 1(·)+2E\|v_ \|^2 1(·); the indicator confines XtkX_t_k to the Euclidean ball of radius R, which is contained in [−R,R]D[-R,R]^D, so Lemma 9 applies on the first event. Inequality (43) uses ∑khk≤1 _kh_k≤ 1 together with DK2η≲D2δ−2((3κ1)∨1)log(1/η)ηDK^2η D^2δ^-2((3 _1) 1) (1/η)\,η and log(1/η)η≤η1/2 (1/η)\,η≤η^1/2 for η≤η0η≤ _0, and 2κ1≤2((3κ1)∨1)2 _1≤ 2 ((3 _1) 1 ). From Lemmata 7 and 8, we observe that, KL(P1−δ∥P^1−δ(v^n))≲ (P_1-δ\| P_1-δ( v_n) ) ∑k=0N−1∫tktk+1at2‖v⋆(Xt,t)−v⋆(Xtk,tk)‖2t+∑k=0N−1∫tktk+1(at−atk)2‖v⋆(Xtk,tk)‖2t _k=0^N-1 _t_k^t_k+1a_t^2\,E\|v_ (X_t,t)-v_ (X_t_k,t_k)\|^2dt+ _k=0^N-1 _t_k^t_k+1(a_t-a_t_k)^2E\|v_ (X_t_k,t_k)\|^2dt +∑k=0N−1hkatk2∥v⋆(Xtk,tk)−v^n(Xtk,tk)∥2 + _k=0^N-1h_ka_t_k^2E\|v_ (X_t_k,t_k)- v_n(X_t_k,t_k)\|^2 ≲ D3δ−c1h+Dδ−2(κ1+κ2)h2+δ−κ1(ϵ2+D2δ−2((3κ1)∨1)η1/2) \ D^3δ^-c_1h+D\,δ^-2( _1+ _2)h^2+δ^- _1 (ε^2+D^2δ^-2 ((3 _1) 1 )η^1/2 ) ≲ D3δ−c1h+δ−κ1(ϵ2+D2δ−2((3κ1)∨1)η1/2), \ D^3δ^-c_1h+δ^- _1 (ε^2+D^2δ^-2 ((3 _1) 1 )η^1/2 ), where the last step using h≤1h≤ 1 and 2(κ1+κ2)≤c12( _1+ _2)≤ c_1, so that the second summand is absorbed by the first. Hence, from Lemma 21, p(Pdata,M(P^1−δ(v^n)))≲ _p(P_data, T_M( P_1-δ( v_n))) p(Pdata,Pnemp)+δκ3+M−12p+M(KL(P1−δ,P^1−δ(v^n)))12p \ W_p(P_data,P^emp_n)+δ _3+M^- 12p+M (KL (P_1-δ, P_1-δ( v_n) ) ) 12p ≲ p(Pdata,Pnemp)+δκ3+M−12p \ W_p(P_data,P^emp_n)+δ _3+M^- 12p +M(D3δ−c1h+δ−κ1(ϵ2+D2δ−2((3κ1)∨1)η1/2))12p +M (D^3δ^-c_1h+δ^- _1 (ε^2+D^2δ^-2 ((3 _1) 1 )η^1/2 ) ) 12p ≤ ≤ p(Pdata,Pnemp)+δκ3+M−12p \ W_p(P_data,P^emp_n)+δ _3+M^- 12p +M(D32pδ−c12ph12p+δ−κ12pϵ1/p+D1pδ−κ12p−(3κ1)∨1pη14p), +M (D 32pδ^- c_12ph 12p+δ^- _12pε^1/p+D 1pδ^- _12p- (3 _1) 1pη 14p ), (44) by the sub-additivity of u↦u1/2pu u^1/2p. We choose M=n2pd,δ=n−1dκ3,h≍n−2p(2p+1d+c12dpκ3)=n−2p(2p+1)d−10κ1+2κ2dκ3,ϵ=n−p(2p+1d+κ12dpκ3),M=n 2pd,\;δ=n^- 1d _3,\;h n^-2p ( 2p+1d+ c_12dp _3 )=n^- 2p(2p+1)d- 10 _1+2 _2d _3,\;ε=n^-p ( 2p+1d+ _12dp _3 ), and η=n−4p(2p+1d+κ1+2((3κ1)∨1)2dpκ3).η=n^-4p ( 2p+1d+ _1+2 ((3 _1) 1 )2dp _3 ). Each of the remaining terms of (44) is then O(n−1/d)O(n^-1/d). Hence, p(Pdata,M(P^1−δ(v^n)))≲p(Pdata,Pnemp)+n−1/d≲n−1/d+n−12plog12p(1/ξ), _p(P_data, T_M( P_1-δ( v_n))) _p(P_data,P^emp_n)+n^-1/d n^-1/d+n^- 12p 12p(1/ξ), with probability at least 1−ξ1-ξ, by Lemma 20. Proof for Estimator v^n,mmc v^mc_n,m: Write ϵc _c for the radius of the ℓ∞ _∞-cover of V appearing in Lemma 14, to distinguish it from the approximation accuracy ϵε of Lemma 9. From Lemma 14 we recall that, with probability at least 1−β1-β, ∑k=0N−1hk‖v^n,mmc(Xtk,tk)−v⋆(Xtk,tk)‖2 _k=0^N-1h_kE\| v^mc_n,m(X_t_k,t_k)-v_ (X_t_k,t_k)\|^2 ≲ ∑k=0N−1hk‖v0(Xtk,tk)−v⋆(Xtk,tk)‖2+DK2Λβm+DϵcK, _k=0^N-1h_kE\|v_0(X_t_k,t_k)-v_ (X_t_k,t_k)\|^2+DK^2 _βm+D\, _cK, where, Λβ=logN(1+(ϵc,,ℓ∞))β. _β= N (1+N ( _c;V, _∞ ) )β. We take m≍ϵc−2K2Λξ′,m _c^-2K^2\, _ξ , which makes the middle term O(DϵcK)O(D\, _cK). This ensures that with probability at least 1−ξ′1-ξ , ∑j=0N−1hj‖v^n,mmc(Xtj,tj)−v⋆(Xtj,tj)‖2≲ _j=0^N-1h_jE \| v^mc_n,m(X_t_j,t_j)-v_ (X_t_j,t_j) \|^2 ∑j=0N−1hj‖v0(Xtj,tj)−v⋆(Xtj,tj)‖2+DϵcK _j=0^N-1h_jE \|v_0(X_t_j,t_j)-v_ (X_t_j,t_j) \|^2+D\, _cK = = ∑j=0N−1hj‖v0(Xtj,tj)−v⋆(Xtj,tj)‖2(‖Xtj‖≤R) _j=0^N-1h_jE \|v_0(X_t_j,t_j)-v_ (X_t_j,t_j) \|^2 1 (\|X_t_j\|≤ R ) +∑j=0N−1hj‖v0(Xtj,tj)−v⋆(Xtj,tj)‖2(∥Xtj∥>R)+DϵcK + _j=0^N-1h_jE \|v_0(X_t_j,t_j)-v_ (X_t_j,t_j) \|^2 1 (\|X_t_j\|>R )+D\, _cK ≲ ∑k=0N−1hkϵ2+∑k=0N−1hk(DK2η+δ−2κ1D2η1/2)+DϵcK _k=0^N-1h_kε^2+ _k=0^N-1h_k (DK^2η+δ^-2 _1D^2η^1/2 )+D\, _cK (45) ≲ ϵ2+DϵcK+D2δ−2((3κ1)∨1)η1/2. \ ε^2+D\, _cK+D^2δ^-2 ((3 _1) 1 )η^1/2. (46) From Lemmata 7 and 8, exactly as above, KL(P1−δ∥P^1−δ(v^n,mmc))≲ (P_1-δ\| P_1-δ( v^mc_n,m) ) D3δ−c1h+δ−κ1(ϵ2+DϵcK+D2δ−2((3κ1)∨1)η1/2). \ D^3δ^-c_1h+δ^- _1 (ε^2+D\, _cK+D^2δ^-2 ((3 _1) 1 )η^1/2 ). Hence, from Lemma 21, p(Pdata,M(P^1−δ(v^n,mmc))) _p(P_data, T_M( P_1-δ( v^mc_n,m))) ≲ p(Pdata,Pnemp)+δκ3+M−12p+M(KL(P1−δ,P^1−δ(v^n,mmc)))12p \ W_p(P_data,P^emp_n)+δ _3+M^- 12p+M (KL (P_1-δ, P_1-δ( v^mc_n,m) ) ) 12p ≤ ≤ p(Pdata,Pnemp)+δκ3+M−12p \ W_p(P_data,P^emp_n)+δ _3+M^- 12p +M(D32pδ−c12ph12p+δ−κ12pϵ1p+D12pδ−κ12pϵc12pK12p+D1pδ−κ12p−(3κ1)∨1pη14p). +M (D 32pδ^- c_12ph 12p+δ^- _12pε 1p+D 12pδ^- _12p _c 12pK 12p+D 1pδ^- _12p- (3 _1) 1pη 14p ). (47) Keeping M,δ,h,ϵM,δ,h,ε and η as above and setting ϵc=ϵ2/K, _c\;=\;ε^2/K, the third term of (47) coincides with the second up to the fixed factor D1/2pD^1/2p, and every term is O(n−1/d)O(n^-1/d) by the computation above. Hence p(Pdata,M(P^1−δ(v^n,mmc)))≲p(Pdata,Pnemp)+n−1/d≲n−1/d+n−12plog12p(1/ξ), _p(P_data, T_M( P_1-δ( v^mc_n,m))) _p(P_data,P^emp_n)+n^-1/d n^-1/d+n^- 12p 12p(1/ξ), with probability at least 1−ξ−ξ′1-ξ-ξ , by Lemma 20. Replacing ξ and ξ′ξ with ξ/2ξ/2 gives the desired result. From Lemma 9, W≲NRDδ−3κ1Dϵ−D(log(1/δ)+log(1/ϵ)+logR)W NR^Dδ^-3 _1Dε^-D ( (1/δ)+ (1/ε)+ R ), with N=h−1N=h^-1. Substituting the choices above gives N=n2p(2p+1)d+20κ1+4κ22dκ3,δ−3κ1D=n3κ1Ddκ3=n6κ1D2dκ3,ϵ−D=nDp(2p+1)d+Dκ12dκ3.N=n 2p(2p+1)d+ 20 _1+4 _22d _3, δ^-3 _1D=n 3 _1Dd _3=n 6 _1D2d _3, ε^-D=n Dp(2p+1)d+ D _12d _3. Since η is polynomial in n, R≍Dlog(1/η)≍DlognR D (1/η) D n, so RD≍(Dlogn)D/2R^D (D n)^D/2 and log(1/δ)+log(1/ϵ)+logR≍logn (1/δ)+ (1/ε)+ R n. Summing exponents, the n-exponents combine as 2p(2p+1)d+Dp(2p+1)d=p(2p+1)(D+2)d 2p(2p+1)d+ Dp(2p+1)d= p(2p+1)(D+2)d and the κ-exponents as (20κ1+4κ2)+6κ1D+Dκ12dκ3=(20+7D)κ1+4κ22dκ3 (20 _1+4 _2)+6 _1D+D _12d _3= (20+7D) _1+4 _22d _3, so W≲np(2p+1)(D+2)d+(20+7D)κ1+4κ22dκ3(Dlogn)D/2logn.W\; \;n p(2p+1)(D+2)d+ (20+7D) _1+4 _22d _3 (D n )^D/2 n. Further, L≲log(1/δ)+log(1/ϵ)+logR≲lognL (1/δ)+ (1/ε)+ R n. From Lemma 9, B≲δ−(3κ1)∨1RD+NB δ^-(3 _1) 1R D+N, where δ−(3κ1)∨1=n(3κ1)∨1dκ3δ^-(3 _1) 1=n (3 _1) 1d _3 and N is as above, so B≲n(2p(2p+1)d+10κ1+2κ2dκ3)∨(3κ1)∨1dκ3Dlogn,B n ( 2p(2p+1)d+ 10 _1+2 _2d _3 ) (3 _1) 1d _3\,D n, and likewise K≲δ−(3κ1)∨1RD≲n(3κ1)∨1dκ3Dlogn,K δ^-(3 _1) 1R D n (3 _1) 1d _3\,D n, the latter unchanged. Finally, by Lemma 23, log(ϵc,,ℓ∞)≤Wlog(2LBL(W+1)L/ϵc)≲Wlog2n ( _c;V, _∞)≤ W (2 L B^L (W+1)^L/ _c ) W ^2n, since logL,logB,logW,log(1/ϵc)≲logn L, B, W, (1/ _c) n and L≲lognL n. With ϵc=ϵ2/K _c=ε^2/K we have ϵc−2K2=ϵ−4K4 _c^-2K^2=ε^-4K^4, whence m m ≍ϵ−4K4(Wlog2n+log(1/ξ)) ε^-4K^4 (W ^2n+ (1/ξ) ) ≍np(2p+1)(D+6)d+(24+7D)κ1+4κ2+4((6κ1)∨2)2dκ3log2n((Dlogn)D/2log3n+log(1/ξ)), n p(2p+1)(D+6)d+ (24+7D) _1+4 _2+4 ((6 _1) 2 )2d _3\; ^2n\; ((D n)^D/2 ^3n+ (1/ξ) ), using ϵ−4=n4p(2p+1)d+4κ12dκ3ε^-4=n 4p(2p+1)d+ 4 _12d _3 and K4≍n8((3κ1)∨1)2dκ3D4log2nK^4 n 8 ((3 _1) 1 )2d _3D^4 ^2n, together with 8((3κ1)∨1)=4((6κ1)∨2)8 ((3 _1) 1 )=4 ((6 _1) 2 ) and 4p(2p+1)d+p(2p+1)(D+2)d=p(2p+1)(D+6)d 4p(2p+1)d+ p(2p+1)(D+2)d= p(2p+1)(D+6)d. ∎ C.2 Supporting Lemmata Lemma 14. Let ϵc∈(0,1) _c∈(0,1), let v0∈v_0 be the network of Lemma 9, and write K for the output bound of ≡ℛ(W,L,B,K)V (W,L,B,K) and Λβ:=log(N(1+(ϵc,,ℓ∞))/β). _β\;:=\; ( N (1 + ( _c; , _∞) ) /β ). If m≥Λβm≥ _β, then with probability at least 1−β1-β over the Monte Carlo sampling, ∑k=0N−1hk‖v^n,mmc−v⋆‖2≲∑k=0N−1hk‖v0−v⋆‖2+DK2Λβm+DϵcK. _k=0^N-1h_kE\| v^mc_n,m-v_ \|^2\; \; _k=0^N-1h_kE\|v_0-v_ \|^2+DK^2 _βm+D\, _cK. (48) Proof. All expectations are conditional on x1,…,xnx_1,…,x_n. Fix k∈0,1,…,N−1k∈\0,1,…,N-1\ and v∈v . Since the Monte Carlo sample is i.i.d., (1m∑j=1m‖v(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2)=‖v(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2 ( 1m _j=1^m \|v(X_t_k^(j),t_k)-( α_t_kX^(j)_1+ β_t_kX_0^(j)) \|^2 )=E \|v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0) \|^2 = = 1n∑i=1nX0‖v(Xi,tk,tk)−(α˙tkxi+β˙tkX0)‖2, 1n _i=1^nE_X_0\|v(X_i,t_k,t_k)-( α_t_kx_i+ β_t_kX_0)\|^2, where Xi,t=αtxi+βtX0X_i,t= _tx_i+ _tX_0 and X0∼P0X_0 P_0. By Bernstein’s inequality for sub-exponential variables (Vershynin, 2026, Theorem 2.9.1), for any fixed k and v, with probability at least 1−β1-β, |1m∑j=1m‖v(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2−‖v(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2| | 1m _j=1^m \|v(X_t_k^(j),t_k)-( α_t_kX^(j)_1+ β_t_kX_0^(j)) \|^2-E\|v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 | ≲ (‖v(Xtk,tk)‖2‖ψ1+‖α˙tkX1+β˙tkX0‖2‖ψ1)(log(1/β)m∨log(1/β)m). ( \|\|v(X_t_k,t_k)\|^2 \|_ _1+ \|\| α_t_kX_1+ β_t_kX_0\|^2 \|_ _1 ) ( (1/β)m\; \; (1/β)m ). (49) Since (v)≤KK(v)≤ K bounds ‖v(⋅,⋅)‖∞\|v(·,·)\|_∞ coordinatewise, ‖v(Xtk,tk)‖2≤DK2\|v(X_t_k,t_k)\|^2≤ DK^2 almost surely, so ‖v(Xtk,tk)‖2‖ψ1≲DK2. \|\|v(X_t_k,t_k)\|^2 \|_ _1\; \;DK^2. Moreover, by D1(a), |α˙t|∨|β˙t|≲δ−κ1| α_t| | β_t| δ^- _1 for t∈[0,1−δ]t∈[0,1-δ], so ‖α˙tkX1+β˙tkX0‖ψ2≲δ−κ1D \|\| α_t_kX_1+ β_t_kX_0\| \|_ _2 δ^- _1 D and hence ‖α˙tkX1+β˙tkX0‖2‖ψ1=‖α˙tkX1+β˙tkX0‖ψ22≲Dδ−2κ1≲DK2, \|\| α_t_kX_1+ β_t_kX_0\|^2 \|_ _1\;=\; \|\| α_t_kX_1+ β_t_kX_0\| \|_ _2^2\; \;Dδ^-2 _1\; \;DK^2, the last step because K≳δ−(3κ1)∨1≥δ−κ1K δ^-(3 _1) 1≥δ^- _1. Hence from (49), for fixed k and v, with probability at least 1−β1-β, |1m∑j=1m‖v(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2−‖v(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2| | 1m _j=1^m \|v(X_t_k^(j),t_k)-( α_t_kX^(j)_1+ β_t_kX_0^(j)) \|^2-E\|v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 | ≲ DK2log(1/β)m. \;DK^2 (1/β)m. (50) Let (ϵc,,ℓ∞)C ( _c;V, _∞ ) be an ϵc _c-cover of V in the ℓ∞ _∞ norm, of cardinality (ϵc,,ℓ∞)N( _c;V, _∞). Applying (50) with β replaced by β/(N(1+(ϵc,,ℓ∞)))β/ (N(1+N( _c;V, _∞)) ) and taking a union bound over the N time indices and the (ϵc,,ℓ∞)+1N( _c;V, _∞)+1 networks in (ϵc,,ℓ∞)∪v0C( _c;V, _∞)∪\v_0\, we obtain: with probability at least 1−β1-β, simultaneously for all v∈(ϵc,,ℓ∞)∪v0v ( _c;V, _∞ )∪\v_0\, ∑k=0N−1hk|1m∑j=1m‖v(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2−‖v(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2| _k=0^N-1h_k | 1m _j=1^m \|v(X_t_k^(j),t_k)-( α_t_kX^(j)_1+ β_t_kX_0^(j)) \|^2-E\|v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 | ≲ DK2Λβm, \;DK^2 _βm, (51) using ∑khk≤1 _kh_k≤ 1. We work on this event for the remainder of the proof. Let v~∈(ϵc,,ℓ∞) v ( _c;V, _∞) be such that ‖v^n,mmc−v~‖∞≤ϵc\| v^mc_n,m- v\|_∞≤ _c. Using ‖a‖2−‖b‖2=⟨a−b,a+b⟩\|a\|^2-\|b\|^2= a-b,a+b , for any k and any x, |‖v^n,mmc(x,tk)−(α˙tkX1+β˙tkX0)‖2−‖v~(x,tk)−(α˙tkX1+β˙tkX0)‖2| |\| v^mc_n,m(x,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2-\| v(x,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 | ≤ ≤\; D‖v^n,mmc−v~‖∞(‖v^n,mmc‖∞+‖v~‖∞+2‖α˙tkX1+β˙tkX0‖) D\,\| v^mc_n,m- v\|_∞ (\| v^mc_n,m\|_∞+\| v\|_∞+2\| α_t_kX_1+ β_t_kX_0\| ) ≤ ≤\; 2Dϵc(DK+‖α˙tkX1+β˙tkX0‖). 2 D\, _c ( DK+\| α_t_kX_1+ β_t_kX_0\| ). (52) Taking expectations in (52) and using ‖α˙tkX1+β˙tkX0‖≲δ−κ1D≲DKE\| α_t_kX_1+ β_t_kX_0\| δ^- _1 D DK gives, for every k, |‖v^n,mmc(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2−‖v~(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2|≲DϵcK, |E\| v^mc_n,m(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2-E\| v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 |\; \;D\, _cK, and the same bound holds with E replaced by the empirical average 1m∑j=1m 1m _j=1^m, since (52) is pointwise in the sample. Chaining these with (51) applied first to v~ v and then to v0v_0, and using that v^n,mmc v^mc_n,m minimises the Monte Carlo objective over ∋v0V v_0, ∑k=0N−1hk‖v^n,mmc(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2 _k=0^N-1h_kE\| v^mc_n,m(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 ≤ ≤\; ∑k=0N−1hk‖v~(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2+(DϵcK) _k=0^N-1h_kE\| v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2+O (D _cK ) ≤ ≤\; 1m∑j=1m∑k=0N−1hk‖v~(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2+(DK2Λβm)+(DϵcK) 1m _j=1^m _k=0^N-1h_k\| v(X_t_k^(j),t_k)-( α_t_kX_1^(j)+ β_t_kX_0^(j))\|^2+O (DK^2 _βm )+O (D _cK ) ≤ ≤\; 1m∑j=1m∑k=0N−1hk‖v^n,mmc(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2+(DK2Λβm)+(DϵcK) 1m _j=1^m _k=0^N-1h_k\| v^mc_n,m(X_t_k^(j),t_k)-( α_t_kX_1^(j)+ β_t_kX_0^(j))\|^2+O (DK^2 _βm )+O (D _cK ) ≤ ≤\; 1m∑j=1m∑k=0N−1hk‖v0(Xtk(j),tk)−(α˙tkX1(j)+β˙tkX0(j))‖2+(DK2Λβm)+(DϵcK) 1m _j=1^m _k=0^N-1h_k\|v_0(X_t_k^(j),t_k)-( α_t_kX_1^(j)+ β_t_kX_0^(j))\|^2+O (DK^2 _βm )+O (D _cK ) ≤ ≤\; ∑k=0N−1hk‖v0(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2+(DK2Λβm)+(DϵcK). _k=0^N-1h_kE\|v_0(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2+O (DK^2 _βm )+O (D _cK ). Finally, conditioning on XtkX_t_k makes the cross term vanish, so that for every v, ∑k=0N−1hk‖v(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2 _k=0^N-1h_kE\|v(X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2 = = ∑k=0N−1hk‖v(Xtk,tk)−v⋆(Xtk,tk)‖2+∑k=0N−1hk‖v⋆(Xtk,tk)−(α˙tkX1+β˙tkX0)‖2, _k=0^N-1h_kE\|v(X_t_k,t_k)-v_ (X_t_k,t_k)\|^2+ _k=0^N-1h_kE\|v_ (X_t_k,t_k)-( α_t_kX_1+ β_t_kX_0)\|^2, this being the discrete analogue of Lemma 11. Applying it to v=v^n,mmcv= v^mc_n,m and to v=v0v=v_0 and subtracting the common second term, which does not depend on v, yields ∑k=0N−1hk‖v^n,mmc−v⋆‖2≲∑k=0N−1hk‖v0−v⋆‖2+DK2Λβm+DϵcK.∎ _k=0^N-1h_kE\| v^mc_n,m-v_ \|^2 _k=0^N-1h_kE\|v_0-v_ \|^2+DK^2 _βm+D _cK. Appendix D Supporting Results D.1 Results for the Velocity Field Lemma 15 (Spatial Derivative of v⋆v_ ). The following holds for any x∈ℝDx∈ R ^D and t>0t>0, v⋆(x,t)= v_ (x,t)= β˙tβtx+α˙tβt−β˙tαtβt(X1|Xt=x)=α˙tαtx+βt2(α˙tαt−β˙tβt)∇logϱt(x), β_t _tx+ α_t _t- β_t _t _tE(X_1|X_t=x)= α_t _tx+ _t^2 ( α_t _t- β_t _t )∇ _t(x), (53) ∇v⋆(x,t)= ∇ v_ (x,t)= αt2βt2(α˙tαt−β˙tβt)Var(X1|Xt=x)+β˙tβtID=12∂t(αt2βt2)Var(X1|Xt=x)+β˙tβtID. _t^2 _t^2 ( α_t _t- β_t _t )Var(X_1|X_t=x)+ β_t _tI_D= 12 _t ( _t^2 _t^2 )Var(X_1|X_t=x)+ β_t _tI_D. (54) Proof. Throughout, fix t∈(0,1)t∈(0,1) with αt≠0,βt>0 _t≠ 0,\ _t>0, and recall the interpolant is given by, Xt=αtX1+βtX0,X1∼Pnemp,X0∼N(0,ID)independent of X1.X_t= _tX_1+ _tX_0,\,X_1 P^emp_n,\,X_0 N(0,I_D)\ independent of X_1. Hence, ϱt _t, the law of XtX_t, is the Gaussian smoothing ϱt(x)=∫ℝDϕ(x,αty,βt2ID)Pnemp(y), _t(x)= _ R ^Dφ(x; _ty, _t^2I_D)\,P^emp_n(dy), where, ϕ(x;μ,Σ):=(2π)−D/2det(Σ)−1/2exp(−12(x−μ)⊤Σ−1(x−μ)).φ(x;μ, ):=(2π)^-D/2 ( )^-1/2 \! (- 12(x-μ) ^-1(x-μ) ). By definition, v⋆(x,t)=[X˙t∣Xt=x]v_ (x,t)=E[ X_t X_t=x]. For notational simplicity, we write m(x,t):=[X1∣Xt=x]∈ℝD,Σ(x,t):=Var(X1∣Xt=x)∈ℝD×D.m(x,t):=E[X_1 X_t=x]∈ R ^D, (x,t):=Var(X_1 X_t=x)∈ R ^D× D. Since βt>0 _t>0, we may write X0=βt−1(Xt−αtX1)X_0= _t^-1(X_t- _tX_1). Hence v⋆(x,t)=[α˙tX1+β˙tX0∣Xt=x]=α˙tm(x,t)+β˙tβt(x−αtm(x,t)),v_ (x,t)=E[ α_tX_1+ β_tX_0 X_t=x]= α_t\,m(x,t)+ β_t _t (x- _tm(x,t) ), using [Xt∣Xt=x]=xE[X_t X_t=x]=x. Collecting terms in m(x,t)m(x,t) gives v⋆(x,t)=β˙tβtx+α˙tβt−β˙tαtβtm(x,t).v_ (x,t)= β_t _t\,x+ α_t _t- β_t _t _t\,m(x,t). (55) To prove equation (54), we note that, ∇xϕ(x,αty,βt2ID)=−βt−2(x−αty)ϕ(x,αty,βt2ID) _xφ(x; _ty, _t^2I_D)=- _t^-2(x- _ty)\,φ(x; _ty, _t^2I_D), differentiation under the integral sign (justified since βt>0 _t>0 and PnempP^emp_n has finite second moments) gives ∇ϱt(x)=−1βt2∫(x−αty)ϕ(x;αty,βt2ID)Pnemp(dy).∇ _t(x)=- 1 _t^2 (x- _ty)\,φ(x; _ty, _t^2I_D)\,P^emp_n(dy). Dividing by ϱt(x)>0 _t(x)>0 and using that ϕ(x,αty,βt2ID)Pnemp(dy)/ϱt(x)φ(x; _ty, _t^2I_D)P^emp_n(dy)/ _t(x) is exactly the conditional law of X1X_1 given Xt=xX_t=x, ∇logϱt(x)=−1βt2[x−αtX1∣Xt=x]=−x−αtm(x,t)βt2.∇ _t(x)=- 1 _t^2\,E[x- _tX_1 X_t=x]=- x- _tm(x,t) _t^2. (56) Solving for m(x,t)m(x,t) (using αt≠0 _t≠ 0): m(x,t)=x+βt2∇logϱt(x)αt.m(x,t)= x+ _t^2∇ _t(x) _t. (57) Substituting (57) into (53), v⋆(x,t)=β˙tβtx+α˙tβt−β˙tαtαtβtx+(α˙tβt−β˙tαt)βtαt∇logϱt(x).v_ (x,t)= β_t _tx+ α_t _t- β_t _t _t _t\,x+ ( α_t _t- β_t _t) _t _t\,∇ _t(x). The coefficient of x simplifies as β˙tβt+α˙tβt−β˙tαtαtβt=β˙tαt+α˙tβt−β˙tαtαtβt=α˙tαt, β_t _t+ α_t _t- β_t _t _t _t= β_t _t+ α_t _t- β_t _t _t _t= α_t _t, and the coefficient of ∇logϱt(x)∇ _t(x) as βt(α˙tβt−β˙tαt)αt=βt2⋅α˙tβt−β˙tαtαtβt=βt2(α˙tαt−β˙tβt). _t( α_t _t- β_t _t) _t= _t^2· α_t _t- β_t _t _t _t= _t^2 ( α_t _t- β_t _t ). Hence v⋆(x,t)=α˙tαtx+βt2(α˙tαt−β˙tβt)∇logϱt(x),v_ (x,t)= α_t _tx+ _t^2 ( α_t _t- β_t _t )∇ _t(x), (58) which together with (53) establishes the first display of the lemma. To complete the proof, we make the following claim, ∇xm(x,t)=αtβt2Σ(x,t). _xm(x,t)= _t _t^2\, (x,t). (59) Write Ni(x):=∫yiϕ(x,αty,βt2ID)Pnemp(y)N_i(x):= y_i\,φ(x; _ty, _t^2I_D)\,P^emp_n(dy), so that mi(x,t)=Ni(x)/ϱt(x)m_i(x,t)=N_i(x)/ _t(x) for i=1,…,Di=1,…,D. Differentiating under the integral sign as before, ∂xjNi(x)= _x_jN_i(x)= −1βt2∫yi(xj−αtyj)ϕ(x;αty,βt2ID)Pnemp(dy) - 1 _t^2 y_i(x_j- _ty_j)\,φ(x; _ty, _t^2I_D)\,P^emp_n(dy) = = −xjβt2Ni(x)+αtβt2ϱt(x)[X1,iX1,j∣Xt=x]. - x_j _t^2N_i(x)+ _t _t^2\, _t(x)\,E[X_1,iX_1,j X_t=x]. Dividing by ϱt(x) _t(x), ∂xjNi(x)ϱt(x)=−xjβt2mi(x,t)+αtβt2[X1,iX1,j∣Xt=x]. _x_jN_i(x) _t(x)=- x_j _t^2\,m_i(x,t)+ _t _t^2\,E[X_1,iX_1,j X_t=x]. By the quotient rule, ∂xjmi(x,t)=∂xjNi(x)ϱt(x)−mi(x,t)∂xjlogϱt(x) _x_jm_i(x,t)= _x_jN_i(x) _t(x)-m_i(x,t)\, _x_j _t(x). Substituting the previous display and (2) (i.e. ∂xjlogϱt(x)=−βt−2(xj−αtmj(x,t)) _x_j _t(x)=- _t^-2(x_j- _tm_j(x,t))), ∂xjmi(x,t)=αtβt2([X1,iX1,j∣Xt=x]−mi(x,t)mj(x,t))=αtβt2Σij(x,t), _x_jm_i(x,t)= _t _t^2 (E[X_1,iX_1,j X_t=x]-m_i(x,t)m_j(x,t) )= _t _t^2\, _ij(x,t), which is (59) in coordinates. Differentiate (53) with respect to x and apply (59): ∇v⋆(x,t)=β˙tβtID+α˙tβt−β˙tαtβt∇xm(x,t)=β˙tβtID+αt(α˙tβt−β˙tαt)βt3Σ(x,t).∇ v_ (x,t)= β_t _tI_D+ α_t _t- β_t _t _t\, _xm(x,t)= β_t _tI_D+ _t( α_t _t- β_t _t) _t^3\, (x,t). Since αt(α˙tβt−β˙tαt)βt3=αt2βt2⋅α˙tβt−β˙tαtαtβt=αt2βt2(α˙tαt−β˙tβt), _t( α_t _t- β_t _t) _t^3= _t^2 _t^2· α_t _t- β_t _t _t _t= _t^2 _t^2 ( α_t _t- β_t _t ), we obtain ∇v⋆(x,t)=αt2βt2(α˙tαt−β˙tβt)Var(X1∣Xt=x)+β˙tβtID.∇ v_ (x,t)= _t^2 _t^2 ( α_t _t- β_t _t )Var(X_1 X_t=x)+ β_t _tI_D. (60) By the quotient rule, ∂t(αt2βt2)=2αtα˙tβt2−2αt2βtβ˙tβt4=2αtβt(α˙tβt−αtβ˙t)βt4=2αt2βt2(α˙tαt−β˙tβt). _t\! ( _t^2 _t^2 )= 2 _t α_t _t^2-2 _t^2 _t β_t _t^4= 2 _t _t( α_t _t- _t β_t) _t^4=2\, _t^2 _t^2 ( α_t _t- β_t _t ). Hence 12∂t(αt2βt2)=αt2βt2(α˙tαt−β˙tβt), 12\, _t\! ( _t^2 _t^2 )= _t^2 _t^2 ( α_t _t- β_t _t ), and substituting into (60) gives ∇v⋆(x,t)=12∂t(αt2βt2)Var(X1∣Xt=x)+β˙tβtID.∇ v_ (x,t)= 12\, _t\! ( _t^2 _t^2 )Var(X_1 X_t=x)+ β_t _tI_D. ∎ Corollary 16. For any t∈[0,1−δ]t∈[0,1-δ], ‖v⋆(x,t)‖≲ \|v_ (x,t)\| δ−κ1‖x‖+δ−κ1D, δ^- _1\|x\|+δ^- _1 D, (61) ‖∇v⋆(x,t)‖≲ \|∇ v_ (x,t)\| Dδ−3κ1. Dδ^-3 _1. (62) Proof. By construction, for any t∈[0,1−δ]t∈[0,1-δ], ‖v⋆(x,t)‖= \|v_ (x,t)\|= ‖β˙tβtx+α˙tβt−β˙tαtβt(X1|Xt=x)‖ \| β_t _tx+ α_t _t- β_t _t _tE(X_1|X_t=x) \| ≲ δ−κ1‖x‖+δ−κ1D. δ^- _1\|x\|+δ^- _1 D. Similarly, ‖∇v⋆(x,t)‖= \|∇ v_ (x,t)\|= ‖αt2βt2(α˙tαt−β˙tβt)Var(X1|Xt=x)+β˙tβtID‖ \| _t^2 _t^2 ( α_t _t- β_t _t )Var(X_1|X_t=x)+ β_t _tI_D \| ≲ Dβt2|α˙tβt−β˙tαtβt|+|β˙tβt| D _t^2 | α_t _t- β_t _t _t |+ | β_t _t | ≲ Dδ−3κ1. Dδ^-3 _1. ∎ Lemma 17. There exists an absolute constant C such that, for every η≤η0η≤ _0, writing R=CDlog(1/η)R=C D (1/η), we have for every t∈[0,1−δ]t∈[0,1-δ] ℙ(‖Xt‖>R)≤ηand‖v⋆(Xt,t)‖2(‖Xt‖>R)≲δ−2κ1D2η1/2.P (\|X_t\|>R )\;≤\;η \|v_ (X_t,t)\|^2 1 (\|X_t\|>R )\; \;δ^-2 _1D^2η^1/2. Proof. Fix t∈[0,1−δ]t∈[0,1-δ]. By Theorem 13, Law(Xt)=ϱt=Law(αtX1+βtX0)Law(X_t)= _t=Law( _tX_1+ _tX_0), a statement about the one-time marginal, which is all that is used below. Since X1∼PnempX_1 P^emp_n is supported in [0,1]D[0,1]^D by Assumption 1 we have ‖X1‖≤D\|X_1\|≤ D, and since αt,βt∈[0,1] _t, _t∈[0,1], ‖Xt‖≤αt‖X1‖+βt‖X0‖≤D+‖X0‖,X0∼(0,ID).\|X_t\|\;≤\; _t\|X_1\|+ _t\|X_0\|\;≤\; D+\|X_0\|, X_0 (0,I_D). (63) Hence ‖Xt‖ψ2≤D+‖X0‖ψ2≲D \|\,\|X_t\|\, \|_ _2≤ D+ \|\,\|X_0\|\, \|_ _2 D, and by Proposition 2.6.1 of Vershynin (2026), ‖Xt‖4≲D2E\|X_t\|^4 D^2. Moreover, by (63) and Gaussian concentration of the norm, ℙ(∥Xt∥≥2D+s)≤ℙ(∥X0∥≥D+s)≤e−s2/2P (\|X_t\|≥ 2 D+s ) (\|X_0\|≥ D+s )≤ e^-s^2/2 for every s≥0s≥ 0; taking s=2log(1/η)s= 2 (1/η) gives the first claim, for a suitable absolute constant C and η≤η0η≤ _0. For the second claim, by Corollary 16 and then the Cauchy–Schwarz inequality, ‖v⋆(Xt,t)‖2(‖Xt‖>R) \|v_ (X_t,t)\|^2 1 (\|X_t\|>R ) ≲ \; δ−2κ1(‖Xt‖2+D) 1(‖Xt‖>R) δ^-2 _1\,E (\|X_t\|^2+D )\, 1 (\|X_t\|>R ) (64) ≤ ≤\; δ−2κ1((‖Xt‖2+D)2)1/2(ℙ(‖Xt‖>R))1/2 δ^-2 _1 (E (\|X_t\|^2+D )^2 )^1/2 (P (\|X_t\|>R ) )^1/2 ≲ \; δ−2κ1(‖Xt‖4+D2)1/2η1/2 δ^-2 _1 (E\|X_t\|^4+D^2 )^1/2\,η^1/2 ≲ \; δ−2κ1Dη1/2≤δ−2κ1D2η1/2. δ^-2 _1D\,η^1/2\;≤\;δ^-2 _1D^2\,η^1/2. In the above calculations, (64) follows from Corollary 16. ∎ Lemma 18 (Temporal velocity of v⋆v_ ). Fix t∈(0,1)t∈(0,1) with αt≠0 _t≠ 0, βt>0 _t>0, and suppose αt,βt _t, _t are twice continuously differentiable in a neighborhood of t. Define ut:=α˙tαt,wt:=β˙tβt,κt:=ut−wt=α˙tαt−β˙tβt.u_t:= α_t _t, w_t:= β_t _t, _t:=u_t-w_t= α_t _t- β_t _t. Then, with E1,t(x)=(X1∣Xt=x)E_1,t(x)=E(X_1 X_t=x), E2,t(x)=(‖X1‖2∣Xt=x)E_2,t(x)=E(\|X_1\|^2 X_t=x), E3,t(x)=(‖X1‖2X1∣Xt=x)E_3,t(x)=E(\|X_1\|^2X_1 X_t=x), and V1,t(x)=Var(X1∣Xt=x)V_1,t(x)=Var(X_1 X_t=x) as in Lemma 15, ∂tv⋆(x,t) _tv_ (x,t) =w˙tx−αt2κt(2wt−ut)βt2V1,t(x)x+αt(u˙t−w˙t+utκt)E1,t(x) = w_t\,x\;-\; _t^2 _t(2w_t-u_t) _t^2\,V_1,t(x)\,x\;+\; _t ( u_t- w_t+u_t _t )\,E_1,t(x) −αt3κt2βt2(E3,t(x)−E1,t(x)E2,t(x)). - _t^3 _t^2 _t^2 (E_3,t(x)-E_1,t(x)E_2,t(x) ). (65) Further, for any t∈[0,1−δ]t∈[0,1-δ] and x∈ℝDx∈ R ^D, ‖∂tv⋆(x,t)‖≲δ−κ1(1+‖x‖)\| _tv_ (x,t)\| δ^- _1(1+\|x\|). Proof. Throughout, write E1,E2,E3,V1E_1,E_2,E_3,V_1 for E1,t(x),E2,t(x),E3,t(x),V1,t(x)E_1,t(x),E_2,t(x),E_3,t(x),V_1,t(x), and g(x):=∇xlogϱt(x)g(x):= _x _t(x). From Lemma 19, taking f(y)=yf(y)=y recovers ∇xE1=αtβt2V1 _xE_1= _t _t^2V_1; taking f(y)=‖y‖2f(y)=\|y\|^2 gives ∇xE2=αtβt2(E3−E2E1) _xE_2= _t _t^2(E_3-E_2E_1). Since tr(V1)=E2−‖E1‖2tr(V_1)=E_2-\|E_1\|^2, ∇xtr(V1,t(x))=∇xE2−2(∇xE1)⊤E1=αtβt2(E3−E2E1−2V1E1). _xtr (V_1,t(x) )= _xE_2-2( _xE_1) E_1= _t _t^2 (E_3-E_2E_1-2V_1E_1 ). (66) The continuity equation ∂tϱt+∇x⋅(ϱtv⋆)=0 _t _t+ _x·( _tv_ )=0 gives, for fixed x, ∂tg(x)=−[∇x2logϱt(x)v⋆(x,t)+∇xv⋆(x,t)g(x)+∇x(∇x⋅v⋆(x,t))]. _tg(x)=- [ _x^2 _t(x)\,v_ (x,t)+ _xv_ (x,t)\,g(x)+ _x ( _x· v_ (x,t) ) ]. (67) Differentiating g(x)=−x−αtE1βt2g(x)=- x- _tE_1 _t^2 in x and using ∇xE1=αtβt2V1 _xE_1= _t _t^2V_1, ∇x2logϱt(x)=αt2βt4V1−1βt2ID=:H. _x^2 _t(x)= _t^2 _t^4V_1- 1 _t^2I_D=:H. (68) Using (58), (54), and (68), Hv⋆(x,t) H\,v_ (x,t) =(αt2βt4V1−1βt2ID)(utx+qtg(x))=αt2utβt4V1x+αt2qtβt4V1g(x)−utβt2x−qtβt2g(x), = ( _t^2 _t^4V_1- 1 _t^2I_D ) (u_tx+q_tg(x) )= _t^2u_t _t^4V_1x+ _t^2q_t _t^4V_1g(x)- u_t _t^2x- q_t _t^2g(x), ∇xv⋆(x,t)g(x) _xv_ (x,t)\,g(x) =(αt2βt2κtV1+wtID)g(x)=αt2κtβt2V1g(x)+wtg(x). = ( _t^2 _t^2 _tV_1+w_tI_D )g(x)= _t^2 _t _t^2V_1g(x)+w_tg(x). Since qt=βt2κtq_t= _t^2 _t, the coefficient of V1g(x)V_1g(x) in the sum is αt2qtβt4+αt2κtβt2=2αt2κtβt2 _t^2q_t _t^4+ _t^2 _t _t^2=2 _t^2 _t _t^2, and the coefficient of g(x)g(x) is wt−qtβt2=wt−κt=2wt−utw_t- q_t _t^2=w_t- _t=2w_t-u_t. Substituting g(x)=−x−αtE1βt2g(x)=- x- _tE_1 _t^2 throughout and collecting terms by type (V1xV_1x, x, V1E1V_1E_1, E1E_1), Hv⋆(x,t)+∇xv⋆(x,t)g(x)=αt2(2wt−ut)βt4V1x−2wtβt2x+2αt3κtβt4V1E1+αt(2wt−ut)βt2E1. Hv_ (x,t)+ _xv_ (x,t)g(x)= _t^2(2w_t-u_t) _t^4V_1x- 2w_t _t^2x+ 2 _t^3 _t _t^4V_1E_1+ _t(2w_t-u_t) _t^2E_1. (69) From (54), ∇x⋅v⋆(x,t)=αt2βt2κttr(V1)+Dwt _x· v_ (x,t)= _t^2 _t^2 _ttr(V_1)+Dw_t, so by (66), ∇x(∇x⋅v⋆(x,t))=αt3κtβt4(E3−E2E1−2V1E1). _x ( _x· v_ (x,t) )= _t^3 _t _t^4 (E_3-E_2E_1-2V_1E_1 ). (70) Adding (70) to (69), the V1E1V_1E_1 contributions cancel exactly (2αt3κtβt4−2αt3κtβt4=0 2 _t^3 _t _t^4- 2 _t^3 _t _t^4=0), leaving ∂tg(x)=−αt2(2wt−ut)βt4V1x+2wtβt2x−αt(2wt−ut)βt2E1−αt3κtβt4(E3−E2E1). _tg(x)=- _t^2(2w_t-u_t) _t^4V_1x+ 2w_t _t^2x- _t(2w_t-u_t) _t^2E_1- _t^3 _t _t^4(E_3-E_2E_1). By the product rule, ∂tv⋆(x,t)=u˙tx+q˙tg(x)+qt∂tg(x),q˙t=βt2(2wtκt+κ˙t),κ˙t=u˙t−w˙t. _tv_ (x,t)= u_t\,x+ q_t\,g(x)+q_t\, _tg(x), q_t= _t^2 (2w_t _t+ κ_t ), κ_t= u_t- w_t. Substituting g(x)=−x−αtE1βt2g(x)=- x- _tE_1 _t^2 and qt∂tg(x)q_t _tg(x), then collecting terms of each type (x; V1xV_1x; E1E_1; E3−E2E1E_3-E_2E_1), we get the final result. Finally, ‖∂tv⋆(x,t)‖≤ \| _tv_ (x,t)\|≤ |w˙t|‖x‖+|αt2κt(2wt−ut)βt2|‖V1,t(x)‖‖x‖+|αt(u˙t−w˙t+utκt)|‖E1,t(x)‖ | w_t|\|x\|+ | _t^2 _t(2w_t-u_t) _t^2 | \|V_1,t(x) \|\|x\|+ | _t ( u_t- w_t+u_t _t ) |\|E_1,t(x)\| +|αt3κt2βt2|‖E3,t(x)−E1,t(x)E2,t(x)‖ + | _t^3 _t^2 _t^2 |\|E_3,t(x)-E_1,t(x)E_2,t(x)\| ≲ δ−κ1(1+‖x‖). δ^- _1(1+\|x\|). ∎ Lemma 19. Fix t∈(0,1)t∈(0,1) with αt≠0 _t≠ 0, βt>0 _t>0. For any measurable f:ℝD→ℝf: R ^D→ R with |f(X1)|<∞E|f(X_1)|<∞, define Ef(x,t):=(f(X1)∣Xt=x)E_f(x,t):=E (f(X_1) X_t=x ). Then, for every x∈ℝDx∈ R ^D, ∇xEf(x,t)=αtβt2Cov(f(X1),X1∣Xt=x). _xE_f(x,t)= _t _t^2\,Cov (f(X_1),X_1 X_t=x ). Proof. Write Nf(x):=ϱt(x)Ef(x,t)=∫f(y)ϕ(x,αty,βt2ID)Pnemp(y)N_f(x):= _t(x)\,E_f(x,t)= f(y)\,φ(x; _ty, _t^2I_D)\,P^emp_n(dy). Differentiating under the integral sign, ∇xNf(x)= _xN_f(x)= ∫f(y)(−x−αtyβt2)ϕ(x,αty,βt2ID)Pnemp(y) f(y) (- x- _ty _t^2 )φ(x; _ty, _t^2I_D)\,P^emp_n(dy) = = −xβt2Nf(x)+αtβt2∫f(y)yϕ(x,αty,βt2ID)Pnemp(y). - x _t^2N_f(x)+ _t _t^2 f(y)\,y\,φ(x; _ty, _t^2I_D)\,P^emp_n(dy). Dividing by ϱt(x) _t(x) and applying the quotient rule together with g(x):=∇xlogϱt(x)=−x−αtE1(x,t)βt2g(x):= _x _t(x)=- x- _tE_1(x,t) _t^2 (Lemma 15), ∇xEf(x,t) _xE_f(x,t) =∇xNf(x)ϱt(x)−Ef(x,t)g(x) = _xN_f(x) _t(x)-E_f(x,t)\,g(x) =−xβt2Ef(x,t)+αtβt2(f(X1)X1∣Xt=x)+x−αtE1(x,t)βt2Ef(x,t) =- x _t^2E_f(x,t)+ _t _t^2E(f(X_1)X_1 X_t=x)+ x- _tE_1(x,t) _t^2E_f(x,t) =αtβt2((f(X1)X1∣Xt=x)−Ef(x,t)E1(x,t)) = _t _t^2 (E(f(X_1)X_1 X_t=x)-E_f(x,t)\,E_1(x,t) ) =αtβt2Cov(f(X1),X1∣Xt=x). = _t _t^2Cov (f(X_1),X_1 X_t=x ). (71) ∎ D.2 Results for the Wasserstein Distance Lemma 20. For any d>dp∗(Pdata)d>d _p (P_data ), with probability at least 1−ξ1-ξ, p(Pnemp,Pdata)≤n−1/d+n−12pD1/2log12p(1/ξ).W_p(P^emp_n,P_data)≤ n^-1/d+n^- 12pD^1/2 12p(1/ξ). Proof. Let fc(y)=supx∈[0,1]D(f(x)−‖x−y‖p)f_c(y)= _x∈[0,1]^D (f(x)-\|x-y\|^p ) denote the c-transform of f. It can be shown that f(x)=infy∈[0,1]D(fc(y)+‖x−y‖p)f(x)= _y∈[0,1]^D (f_c(y)+\|x-y\|^p ). From Kantorovich duality, it is well known that, for any p≥1p≥ 1, p(μ,ν)=supf∈b([0,1]D)(∫fμ−∫fcν), _p^p(μ,ν)= _f _b ([0,1]^D ) ( fdμ- f_cdν ), where bC_b denotes the set of all bounded continuous functions on [0,1]D[0,1]^D. By noting that (f+a)c=fc+a(f+a)_c=f_c+a for any constant a, we can restrict this function class by shifting f such that infx∈[0,1]Df(x)=0 _x∈[0,1]^Df(x)=0. Since f satisfies f(x)−f(y)≤‖x−y‖p≤Dp/2f(x)-f(y)≤\|x-y\|^p≤ D^p/2 under the optimal transform, this shift guarantees 0≤f(x)≤Dp/20≤ f(x)≤ D^p/2. Consequently, for any y∈[0,1]Dy∈[0,1]^D, we have fc(y)=supx∈[0,1]D(f(x)−‖x−y‖p)≥f(y)−‖y−y‖p=f(y)≥0. f_c(y)= _x∈[0,1]^D (f(x)-\|x-y\|^p )≥ f(y)-\|y-y\|^p=f(y)≥ 0. Hence, we can take the supremum on the set of all bounded continuous functions from [0,1]D→[0,Dp/2][0,1]^D→[0,D^p/2]. Thus, p(μ,ν)=supf∈b([0,1]D)∩0≤f≤Dp/2(∫fdμ−∫fcdν). _p^p(μ,ν)= _f _b ([0,1]^D )∩\0≤ f≤ D^p/2\ ( fdμ- f_cdν ). Let W(x1,…,xn)= W(x_1,…,x_n)= p(Pdata,Pnemp) _p^p(P_data,P^emp_n) = = supf∈b([0,1]D)∩0≤f≤Dp/2(1n∑i=1nf(xi)−∫fcdPdata). _f _b ([0,1]^D )∩\0≤ f≤ D^p/2\ ( 1n _i=1^nf(x_i)- f_cdP_data ). Hence, W(x1,…,xi,…,xn)−W(x1,…,xi′,…,xn) W(x_1,…,x_i,…,x_n)-W(x_1,…,x_i ,…,x_n) ≤ ≤ supf∈b([0,1]D)∩0≤f≤Dp/2(1n∑j=1nf(xj)−∫fcdPdata) _f _b ([0,1]^D )∩\0≤ f≤ D^p/2\ ( 1n _j=1^nf(x_j)- f_cdP_data ) −supf∈b([0,1]D)∩0≤f≤Dp/2(1n∑j≠if(xj)+f(xi′)n−∫fcdPdata) - _f _b ([0,1]^D )∩\0≤ f≤ D^p/2\ ( 1n _j≠ if(x_j)+ f(x_i )n- f_cdP_data ) ≤ ≤ 1nsup0≤f≤Dp/2(f(xi)−f(xi′)) 1n _0≤ f≤ D^p/2 (f(x_i)-f(x_i ) ) ≤ ≤ Dp/2n. D^p/2n. Applying the bounded difference inequality, we get that with probability at least 1−ξ1-ξ, p(Pnemp,Pdata)≤ _p^p(P^emp_n,P_data)≤ p(Pnemp,Pdata)+Dplog(1/ξ)2n \ EW_p^p(P^emp_n,P_data)+ D^p (1/ξ)2n ≤ ≤ n−p/d+Dplog(1/ξ)2n, \ n^-p/d+ D^p (1/ξ)2n, where the last inequality follows from Chakraborty et al. (2026, Theorem 10). The desired inequality now follows from taking the p-th root of both sides and applying the subadditivity inequality (a+b)1/p≤a1/p+b1/p(a+b)^1/p≤ a^1/p+b^1/p (valid for a,b≥0a,b≥ 0 and p≥1p≥ 1). ∎ Lemma 21. For any M>0M>0 and any estimator v v of v⋆v_ , p(Pdata,M(P^1−δ(v^)))≲ _p(P_data, T_M( P_1-δ( v))) p(Pdata,Pnemp)+δκ3+M−12p+M(KL(P1−δ,P^1−δ(v^)))12p. _p(P_data,P^emp_n)+δ _3+M^- 12p+M (KL (P_1-δ, P_1-δ( v) ) ) 12p. Proof. We first observe the following: p(Pdata,M(P^1−δ(v^))) _p(P_data, T_M( P_1-δ( v))) ≤ ≤ p(Pdata,M(P1−δ))+p(M(P1−δ),M(P^1−δ(v^))) _p(P_data, T_M(P_1-δ))+W_p( T_M(P_1-δ), T_M( P_1-δ( v))) ≤ ≤ p(Pdata,M(P1−δ))+2M(TV(M(P1−δ),M(P^1−δ(v^))))1/p _p(P_data, T_M(P_1-δ))+2M (TV ( T_M(P_1-δ), T_M( P_1-δ( v)) ) )^1/p ≤ ≤ p(Pdata,M(P1−δ))+2M(TV(P1−δ,P^1−δ(v^)))1/p _p(P_data, T_M(P_1-δ))+2M (TV (P_1-δ, P_1-δ( v) ) )^1/p ≲ p(Pdata,M(P1−δ))+M(KL(P1−δ,P^1−δ(v^)))12p _p(P_data, T_M(P_1-δ))+M (KL (P_1-δ, P_1-δ( v) ) ) 12p ≤ ≤ p(Pdata,Pnemp)+p(Pnemp,P1−δ)+p(P1−δ,M(P1−δ))+M(KL(P1−δ,P^1−δ(v^)))12p. _p(P_data,P^emp_n)+W_p(P^emp_n,P_1-δ)+W_p(P_1-δ, T_M(P_1-δ))+M (KL (P_1-δ, P_1-δ( v) ) ) 12p. We note that, p(Pnemp,P1−δ)≤ _p^p(P^emp_n,P_1-δ)≤ X1∼Pnemp,X0∼(0,ID)‖X1−(α1−δX1+β1−δX0)‖p _X_1 P^emp_n,X_0 (0,I_D)\|X_1-( _1-δX_1+ _1-δX_0)\|^p ≲ (1−α1−δ)pX1∼Pnemp‖X1‖p+β1−δpX1∼Pnemp,X0∼(0,ID)‖X0‖p (1- _1-δ)^pE_X_1 P^emp_n\|X_1\|^p+ _1-δ^pE_X_1 P^emp_n,X_0 (0,I_D)\|X_0\|^p ≲ δpκ3. \ δ^p _3. Taking the p-th root, we get, p(Pnemp,P1−δ)≲δκ3W_p(P^emp_n,P_1-δ) δ _3. Further, p(P1−δ,M(P1−δ))≤ _p^p(P_1-δ, T_M(P_1-δ))≤ ∥X1−δ∥p∥X1−δ∥2≥M \|X_1-δ\|^p 1\\|X_1-δ\|_2≥ M\ ≤ ≤ ‖X1−δ‖2pℙ(‖X1−δ‖2≥M) E\|X_1-δ\|^2p P (\|X_1-δ\|_2≥ M ) ≲ ‖X1−δ‖2M E\|X_1-δ\|_2M ≲ M−1/2, M^-1/2, so that, taking p-th roots, p(P1−δ,M(P1−δ))≲M−1/(2p)W_p(P_1-δ, T_M(P_1-δ)) M^-1/(2p). Hence, p(Pdata,M(P^1−δ(v^)))≲ _p(P_data, T_M( P_1-δ( v))) p(Pdata,Pnemp)+δκ3+M−12p+M(KL(P1−δ,P^1−δ(v^)))12p. _p(P_data,P^emp_n)+δ _3+M^- 12p+M (KL (P_1-δ, P_1-δ( v) ) ) 12p. ∎ Appendix E Supporting Results from the Literature Lemma 22. For any ε∈(0,1) ∈(0,1), there exists a ReLU network architecture ℛ(W,L,B,K),RN\! (W,L,B,K ), such that for any target 11-Lipschitz function f defined on [0,1]D[0,1]^D with f(0)=0f(0)=0, the architecture yields an approximation f f satisfying ‖f−f^‖∞≤ε.\|f- f\|_∞≤ . Moreover, the Lipschitz continuity of f f is bounded by ‖f^(x)−f^(y)‖≤10D‖x−y‖1,∀x,y∈[0,1]D.\| f(x)- f(y)\|≤ 10D\,\|x-y\|_1, ∀\,x,y∈[0,1]^D. The configuration of the network architecture is K=D,B=O(1),L=O(log(1/ε)+D),K= D,B=O(1),L=O\! ( (1/ )+D ), and W=O(ε−D(log(1/ε)+D)).W=O\! ( ^-D ( (1/ )+D ) ). Lemma 23 (Lemma 21 of Nakada and Imaizumi (2020)). Let ℱ=ℛ(W,L,B)F=RN(W,L,B) be a space of ReLU networks with the number of weights, the number of layers, and the maximum absolute value of weights bounded by W, L, and B, respectively. Then, log(ϵ,ℱ,ℓ∞)≤Wlog(2LBL(W+1)L1ϵ). (ε;F, _∞)≤ W (2LB^L(W+1)^L 1ε ).