Paper deep dive
From Distributions to Stochastic Processes: Neural Approximation of Measure-Valued Maps
Yichen Wang, Ziyi Wang, Wenlian Lu, Chenghuang Shen, Jianfeng Liu, Zhengdong Xiao, Longjiu Luo, Qianrong Wang
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Learning mappings between probability distributions arises naturally when inputs and outputs are represented by populations of samples rather than individual observations. We develop an approximation-theoretic framework for distribution-to-distribution learning and extend it to mappings between stochastic processes. For continuous operators on $W_2$-compact families of finite-dimensional probability laws, we establish uniform neural approximation in the 2-Wasserstein metric using finite law statistics, a simplex-valued neural map, and a shared atomic output support that guarantees valid probability measures. We further extend this principle to probability laws on separable Hilbert spaces through finite-rank orthogonal projections. These results establish the representational feasibility of learning transformations between probability laws rather than deterministic vectors or functions. To demonstrate practical relevance, we study two problems naturally defined at the distribution level: prediction of first-passage-time distributions for an Ornstein--Uhlenbeck process and nonlinear response-path laws of a Duffing oscillator. Because the theory is model-agnostic and broader than any single practical architecture, the experiments use task-adapted neural models rather than reproducing the theoretical construction exactly. In both problems, the proposed models outperform a fixed-feature MLP baseline and distribution-space kernel regression. These experiments complement the theory by demonstrating the practical learnability of distribution-to-distribution transformations in random systems.
Tags
Links
- Source: https://arxiv.org/abs/2609.33649v1
- Canonical: https://arxiv.org/abs/2609.33649v1
Trouble viewing inline? Open PDF directly →
Full Text
108,651 characters extracted from source content.
Expand or collapse full text
From Distributions to Stochastic Processes: Neural Approximation of Measure-Valued Maps Yichen Wang Affiliation: Shanghai Center for Mathematical Sciences, Shanghai, China Email: yichenwang23@m.fudan.edu.cn Ziyi Wang Affiliation: Department of Mathematics, Fudan University, Shanghai, China Email: 21110180058@m.fudan.edu.cn Wenlian Lu Affiliation: Center for Applied Mathematics, Fudan University, Shanghai, China Email: wenlian@fudan.edu.cn Chenghuang Shen Affiliation: Alibaba Group, Hangzhou, China Affiliation: Shanghai Center for Mathematical Sciences, Shanghai, China Email: austinshen9618@hotmail.com Jianfeng Liu Affiliation: Alibaba Group Email: ljf7451320@hotmail.com Zhengdong Xiao Affiliation: Alibaba Group Email: tianqiao.xzd@alibaba-inc.com Longjiu Luo Affiliation: Alibaba Group Email: luolongjiu@gmail.com Qianrong Wang Affiliation: Alibaba Group Email: qianrong.wqr@alibaba-inc.com Abstract Learning mappings between probability distributions arises naturally in settings where both inputs and outputs are represented by populations of samples rather than individual observations. This paper develops an approximation-theoretic framework for such distribution-to-distribution learning problems and extends it to mappings between stochastic processes. For continuous operators acting on W2W_2 -compact families of finite-dimensional input laws, we establish uniform neural approximation in the 2-Wasserstein metric. Our construction relies on finite law statistics as input representations, a simplex-valued neural map, and a shared atomic output support that guarantees the output is always a valid probability measure. We then extend this approximation principle to probability laws on separable Hilbert spaces via fixed finite-rank orthogonal projections. Together, these results establish the representational feasibility of learning transformations whose inputs and outputs are probability laws rather than deterministic vectors or functions. To assess practical relevance, we construct numerical problems whose target transformations are naturally defined at the level of distributions or probability laws: predicting the first-passage-time distribution of an Ornstein–Uhlenbeck process and the nonlinear response-path laws of a Duffing oscillator. Since the theoretical framework is deliberately model-agnostic and covers a substantially broader class of operators than can be efficiently represented by any single practical architecture, our experiments instantiate the framework using task-adapted neural models rather than reproducing the theoretical construction verbatim. Both models outperform a fixed-feature MLP baseline and distribution-space kernel regression. These proof-of-concept demonstrations complement the approximation theory by supporting the learnability and practical relevance of distribution-to-distribution transformations for random systems. 1 Introduction Motivation and problem setup. Many random systems are naturally described at the population level. In uncertainty propagation, mean-field dynamics, ensemble forecasting, stochastic PDEs, and Bayesian inverse problems, the desired prediction is an output law rather than the response of one prescribed input realization (Lord et al., 2014; Walsh, 1986; Sznitman, 1991; Carmona and Delarue, 2018; Stuart, 2010; Evensen, 2003). When the observations are paths or fields, the same question concerns laws on function spaces. In both cases, probability distributions are the mathematical objects of interest, while finite ensembles are observations of those objects. We formulate this task as learning a distributional operator :⊂2()⟶2(),μ⟼(μ),G:K _2(X) _2(Y), μ (μ), where 2P_2 denotes probability measures with finite second moment. Training observations are pairs of empirical laws (μ^i,ν^i)i=1n\( μ_i, ν_i)\_i=1^n. The input and output ensembles within a pair may be sampled independently, so there need not be a meaningful map between individual particles. The target is instead the law-level relation νi=(μi) _i=G( _i). We measure approximation in the 2-Wasserstein metric W2W_2. The finite-dimensional setting takes =ℝdxX=R^d_x and =ℝdyY=R^d_y; the process setting takes X and Y to be separable Hilbert spaces of paths or fields. Related work and the theoretical gap. Classical universal approximation theory treats continuous maps between finite-dimensional vectors (Cybenko, 1989). Neural-operator theory extends the encoding–mapping–reconstruction principle to deterministic functions, with DeepONet and Fourier neural operators providing prominent examples (Lu et al., 2021; Li et al., 2021; Kovachki et al., 2023). Our input and output, however, are probability laws. A deterministic operator between functions does not establish that an input law can be encoded by finitely many observable statistics, transformed continuously, and decoded into a valid output law with uniform error control. Distribution regression and invariant set networks address complementary parts of the problem. Distribution regression maps an empirical law to a scalar or vector, while distribution-to-distribution regression permits a distributional response (Szabó et al., 2016; Oliva et al., 2013). DeepSets and Set Transformers provide practical permutation-invariant encoders for unordered samples (Zaheer et al., 2017; Lee et al., 2019), but an architecture for finite sets alone does not give a uniform approximation theorem on a space of laws. More specifically, Oliva et al. (2013) assume densities on compact Euclidean domains and derive an expected L2L_2 density-risk bound for a kernel-smoothed orthogonal-series estimator. Their setting does not yield uniform W2W_2 approximation or guarantee that a truncated output series is a probability density. Pham and Warin (2023) approximate maps of the form (x,μ)↦V(x,μ)(x,μ) V(x,μ), whose outputs are vectors or functions rather than probability laws. Measure-valued approximation is closer to our objective. The probabilistic transformer maps finite-dimensional features to mixtures of fixed measures and proves compact-uniform approximation for 1(ℝD)P_1(R^D)-valued maps in W1W_1 (Kratsios, 2021); the general framework of Kratsios et al. (2023) covers maps between broad classes of metric spaces, including probability-valued targets. Transformer dynamics can also interpolate finitely many prescribed input–target measure pairs (Geshkovski et al., 2024). These results establish important general measure-valued approximation principles. What remains unresolved for our setting is their combination with a law-valued input accessed through an empirical ensemble: a finite encoder that captures W2W_2 rather than only weak behavior, a decoder that always returns a probability measure, a uniform bound for replacing feature integrals by sample averages, and a lifting from finite-dimensional laws to laws on path spaces. The gap is particularly delicate because W2W_2 convergence requires both weak convergence and control of second moments, and path-space reduction adds input- and output-projection errors. Our theory. For a W2W_2-compact family ⊂2(ℝdx)K _2(R^d_x) and a W2W_2-continuous operator :→2(ℝdy)G:K _2(R^d_y), Theorem 2.1 constructs a finite law encoder from bounded test-function expectations and a truncated second-moment feature. A neural map with a softmax output assigns weights to a common finite set of output atoms. The bounded tests control weak behavior, the moment feature supplies the additional information required by W2W_2, and the simplex-valued decoder guarantees nonnegative weights with unit mass. Compactness and continuity then give uniform W2W_2 approximation of the target operator. For empirical inputs, Proposition B.5 replaces population feature integrals by independent sample averages and separates the population approximation error from a uniform high-probability query-sampling error. Theorem 2.2 extends the result to laws on separable Hilbert spaces. Fixed finite-rank orthogonal maps encode paths into coordinate vectors, the finite-dimensional distributional operator acts on the induced coordinate laws, and synthesis maps reconstruct output paths. Uniform projection on W2W_2-compact law families decomposes the total error into input projection, finite-dimensional neural approximation, and output reconstruction. The argument remains entirely at the level of probability laws and therefore does not require eventwise pairing of input and output paths. Together, the two theorems provide the missing W2W_2-specific link between empirical law encoding, valid probability-valued decoding, and process-law approximation. They establish representational feasibility rather than statistical training rates or a mandatory practical architecture. Contributions. Our contributions are: • We formulate distribution-to-distribution learning from finite, unpaired input and output ensembles as approximation of a continuous operator between probability-law spaces. • We prove uniform neural approximation in W2W_2 on compact families of finite-dimensional input laws using finitely many law statistics, neural simplex weights, and a valid common-support output law. We also give a uniform query-sample bound for empirical feature averages. • We lift the approximation result to process and random-field laws on separable Hilbert spaces through fixed finite-rank projections, with explicit control of input projection, finite-dimensional approximation, and output reconstruction. • We demonstrate the framework on a censored Ornstein–Uhlenbeck first-passage law and nonlinear Duffing response-path laws. Task-specific categorical and generative decoders assess practical learnability against fixed-feature and distribution-space kernel baselines without identifying the existence construction with a unique architecture. 2 Universal Approximation Theory for Operators between Probability Laws We establish two universal approximation results: one for operators between finite-dimensional probability laws and one for laws on separable Hilbert spaces. Both use encoding, neural mapping, and valid probability-law reconstruction. Supporting constructions, empirical-feature bounds, and projection lemmas are summarized below; their formal statements and all proofs appear in Appendix B. 2.1 Problem setup and assumptions We begin with the standard probability-space notation used throughout. If (,d)(X,d_X) is a Polish space, ()P(X) denotes its Borel probability measures and, for 1≤p<∞1≤ p<∞, p()=μ∈():∫d(x,x0)pμ(x)<∞,P_p(X)= \μ (X): _Xd_X(x,x_0)^p\,dμ(x)<∞ \, (1) where the definition does not depend on the fixed point x0∈x_0 . For μ,ν∈p()μ,ν _p(X), let Π(μ,ν) (μ,ν) be their set of couplings and define Wp(μ,ν)=infπ∈Π(μ,ν)(∫×d(x,y)pπ(x,y))1/p.W_p^X(μ,ν)= _π∈ (μ,ν) ( _X×Xd_X(x,y)^p\,dπ(x,y) )^1/p. (2) When the underlying space is clear, we simply write WpW_p. For a measurable map F:→F:X , F#μF_\#μ denotes the pushforward law, and for i.i.d. samples X1,…,XNin∼μX_1,…,X_N_in μ we write μ^Nin=1Nin∑b=1NinδXb μ_N_in= 1N_in _b=1^N_in _X_b (3) for the empirical measure. Let X and Y be Polish spaces. A distributional operator is a map :⊂2()⟶2().G:K _2(X) _2(Y). Its uniform approximation error on K is measured by supμ∈W2((μ),(μ)) _μ W_2^Y(N(μ),G(μ)). The finite-dimensional theory takes =ℝdxX=R^d_x and =ℝdyY=R^d_y and uses the following assumptions throughout: (A1) ⊂2(ℝdx)K _2(R^d_x) is W2W_2-compact. Consequently, the family is tight and its second moments are uniformly integrable: limR→∞supμ∈∫∥x∥2>R∥x∥22dμ(x)=0. _R→∞ _μ _\\|x\|_2>R\\|x\|_2^2\,dμ(x)=0. (4) (A2) :→2(ℝdy)G:K _2(R^d_y) is W2W_2-continuous. Hence G is uniformly continuous on K and ()G(K) is W2W_2-compact. (A3) The chosen continuous, nonpolynomial activation gives uniform neural-network approximation on compact Euclidean sets, as in the classical universal approximation theorem. Bounded clipping, or an equivalent bounded output activation, is available for the feature networks. Two kinds of input statistics are needed. Choose bounded weak-convergence tests ψ1,…,ψM∈Cb(ℝdx) _1,…, _M∈ C_b(R^d_x) and retain the second moment as a separate coordinate: Ewk(μ) E_ wk(μ) =(∫ψ1μ,…,∫ψMμ), = ( _1\,dμ,…, _M\,dμ ), m2(μ) m_2(μ) =∫‖x‖22μ(x), = \|x\|_2^2\,dμ(x), E(μ) E(μ) =(Ewk(μ),m2(μ))∈ℝM+1. =(E_ wk(μ),m_2(μ)) ^M+1. (5) The bounded coordinates determine weak behavior, while convergence of the second moments upgrades weak convergence to W2W_2 convergence. Compactness then permits a finite choice of tests at any prescribed accuracy. The construction has three layers. First, E represents the input law to the accuracy required by the uniform continuity of G. Second, a continuous simplex-valued map assigns weights to a common finite output support. Third, ordinary neural networks approximate the feature functions and simplex map. The Tietze–Dugundji extension theorem supplies the simplex-valued extension needed in the second layer, while tightness and (4) control neural approximation of the feature integrals. 2.2 Universal approximation of finite-dimensional laws The first theorem combines finite law representation with neural realization. The encoder integrates bounded neural features, including a truncated approximation of the second moment, and a softmax network predicts weights on a common finite output support. Theorem 2.1 (Universal approximation of distributional operators). Assume (A1)–(A3). Let ΔN=w∈[0,1]N:∑j=1Nwj=1 _N=\w∈[0,1]^N: _j=1^Nw_j=1\. For every ε>0 >0, there exist M,N∈ℕM,N , common output atoms a1,…,aN∈ℝdya_1,…,a_N ^d_y, bounded scalar feature networks ψ^i ψ_i, a bounded neural approximation q^R q_R of q(x)=‖x‖22q(x)=\|x\|_2^2 obtained by truncation at some radius R, and a softmax network p^:ℝM+1→ΔN p:R^M+1→ _N. Define E^wk,R(μ) E_ wk,R(μ) =(∫ψ^1μ,…,∫ψ^Mμ), = ( ψ_1\,dμ,…, ψ_M\,dμ ), E^R(μ) E_R(μ) =(E^wk,R(μ),∫q^Rμ), = ( E_ wk,R(μ), q_R\,dμ ), (6) and ^θ(μ)=∑j=1Np^j(E^R(μ))δaj. N_θ(μ)= _j=1^N p_j( E_R(μ)) _a_j. (7) Then the networks can be chosen so that supμ∈W2(^θ(μ),(μ))<ε. _μ W_2( N_θ(μ),G(μ))< . Construction and proof outline. Compactness permits finitely many bounded tests and a second-moment coordinate to distinguish input laws to any required W2W_2 tolerance (Lemma B.1). Uniform continuity of G transfers this tolerance to the outputs. Local finite-support output approximations are collected on common atoms, and a partition of unity gives continuous simplex weights. The formal continuous construction is Theorem B.2 in Appendix B.1. Bounded tests are then approximated by clipped networks on large balls; tightness controls the complement, and uniform integrability controls truncation of the quadratic feature (Lemma B.3). Regularizing the weights, approximating their logits, and applying softmax preserves nonnegativity and unit mass. Common-support stability converts weight error to W2W_2 error (Lemma B.4). The full proof of Theorem 2.1 is in Appendix B.2. The fixed atoms are an existence device shared by all input laws. Practical models may use richer structured decoders, but their approximation properties require separate justification. Access through empirical ensembles. For an already realized neural operator, the population feature integrals can be replaced by averages of independent query samples. Hoeffding concentration and the output network’s modulus of continuity give a high-probability bound uniform over the queried law in K. This separates population approximation error from empirical feature-estimation error. With a Lipschitz weight network, the latter contributes an output W2W_2 bound of order Nin−1/4N_in^-1/4 for fixed confidence and realized architecture. The formal bound and its constants are in Proposition B.5, Appendix B.3. This is a query-sample guarantee, not a training, optimization, or network-width rate. Boundedness includes the truncated quadratic feature; the same concentration conclusion need not hold for an untruncated second moment. 2.3 Universal approximation of Hilbert-space laws For process and random-field laws, finite-dimensionalization acts on each path and therefore produces a coordinate law, rather than a single feature vector representing the whole input distribution. Let X and Y be separable Hilbert spaces. Fix increasing finite-dimensional subspaces Xn⊂X_n and Ym⊂Y_m with dense unions, with dn=dimXnd_n= X_n and qm=dimYmq_m= Y_m. Orthonormal coordinates define analysis maps ℰn:→ℝdnE_X^n:X ^d_n and ℰm:→ℝqmE_Y^m:Y ^q_m and synthesis maps n:ℝdn→D_X^n:R^d_n and m:ℝqm→D_Y^m:R^q_m . Their compositions Pn=n∘ℰn,Pm=m∘ℰmP_X^n=D_X^n _X^n, P_Y^m=D_Y^m _Y^m are orthogonal projections. Their pushforwards approximate every W2W_2-compact law family uniformly; the formal reconstruction and stability results are Lemmas B.7 and B.6 in Appendix B.4. Theorem 2.2 (Universal approximation on separable Hilbert-space laws). Assume the activation condition (A3). Let X and Y be separable Hilbert spaces equipped with the orthogonal coordinate systems above. Let ⊂2()K _2(X) be W2W_2^X-compact, and let :⟶2()G:U _2(Y) be W2W_2-continuous on an open neighborhood U of K. For every ε>0 >0, there exist n,m∈ℕn,m and a finite-dimensional neural distributional operator ^θ:2(ℝdn)→2(ℝqm) N_θ:P_2(R^d_n) _2(R^q_m) of the form in Theorem 2.1 such that θ,n,m=(m)#∘^θ∘(ℰn)#N_θ,n,m=(D_Y^m)_\# N_θ (E_X^n)_\# (8) satisfies supμ∈W2(θ,n,m(μ),(μ))<ε. _μ W_2^Y(N_θ,n,m(μ),G(μ))< . Proof outline and interpretation. Uniform projection accuracy first makes (Pn)#μ(P_X^n)_\#μ close to μ over K. The open neighborhood ensures that projected input laws remain in the domain of G, and continuity controls the induced output perturbation. Uniform output projection reduces the target to a finite-dimensional coordinate law. Theorem 2.1 approximates the resulting law operator, and synthesis is an isometry. The triangle inequality thus separates input projection, finite-dimensional neural approximation, and output reconstruction errors. The complete proof is in Appendix B.4; further path-law construction details are in Appendix B.5. The result requires no eventwise coupling of input and output paths. Its orthogonal coordinate systems are fixed independently of the input law; a distribution-specific PCA or Karhunen–Loève basis, or a learned encoder, requires additional stability and uniform-reconstruction analysis. 3 Experiments Theoretical results are typically established through general constructions to ensure broad applicability, but such constructions are often not directly implementable in practice. We therefore design concrete and computationally feasible distributional neural operators for controlled finite-dimensional Gaussian distributions, first-passage times of Ornstein–Uhlenbeck (OU) processes, and nonlinear Duffing oscillators. Rather than mechanically reproducing the theoretical constructions, we follow their underlying principles and tailor each model to the structure of its problem. The controlled Gaussian experiment acts on finite-dimensional distributions, whereas the OU and Duffing tasks use empirical laws of stochastic-process paths. Owing to space limitations, the main text focuses on the OU and Duffing tasks; complete results for the controlled Gaussian task are provided in Appendix E, and the corresponding theoretical extensions appear in Appendix C. Figure 1: A general example of the model pipeline. The model used for a specific task may differ in certain details, but broadly follows this design rationale. For task-specific model structures, see Figures 5 and 6 in the appendix. To the best of our knowledge, these tasks have no established benchmark or off-the-shelf baseline. We therefore compare against task-appropriate fixed-feature and kernel-regression alternatives. Implementation details for data generation and the comparison methods are provided in Appendix D. Unless otherwise specified, input and output ensembles are sampled independently within each law pair, so the models cannot exploit trajectory-level correspondence. Every experiment uses 1,000 training law pairs and 200 test law pairs, which exhaust the 1,200-law dataset. The test set is used only to evaluate final model performance. Neural experiments are repeated with five random seeds. The default neural training configuration is AdamW for 1,000 epochs with initial learning rate 10−310^-3 and weight decay 10−210^-2. All computations are performed on one NVIDIA GeForce RTX 4090 GPU with 24 GB of memory. Source code is available at https://anonymous.4open.science/r/distributional_operator-3E6C. 3.1 Experiment I: OU First-Passage Distribution 3.1.1 Problem The OU first-passage problem has a standard interpretation as a leaky integrate-and-fire (LIF) neuron. In this model, the membrane potential evolves until it first reaches a firing threshold, at which point a spike is recorded and the potential is reset (Burkitt, 2006). Under the classical diffusion approximation with Gaussian white-noise input, the subthreshold membrane potential is an OU process. Consequently, the elapsed time from a reset to the next threshold crossing is one inter-spike interval (ISI), whose law is the OU first-passage-time law (Ricciardi and Sato, 1988; Alili et al., 2005). Specifically, the membrane potential follows dVt=(−Vtγ+m)dt+qdWt,V(0)=v0,dV_t= (- V_tγ+m )dt+ q\,dW_t, V(0)=v_0, (9) with threshold VthV_th and first-passage time T=inft≥0:V(0)=v0,Vt≥Vth.T= \t≥ 0:V(0)=v_0,V_t≥ V_th\. (10) The model does not observe (m,q)(m,q) directly. Instead, the input law is that of the drifted Brownian process dIt=mdt+qdWt,I(0)=0,dI_t=m\,dt+ q\,dW_t, I(0)=0, (11) with the same (m,q)(m,q). The target distributional operator is :ℒ(I)⟼ℒ(T),G:L(I) (T), where the first-passage law is represented by 48 finite-time categories and one censoring category. The dataset contains 1,200 input–output law pairs. Each input law is represented by NinN_in paths on 256 observation times, and each target law by NoutN_out independently simulated reset trials. The ensembles share the same law parameters but have no samplewise pairing. The parameter family, simulation scheme, censoring rule, binning, and split are detailed in Appendix D.1.1. 3.1.2 Model Instantiation Following the process pipeline, we use a training-fitted PCA encoder ℰIE_I to preprocess the paths, selecting the truncation dimension dXd_X so that the cumulative explained-variance ratio exceeds 0.990.99. For the iith input ensemble, the Distributional Operator computes αij _ij =ℰI(Iij)∈ℝdX, =E_I(I_ij) ^d_X, fij f_ij =ϕpath(αij), = _path( _ij), hi h_i =1Nin∑j=1Ninϕset(fij), = 1N_in _j=1^N_in _set(f_ij), ℓi _i =ρθ(hi),p^i=softmax(ℓi)∈Δ49. = _θ(h_i), p_i=softmax( _i)∈ _49. The shared path and set encoders followed by mean aggregation make the model invariant to the ordering of the input paths. The first 48 entries of p^i p_i parameterize the finite first-passage intervals and the last entry represents censoring. Training minimizes the count-weighted categorical negative log-likelihood. PCA fitting, layer widths, parameter counts, optimization, checkpointing, and complete baseline configurations are reported in Appendix D.1.2. 3.1.3 Evaluation Metrics and Results Negative log-likelihood is the categorical training objective and measures predictive fit. Hellinger distance and KL divergence compare the full 49-category distributions from complementary symmetric and directed perspectives. Finite-bin W2W_2 measures displacement of probability mass over passage time after removing and renormalizing the censoring category, while tail error measures the absolute error in the probability of no crossing. Lower is better for every metric; formal definitions are in Appendix D.1.3. Model NLL Hellinger KL divergence W2,finW_2,fin Tail error Distributional Operator 3.11176±0.000853.11176± 0.00085 0.18055±0.000480.18055± 0.00048 0.10636±0.000850.10636± 0.00085 0.16634±0.002860.16634± 0.00286 0.00373±0.000180.00373± 0.00018 Fixed-Feature MLP 3.11680±0.000113.11680± 0.00011 0.18592±0.000130.18592± 0.00013 0.11140±0.000110.11140± 0.00011 0.20957±0.000260.20957± 0.00026 0.00399±0.000030.00399± 0.00003 Kernel regression 3.152413.15241 0.216160.21616 0.147010.14701 0.345670.34567 0.006390.00639 Table 1: Test results for input process laws mapped to OU first-passage-time laws. Entries are the mean ± standard deviation over five runs with random initializations. Kernel regression is deterministic and therefore has zero standard deviation. Lower is better for every metric. For a qualitative theoretical reference, we evaluate the OU first-passage-time density using its Bessel-bridge representation (Alili et al., 2005). This reference density is independent of the empirical target histograms used for training and lets us compare the shapes predicted by the learned and kernel-based models against a continuous first-passage-time calculation. Figure 2: OU first-passage-time probability densities for three representative test laws. The continuous blue curve is the theoretical reference calculated from the Bessel-bridge representation. The step curves show the densities predicted by the Distributional Operator, kernel regression, and the Fixed-Feature MLP on the 48 finite time bins. Agreement with the reference curve compares how well the models recover the location, height, and tail of the first-passage-time density. The Distributional Operator gives the lowest mean NLL, Hellinger distance, KL divergence, finite-bin W2W_2, and tail error. Kernel regression is worse on every metric, indicating that local averaging between training laws does not capture the map as accurately. The density comparisons in Figure 2 show that the learned operator recovers the location, peak, and tail of the theoretical first-passage-time density across representative test laws. 3.2 Experiment I: Duffing Uncertainty Propagation 3.2.1 Problem We study a distribution-to-distribution map induced by the fixed hardening Duffing oscillator Y¨(t)+0.2Y˙(t)+Y(t)+Y(t)3=X(t),Y(0)=Y˙(0)=0, Y(t)+0.2 Y(t)+Y(t)+Y(t)^3=X(t), Y(0)= Y(0)=0, on t∈[0,20]t∈[0,20]. Given a probability law μ of stochastic forcing paths, the target is the response law under the Duffing solution operator S, namely :μ⟼#μ.G:μ _\#μ. This is a distributional uncertainty-propagation problem: the model must characterize an ensemble of nonlinear responses rather than predict the response to one prescribed forcing trajectory. Such laws arise in nonlinear vibration and energy-harvesting models under stochastic or broadband excitation (Daqaq, 2010; Jia, 2020). The dataset contains 1,200 input–output law pairs. Each law is represented by NinN_in forcing paths and NoutN_out independently sampled response paths on 256 observation times. The ensembles share the same forcing-law parameters but have no samplewise pairing. The forcing family, numerical solver, and measure-level construction are detailed in Appendix D.2. 3.2.2 Model Instantiation Following the process pipeline, we use training-fitted PCA encoders ℰXE_X and ℰYE_Y to preprocess the forcing and response paths. Their truncation dimensions dXd_X and dYd_Y are selected so that each cumulative explained-variance ratio exceeds 0.990.99. For the iith input ensemble, αij _ij =ℰX(Xij)∈ℝdX, =E_X(X_ij) ^d_X, ci c_i =ρθ(1Nin∑j=1Ninϕθ(αij)), = _θ ( 1N_in _j=1^N_in _θ( _ij) ), zir z_ir ∼(0,IdY),βir=Fθ(zir,ci), (0,I_d_Y), _ir=F_θ(z_ir;c_i), Y^ir Y_ir =Y(βir),r=1,…,Nout, =D_Y( _ir), r=1,…,N_out, where FθF_θ is a conditional RealNVP and YD_Y is the inverse output-PCA reconstruction map (Dinh et al., 2017). Thus the empirical forcing-process law is projected, encoded by expectation features, combined with independent latent draws, and decoded into samples from the predicted response-process law. The conditional generator uses six affine coupling transformations. For a masked partition h=(hℐ,h)h=(h_I,h_J), each transformation has the form hℐ′ h _I =hℐ, =h_I, h′ h _J =h⊙exp(s~ℓ(hℐ,ci))+tℓ(hℐ,ci),s~ℓ=2tanh(sℓ/2). =h_J ( s_ (h_I,c_i))+t_ (h_I,c_i), s_ =2 (s_ /2). Training minimizes the debiased Sinkhorn divergence between NoutN_out generated and NoutN_out target response paths. Layer widths, parameter counts, PCA fitting, coupling masks, optimization schedules, and complete baseline configurations are reported in Appendices D.2.2 and D.2.3. 3.2.3 Evaluation Metrics and Results Sinkhorn divergence is the training objective and measures regularized transport discrepancy between response ensembles. RBF MMD compares kernel mean embeddings and is sensitive to broad distributional differences. Sliced W2W_2 averages one-dimensional transport discrepancies over random path-space projections, while energy distance compares within- and between-ensemble Euclidean distances. Formal definitions and evaluation settings are in Appendix D.2.4. We compared the performance of distributional operator (proposed model) with that of two baseline models on the test data in Table 2.Lower values indicate closer agreement for all four metrics. To mitigate the effects of randomness, the evaluation results were computed using five different random seeds and are reported as the mean ± standard deviation of the returns, except for the kernel regression model, which is deterministic and therefore reported as a single value. To intuitively illustrate the predictive performance, we present 3the results of the models on three randomly selected test datasets. Model Sinkhorn MMD SW2SW_2 Energy Distributional Operator 8.2050±0.27868.2050± 0.2786 0.0760±0.00100.0760± 0.0010 0.1639±0.00360.1639± 0.0036 0.3407±0.01040.3407± 0.0104 Fixed-Feature MLP 9.9426±0.67959.9426± 0.6795 0.0777±0.00170.0777± 0.0017 0.1718±0.00360.1718± 0.0036 0.3890±0.01730.3890± 0.0173 Kernel regression 18.629118.6291 0.12550.1255 0.24980.2498 0.86500.8650 Table 2: Test performance on the Duffing dataset. Neural-model entries are the mean ± sample standard deviation on the 200-law test set. Kernel regression is a single fit with a fixed bandwidth. Sinkhorn denotes the Sinkhorn divergence, MMD the maximum mean discrepancy, SW2SW_2 the sliced 2-Wasserstein distance, and Energy the energy distance; lower values indicate better agreement. Figure 3: Conditional distribution prediction for the Duffing oscillator. Target and predicted response distributions are shown for three randomly chosen test forcing measures. Solid blue lines and shaded regions denote the reference mean and 5th–95th percentile envelope, while dashed orange lines and hatched regions show the corresponding model predictions. Agreement in both statistics demonstrates the model’s accuracy in approximating the nonlinear pushforward from forcing to response distributions. 4 Conclusion We developed a constructive approximation theory for continuous distribution-to-distribution maps in the 2-Wasserstein metric. On W2W_2-compact input families, finitely many law statistics, a neural simplex-valued map, and a shared atomic support suffice for uniform approximation by valid output measures. We also quantified the error incurred when those statistics are estimated from a finite unordered ensemble. By combining finite-rank analysis and synthesis maps with the finite-dimensional construction, the same principle extends to probability laws on separable Hilbert spaces and therefore to distributional operators acting on stochastic processes. The OU first-passage and Duffing uncertainty-propagation experiments show that this law-level formulation is computationally meaningful and learnable from unpaired ensembles. The task-adapted categorical and conditional-flow models outperform the fixed-feature and kernel-regression comparators on the reported metrics. These architectures instantiate the framework’s central principles but do not reproduce the proof’s finite-support construction, and the reported comparisons are not evidence of architectural optimality: the baseline set is necessarily limited, stronger designs may exist, and the remaining predictive errors leave substantial room for improvement. Future work should tighten the connection between constructive theory and practical architectures, weaken the compactness and neighborhood assumptions, obtain sharper statistical and optimization rates, and evaluate the framework on broader and more demanding stochastic systems. References Alili et al. (2005) L. Alili, P. Patie, and J. L. Pedersen Representations of the first hitting time density of an ornstein–uhlenbeck process. Stochastic Models 21 (4), p. 967–980. External Links: Document, Link Cited by: §3.1.1, §3.1.3. Burkitt (2006) A. N. Burkitt A review of the integrate-and-fire neuron model: i. homogeneous synaptic input. Biological Cybernetics 95 (1), p. 1–19. External Links: Document, Link Cited by: §3.1.1. Carmona and Delarue (2018) R. Carmona and F. Delarue Probabilistic theory of mean field games with applications. Springer. Cited by: §1. Cybenko (1989) G. Cybenko Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems 2 (4), p. 303–314. Cited by: §1. Daqaq (2010) M. F. Daqaq Response of uni-modal duffing-type harvesters to random forced excitations. Journal of Sound and Vibration 329 (18), p. 3621–3631. External Links: Document Cited by: §3.2.1. Dinh et al. (2017) L. Dinh, J. Sohl-Dickstein, and S. Bengio Density estimation using real NVP. In International Conference on Learning Representations, Cited by: §3.2.2. Evensen (2003) G. Evensen The ensemble kalman filter: theoretical formulation and practical implementation. Ocean Dynamics 53 (4), p. 343–367. Cited by: §1. Geshkovski et al. (2024) B. Geshkovski, P. Rigollet, and D. Ruiz-Balet Measure-to-measure interpolation using transformers. arXiv preprint arXiv:2411.04551. Cited by: §1. Jia (2020) Y. Jia Review of nonlinear vibration energy harvesting: duffing, bistability, parametric, stochastic and others. Journal of Intelligent Material Systems and Structures 31 (7), p. 921–944. External Links: Document Cited by: §3.2.1. Kovachki et al. (2023) N. B. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar Neural operator: learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research 24 (89), p. 1–97. Cited by: §1. Kratsios et al. (2023) A. Kratsios, C. Liu, M. Lassas, M. V. de Hoop, and I. Dokmanić An approximation theory for metric space-valued functions with a view towards deep learning. arXiv preprint arXiv:2304.12231. Cited by: §1. Kratsios (2021) A. Kratsios Universal regular conditional distributions. arXiv preprint arXiv:2105.07743. Cited by: §1. Lee et al. (2019) J. Lee, Y. Lee, J. Kim, A. R. Kosiorek, S. Choi, and Y. W. Teh Set transformer: a framework for attention-based permutation-invariant neural networks. In International Conference on Machine Learning, p. 3744–3753. Cited by: §1. Li et al. (2021) Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, Cited by: §1. Lord et al. (2014) G. J. Lord, C. E. Powell, and T. Shardlow An introduction to computational stochastic PDEs. Cambridge University Press. Cited by: §1. Lu et al. (2021) L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis DeepONet: learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. Nature Machine Intelligence 3 (3), p. 218–229. Cited by: §1. Oliva et al. (2013) J. Oliva, B. Póczos, and J. Schneider Distribution to distribution regression. In International Conference on Machine Learning, p. 1049–1057. Cited by: §1. Pham and Warin (2023) H. Pham and X. Warin Mean-field neural networks: learning mappings on wasserstein space. Neural Networks 168, p. 380–393. Cited by: §1. Ricciardi and Sato (1988) L. M. Ricciardi and S. Sato First-passage-time density and moments of the ornstein–uhlenbeck process. Journal of Applied Probability 25 (1), p. 43–57. External Links: Document, Link Cited by: §3.1.1. Stuart (2010) A. M. Stuart Inverse problems: a bayesian perspective. Acta Numerica 19, p. 451–559. Cited by: §1. Szabó et al. (2016) Z. Szabó, B. K. Sriperumbudur, B. Póczos, and A. Gretton Learning theory for distribution regression. Journal of Machine Learning Research 17 (152), p. 1–40. Cited by: §1. Sznitman (1991) A. Sznitman Topics in propagation of chaos. In Ecole d’Eté de Probabilités de Saint-Flour XIX–1989, p. 165–251. Cited by: §1. Walsh (1986) J. B. Walsh An introduction to stochastic partial differential equations. In Ecole d’Eté de Probabilités de Saint-Flour XIV–1984, p. 265–439. Cited by: §1. Zaheer et al. (2017) M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Póczos, R. Salakhutdinov, and A. J. Smola Deep sets. In Advances in Neural Information Processing Systems, Vol. 30. Cited by: §1. AI Use Disclosure Generative AI tools were used to polish the writing and improve readability; to search for and retrieve relevant literature and technical information; to prepare initial drafts of portions of the manuscript; and to assist in generating the synthetic datasets used in the experiments. The authors independently checked the retrieved sources and reviewed, revised, and validated all AI-assisted text, citations, data-generation procedures, and resulting datasets. The authors take responsibility for the final content of this work, including all text, claims, and artifacts produced with the aid of generative AI. Appendix A Table of Notations Notation Meaning ℝR, ℕN Real numbers and positive integers. X, Y Input and output state spaces; in the process setting these are separable Hilbert spaces. 2()P_2(X) Borel probability measures on X with finite second moment. W2W_2, W2W_2^X Quadratic Wasserstein distance, with the superscript indicating the underlying state space when needed. K A W2W_2-compact family of admissible input laws. μ, ν Generic input and output probability laws. G Target distributional operator, mapping an input law to an output law. N, ^θ N_θ Continuous finite-support approximation and its neural realization. δx _x Dirac probability measure concentrated at x. F#F_\#, F#μF_\#μ Pushforward of a measure through a measurable map F. E(μ)E(μ), E^(μ) E(μ) Exact finite-dimensional law encoder and its neural approximation. ψi _i, ψ^i ψ_i Bounded test functions and their neural approximations. m2(μ)m_2(μ) Second moment ∫∥x∥22μ(x) x _2^2\,dμ(x). aja_j, pjp_j Shared output-support points and their simplex-valued weights. ΔN _N Probability simplex q∈[0,1]N:∑j=1Nqj=1\q∈[0,1]^N: _j=1^Nq_j=1\. M, N Number of input features and output support points in the theoretical construction. NinN_in Number of samples representing an input distribution in a theoretical query or experiment. NoutN_out Number of samples representing an output distribution in an experiment. NbatchN_batch Number of law-level examples in an optimization mini-batch. ℰnE_X^n, nD_X^n Finite-rank encoder and decoder for the Hilbert-space state X. PnP_X^n Orthogonal projection n∘ℰnD_X^n _X^n. E, Pr Expectation and probability. suppsupp, LipLip, diamdiam Support, Lipschitz constant, and diameter. ε , δ Approximation tolerance and confidence level. Appendix B Definitions and Theory Details We collect the definitions, supporting results, and proofs for the main-text universal approximation theorems. B.1 Continuous finite-support construction Lemma B.1 (Finite separation in W2W_2). Let ⊂2(ℝd)K _2(R^d) be W2W_2-compact and let η>0η>0. Then there are M∈ℕM , functions ψ1,…,ψM∈Cb(ℝd) _1,…, _M∈ C_b(R^d), and r>0r>0 such that, for all μ,ν∈μ,ν , ‖Ewk(μ)−Ewk(ν)‖2<r, \|E_ wk(μ)-E_ wk(ν)\|_2<r, |m2(μ)−m2(ν)|<r⟹W2(μ,ν)<η. |m_2(μ)-m_2(ν)|<r W_2(μ,ν)<η. Proof. Choose a countable convergence-determining family in Cb(ℝd)C_b(R^d). If no finite subfamily worked, compactness would give sequences μk,νk _k, _k with all first k test integrals and their second moments approaching one another, but with W2(μk,νk)≥ηW_2( _k, _k)≥η. Passing to W2W_2-convergent subsequences gives limits with identical bounded-test integrals and hence the same weak limit. Convergence of the second moments then characterizes W2W_2 convergence, a contradiction. ∎ Theorem B.2 (Continuous finite-support approximation). Let ⊂2(ℝdx)K _2(R^d_x) be W2W_2-compact and let :→2(ℝdy)G:K _2(R^d_y) be W2W_2-continuous. For every ε>0 >0, there exist M,N∈ℕM,N , functions ψ1,…,ψM∈Cb(ℝdx) _1,…, _M∈ C_b(R^d_x), fixed points a1,…,aN∈ℝdya_1,…,a_N ^d_y, and a continuous map p=(p1,…,pN):ℝM+1→ΔNp=(p_1,…,p_N):R^M+1→ _N, where ΔN=q∈[0,1]N:∑j=1Nqj=1 _N=\q∈[0,1]^N: _j=1^Nq_j=1\, such that 0(μ)=∑j=1Npj(E(μ))δajN_0(μ)= _j=1^Np_j(E(μ)) _a_j (12) satisfies supμ∈W2(0(μ),(μ))<ε, _μ W_2(N_0(μ),G(μ))< , with E defined in (5). Proof of Theorem B.2. By uniform continuity of G and Lemma B.1, choose a finite feature map E and radius r so that nearby features produce outputs within ε/4 /4. Cover the compact set E()E(K) by finitely many balls centered at E(μℓ)E( _ ). A continuous partition of unity αℓ _ subordinate to this cover may be obtained from normalized tent functions. Finitely supported laws are dense in 2(ℝdy)P_2(R^d_y), so choose ν^ℓ ν_ with W2(ν^ℓ,(μℓ))<ε/4W_2( ν_ ,G( _ ))< /4. Collect the finitely many atoms of all ν^ℓ ν_ into a common list a1,…,aNa_1,…,a_N and write their weights as βℓj _ j. Define on E()E(K) qj(z)=∑ℓαℓ(z)βℓj,0(μ)=∑jqj(E(μ))δaj.q_j(z)= _ _ (z) _ j, _0(μ)= _jq_j(E(μ)) _a_j. Mixing couplings for the active ℓ gives W22((μ),0(μ))≤∑ℓαℓ(E(μ))W22((μ),ν^ℓ)<ε2.W_2^2(G(μ),N_0(μ))≤ _ _ (E(μ))W_2^2(G(μ), ν_ )< ^2. Finally extend the continuous map q:E()→ΔNq:E(K)→ _N to all of ℝM+1R^M+1 using the Dugundji extension theorem (the simplex is convex). ∎ B.2 Neural realization and proof of Theorem 2.1 Lemma B.3 (Uniform neural approximation of the features). Let ⊂2(ℝd) K _2(R^d) be tight and satisfy (4). For finitely many ψi∈Cb(ℝd) _i∈ C_b(R^d) and every ρ>0ρ>0, bounded neural networks ψ^i ψ_i and a bounded neural network q^R q_R exist such that supμ∈|∫(ψ^i−ψi)μ| _μ∈ K | ( ψ_i- _i)dμ | <ρ, <ρ, supμ∈|∫q^Rμ−m2(μ)| _μ∈ K | q_Rdμ-m_2(μ) | <ρ. <ρ. Proof. For a bounded test, choose a ball carrying uniformly high probability, approximate the test uniformly there, and clip the network to a fixed bound; tightness controls the complement. For q(x)=‖x‖22q(x)=\|x\|_2^2, choose R so that the tail integral in (4) is small, multiply q by a continuous cutoff supported in the ball of radius 2R2R, and uniformly approximate and clip that compactly supported function. Its integral differs from m2m_2 only by the approximation error and the uniformly small tail. ∎ Lemma B.4 (Common-support stability). For a1,…,aN∈ℝda_1,…,a_N ^d, w,v∈ΔNw,v∈ _N, and D=maxj,k‖aj−ak‖2D= _j,k\|a_j-a_k\|_2, W22(∑jwjδaj,∑jvjδaj)≤D22‖w−v‖1.W_2^2\! ( _jw_j _a_j, _jv_j _a_j )≤ D^22\|w-v\|_1. Proof. Match the mass common to wjw_j and vjv_j at zero cost and transport the remaining mass, whose total is ‖w−v‖1/2\|w-v\|_1/2, across distances at most D. ∎ The same neural realization argument applies to a tight input family with uniformly integrable second moments whenever the continuous finite-support construction is uniform on that family. Proof of Theorem 2.1. Apply Theorem B.2 with error ε/2 /2. Lemma B.3 uniformly approximates its feature vector. The resulting feature vectors lie in a compact neighborhood of E()E(K). Uniform continuity of p makes the induced weight error small. To realize p by softmax, use the positive regularization pjτ=pj+τ1+Nτ.p_j^τ= p_j+τ1+Nτ. Approximate the continuous logits logpjτ p_j^τ uniformly on that compact neighborhood, and apply softmax. Lemma B.4 converts the final weight error to W2W_2 error. ∎ B.3 Uniform empirical-feature bound Proposition B.5 (Uniform empirical-feature error bound). Fix any neural realization from Theorem 2.1, and write fi=ψ^if_i= ψ_i for i≤Mi≤ M and fM+1=q^Rf_M+1= q_R. Define the realized feature ranges ℓi _i =infx∈ℝdxfi(x), = _x ^d_xf_i(x), ui u_i =supx∈ℝdxfi(x), = _x ^d_xf_i(x), si s_i =ui−ℓi, =u_i- _i, the feature box =∏i=1M+1[ℓi,ui]Z= _i=1^M+1[ _i,u_i], and the modulus of continuity of the output network on this box by ωp^,(t)=supz,z′∈‖z−z′‖∞≤t‖p^(z)−p^(z′)‖1. _ p,Z(t)= _ subarraycz,z \\ \|z-z \|_∞≤ t subarray\| p(z)- p(z )\|_1. Let ηθ=supμ∈W2(^θ(μ),(μ)),D=maxj,k‖aj−ak‖2. _θ= _μ W_2( N_θ(μ),G(μ)), D= _j,k\|a_j-a_k\|_2. If X1,…,XNin∼i.i.d.μX_1,…,X_N_in .i.d. μ and the feature integrals in (6) are replaced by empirical averages, then, for every t>0t>0, supμ∈Pr(W2(^θ,Nin(μ),(μ))>ηθ+Dωp^,(t)2)≤2∑1≤i≤M+1si>0exp(−2Nint2si2). _μ \! (W_2( N_θ,N_in(μ),G(μ))> _θ+D _ p,Z(t)2 )≤ 2 _ subarrayc1≤ i≤ M+1\\ s_i>0 subarray \! (- 2N_int^2s_i^2 ). (13) In particular, with smax=maxisis_ = _is_i and tNin,δ=smaxlog(2(M+1)/δ)2Nin,t_N_in,δ=s_ (2(M+1)/δ)2N_in, the following bound holds with probability at least 1−δ1-δ, uniformly in μ∈μ : supμ∈Pr(W2(^θ,Nin(μ),(μ))>ηθ+Dωp^,(tNin,δ)2)≤δ. _μ \! (W_2( N_θ,N_in(μ),G(μ))> _θ+D _ p,Z(t_N_in,δ)2 )≤δ. (14) If p p is LpL_p-Lipschitz from ℓ∞ _∞ to ℓ1 _1 on Z, the second term in (14) is at most DLpsmax2(log(2(M+1)/δ)2Nin)1/4.D L_ps_ 2 ( (2(M+1)/δ)2N_in )^1/4. Consequently, for every ε>0 >0 and δ∈(0,1)δ∈(0,1) the networks and NinN_in can be chosen so that supμ∈Pr(W2(^θ,Nin(μ),(μ))>ε)<δ. _μ \! (W_2( N_θ,N_in(μ),G(μ))> )<δ. Proof of Proposition B.5. For μ∈μ , set z(μ)=(∫f1μ,…,∫fM+1μ),zNin(μ)=1Nin∑b=1Nin(f1(Xb),…,fM+1(Xb)).z(μ)= ( f_1\,dμ,…, f_M+1\,dμ ), z_N_in(μ)= 1N_in _b=1^N_in (f_1(X_b),…,f_M+1(X_b) ). Both vectors belong to Z. For every coordinate with si>0s_i>0, Hoeffding’s inequality gives Pr(|zNin,i(μ)−zi(μ)|>t)≤2exp(−2Nint2si2). \! (|z_N_in,i(μ)-z_i(μ)|>t )≤ 2 \! (- 2N_int^2s_i^2 ). Coordinates with si=0s_i=0 have zero empirical error. Hence a union bound, whose right-hand side is independent of μ, yields supμ∈Pr(‖zNin(μ)−z(μ)‖∞>t)≤2∑1≤i≤M+1si>0exp(−2Nint2si2). _μ \! (\|z_N_in(μ)-z(μ)\|_∞>t )≤ 2 _ subarrayc1≤ i≤ M+1\\ s_i>0 subarray \! (- 2N_int^2s_i^2 ). (15) On the complementary event, the definition of the modulus gives ‖p^(zNin(μ))−p^(z(μ))‖1≤ωp^,(t).\| p(z_N_in(μ))- p(z(μ))\|_1≤ _ p,Z(t). Applying Lemma B.4 to the empirical- and population-feature outputs therefore gives W2(^θ,Nin(μ),^θ(μ))≤Dωp^,(t)2.W_2( N_θ,N_in(μ), N_θ(μ))≤ D _ p,Z(t)2. The triangle inequality and the definition of ηθ _θ, combined with (15), prove (13). For t=tNin,δt=t_N_in,δ, each nonzero-range summand on the right of (13) is at most δ/(M+1)δ/(M+1) after including its leading factor 22, so their sum is at most δ. This proves (14). If p p is LpL_p-Lipschitz, then ωp^,(t)≤Lpt _ p,Z(t)≤ L_pt, which gives the displayed Lipschitz specialization. Finally, apply (14) with confidence level δ/2δ/2. Theorem 2.1 permits choosing ηθ<ε/2 _θ< /2. Since Z is compact and p p is continuous, ωp^,(t)→0 _ p,Z(t)→ 0 as t↓0t 0. Taking NinN_in sufficiently large makes the second term in (14) smaller than ε/2 /2, proving the last assertion. ∎ B.4 Hilbert-space reconstruction and proof of Theorem 2.2 Lemma B.6 (Lipschitz pushforward and reconstruction). If F:→F:X is Lipschitz, then W2(F#μ,F#ν)≤Lip(F)W2(μ,ν).W_2^Y(F_\#μ,F_\#ν) (F)W_2^X(μ,ν). Moreover, for a measurable P:→P:X , W2(P#μ,μ)≤(∫‖Px−x‖2μ(x))1/2.W_2^X(P_\#μ,μ)≤ ( \|Px-x\|_X^2dμ(x) )^1/2. Proof. Push a coupling of μ,νμ,ν through (F,F)(F,F) for the first inequality. For the second, push μ through x↦(Px,x)x (Px,x). ∎ Lemma B.7 (Uniform Wasserstein autoencoding in a separable Hilbert space). Let X be a separable Hilbert space and let X1⊂X2⊂⋯⊂X_1⊂ X_2⊂·s be finite-dimensional subspaces whose union is dense. Write dn=dimXnd_n= X_n, choose an orthonormal basis e1n,…,ednne_1^n,…,e_d_n^n of XnX_n, and define ℰnx _X^nx =(⟨x,e1n⟩,…,⟨x,ednn⟩), = ( x,e_1^n _X,…, x,e_d_n^n _X ), nz _X^nz =∑j=1dnzjejn,Pn=n∘ℰn. = _j=1^d_nz_je_j^n, P_X^n=D_X^n _X^n. Then ℰnE_X^n is a contraction, nD_X^n is an isometry, and PnP_X^n is the orthogonal projection onto XnX_n. Moreover, every W2W_2^X-compact set ⊂2()K _2(X) satisfies limn→∞supμ∈W2((Pn)#μ,μ)=0. _n→∞ _μ W_2^X ((P_X^n)_\#μ,μ )=0. (16) Proof. The norm claims and the identity of PnP_X^n follow from orthonormality. For μ∈2()μ _2(X), set rn(μ)=(∫‖(I−Pn)x‖2μ(x))1/2.r_n(μ)= ( _X\|(I-P_X^n)x\|_X^2\,dμ(x) )^1/2. The coupling x↦(Pnx,x)x (P_X^nx,x) and Lemma B.6 give W2((Pn)#μ,μ)≤rn(μ)W_2^X((P_X^n)_\#μ,μ)≤ r_n(μ). Because Pnx→xP_X^nx→ x and ‖(I−Pn)x‖≤‖x‖\|(I-P_X^n)x\|_X≤\|x\|_X, dominated convergence gives rn(μ)→0r_n(μ)→ 0 for each μ. The convergence is uniform on K. Indeed, for any coupling π of μ and ν, the reverse triangle inequality in L2(π,)L^2(π;X) and the contractivity of I−PnI-P_X^n yield |rn(μ)−rn(ν)|≤(∫‖(I−Pn)(x−y)‖2π(x,y))1/2≤(∫‖x−y‖2π(x,y))1/2.|r_n(μ)-r_n(ν)|≤ ( \|(I-P_X^n)(x-y)\|_X^2\,dπ(x,y) )^1/2≤ ( \|x-y\|_X^2\,dπ(x,y) )^1/2. Taking the infimum over π shows that every rnr_n is 11-Lipschitz on 2()P_2(X). Pointwise convergence together with this common Lipschitz bound and compactness of K implies supμ∈rn(μ)→0 _μ r_n(μ)→ 0, proving (16). ∎ Proof of Theorem 2.2. Write Pn=n∘ℰnP_X^n=D_X^n _X^n and define PmP_Y^m similarly. Lemma B.7 gives supμ∈W2((Pn)#μ,μ)⟶0. _μ W_2^X ((P_X^n)_\#μ,μ ) 0. (17) Since K is compact and U is an open neighborhood of K, (17) implies that (Pn)#μ∈(P_X^n)_\#μ for every μ∈μ and all sufficiently large n. Continuity of G then implies supμ∈W2(((Pn)#μ),(μ))⟶0. _μ W_2^Y (G((P_X^n)_\#μ),G(μ) ) 0. (18) To see the uniformity in (18), otherwise choose nk→∞n_k→∞ and μk∈ _k violating it. After passing to a subsequence, μk→μ∈ _k→μ in W2W_2^X; then (17) also gives (Pnk)#μk→μ(P_X^n_k)_\# _k→μ, contradicting continuity of G at μ. Fix n so that the left-hand side of (18) is less than ε/3 /3. The family n=((Pn)#μ):μ∈⊂2() L_n=\G((P_X^n)_\#μ):μ \ _2(Y) is compact. Applying Lemma B.7 in Y, choose m such that supν∈nW2((Pm)#ν,ν)<ε/3. _ν∈ L_nW_2^Y((P_Y^m)_\#ν,ν)< /3. (19) The encoded input family n=(ℰn)#K_n=(E_X^n)_\#K is compact. On nK_n define n,m(λ)=(ℰm)#((n)#λ).S_n,m(λ)=(E_Y^m)_\#G((D_X^n)_\#λ). Both pushforwards are Lipschitz by Lemma B.6, so n,mS_n,m is a continuous operator between finite-dimensional law spaces. Theorem 2.1 therefore provides ^θ N_θ such that supλ∈nW2(^θ(λ),n,m(λ))<ε/3. _λ _nW_2 ( N_θ(λ),S_n,m(λ) )< /3. (20) For μ∈μ , put λ=(ℰn)#μλ=(E_X^n)_\#μ and ν=((Pn)#μ)ν=G((P_X^n)_\#μ). Since mD_Y^m is an isometry, Lemma B.6, (20), (19), and (18) give W2(θ,n,m(μ),(μ)) W_2^Y(N_θ,n,m(μ),G(μ)) ≤W2((m)#^θ(λ),(Pm)#ν) ≤ W_2^Y ((D_Y^m)_\# N_θ(λ),(P_Y^m)_\#ν ) +W2((Pm)#ν,ν)+W2(ν,(μ))<ε. +W_2^Y((P_Y^m)_\#ν,ν)+W_2^Y(ν,G(μ))< . Taking the supremum over μ∈μ proves the claim. ∎ B.5 Additional path-law construction details B.5.1 Distributional Approximation in Wasserstein Space From a distributional perspective, a random-field operator does not act directly on a representation of random variables defined on a fixed probability space; instead, it acts on probability measures over path spaces. Following the notation introduced earlier, let 2()P_2(X) denote the space of probability measures on X with finite second moments, and let W2W_2^X denote the Wasserstein distance induced by the norm on X. For a measurable map F:→F:X , write F#μF_\#μ for the pushforward measure. In the distributional formulation, a random-field operator is expressed as a map between probability measures on path spaces. Let X and Y be the path spaces of the input and output random fields, respectively. We consider a continuous operator of the form :⊂2()⟶2(),G:U _2(X) _2(Y), (21) where continuity is understood with respect to W2W_2^X and W2W_2^Y on the domain and codomain, respectively. Given a W2W_2^X-compact set ⊂K , our goal is to construct a random-field distributional operator approximator θ,n,mN_θ,n,m such that, for every ε>0 >0, supμ∈W2(θ,n,m(μ),(μ))<ε. _μ W_2^Y\! (N_θ,n,m(μ),G(μ) )< . (22) The domain is denoted by U because the subsequent proof must compare not only (μ)G(μ) with the output of the approximating operator, but also (μ)G(μ) with the output associated with the input distribution reconstructed by the input autoencoder. The target operator must therefore be defined and continuous on a set containing both the original family of input distributions and its family of reconstructed distributions. To apply the finite-dimensional distributional operator approximation result to random-field distributions on path spaces, we use an encoding–finite-dimensional distributional operator–decoding architecture. Let ℰn _X^n :→ℝn, :X ^n, n _X^n :ℝn→, :R^n , (23) ℰm _Y^m :→ℝm, :Y ^m, m _Y^m :ℝm→ :R^m (24) be the encoders and decoders for the input and output path spaces, respectively. Define the corresponding path-space reconstruction operators by Pn=n∘ℰn,Pm=m∘ℰm.P_X^n=D_X^n _X^n, P_Y^m=D_Y^m _Y^m. (25) If X is an X-valued random field with μ=ℒ(X)μ=L(X), then the encoded random variable is ℰn(X)∈ℝnE_X^n(X) ^n, whose law is (ℰn)#μ=ℒ(ℰn(X))∈2(ℝn).(E_X^n)_\#μ=L\! (E_X^n(X) ) _2(R^n). (26) Thus, finite-dimensional encoding does not turn the random-field distribution μ into a single deterministic vector. Rather, it pushes the path-space random variable forward to a finite-dimensional random vector while preserving its probability distribution. For fixed n and m, the finite-dimensional distributional operator approximation theory developed above provides an abstract finite-dimensional neural operator module ^θ:2(ℝn)⟶2(ℝm). N_θ:P_2(R^n) _2(R^m). (27) The random-field distributional operator approximator induced by this module is defined by θ,n,m=(m)#∘^θ∘(ℰn)#.N_θ,n,m=(D_Y^m)_\# N_θ (E_X^n)_\#. (28) This architecture can be represented by the following commutative diagram: 2()→2()(ℰn)#↓↑(m)#2(ℝn)→^θ2(ℝm). array[]cP_2(X)& \ G\ &P_2(Y)\\[3.99994pt] (E_X^n)_\# && (D_Y^m)_\#\\[3.99994pt] P_2(R^n)& \ N_θ\ &P_2(R^m). array (29) In diagram (29), the upper arrow is the target random-field distributional operator G; the left arrow pushes probability measures on the input path space forward to a finite-dimensional distribution space; the lower arrow is the neural operator module between finite-dimensional distribution spaces; and the right arrow decodes the finite-dimensional output distribution back to the output path space. The central question can therefore be stated as follows: under what conditions does the approximating path formed by the left, lower, and right arrows approximate the target operator G uniformly in W2W_2^Y? The proof strategy can be understood by examining the sources of error in the commutative diagram. First, on the input side, we must control the information loss caused by finite-dimensional encoding. In other words, after the original input distribution μ has been encoded and decoded, the resulting reconstructed distribution should be close to the original distribution in W2W_2^X. If this error can be controlled uniformly over the family of input distributions under consideration, and if the target operator G is continuous with respect to Wasserstein distance, then the finite-dimensional reduction on the input side cannot produce an uncontrollable error in the output distribution. Second, once the input has been encoded, the problem becomes one of operator approximation between finite-dimensional distribution spaces. This part is not developed again here; it is supplied by the finite-dimensional distributional operator approximation theory established above and is denoted uniformly by ^θ N_θ. Its role is to approximate the finite-dimensional distributional map induced by the original operator G between the input and output encoding spaces, namely, the lower arrow in the commutative diagram. Finally, on the output side, we must control the error introduced when the finite-dimensional output distribution is decoded back to path space. Even if the finite-dimensional distributional module gives a good approximation in 2(ℝm)P_2(R^m), the output decoder must still recover the target output distribution stably in W2W_2^Y. The output path space must therefore possess a corresponding uniform autoencoding property in the Wasserstein sense. In summary, the approximation error of a random-field distributional operator can be divided into three components: the input-side autoencoding error, the finite-dimensional distributional operator approximation error, and the output-side autoencoding error. The main theorem that follows is organized precisely around this decomposition. It first controls the reconstruction errors on the input and output sides through finite-dimensional autoencoding structures on the path spaces, then invokes the finite-dimensional distributional operator approximation result established above to control the error along the lower arrow of the commutative diagram, and finally obtains uniform Wasserstein approximation of G by θ,n,mN_θ,n,m. It follows that the key to the subsequent theory is not to reconstruct the finite-dimensional distributional operator, but to prove that the three errors above can be controlled simultaneously in the Wasserstein sense. This formulation does not require the input and output random fields to be paired pointwise on the same event space. It is therefore better suited to learning problems in which random-field distributions are accessed only through sample sets or empirical distributions. The following theory will make this error decomposition based on diagram (29) rigorous. Appendix C Supplementary Theory Proof For Experiments C.1 Structured Output Representations and Neural Realizations General finite-support output measures are expressive, but their support size can grow rapidly with the covering number of the input set and the complexity of local output laws. Practical problems often provide additional output structure. We consider common compact support, finite-dimensional parametric families, generative pushforwards, and conditional densities. Lemma C.1 (Finite-dimensional factorization). Let ⊂2(ℝdx)K _2(R^d_x) be W2W_2-compact, let Y be a Banach space, and let F:→F:K be continuous. For every ε>0 >0, there are a finite-dimensional encoder E:→ℝME:K ^M and a continuous map H:ℝM→H:R^M such that supμ∈‖H(E(μ))−F(μ)‖<ε. _μ \|H(E(μ))-F(μ)\|_Y< . If F()F(K) lies in a convex set 0Y_0, then H(E())H(E(K)) may be required to lie in 0Y_0. Proof. Use Lemma B.1 and a finite cover of E()E(K) with a partition of unity αℓ _ . For centers μℓ _ , define H(z)=∑ℓ=1Lαℓ(z)F(μℓ).H(z)= _ =1^L _ (z)F( _ ). Uniform continuity controls the error. This is a convex combination, so it preserves 0Y_0. Extend H from the compact feature set by the Tietze–Dugundji theorem. ∎ C.1.1 Common Support on a Compact Set Proposition C.2 (Shared support points). Suppose the assumptions of Theorem B.2 hold and there is a compact D⊂ℝdyD ^d_y such that supp(μ)⊂DsuppG(μ)⊂ D for every μ∈μ . For every ε>0 >0, there are common points a1,…,aN∈Da_1,…,a_N∈ D, an encoder E, and continuous simplex-valued weights pjp_j such that (μ)=∑j=1Npj(E(μ))δaj,supμ∈W2((μ),(μ))<ε.N(μ)= _j=1^Np_j(E(μ)) _a_j, _μ W_2(N(μ),G(μ))< . The support size N may be bounded by the covering number of D at a scale proportional to ε . Proof. Choose an r-net ajj=1N\a_j\_j=1^N of D and a continuous partition of unity λj\ _j\ subordinate to it. For any law ν on D, set Q(ν)=∑j=1N(∫Dλj(y)ν(y))δaj.Q(ν)= _j=1^N ( _D _j(y)dν(y) ) _a_j. The coupling λj(y)dν(y) _j(y)dν(y) gives W2(ν,Q(ν))≤rW_2(ν,Q(ν))≤ r. The map q(μ)=(∫λ1(μ),…,∫λN(μ))q(μ)=( _1dG(μ),…, _NdG(μ)) is continuous. Apply Lemma C.1 to factor q through E and choose the weight error so that the common-support coupling bound is below ε−r -r. ∎ C.1.2 Finite-Dimensional Parametric Output Families Let =νa:a∈Θ⊂2(ℝdy) F=\ _a:a∈ \ _2(R^d_y) be a parametric family. Suppose (μ)=νA(μ)G(μ)= _A(μ) for a continuous parameter map A:→ΘA:K→ . Proposition C.3 (Parametric output family). Assume a↦νa _a is W2W_2-continuous on a compact convex set KΘ⊂ΘK_ ⊂ containing A()A(K). For every ε>0 >0, there are E and continuous H:E()→KΘH:E(K)→ K_ such that supμ∈W2(νH(E(μ)),(μ))<ε. _μ W_2( _H(E(μ)),G(μ))< . Proof. Uniform continuity of a↦νa _a converts a sufficiently small parameter error into a W2W_2 error below ε . Apply Lemma C.1 to A while preserving the convex parameter set. ∎ Continuity of A is an additional assumption: continuity of G does not by itself imply a continuous parameter selection. Nonidentifiability, label switching in mixture models, or a discontinuous inverse parameterization can prevent such a selection. Legal parameters must also be enforced. For Gaussian outputs, a Cholesky factor or a positive transform of eigenvalues guarantees a positive semidefinite covariance. Corollary C.4 (Neural realization of a parametric output family). Under Proposition C.3, suppose there is r>0r>0 such that Cr=a∈ℝp:dist(a,A())≤r⊂KΘ.C_r=\a ^p:dist(a,A(K))≤ r\⊂ K_ . For every ε>0 >0, there are input-feature networks and a finite-dimensional network H^θ:ℝM→ℝp H_θ:R^M ^p such that ^θ(μ)=νH^θ(E^(μ)),supμ∈W2(^θ(μ),(μ))<ε. N_θ(μ)= _ H_θ( E(μ)), _μ W_2( N_θ(μ),G(μ))< . C.1.3 Generative Pushforward Representation Let ρ∈2(ℝdz)ρ _2(R^d_z) be a fixed latent law. Suppose there is a measurable map G:×ℝdz→ℝdyG:K×R^d_z ^d_y with (μ)=G(μ,⋅)#ρG(μ)=G(μ,·)_\#ρ. Proposition C.5 (Generative pushforward). If μ↦G(μ,⋅)μ G(μ,·) is continuous from (,W2)(K,W_2) to L2(ρ,ℝdy)L^2(ρ;R^d_y), then for every ε>0 >0 there are an encoder E and continuous H:E()→L2(ρ,ℝdy)H:E(K)→ L^2(ρ;R^d_y) such that H(μ)=H(E(μ))#ρ,supμ∈W2(H(μ),(μ))<ε.N_H(μ)=H(E(μ))_\#ρ, _μ W_2(N_H(μ),G(μ))< . Proof. Apply Lemma C.1 in L2(ρ,ℝdy)L^2(ρ;R^d_y) with error ε . Coupling the two outputs by the same Z∼ρZ ρ gives W22(H(μ),(μ))≤‖H(E(μ))−G(μ,⋅)‖L2(ρ)2.W_2^2(N_H(μ),G(μ))≤\|H(E(μ))-G(μ,·)\|_L^2(ρ)^2. ∎ The existence of a generator for each individual law does not imply a continuous choice of generators across a family. The pushforward map is many to one, and mode permutations or phase choices can make a generator selection discontinuous even when the laws vary continuously. Polynomial chaos is one example. If Z∼ρZ ρ and Y=f(Z)=∑i=0∞aiψi(Z)Y=f(Z)= _i=0^∞a_i _i(Z) in L2(ρ)L^2(ρ), then the truncation fM(Z)=∑i=0Maiψi(Z)f_M(Z)= _i=0^Ma_i _i(Z) satisfies W22(f#ρ,fM#ρ)≤‖f−fM‖L2(ρ)2.W_2^2(f_\#ρ,f_M_\#ρ)≤\|f-f_M\|_L^2(ρ)^2. For standard Gaussian ρ, the basis may be the Hermite polynomials. If the coefficients depend continuously on the input law, a network can predict them from finite distributional features. Corollary C.6 (Neural-operator realization of a generator). Under Proposition C.5, for every ε>0 >0 there are input-feature networks and a deterministic neural operator H^θ:ℝM→L2(ρ,ℝdy) H_θ:R^M→ L^2(ρ;R^d_y) such that ^θ(μ)=H^θ(E^(μ))#ρ,supμ∈W2(^θ(μ),(μ))<ε. N_θ(μ)= H_θ( E(μ))_\#ρ, _μ W_2( N_θ(μ),G(μ))< . C.1.4 Conditional Density Representation Proposition C.7 (Existence of a conditional-density representation). Let D⊂ℝdyD ^d_y be compact with positive Lebesgue measure. Suppose every (μ)G(μ) is supported on D and has density pμp_μ, and the map P:→L1(D)P:K→ L^1(D), P(μ)=pμP(μ)=p_μ, is continuous. For every ε>0 >0, there are E and a continuous map H:E()→L1(D)H:E(K)→ L^1(D) such that every H(z)H(z) is a probability density and H(μ)(dy)=H(E(μ))(y)dy,supμ∈W2(H(μ),(μ))<ε.N_H(μ)(dy)=H(E(μ))(y)dy, _μ W_2(N_H(μ),G(μ))< . Proof. The density set (D)=p∈L1(D):p≥0 a.e.,∫Dp(y)dy=1A(D)= \p∈ L^1(D):p≥ 0 a.e.,\ _Dp(y)dy=1 \ is closed and convex. Apply Lemma C.1 to P in L1(D)L^1(D) while preserving (D)A(D). For laws supported on D, W22(H(μ),(μ))≤diam(D)22‖H(E(μ))−pμ‖L1(D).W_2^2(N_H(μ),G(μ))≤ diam(D)^22\|H(E(μ))-p_μ\|_L^1(D). ∎ A normalized exponential preserves positivity and unit mass: (h)(y)=exp(h(y))∫Dexp(h(ζ))ζ. S(h)(y)= (h(y)) _D (h(ζ))dζ. Alternatively, for a measurable partition Qjj=1J\Q_j\_j=1^J and simplex weights αj _j, the bin density pα(y)=∑j=1Jαj|Qj|Qj(y)p_α(y)= _j=1^J _j|Q_j|1_Q_j(y) is automatically valid and converges in L1L^1 as the partition is refined. Lemma C.8 (Stability of the normalized exponential). For compact D of positive Lebesgue measure and h1,h2∈C(D)h_1,h_2∈ C(D), ‖(h1)−(h2)‖L1(D)≤2‖h1−h2‖C(D).\| S(h_1)- S(h_2)\|_L^1(D)≤ 2\|h_1-h_2\|_C(D). Proof. Set r=h1−h2r=h_1-h_2, ht=h2+trh_t=h_2+tr, and pt=(ht)p_t= S(h_t). Then ∂tpt(y)=pt(y)(r(y)−∫Dr(ζ)pt(ζ)ζ), _tp_t(y)=p_t(y) (r(y)- _Dr(ζ)p_t(ζ)dζ ), so ‖∂tpt‖L1(D)≤2‖r‖C(D)\| _tp_t\|_L^1(D)≤ 2\|r\|_C(D). Integrating over t∈[0,1]t∈[0,1] proves the result. ∎ Corollary C.9 (Neural realization of normalized-exponential densities). Under Proposition C.7, suppose a continuous energy map A:→C(D)A:K→ C(D) satisfies pμ(y)=exp(A(μ)(y))∫Dexp(A(μ)(ζ))ζ.p_μ(y)= (A(μ)(y)) _D (A(μ)(ζ))dζ. For every ε>0 >0, there are feature networks and a neural operator H^θ:ℝM→C(D) H_θ:R^M→ C(D) such that p^θ(y∣μ) p_θ(y μ) =exp(H^θ(E^(μ))(y))∫Dexp(H^θ(E^(μ))(ζ))ζ, = ( H_θ( E(μ))(y)) _D ( H_θ( E(μ))(ζ))dζ, ^θ(μ)(dy) N_θ(μ)(dy) =p^θ(y∣μ)dy, = p_θ(y μ)dy, and supμ∈W2(^θ(μ),(μ))<ε. _μ W_2( N_θ(μ),G(μ))< . Proof. Approximate A after finite-dimensional factorization and input-feature approximation, obtaining supμ∈‖H^θ(E^(μ))−A(μ)‖C(D)<γ. _μ \| H_θ( E(μ))-A(μ)\|_C(D)<γ. Lemma C.8 gives an L1L^1 density error below 2γ2γ, and therefore W22(^θ(μ),(μ))≤diam(D)2γ.W_2^2( N_θ(μ),G(μ)) (D)^2γ. Choose γ<ε2/diam(D)2γ< ^2/diam(D)^2. ∎ If pμp_μ is strictly positive and μ↦logpμ p_μ is continuous into C(D)C(D), one may take A(μ)=logpμA(μ)= p_μ. Densities with zeros require a positive regularization or bin-density approximation first. Appendix D Data Generation and Experiments Protocol D.1 OU First-Passage Time D.1.1 Data Generation The 1,200 laws are divided equally among three drift regimes: subthreshold:m∼Unif(0.75,0.95),400 laws,balanced:m∼Unif(0.95,1.05),400 laws,suprathreshold:m∼Unif(1.05,1.25),400 laws. array[]lsubthreshold:&m (0.75,0.95),&400 laws,\\ balanced:&m (0.95,1.05),&400 laws,\\ suprathreshold:&m (1.05,1.25),&400 laws. array The master seed is zero. Parameter sampling uses seed 101, and q is sampled independently with logq∼Unif(log0.1,log0.35). q ( 0.1, 0.35). Each law is assigned a unique identifier from 0 through 1199. For the target, 200 independent reset trials follow dVt dV_t =(−Vt+m)dt+qdWt, =(-V_t+m)dt+ q\,dW_t, V0 V_0 =0, =0, using random seed 303, step size 0.0020.002, threshold 1, and terminal time 8. A crossing within a step is linearly interpolated; trials without a crossing by time 8 are censored. Uncensored times are assigned to 48 equal-width bins with edges tb=8b/48t_b=8b/48, and censored trials form category 49. Dividing the 49 category counts by 200 gives the target mass. The input process file is generated from the same parameter table. For every law, 200 independent paths satisfy dXt=mdt+qdWt,X0=0,dX_t=m\,dt+ q\,dW_t, X_0=0, and are evaluated at tk=8k/255t_k=8k/255, k=0,…,255k=0,…,255. Independent Gaussian increments with variance 8/2558/255 give the exact transition law on this grid, and cumulative summation preserves the temporal correlation within each path. A deterministic stream generated by SeedSequence([0, law_id]) makes every law reproducible independently of generation order. The resulting array has shape 1200×200×2561200× 200× 256 and is stored in single precision. The input and target files are joined by their unique law identifiers rather than by row position, with matched parameters and regime labels verified before use. A deterministic regime-stratified split with seed zero assigns 1,000 complete laws to training and 200 to testing; all paths from a law remain in the same split. There is no validation split, and the test laws are excluded from PCA fitting, optimization, scheduling, early stopping, and checkpoint selection. Figure 4: Predicted and empirical interspike-interval distributions. Ground-truth probability masses (left) and those predicted by the Process Distributional Operator (right) are shown for 50 test examples selected uniformly at random and ordered by increasing empirical mean interspike interval (ISI). Each heatmap row represents one input law, and the columns span the 48 finite ISI bins from 0 to 8. The narrow column beside each heatmap shows the probability of no threshold crossing before the end of the observation window. All panels use the same color scale. D.1.2 Model and Baseline Implementations For the Distributional Operator, one global PCA is fitted to the 200,000 training paths after flattening the law and path axes. The 256 time coordinates are centered using their training means but are not standardized coordinatewise. The mean and covariance are accumulated in double precision, and a symmetric eigendecomposition orders the components by decreasing eigenvalue. The smallest dimension whose cumulative explained-variance ratio is strictly greater than 0.990.99 is 14, with achieved ratio 0.990150110.99015011. The fitted mean and components are then applied unchanged to both splits. The PCA state, threshold, split seed, and training-law identifiers are stored with the checkpoint so evaluation never refits the transform. The resulting model input has shape Nbatch×Nin×14N_batch× N_in× 14. Each compressed path passes through a shared 14→32→3214→ 32→ 32 SiLU path MLP. A DeepSets element network then maps 32→32→3232→ 32→ 32 with a SiLU hidden activation and no dropout, and averages the resulting 32-dimensional embeddings over the input-sample axis. Finally, a 32→32→4932→ 32→ 49 outer network with a GELU hidden activation returns the logits. The complete sequence of tensor shapes is (Nbatch,Nin,256) (N_batch,N_in,256) →PCA(Nbatch,Nin,14) PCA(N_batch,N_in,14) →ϕpath(Nbatch,Nin,32) _path(N_batch,N_in,32) →DeepSetsmean(Nbatch,32)→(Nbatch,49). DeepSets\ mean(N_batch,32) ρ(N_batch,49). The neural model has 5,265 trainable parameters. Softmax returns the predicted bin masses and cumulative summation returns the predicted CDF. It is trained for 1,000 epochs with AdamW, learning rate 10−310^-3, weight decay 10−210^-2, mini-batches of 64 laws, no learning-rate scheduler, and no early stopping; the saved final checkpoint is evaluated once on the untouched test split. Figure 5: Distributional Operator used for the OU first-passage task. An empirical law of PCA-compressed input paths is encoded by shared path and DeepSets networks, mapped to categorical output-law parameters, and sampled when realizations of the predicted first-passage distribution are required. The kernel baseline converts the input particles into a normalized 32-category histogram over the global training-particle range. For query X and training law XiX_i, let d(X,Xi)d(X,X_i) be the exact one-dimensional quadratic Wasserstein distance between the piecewise-uniform histograms. It predicts p^(X)=∑iwi(X)pi,wi(X)=exp[−d(X,Xi)2/(2h2)]∑jexp[−d(X,Xj)2/(2h2)],h=0.15. p(X)= _iw_i(X)p_i, w_i(X)= [-d(X,X_i)^2/(2h^2)] _j [-d(X,X_j)^2/(2h^2)], h=0.15. Self-matches are excluded for training-law evaluation, and distances are computed in blocks of 256 training laws. The fixed-feature MLP replaces learned particle aggregation by the sample mean, population variance, and eight empirical sine and cosine moments at frequencies drawn once from (0,1)N(0,1). The resulting 18 features pass through a linear 18→3218→ 32 map, GELU, and a 32→4932→ 49 logit map. It is trained for 1,000 epochs with mini-batches of 64 laws using the count-weighted categorical NLL objective. D.1.3 Evaluation Metrics Five runs use mini-batches of 64 under the default training configuration. For normalized target histogram pip_i and prediction p^i p_i, we report count-normalized NLL, Hellinger distance, empirical-to-prediction KL divergence, tail error, and finite-support W2W_2: NLLi _i =−∑k=149piklogp^ik, =- _k=1^49p_ik p_ik, H(pi,p^i) H(p_i, p_i) =[12∑k=149(pik−p^ik)2]1/2, = [ 12 _k=1^49( p_ik- p_ik)^2 ]^1/2, DKL(pi∥p^i) D_KL(p_i\| p_i) =∑k:pik>0piklogpikp^ik, = _k:p_ik>0p_ik p_ik p_ik, Etail(pi,p^i) E_tail(p_i, p_i) =|pi,49−p^i,49|. =|p_i,49- p_i,49|. For finite-support W2W_2, category 49 is removed, the first 48 masses are renormalized, and the exact one-dimensional quadratic Wasserstein distance is computed between the associated piecewise-uniform laws. Metrics are averaged over the 200 test laws and then summarized across the five runs. D.2 Nonlinear Duffing Oscillator D.2.1 Data Generation The Duffing coefficients are fixed throughout the dataset at linear stiffness a=1a=1, cubic stiffness b=1b=1, and damping c=0.2c=0.2. A forcing-law parameter is θ=(A,f,ση,ℓη)∈[0.5,2.0]×[0.10,0.30]×[0.05,0.40]×[0.10,1.00].θ=(A,f, _η, _η)∈[0.5,2.0]×[0.10,0.30]×[0.05,0.40]×[0.10,1.00]. For each realization, independently sample Φ∼Unif[0,2π] [0,2π] and a zero-mean Gaussian process η with squared-exponential covariance Kθ(s,t)=ση2exp(−(s−t)22ℓη2).K_θ(s,t)= _η^2 \! (- (s-t)^22 _η^2 ). The resulting stochastic forcing path is X(t)=Asin(2πft+Φ)+η(t).X(t)=A (2π ft+ )+η(t). Thus (A,f,ση,ℓη)(A,f, _η, _η) identifies a probability law, whereas the phase and Gaussian-process path are realization-level random variables and are resampled for every particle. We use a Latin-hypercube design to choose 1,200 parameter vectors over the four-dimensional box above, providing space-filling coverage near potentially sharp response transitions. The split is performed at the measure level: 1,000 settings are used for training and 200 for testing. No realizations from one forcing law appear in multiple splits. An optional extrapolation benchmark may instead reserve boundary bands in A or f, but it is kept separate from the primary interpolation split. Both forcing and displacement paths are stored at Nt=256N_t=256 equally spaced observation times tk=20k255,k=0,…,255.t_k= 20k255, k=0,…,255. The ODE is not solved directly on this coarse grid. Each forcing realization is generated on a fine grid with step size Δtsolve=0.005 t_solve=0.005 and interpolated inside the first-order system y˙1=y2,y˙2=X(t)−0.2y2−y1−y13. y_1=y_2, y_2=X(t)-0.2y_2-y_1-y_1^3. We integrate from zero initial conditions with adaptive Dormand–Prince RK45, using relative tolerance 10−810^-8 and absolute tolerance 10−1010^-10, and then sample the solution on the observation grid. Before full generation, a convergence study on representative parameter settings compares Δt∈0.02,0.01,0.005,0.0025 t∈\0.02,0.01,0.005,0.0025\ and requires relative trajectory error below 10−410^-4 against the finest reference. For each parameter vector θi _i, let μi _i be the corresponding forcing-process law. Draw two independent ensembles, Xi,j∼iidμi,Xi,j′∼iidμi,j=1,…,128,X_i,j iid _i, X _i,j iid _i, j=1,…,128, and propagate only the second ensemble through the fixed Duffing solution operator, Yi,j=(Xi,j′)Y_i,j=S(X _i,j). The empirical input and target laws are therefore μ^X(i) μ_X^(i) =1128∑j=1128δXi,j, = 1128 _j=1^128 _X_i,j, μ^Y(i) μ_Y^(i) =1128∑j=1128δYi,j. = 1128 _j=1^128 _Y_i,j. This deliberately removes trajectory-level correspondence while retaining the shared measure parameter. The stored arrays X and Y each have shape 1,200×128×2561,200× 128× 256 in measure,sample,time order, and a 1,200×41,200× 4 metadata array stores the parameter vectors for diagnostics and split construction, not as model input. For the neural configurations, a fixed random permutation assigns the first 1,000 law indices to training and the remaining 200 to testing. The training loader is shuffled; the test loader is not shuffled. Mini-batches contain eight laws, use no worker processes, and contain all 128 particles from each selected law. D.2.2 Distributional Operator Implementation Figure 6: Process-valued Distributional Operator used for the Duffing task. Top: the input-process ensemble is projected to finite-dimensional coordinate laws and encoded into a condition through expectation features. Bottom: the condition and independent latent draws drive a conditional generator, whose projected output samples are reconstructed as response-process paths. Separate PCA transforms are fitted to the pooled training forcing and response trajectories, never to test data. The 256-dimensional forcing and response grids are truncated to 24 and 13 principal coordinates, respectively. PCA fitting accumulates the training mean and scatter matrix in 64-law blocks in double precision and uses a symmetric eigendecomposition. Each retained coordinate is then standardized using its population mean and standard deviation over all trajectories in the training laws; test data reuse these stored transformations. Near-constant coordinates use unit scale. Inputs to the learned encoder are the normalized 24-dimensional forcing coefficients. Generated normalized 13-dimensional response coefficients are unstandardized and reconstructed on the original 256-point grid before the training loss and reported metrics are evaluated. The DeepSets particle map and outer map each have two width-256 SiLU hidden layers. The particle embedding and law context both have dimension 128, and the aggregation is the sample mean. The conditional RealNVP has six affine coupling layers with alternating even–odd masks. Each scale–translation conditioner takes the masked 13-dimensional coefficient vector together with the 128-dimensional context, uses two width-256 SiLU hidden layers, and has a zero-initialized output layer. Its log-scale is bounded as s~=2tanh(s/2) s=2 (s/2). The base law is (0,I13)N(0,I_13), temperature is one, and 128 particles are generated for each input law. Training minimizes the mean debiased GeomLoss Sinkhorn divergence with p=2p=2, blur 0.010.01, scaling 0.50.5, automatic backend selection, and no auxiliary mean loss. AdamW uses learning rate 10−310^-3 and weight decay 10−310^-3; the global gradient norm is clipped to one. Training runs for at most 1,000 epochs. D.2.3 Baseline Implementations The Fixed-Feature MLP operates on the raw 256-point paths without PCA. A two-hidden-layer width-64 SiLU network encodes each forcing path into 64 features. A DeepSets module with two width-64 SiLU layers maps these features to 128-dimensional sample embeddings, concatenates their population mean and standard deviation (with variance stabilization 10−610^-6), and applies two more width-64 SiLU layers to produce a 64-dimensional law context. Given this context and 32-dimensional standard Gaussian noise at temperature one, a three-hidden-layer width-64 SiLU MLP generates a 32-dimensional output code; a two-hidden-layer width-64 SiLU decoder maps each code directly to a 256-point response path. It generates 128 paths per law and trains for 200 epochs with batches of eight laws. Its debiased Sinkhorn loss uses p=2p=2, blur 0.050.05, scaling 0.50.5, and automatic backend selection. AdamW uses learning rate 10−310^-3 and weight decay 10−210^-2, with global gradient norm clipped to one. The Nadaraya–Watson baseline computes Sinkhorn divergences dkd_k between a query forcing law and each of the 1,000 training laws, then uses wk∝exp[−dk/(2h2)]w_k [-d_k/(2h^2)]. Its 64 predicted response particles are the elementwise weighted averages Y^j=∑k=11000wkYk,j,j=1,…,64, Y_j= _k=1^1000w_kY_k,j, j=1,…,64, of the stored target-particle arrays, rather than samples from a discrete mixture over training laws. The common train–test split and a fixed subset of 64 particles are used for every law. Input-law distances are debiased GeomLoss Sinkhorn divergences with p=2p=2, blur 0.050.05, scaling 0.50.5, the tensorized backend, four reference laws per distance chunk, and two query laws per batch. The bandwidth is fixed at h=0.005h=0.005, and the reported baseline is a single deterministic fit. D.2.4 Evaluation Metrics For a target empirical law ν=n−1∑i=1nδyiν=n^-1 _i=1^n _y_i and a predicted empirical law ν^=m−1∑j=1mδy^j ν=m^-1 _j=1^m _ y_j, let Π(ν^,ν) ( ν,ν) denote their couplings. With ground cost c and regularization ε>0 >0, define OTε(ν^,ν) _ ( ν,ν) =infπ∈Π(ν^,ν)∫c(x,y)dπ(x,y)+εKL(π∥ν^⊗ν), = _π∈ ( ν,ν) \ c(x,y)dπ(x,y)+ (π\| ν ν) \, Sε(ν^,ν) S_ ( ν,ν) =OTε(ν^,ν)−12OTε(ν^,ν^)−12OTε(ν,ν). =OT_ ( ν,ν)- 12OT_ ( ν, ν)- 12OT_ (ν,ν). For a positive-definite kernel k, the biased empirical MMD and empirical energy distance are MMD^k MMD_k =[1m2∑j,j′k(y^j,y^j′)+1n2∑i,i′k(yi,yi′)−2mn∑j,ik(y^j,yi)]1/2, = [ 1m^2 _j,j k( y_j, y_j )+ 1n^2 _i,i k(y_i,y_i )- 2mn _j,ik( y_j,y_i) ]^1/2, ED ED =2mn∑j,i‖y^j−yi‖2−1m2∑j,j′‖y^j−y^j′‖2−1n2∑i,i′‖yi−yi′‖2. = 2mn _j,i\| y_j-y_i\|_2- 1m^2 _j,j \| y_j- y_j \|_2- 1n^2 _i,i \|y_i-y_i \|_2. For Pθ(z)=⟨θ,z⟩P_θ(z)= θ,z and uniform measure σ on the unit sphere, sliced 2-Wasserstein distance is SW22(ν^,ν)=∫W22((Pθ)#ν^,(Pθ)#ν)σ(θ).SW_2^2( ν,ν)= W_2^2\! ((P_θ)_\# ν,(P_θ)_\#ν )dσ(θ). The model-specific training losses and schedules are given in Appendices D.2.2 and D.2.3. Final evaluation generates 128 predicted paths for each target law. We report the model’s configured debiased Sinkhorn divergence, biased RBF MMD, sliced W2W_2, and empirical energy distance on the reconstructed 256-point paths. For MMD, a separate median-heuristic bandwidth is computed for each law from all prediction–target squared distances and the nonzero within-sample squared distances; the reported value is MMD, not MMD squared. Sliced W2W_2 uses 128 fixed unit-normalized Gaussian projections. Metrics are averaged over laws. The conditional-flow values in Table 2 are test-set mean ± sample standard deviation; kernel regression is evaluated once. Appendix E Supplementary Experiments The controlled Gaussian study below is a purely supplementary experiment. Its complete setup, results, and implementation details are reported here. E.1 Controlled Gaussian Distributional Mapping E.1.1 Problem Setup This controlled task learns a map from four-dimensional Gaussian-mixture laws to four-dimensional Gaussian laws from unpaired ensembles of 200 input and 200 output particles. The target mean and covariance are nonlinear functions of analytic input-law moments and Fourier features, while the Distributional Operator sees only the empirical input law and predicts a full-covariance Gaussian. Complete data-generation, baseline, architecture, optimization, and metric definitions are given in Appendix E.1.4. E.1.2 Model Instantiation The Distributional Operator applies a shared particle encoder ϕθ:ℝ4→ℝ64 _θ:R^4 ^64 to each input sample and averages the resulting features. The pooled representation is processed by an outer network ρθ _θ and a Gaussian output head: h(b) h^(b) =1200∑j=1200ϕθ(xj(b)), = 1200 _j=1^200 _θ(x_j^(b)), (μ^(b),u^(b)) ( μ^(b), u^(b)) =ρθ(h(b)), = _θ(h^(b)), L^(b) L^(b) =ℛltr(u^(b)),Σ^(b)=L^(b)L^(b)⊤. =R_ltr( u^(b)), ^(b)= L^(b) L^(b) . The particle encoder has two width-64 GELU hidden layers, and the outer network has three width-128 GELU hidden layers. The Gaussian head outputs four mean coordinates and ten entries of a lower-triangular Cholesky factor. Softplus diagonal entries with an offset of 10−510^-5 guarantee a positive-definite covariance. Thus the predicted law is (μ^(b),Σ^(b))N( μ^(b), ^(b)) and is invariant to the ordering of the input particles. E.1.3 Evaluation Training minimizes empirical output-sample NLL. We additionally report the closed-form Gaussian W2W_2 distance, target-to-prediction KL divergence, and Hellinger distance, averaged over the 200 test laws. Full metric definitions and optimization details are in Appendix E.1.4. We compare against the Fixed-Feature MLP and distribution-space kernel regression; see Appendix E.1.4 for their implementations. Model NLL W2W_2 KL divergence Hellinger distance Distributional Operator 4.016462±0.002379 4.016462± 0.002379 0.162049±0.004262 0.162049± 0.004262 0.048246±0.001804 0.048246± 0.001804 0.099313±0.002195 0.099313± 0.002195 Fixed-Feature MLP 4.033581±0.0026874.033581± 0.002687 0.187644±0.0035890.187644± 0.003589 0.066297±0.0023620.066297± 0.002362 0.112535±0.0023480.112535± 0.002348 Kernel regression 4.776576±0.0000004.776576± 0.000000 0.741837±0.0000000.741837± 0.000000 0.806815±0.0000000.806815± 0.000000 0.396992±0.0000000.396992± 0.000000 Table 3: Test results for the controlled Gaussian mapping. Entries are the mean ± standard deviation over five runs. Kernel regression is deterministic and therefore has zero standard deviation. Lower is better for every metric. Figure 7: Predicted versus true output means on the test set for the trained Distributional Operator. Each panel shows one of the four mean components, with one point per test distribution; the dashed diagonal y=xy=x denotes perfect prediction. Figure 8: Elementwise mean absolute error of the predicted output covariance matrices on the test set for the trained Distributional Operator. For covariance entry (i,j)(i,j), the displayed value is Ntest−1∑k=1Ntest|Σ^ij(k)−Σij(k)|N_test^-1 _k=1^N_test| _ij^(k)- _ij^(k)|; darker errors identify components that are systematically more difficult to predict. Figure 7 compares the predicted and true mean components directly, while Figure 8 summarizes covariance error across all test distributions. E.1.4 Additional Details Data Generation For each of the 1,200 laws, the number of mixture components CkC_k is sampled uniformly from 1,2,3\1,2,3\. Independent unit-rate exponential variables are normalized to give mixture weights πkc=ukc∑c′=1Ckukc′. _kc= u_kc _c =1^C_ku_kc . Component means and coordinatewise standard deviations are sampled independently as mkcj∼Unif[−2,2],skcj∼Unif[0.15,0.80],m_kcj [-2,2], s_kcj [0.15,0.80], giving the input law Pk=∑c=1Ckπkc(mkc,diag(skc2)).P_k= _c=1^C_k _kcN\! (m_kc,diag(s_kc^2) ). We draw S=200S=200 independent input particles from each PkP_k. The deterministic law feature zk∈ℝ24z_k ^24 concatenates the analytic mixture mean, elementwise variance, and eight sine and eight cosine Fourier features. Eight frequency vectors ωr∼(0,I4) _r (0,I_4) are sampled once and shared by the dataset. Specifically, m¯k m_k =∑cπkcmkc, = _c _kcm_kc, vk v_k =∑cπkc(skc2+mkc2)−m¯k2, = _c _kc(s_kc^2+m_kc^2)- m_k^2, akr a_kr =∑cπkcexp(−12∑jskcj2ωrj2)sin(ωr⊤mkc), = _c _kc \! (- 12 _js_kcj^2 _rj^2 ) ( _r m_kc), bkr b_kr =∑cπkcexp(−12∑jskcj2ωrj2)cos(ωr⊤mkc), = _c _kc \! (- 12 _js_kcj^2 _rj^2 ) ( _r m_kc), zk z_k =[m¯k,vk,ak,bk]. =[ m_k,v_k,a_k,b_k]. One fixed random map is shared by all laws. Every entry of Wμ,Wd,WoW_μ,W_d,W_o is sampled independently from (0,1/24)N(0,1/24). The output mean is μk=Wμzk _k=W_μz_k. A lower-triangular matrix LkL_k is constructed by (Lk)ii=0.25+0.751+exp[−(Wdzk)i]−1,(L_k)_i=0.25+0.75\1+ [-(W_dz_k)_i]\^-1, and (Lk)ij=0.2tanh((Wozk)ij)(L_k)_ij=0.2 ((W_oz_k)_ij) for i>ji>j. Thus Σk=LkLk⊤+10−4I4,Qk=(μk,Σk). _k=L_kL_k +10^-4I_4, Q_k=N( _k, _k). Finally, 200 independent output particles are drawn from QkQ_k. The analytic targets support evaluation but are not used by the principal training loss. Baseline Implementations The Fixed-Feature MLP computes the empirical mean and unbiased empirical covariance of each input ensemble, concatenates all four mean and sixteen covariance entries, and maps the resulting 20 features through four width-128 GELU layers. The same full-covariance Gaussian head used by the Distributional Operator returns the predicted mean and covariance. It uses the same split, objective, optimizer, scheduler, and evaluation procedure as the main model. The parameter-free kernel estimator stores empirical input and output means and covariances for all 1,000 training laws. For a query, it computes Gaussian W2W_2 distances dkd_k between its empirical input moments and those of each training law, then assigns radial-basis weights αk=exp[−dk2/(2h2)]∑ℓ=11000exp[−dℓ2/(2h2)],h=1. _k= [-d_k^2/(2h^2)] _ =1^1000 [-d_ ^2/(2h^2)], h=1. It interpolates empirical output moments as μ^=∑kαkm~ky,Σ^=∑kαkC~ky+10−5I4. μ= _k _k m_k^y, = _k _k C_k^y+10^-5I_4. The covariance is symmetrized before stabilization and likelihood evaluation. Evaluation Metrics A fixed random permutation defines the common train–test split. Five runs use mini-batches of 64 under the default training configuration. For target p=(μ,Σ)p=N(μ, ) and prediction q=(μ^,Σ^)q=N( μ, ), we report empirical output-sample NLL, Gaussian W2W_2, target-to-prediction KL divergence, and Hellinger distance. The latter two are DKL(p∥q) D_KL(p\|q) =12[tr(Σ^−1Σ)+(μ^−μ)⊤Σ^−1(μ^−μ)−4+logdetΣ^detΣ], = 12 [tr( ^-1 )+( μ-μ) ^-1( μ-μ)-4+ ], H(p,q) H(p,q) =1−BC(p,q), = 1-BC(p,q), BC(p,q) (p,q) =(detΣ)1/4(detΣ^)1/4det((Σ+Σ^)/2)1/2exp[−18(μ−μ^)⊤(Σ+Σ^2)−1(μ−μ^)]. = ( )^1/4( )^1/4 (( + )/2)^1/2 \! [- 18(μ- μ) ( + 2 )^-1(μ- μ) ]. Each metric is averaged over the 200 test laws; NLL is additionally averaged over the 200 observed output particles per law.