Stochastic Process Learning via Operator Flow Matching
Authors: Yaozhong Shi
Zachary E Ross
Domniki Asimaki
Kamyar Azizzadenesheli
Venue: NeurIPS 2025 spotlight ()
Keywords: Operator Learning
Fourier Neural Operator
Bayesian Methods
Gaussian Process
AI for Science
Primary area: machine_learning_for_sciences

Abstract: Expanding on neural operators, we propose a novel framework for stochastic process learning across arbitrary domains. In particular, we develop operator flow matching (OFM) for learning stochastic process priors on function spaces. OFM provides the probability density of the values of any collection of points and enables mathematically tractable functional regression at new points with mean and density estimation. Our method outperforms state-of-the-art models in stochastic process learning, functional regression, and prior learning.

Stochastic Process Learning via Operator Flow Matching
Yaozhong Shi California Institute of Technology yshi5@caltech. edu
Zachary E. Ross
California Institute of Technology zrossêcaltech, edu
Domniki Asimaki
Kamyar Azizzadenesheli
California Institute of Technology dommiki@caltech. edu
NVIDIA Corporation kaazizzad@gmail. com
Abstract
Expanding on neural operators, we propose a novel framework for stochastic process learning across arbitrary domains. In particular, we develop operator flow matching (OFM) for learning stochastic process priors on function spaces. OFM provides the probability density of the values of any collection of points and enables mathematically tractable functional regression at new points with mean and density estimation. Our method outperforms state-of-the-are models in stochastic process learning, functional regression, and prior learning.
Introduction
Stochastic processes are foundational many domains, from functional regression and physics- based data assimilation, to financial markets, geophysics, and black box optimization. Stochastic processes can serve as prior distributions over functions and can provide the density of any finite collection of points. Conventionally, these priors are designed by hand from predefined Gaussian processes (GP) and their variants, and therefore assume that they adequately describe the phenomena of interest. However, many processes modeled in the natural world not well described by GP, Fig 2. Such models limit the ficxibility and generalizability of these stochastic processes in world applications, leaving behind significant challenges for more general stochastic process learning (SPL).
In SPL. the prior over the stochastic process is learned from data, i.e., set of historical point evaluations in past experiments. Learning the prior over the process is crucial for universal functional regression (UFR), which recently proposed Bayesian scheme for functional regression and takes GP-regression special case when the prior is Gaussian [1]. UFR is important to scientific and engineering domains, including reanalysis, data completion-assimilatioon uncertainty quantification, and black box optimization.
In this paper, we introduce a novel operator learning framework termed operator flow matching (OFM) for learning priors over stochastic processes through the joint distribution of any collection of points. To achieve this, we theoretically and empirically generalize marginal optimal transport flow matching [2] to infinite-dimensin
function spaces where we map GP into prior over function spaces through flow differential equation. We then derive SPL from the function space derivation and learn prior over arbitrary sets of points. For SPL, we map any collection of pointwise evaluations of GP to pointwise evaluations of target functions. This allows us to learn prior distributions over more general stochastic processes, hence enabling sampling values of any collection of points with their associated density and facilitating efficient UFR. We leverage this capability by extending neural operators [3] dd]-des initially to map functions between
spaces-to to maps
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

between collections of points deploying their functional convergence properties. This serves as the essential architecture block OFM.
After learning the prior and having access to the densities, OFM can be used for UFR, where given any collection of points of the underlying function, estimate the posterior
value of any new collection of points and efficiently sample from their posterior values using stochastic gradient Langevin dynamics (SGLD) [4], i.e., Gaussian sampling the input GP space (Fig 1). We show that OFM significantly outperforms state-of-the-art (SOTA) methods, including deep GPs, conditional models, and operator flows (OrFLow) [1, 5-10].
Generalizing Gerressssson to regression over general stochastic process in a practical and imple- mentable way demands unified framework that spans many key fields. Because readers will have diverse backgrounds and expertise, we provide set potential questions and answers in Appendix A further clarify the development and highlight contributions.
Prior Learning over stochastic process
Posterior sampling given partial observations
(
log (u(x)) K((x))
Partial noley observations
OFM
posteror distribution
U
0 U1
Bijective
OFM
U
Ug
Postertor samples
GP semple
data sampte
trainable
prior "Un: space UI space
Two-phase strategy for prior learning and posterior sampling. In the prior-learning phase, OFM leverages the marginal optimal-transport path to learn bijective mapping between predefined GP and the unknown stochastic process that generates the training data. In the posterior-sampling phase, the learned prior is frozen; given noisy, partial observations, the exact posterior is obtained Bayes' theorem. SGLD, aided by the Hutchinson trace estimator, then enables efficient and robust sampling.
Figure I:
Lastly, our contributions are summarized as follows:
First work to extend flow matching / stochastic interpolants / rectified flow [11-13] to stochastic processes via operator learning. The formulation enables likelihood estimation for function values at any collection of points for the target stochastic process. We also contribute to the development marginal (dynamic) optimal transport in infinite-dimensional flow matching through optimal coupling and dynamic Kantorovich formulations.
First integration of flow matching with functional regression yields a unified framework for prior and posterior sampling that is applicable both generation and regression tasks.
A practical generalization of GP regression that provides the exact prior and posterior density over an unknown stochastic process (whether or non-GP). contrast, previous methods such as deep GPs and conditional models work with approximate posteriors. Our method achieves SOTA performance in all challenging functional regression tasks.
This work provides unified perspective bridging several important fields, which opens new research directions for problems in science and engineering. Additionally, we present extensive ablation and scaling studies that demonstrate the effectiveness of each component in our framework.
Related Work
Neural operators. Neural operators constitute paradigm in machine learning for learning maps between function spaces, generalization of conventional neural networks that map between finite di- mensional spaces [14, 15]. Among neural operator architectures, Fourier neural operators (FNO) [14] enable convolution the spectral domain and are effective for operator learning [16-21]. In this work, use this as neural operator architecture.

---

15
Observations
(b) Ground truth
OFM mean
(d) OFM sample
(e) GP sample
(I) ConvCNP mean
(g) CN sample
Figure 2: Operator Flow Matching (OFM) regression on Navier-Stokes functional data with resolution
64 x 64. (a) 32 random observations (only 0.7%). (b) Ground truth sample. (c) Predicted mean from OFM. (d) One posterior sample from OFM. (e) One posterior sample from the best fitted GP. (f) Predicted mean from ConvCNP. (g) One posterior sample from ConvCNP.
Direct function samples. There is body of work on generative models dedicated to learning distributions over functions, such that direct sampling on the function space is possible. For example, generative adversarial neural operators (GANO) generalize generative adversarial nets on finite dimensional spaces to function spaces [22, 23), yielding a neural operator generative model that maps GP to data functions [3]. Other works in this area have followed the success of diffusion models [24, 25] in finite dimensional spaces, e.g., denoising diffusion operators generalize diffusion models to function spaces by using GP means to add noise and neural operators to learn the score operator on function-valuee data [26-28]. Moreover, the same principle has been deployed work. to generalize flow matching [11] |] to functional spaces [29], approach closely related However, these works on learning generative models on function spaces do not support UFR the way GP-regression does because they (i) focus solely generating function samples, (ii) do not clarify how to model a stochastic process on point value sequential generation, and (iii) do not provide point evaluation of probability density. In contrast, OFM enables exact density calculation and conditional posterior sampling.
Stochastic processes. Earlier works SPL have focused on hand-tunce tenettundd methods in the style of Ge-regression. these cases, an expert tunes the GP parameters given set of experimental samples. More advanced methods rely deep GPs, in which network of GPs is stacked on top of each other. The parameters deep GPs are commonly optimized by minimizing the variational free energy, which serves bound on the negative log marginal likelihood. [30, 31]. Deep GPs have limitations terms of learnability, expressivity, and computational complexity. Warped GPs [32] and transforming GP [33] methods use historical data to learn pointwise transformation of GP values and achieve on par performance compared to deep GP type methods. The pointwise nature of such approaches limits their generality.
Another line of work proposes learning the conditional distribution of point values, inspired by variational inference and designed for sampling from function spaces; we refer to this line of models as conditional models' [7]. This method trains a model to map any collection points and their values to vector, used as an input decoder that maps any collection of points to their values. The architectures used in these models (including the decoder) are not mathematically consistent as the number of points grows, limiting the approach to finite dimensions. In particular, Convolutional Conditional Neural Processes (ConvCNPs) [6] use convolutional architectures to capture local, translation-invarian structure, but they assume stationarity and struggle with long-range dependencies or highly non-uniform data. The diffusion based variants [34] also use uncorrelated Gaussian noise, and the results do not exist in function spaces [22, 26]. Furthermore, methods based conditional models are unable to provide density estimation for collections of points, as needed for UFR and SPL in general.
Separately, OPFLOW introduced invertible neural operators that are trained to map any collection of points sampled from GP to new collection of points in the data space [1], using the maximum likelihood principle. This method mathematically consistent as the resolution grows, captures
Neural process (NP) is prominent example of conditional models

---

the likelihood of any collection of points, and allows for UFR using SGLD. However, similar to normalizing flow [35] methods in finite dimensional domains, the use of invertible deep learning models makes their training challenge, particularly with regard to expressiveness, as also described in the original OPFr.ow work. Finally, we strongly encourage readers to consult Appendix and Q for an in-depth comparison of stochastic process learning and other generative frameworks.
3 Operator Flow Matching
Here, we introduce the problem setting and notations used for OFM in function space. We recommend that readers consult Appendix B-E for a foundational overview of SPL, UFR and related background. Subsequently, present the framework of OFM, which extends marginal optimal transport flow matching [2] to dufinite-dimensin
spaces. We further demonstrate the generalization of flow matching to stochastic processes as it is induced from OFM function spaces. Finally, illustrate how to evaluate exact and tractable likelihoods for any point evaluation of functions using OFM, making it applicable in the UFR setting.
For real separable Hilbert space (H, (t--), II-II): equipped with the Borel o algebra of measurable sets denoted by B(H), we introduce two measures on B(H), as the reference measure and #1 the data measure. Consider function ho sampled from vo, such that ho ~ vo. For a smooth time-varying infinite dimensional vector field G. H + H. we define an ordinary differential equation (ODE)
00, Ge (ho))
(1)
Ōt
with initial condition Po (ho)) = hyg where h = h4 = 2 (ho) represents a function hg transported along vector field from time 0 to time t. The diffcomorphism &, induces a pushforward measure H4 := (Ho)) with (4g = Mo. and we refer to H4 as the path of probability measure. The goal is to construct path of probability measure such that at t = /a 11- The dynamic relationship between the time varying measure 14 and vector field G. can be characterized by the continuity equation:
(prGr)
Ōt
In practice we use Eq. 2 in its weak form [29, 36] to check whether given vector field 9, generates the target μt:
J. J @x(h, Øl + (Go(h), (n.p(h,
Where yp E Cyl(H x [0, 1]), and Cyl(H x [o, 1I) is the space of smooth cylindrical test functions. Suppose that the varying vector field G. and the induced Me satisfying Eq. are known. We parameterize G. with a neural operator Go : [0, 1] x H + H and regress Go to target G. through flow matching objective.
CtM = 10.1 1| (t. h) - (h)l2
(4)
However, similar to its finieedimeensionl counterpart, G, is typically unknown. Moreover, there are infinitely many paths of probability measures that satisfy Eq. and ensure #1 v1. Therefore, is necessary to specify a path of probability measures to effectively guide the learning of Ga.
By constructing the appropriate Gaussian and conditional probability measures and leveraging optimal coupling together with the dynamic Kantorovich formulation [37], we propose marginal (dynamic) optimal- optimal-transport flow matching in function space; detailed theory development and proofs provided in Appendix and G. In the next subsection, we show how to generalize flow matching to stochastic processes, which is induced by the function space derivation.
3.1 Generalizing Flow Matching to Stochastic Processes
Stochastic processes are inherently of points (Brémaud [38], Chapter 1). We generalize the above marginal optimal transport flow matching on function spaces to stochastic processes by defining the transport map on any collection of points. We then show that, as the collection of points covers the space in the limit, this generalization recovers dimente-dimeniooall flow matching implemented with neural operators.
and define distributions over any collection

---

For any n and points {I1, D2........... Fn), consider an ODE system in which a vector of random variables un E R" is gradually transformed into u, E IR", for which, the ith entry is equal u(x,). via a smooth, time-varying vector field, denoted by 9. with abuse of notation.
U4 == 4; = uo + L. Go(u))ds
(5)
Here, the neural operator is applied a collection of point evaluations. Given the set of points and the density of Po := Po (fuo(r1), uo uo (Ix,,) }}) = N (0, (2 [nn)), }),, where K n n covariance matrix with entries described by kernel function k(z;, and uo po. the time-varying density Pt induced by the diffeomorphism D, or Gj can be computed, extending the transport equation Eq. collections of points,
Opr(ue) (V (Ge))))(u)) Öt
Eq. shows that constructing Pt is equivalent to constructing G. for finite entries for which the analysis carries to finite collections of random variables. In the following, we refer Pa as the marginal probability path induced by G. for the given collection of points. From Eq. 6, the log density can be computed through integration,
log pa (w,) = log Po(un) - f((v (u)(u,)s
(7)
In this formulation, we are seeking a specific vector field that transports density 4o to target density qi for any and any collection of points (x1, ,2,,.., with boundary conditions po 2o,Pi Q!- Extending optimal transport flow matching to stochastic processes, we parameterize the vector field G. with neural operator Ge. which is optimized through the flow matching objective for SPL,
CFM := sus
sup (21....fa)
100))(0, (t un) - (2))))
(8)
Note that Pt and uy depend on the point collocations. In the above equation, the suprema are intractable and we replace them with an expectation soft approximation (see Appendix A. Q5 for a detailed discussion of the approximation). Moreover, the true G. is usually unknown and to address it, we derive a probability path conditioned on latent variable of the same alphabet size as the collection. Consequently, the marginal probability path pe (u,) is a mixture of conditional probability paths (uel))).
pr(u)) = / pc(uelz)((z)dz
(9)
Ge(we) = Ba( Ge(uel=)pn(uu))
(10)
(wt)
Given Eq. 10, the conditional flow matching (CFM) objective is defined as,
LCFM == E,,,,,,.,,, Be.ggg....,...gg..ppp....... (o ux) - Gr(ue|=)11
(11)
Eqs. 11 and 10 have an identical gradient for 8, i.c. (O = (0). Inspired by the finite dimensional developments [2], the variable is chosen couple (uo, u1) from the coupling π (ug. un) = q(2), which is achieved by minimizing the dynamic 2- Wasserstein distance,
Wayn(yo, 0)) = inf /
(I.. dudd
(12)
pr.9. 22 Jo
Under mild conditions on R", this is equivalent to the statie 2- Wasserstein distance,
Wan (go, 41)2 = inf / well
-u j ||u - wall 2 (eddduuo u1).
(13)
Considering the class of Gaussian conditional probability paths pr (z) ((uglm(()) (z)) (f21, .2,... 22 ,,, x)), with conditional flow (uull) = T + m
Specially, we choose m. = luj + (1 - 1)ul and a. = o, where a > is a small constant. Then we have the following closed-form expression for the correspondngg vector field (full derivation provided in Appendix H)
( d1 uo
(14)
5

---

With the aforementioned developments, for any collection of points, we transport Gaussian distribu- tion a target distribution. The Gaussian distribution drawn from a GP. with its covariance matrix K (a , .Fn) determined by the kernel function k(zi, xi)) of the GP. According to Kolmogorov extension theorem (KET) [39], there exists a valid stochastic process whose lnite-dimensional marginal is the pushforward distribution under Gg. This demonstrates that the generalization of flow matching to dimesio-dimeniionll spaces with neural operators naturally induces the generalization of flow matching to stochastic processes. In the scenario where the limit of points covers the space, these two become equivalent. For a detailed explanation and proof, please refer to Appendix C and D.
Our framework extends seamlessly to alternative probability paths, such as those in stochastic interpolants and rectified flow [12, 13]. because the generalization to stochastic processes is decoupled from any specific path. However, we focus on the OT path in this work for two reasons: (1) ablation and scaling studies show that it accelerates inference and enables more accurate prior learning compared to independent coupling (Table 7, and Appendix P); (2) theoretically, the OT path (straight line) simplifies the Jacobian evaluation required for likelihood estimation (discussed in the next subsection), yieldiin greater stability and speed than arbitrary paths.
3.2 Likelihood Estimation and Bayesian Universal Functional Regression
We parameterize Go with FNO [14] ensure our model is resolution agnostic, and assume Go learns a map from P/ to 11, which serves as the prior. In practice, we deal with discretized evaluations of functions that may have different sampling rate and resolution. For instance, consider function sampled from V1. observed on collection of points i1 := {u(z1), u(z2)..... .... ul u(=m) }:: thus we have a density function P(u1) defined on collection of points {-1, Z2, ..... Im ), where P(u1) is derived from measure V/1 This is similar to how a multivariate Gaussian distribution can be derived from a Gaussian measure characterized by a Gaussian process. Therefore, we can rewrite Eq. 7 as:
log P(un) log P(wo) [' ['(V. -o))(u,)
(15)
where uo is drawn from the reference Gaussian measure Vy, which also defined on the collection of points {æ,, I2, ..., Fm h. Thus, P(uo) is a multivariate Gaussian with tractable density function. Furthermore, with the probability density function PP(u)), we can evaluate the precise likelihood of any un from PP( via Eq. 15. However, following a similar argument to Grathwohl et al. [40], the computation of V GG((u) incurs a cost of O(m²) where m is the cardinality of (#1, 32, X2,,..: F. This quadratic time complexity renders the likelihood calculation prohibitively expensive. To address this issue, we adopt the strategy proposed in Grathwohl et al. [40]. using the unbiased Skilling-Hutching trace estimator [41, 42] to approximate the divergence term. This technique reduces the computation cost to O(m). which is the same as the cost of inference, thereby streamlining the evaluation process. The estimator is implemented as follows:
3go (Go(u. t)
V Go(ue) = Excesle"
(16)
Ou
In the unbiased trace estimator, the random variable & is characterized by RE()) = 0 and Cov(e) T. This estimator benefits from the optimal transport nature of the map which gives rise to a direct line. The gradient computation in Eq. 16 can be efficiently handled with reverse-mode automatic differentiation, allowing for precise estimation with arbitrary error by averaging over a sufficient number of runs, which can benefit from parallel computing of GPUs.
With the efficient tool established for estimating the likelihood of any discretized function samples, we now turn our attention to UFR, i.e., Bayesian functional regression with learned prior. Consider a collection of pointwise observations of the underlying unknown function drawn from V/, corrupted with Gaussian noise, denoted as (ū(x1), ((22)..... ... ū(x,)) or (ū(x:)) = We specifically focus on Gaussian white noise characterized by NO0 "),, such that ((;;) = m(x;) +6; for E (1..... n) (depending on the nature of the problem, the noise may also be considered as correlated GP noise). In UFR setting, we are interested in the posterior distribution over new m 2 n points that include the observation points.
Proposition 3.1. Given noisy observations (((x)))] the posterior distribution is
log I. ( (fu(z;)}m ( -. [x:)) u( log P (fu(x;) )m +C (17)
20

---

Where the constant C' = log (2go²) log ([ū(2;))
Proof. This is derived from Bayes' theorem, along with the translation invariance property of the Gaussian distribution, see the full proof in Appendix
Given this form of the posterior distribution, we adopt SGLD [4] to efficiently sample from it, and then derive statistical features of interest, e.g. mean, maximum a posteriori, and posterior uncertainty, ie., variance, from the posterior samples. More specifically, we follow the posterior sampling strategy developed by Shi et al. [1], which, given an invertible framework, suggests SGLD sampling within the imput GP space where the Gaussian measure My is defined and Langevin dynamics is native to, and then mapping to the data function space (where data measure 1/1 is defined). We also present a scaling study in Appendix P analyzing the variance introduced by the Hutchinson trace estimator, demonstrating its robustness and effectiveness when paired with SGLD sampling. Last, Eq 17 offers unified view of prior and posterior sampling with flow matching models by showing that when no observations are present, the posterior collapses to the prior (up to a constant), making the posterior sampling process identical to that of the prior.
Finally, we provide a plain-language summary of the framework to help the reader better understand the paper (with detailed discussion available in Appendix Q).
OFM is an expressive, flow-based generative model that learns a prior over functions (stochastic process): neural operator parameterizes continuous probability flow that transports samples from a simple reference Gaussian process to data-like functions, yielding explicit prior with tractable density.
It excels at funetional regression by treating functions first-class objects rather than mere pointwise values (unlike NPs), yielding predictions that are consistent across resolutions and arbitrary query sets.
The learned flow is invertible, enabling change-of- of-varibaless for likelihoods and principled Bayesian regression, yielding calibrated uncertainty from few observations.
Experiments
In this section, we demonstrate the superior regression performance compared several baselines across a variety of function datasets, including both Gaussian and highly non- non-Gaussian Processes. As baselines, we employ standard Regression [8], Deep GPs [9, 10), Conditional models [5-7], and OrFLow [1].
For our function datasets, we analyze: (1) Gaussian and non-Gaussian with known posterior, including ID GPs, 2D GPs, and ID Truncated GPs, (TGP). (2) Highly non-GPs, datasets with unknown posterior, such those derived from Navier-Stokes equations [14], black hole dataset from expensive Monte Carlo simulation, and 2D Signed Distance Functions extracted from MNIST digits (MNIST-SDF) [43]. During regression, we assume that the prior Go always successfully trained and remains frozen. Details about the learning process for priors and experimental setup for regression are provided in the Appendix M, O.
65+
85
mein
moin
mean mean
1.
60
Osservatons
samples
(a)
Ground truth
(b) OFM
(c) Uncertainty comparison
Figure 3: OFM regression on GP data. (a) Ground truth GP regression with observed data and predicted samples. (b) OFM regression with observed data and predicted samples. (c) Standard deviation comparison between true GP and OFM predictions.
ID GP data. This experiment replicates the results of classical GP regression, wherein the posterior distributions are preciscly known in closed form. The process involves generating single new
7

---

GieeNatees
EB Soud F74t
68
31
(a) Ground truth
(b) OFM
(c) Uncertainty of OFM
(d) GP prior
(c) Uncertainty of GP
Figure 4: OFM regression TGP data.
realization from the data measure V1. We then select observations at n = 6 randomly chosen positions, incorporating predefined noise level. The posterior is inferred across
128 positions, which includes estimating noise-free values at the observation points. We evaluate our results with two commonly used quantities in the GP literature (1) Standardized Mean Squared Error (SMSE) that normalizes the mean squared error by the variance of the ground truth; and (2) Mean Standardized Log Loss (MSLL), originally introduced by Williams and Rasmussen [8], defined as:
log(2xo²)) + (fu(x))) [2 {ū(z:))?
- -oogp(tu22))
(18)
20
where {ū(x;)]=_ represents observations, { (u(x))) indicate the new positions queried, and the test data (true posterior samples). Meanwhile, [ū(x;))?= )== and o² are predicted mean and variances from the model. We average out SMSE and MSLL over test dataset containing 1000 true posterior samples for models. The performance of each model is detailed Table 1. From Fig. 3, the regression with OFM matches the analytical solution well and provides realistic posterior samples. We further include GP regression tasks using more complex kernels (Gibbs and Rational Quadratic kernel), as shown in Appendix N, OFM consistently outperforms all comparative methods across all metries.
(a) Observations
(b) Ground truth
(c) OFM mean
(d) OFM sample
231
230
338
318
405
111
(e) GP sample
(I) ConvCNP mean (g) CNnVCNP sample
Figure 5: OFM regression black hole data with resolution 64 x 64. 32 random observations. (b) Ground truth sample. (c) Predicted mean from OFM. (d) One posterior sample from OFM. (e) One posterior sample from best fitted GP. (f) Predicted mean from ConvCNP, (g) One posterior sample from ConvCNP
Truncated GP data. In this experiment, we analyze the regression performance of OFM for tractable non-GP. Specifically, we work on truncated [1, 44], which constrains the function amplitude within a specified range. This is achieved by applying sampling-rejection strategy on samples from GP prior. We set the bounds of our TGP to 1-1.2 1.2) and perform regression using observations only at three points, while estimating the posterior across = 128 points. Subsequently, we sample 1000 true TGP posteriors from the GP prior to calculate the mean and standard deviation. Traditional metrics like MSLL and SMSE, which assume Gaussian posterior, not suitable for TGP. We

---

Table 1: Comparison of OFM with baseline models: GP regression; OpFlow [1]; Conditional models:
NP ( [7): Attentive NP ( [5], ANP); Convolutional Conditional NP ( [6), ConvCNP); Deep variational GP ( [9], DGP); Deep Sigma Point Process [10], DSSP); Metries SMSE and MSLL used for ID and 2D GP example. Mean squared error for the predicted mean (μ) and standard deviation (σ) are
used for TGP example. Performance of GP regression for ID and 2D GP removed (marked with
-"), which are taken as the ground truth. Best performance in bold.
Dataset
ID GP
2D GP
ID TGP
Algorithm J. Metric
SMSE
MSLL
SMSE
MSLL
prior
6 14-100
1.6
NP
6.1-10 10-1
4.5- 10
1.7 10
2.1-10
1.0-10-
1.9100-2
ANP
5.1-10-1
9.8. 10
1.6-10"
1.1-10
1.4-0-1
1.7-10-
ConvCNP
5.6-10-1
207-10
1.7-10-1
4.5 4.
16-10-2
2.1-10-3
DGP
4.1-10-1
6.--
1.8-10
4.2-10
4.9-10-1
1.4-10-2
DSPP
4.7-10-1
6.5
1.9-10
6.6-10
1.1-10-2
1.3-10-2
OpFlow
5.0-10-1
2.0-10-1
1.4
1.1.10-1
1.3
3.9-10-3
(Fr(Ours)
4.1-10
5.5
5.5
10
1.3 10
1.6-
5.2 10
9.5
10
evaluate the performance using the mean squared error for both the predicted mean and standard deviation. The results are reported in Table. 1, and illustrated in Fig. 4. OFM accurately learns the specified bounds and provides accurate estimations of mean and standard deviation, along with realistic posterior samples. In contrast, directly applying GP regression exceeds the bounds and yields unrealistic posterior samples.
2D GP data. Similar to the ID GP example, extend our regression analysis to 2D GP data. As shown in Fig. 6 and detailed in Table I, OFM provide accurate posterior estimation. The relative error shown in Fig. 6 is the absolute error normalized by the maximum absolute value of the mean prediction derived from the ground truth GP regression.
Navier-Stokes, Black hole and MNIST-SDE datasets. We collected 2D Navier-Stokes, dataset and applied OFM for the regression. Unlike the GP experiments, where MSLL and SMSE score serve as standard benchmarks, evaluating the performance of models on general non-GPs presents a significant challenge due to the difficulty of determining the true posterior and lack of benchmarks. Moreover, the evidence term in the posterior (Eq. 41) is intractable, so the likelihood cannot serve as a meaningful evaluation metric our setting. A detailed discussion is provided in Appendix A.
We present the predicted mean and a posterior sample in Fig 2 for visual comparison with the ground truth. The predicted mean and the posterior sample are closely aligned with the ground truth. In contrast, traditional GP regression and NP models failed to accurately capture the dynamics of the Navier-Stokes data. In Fig. 5, we conduct similar analysis using a simulated black hole dataset. Here, OFM provides a much more realistic mean and posterior sample that capture the density and swirling patterns of the black hole. Once again, GP regression and NP fail to capture these key statistics. We observe similar outcomes when applying OPM to the MNIST-SDF example (Fig 8), where OFM correctly recognizes the number "7" while GP regression does not.
5 Conclusion
In this paper, we proposed Operator Flow Matching (OFM) for stochastic process learning, which generalizes flow matching models to anfnite-dimensiooall space and stochastic process with optimal transport path. OFM efficiently computes the probability density for any finite collection of points and supports mathematically tractable functional regression. We extensively tested OFM across a diverse range of datasets, including those with closed-form GP and non-GP data, as well as highly non-GP such as Navier-Stokes and black hole data. In comparative evaluations, OFM consistently outperformed all baseline models, establishing new standards in stochastic process learning and regression.
Despite SOTA accuracy, our method is presently limited to low-dimensional domains and demands larger datasets and more compute than GP-based baselines. see Appendix A, and Q for details. Python code available at https: : //githtb.. wwwmysssii//SL..0.

---

Acknowledgments
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Science Foundations for Energy Earthshot under Award Number DE-SC0024705. ZER is supported by fellowship from the David and Lucile Packard Foundation. We also thank Charles Gammie, Ben Prather, Abhishek Joshi, Vedant Dhruv, and Chan for providing the black hole simulations.
References
[1] Yaozhong Shi, Angela Gao, Zachary F. Ross, and Kamyar Azizzadenesheli Universal Functional Regression with Neural Operator Flows, November 2024. URL http://arxiv. org/abs/2404. 02986. arXiv:2404.02986 [cs] version: 3.
[2] Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector- Brooks, Guy Wolf, and Yoshua Bengio. Improving and generalizing flow-based generative models with minibatch optimal transport, March 2024. URL http: /ttp:://a/xivvioo / 00482. arXiv:2302. 20022322.00082 [es].
[3] Kamyar Azizzadeneshelii Nikola Kovachki, Zongyi Li, Miguel Liu-Schiaffini, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and de- sign. Nature Reviews Physics, 6(5):320-328, May 2024. ISSN 2522-5820. doi: 10. 1038/ $422544-04-007712-3 URL https: htopps://wwwt http::://www..atur..m. naw.re. wwtpaaarceeeeeeee-224444444 Publisher: Nature Publishing Group.
[4] Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML' 11, pages 681 681-688, Madison, WI, USA, June 2011. Omnipress. ISBN 978-1- 4503-0619-5.
[5] Hyunjik Kim, Andriy Mnih, Jonathan Schwarz, Marta Garnelo, Ali Eslami, Dan Rosenbaum, Oriol Vinyals, and Yee Whye Teh. Attentive Neural Processes, July 2019. URL http: [cs, stat].
/ //arxiv. org/abs/1901 05761. arXiv:1901.0555 .05761
[6] Jonathan Gordon, Wessel Bruinsma, Andrew K. Foong, James Requeima, Yann Dubois, and Richard E. Turner. Convolutional Conditional Neural Processes, June 2020. URL http: //arxiv. wwt//ogg///////ee/eee////eeeee/e htsp//bb//99..... 13556. arXiv:1910.13555 [stat].
[7] Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo J. Rezende, S. M. Ali Eslami, and Yee Whye Teh. Neural Processes, July 2018. URL http://arxivv wrgg abs/ 1807. 01622. arXiv: 807.00662 [cs, stat].
[8] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processss for machine learning. volume 2. MIT press Cambridge, MA, 2006. Issue: 3.
[9] Hugh Salimheni and Marc Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. Advances in neural information processing systems, 30, 2017.
[10] Martin Jankowiak, Geoff Pleiss, and Jacob Gardner. Deep Sigma Point Processes, December
2020. URL http: arxiv. org/ ab////b///222 abs/2002.09112. arXiv:2002.09 [cs, stat].
[11] Yaron Lipman, Ricky T. Q. Chen, Heli Ben- Ken-Hamu, Maximilian Nickel, and Matt Le. Flow Matching for Generative Modeling. February 2023. URL http::/arxti http//awwwi.rraaraa/////aa/ a2g/abs/2200. 02747. arXiv:2210.02747 [cs, stat].
[12] Michael S. Albergo and Eric Vanden-Eijnden Building Normalizing Flows with Stochastic Interpolants, March 2023. URL. http:///a/// wrg. abs /aa.b///229... 15571. arXiv:2209.1557 (cs, stat].
[13] Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow, September 2022. URI. http: http://axxxv.r//eee / 2209. 03003, arXiv:2209.0320 [cs].
10

---

[14] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differ- ential Equations, May 2021. URL. http://arxiv.... org/ http//abs2200...0000 arXiv:2010.0889, [es, math].
[15] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs. Journal of Machine Learning Research, 24(89):1-97, 2023.
ISSN 1533-7928. URL h jmlr org/ www.oogg.....eee..o https/ppprs//vvv/ 11//2///11110 html.
[16] Yan Yang, Angela Gao, Jorge C. Castellanos, Zachary E. Ross, Kamyar Azizzadenesheli, and Robert W. Clayton. Seismic wave propagation and inversion with Neural Operators, October 2021. URL httpp///aaxxx...//aabseas/ee/e// 054211111.. arXiv:2108.0555 0..411
[17] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li. Kamyar Azizzadeneshelii Pedram Hassanzadeh, Karthik Kashinath, and Animashree Anandkumar. FourCastNet: Global Data- driven High-resolution Weather Model using Adaptive Fourier Neural Operators, February 2022.
URL attp::/arrxe wtt.aaiiiiigg...eaeeeeoeeee www.orgg/////oooooooooooooo/.oo/ 2bb/ 11214. arXiv:2202.112 11214
[18] Gege Wen, Zongyi Li, Qirui Long, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M. Benson. Real-time high-resolution CO2 geological storage prediction using nested Fourier neural operators. Energy & Environmental Science, 16(4):1732-1741, April 2023. ISSN 1754-5706. doi: 10. 00.10399D2FF04244F.. URL. https: / rrcbbb..s.... rsc. org...rg/eee wrggen////ee httg/////ooeeet content/ articlelanding/205 /2222//444.2
Publisher: The Royal Society of Chemistry.
[19] Yan Yang, Angela F. Gao, Kamyar Azizzadenesheli, Robert W. Clayton, and Zachary E. Ross. Rapid Seismie Waveform Modeling and Inversion with Neural Operators, April 2023. URL
http:///rxxv.orggg
11955. arXiv:2209.1 1955.
[20] Hongyu Sun, Zachary E. Ross, Weiqiang Zhu, and Kamyar Azizzadenesheli. Phase Neu- ral Operator for Multi-Station Picking of Seismic Arrivals. Geophysical Research Let- ters, 50024))22233LL1064344 2023. ISSN 1944-8007. doi: 10.1029/2023GL1164444 1029 0029/22236L104434, _eprint:
URL / wil wiley. com/ cttp:////.willlllcc/////////////
[21] Zongyi Li, Nikola Borislavov Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Prakash Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, and Anima Anandkumar. -eometry-Informed Neural Operator for Large-Scale 3D PDEs, September 2023. URL http: http:p////xxvvoiieeee org/ baarg//bb//// 2209 .0089.. arXiv:2309.005:
[22] Md Ashiqur Rahman, Manuel A. Florez, Anima Anandkumar, Zachary E. Ross, and Kamyar Azizzadenesheli. Generative Adversarial Neural Operators, October 2022. URL http:// arxiv. www.g//g/../...../////// 03017. arXiv:2205.03017 [cs, math].
[23] Yaozhong Shi, Grigorios Lavrentiadis, Domniki Asimaki, Zachary E. Ross, and Kamyar Azizzadenesheli. Broadband Ground-Motion Synthesis via Generative Adversarial Neu- ral Operators: Development and Validation. Bulletin of the Seismological Society of America, 114(4):2151-117 August 2024. ISSN 0037-1106, 1943-3573, doi: 10.1785/ 0120230207. URL
[24] Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-Based Generative Modeling through Stochastic Differential Equations, February 2021. /
1345616111111..11111 arXiv:201 1.13456 [cs, stat).
[25] Jonathan Ho, Ajay Jain, and Pieter Abheel. Denoising Diffusion Probabilistic Models, December 2020. URL http: / arwiaaivvvvv../.aaaaaaaaa wttp://....//b/bb/aa http://bb//0000 1 11239. arXiv:2006.1 1239 [cs, stat].
[26] Jae Hyun Lim, Nikola B. Kovachki, Ricardo Baptista, Christopher Beckham, Kamyar Aziz- zadenesheli, Jean Kossaifi, Vikram Voleti, Jiaming Song, Karsten Kreis, Jan Kautz, Christopher Pal, Arash Vahdat, and Anima Anandkumar. Score-based Diffusion Models in Function Space, November 2023. URL http: httt::///xxiv www.:oorg/////eeeeeeoee 07400. 000000002.00 [es, math, stat].

---

[27] Jakiw Pidstrigach, Youssef Marzouk, Sebastian Reich, and Sven Wang. Infinite- Dimensional Diffusion Models, October 2023. URL httphtt://aree org abs/ 2302. 10130. arXiv:2302.1011 10130 [es, math, stat].
[28] Gavin Kerrigan, Justin Ley, and Padhraic Smyth. Diffusion Generative Models in Infinite Dimensions, February 2023. URL http://arziv wwwpeiiiii....vaaaeeee wwwggaa/////aree 00886. arXiv:2212.00886
[cs, stat].
[29] Gavin Kerrigan, Giosue Migliorini, and Padhraic Smyth. Functional Flow Matching, December
2023. URL http://arxiv,. org/ /aa///bbs////.../ 17209. urXiv:2305.17200 [cs, stat].
[30] Andreas Damianou and Neil D Lawrence. Deep gaussian processes. pages 207-215, PMLR, 2013.
[31] Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai. When Gaussian process meets big data: A review of scalable GPs. IEEE transactions on neural networks and learning systems, 31 (11):4405 44233 2020. ISSN 2162-237X. Publisher: IEEE.
[32] Peng Kou, Feng Gao, and Xiaohong Guan. Sparse online warped Gaussian process for wind power probabilistic forecasting. Applied energy, 108:410-428, 2013. ISSN 0306-2619, Pub- lisher: Elsevier.
[33] Juan Maroñas, Oliver Hamelijnck, Jeremias Knoblauch, and Theodoros Damoulas. Trans- forming Gaussian processes with normalizing flows. pages 1081-1089, PMLR, 2021. ISBN 2640-3498.
[34] Vincent Dutordoir, Alan Saul, Zoubin Ghahramani, and Fergus Simpson. Neural Diffusion Processes, June 2023. URL http://arxivv http//aaxxvvv..aaaaaae wwwgg//s////...../////g// 03992. arXiv:2206.03992 [es, stat].
[35] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-66) -6, 2021. ISSN 1533-7928.
[36] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient in metric spaces and in the space of probability measures. Springer Science & Business Media. 2008. ISBN 3-7643-8722-X.
[37] Lénaïc Chizat, Gabriel Peyré, Bernhard Schmitzer, and François- Xavier Vialard. Unbalanced optimal transport: Dynamic and Kantorovich formulations. Journal Functional Analysis, 27401133000-3123, June 2018. ISSN 0022-1236. doi: 10. if.1 201..03.0800.0 URL
[38] Pierre Brémaud. Probability Theory and Stochastic Processes. Universitext. Spriiee Interna- tional Publishing, Cham, 2020. ISBN /ttp:///ce//aa.caaaa//////////20 978- 97833-030-40183-
http:///wwww.cceeceereeee...
970-000001-22---3330-00--1---213 doi: 10.1007/
htthtt////liie springer. www /ttp:////o..1111////-----------0
[39] Andreï Nikolaevich Kolmogorov and Albert T Bharucha-Reid. Foundations of the theory of probability: Second English Edition. Courier Dover Publications, 2018. ISBN 0-486-82159-5.
[40] Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud. FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models, October 2018. URL http://arxiv wwwgg//bs//1a/ 18800. 01367. urXiv:1810.011
[cs, stat].
[41] Michael Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in -Simsiici-Smmaiiionn and Computation, 18(3): 1059- 1076, 1989. ISSN 0361-0918. Publisher: Taylor & Francis.
[42] John Skilling. The cigenvalues of mega-dimensional matrices. Maximum Entropy and Bayesian Methods: Cambridge, England, 1988, pages 455-466, 1989. ISSN 14484. Publisher: Springer.
[43] Vincent Sitzmann, Eric Chan, Richard Tucker, Noah Snavely, and Gordon Wetzstein. Metasdf: Meta-learning signed distance functions. Advances in Neural Information Processing Systems, 33:10136-10144 2020.
12

---

[44] Laura P Swiler, Mamikon Gulian, Ari L Frankel, Cosmin Safta, and John D Jakeman, A survey of constrained Gaussian process regression: Approaches and implementation challenges. Journal of Machine Learning for Modeling and Computing, 1(2), 2020. ISSN 2689-3967.
Publisher: Begel House Inc.
[45] Gavin Kerrigan, Giosue Migliorini, and Padhraic Smyth. Dynamic Conditional Optimal Transport through Simulation Crmue Flows, May 2024. URL http: /arxiv///xivv///aa www.oosg///////o.oooooooo ab20444.
04240. utXiv:2404044200
[46] Yang Jin, Zhicheng Sun, Ningyuan Li, Kun Xu, Kun Xu, Hao Jiang, Nan Zhuang, Quzhe Huang, Yang Song, Yadong Mu, and Zhouchen Lin. Pyramidal Flow Matching for Effi- cient Video Generative Modeling, March 2025. URL http:///://azi http::///aaiii..oaab//////ss arg/abb//4410... 24100.054...... arXiv:2410.0595. [es].
[47] Qihao Liu, Xi Yin, Alan Yuille, Andrew Brown, and Mannat Singh. Flowing from Words to Pixels: A Noise-Free Framework for Cross-Modalls Evolution, March 2025. URL http: //arxiv.org wtgg/bs/////11// 15213. arXiv:2412.15211 [cs].
[48] Cédric Villani. Optimal Transport, volume 338 of Grundlehren der mathematischen Wis- senschaften. Springer, Berlin, Heidelberg, 2009, ISBN 978-3-540-71049-3 978-3-540- 71050-9. doi: 11111007/978--550-7000000 URL http://link. springer www/co.mmm.c../m.c 1007/ 978- 103--540-715000
[49] Vladimir Igorevich Bogachev. Gaussian measures. American Mathematical Soc., 1998. ISBN 0-8218-1054-5. Issue: 62.
[50] Leonid Vasilevich Kantorovich and SG Rubinshtein. On space of totally additive functions. Vestnik of the St. Petersburg University: Mathematics, 13(7):52-59, 1958. ISSN 1063-4541. Publisher: Allerton Press, Inc.
[51] Chunyuan Li, Changyou Chen, David Carlson, and Lawrence Carin. Preconditioned Stochastic Gradient Langevin Dynamics for Deep Neural Networks, December 2015. URL http:// arxiv. Dexiv.org/are orgooog////oo hts/bb//1512 07666. arXiv:1512.07666.
[52] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural Ordi- nary Differential Equations, December 2019. URL http:://arast wwwaaaawiivrr.orraeee orgo/oog////.oooooooo aa//a//s////0// 07366. arXiv:1806.07366 [es, stat].
[53] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP, February 2017. URL http:
brgg/bss/600 0tte//665..0800 08803. utXiv:1605.088033
[54] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to
the Monge-Kantorovve mass transfer problem. Numerische Mathematik, 84(3):375-393, 2000. URL https: iap. wt//at//aaaaaa senttuulliess/ee tetiittteelllaalaeeeeeeee http:///mm//22//raarrrrr00// Transsprttpptimmlla 20 (01.6777_vvc--00 WWW_VVCCXXDDD203-22222 pdf. Publisher: Springer- Verlag Berlin/HHeidelherg.
[55] Emilien Dupont, Yee Whye Teh, and Arnaud Doucet. Generative Models as Distributions of
Functions, February 2022. URL http://arxivv
http//bs//00..... 04776. arXiv:2102.04776.
[56] Yaozhong Shi, Zachary E. Ross, Domniki Asimaki, and Kamyar Azizzadenesheli. Mesh- Informed Neural Operator: A Transformer Generative Approach, June 2025. URL http:
//arxiv. org barg//b///e abs/ 16656. arXiv:2506.1666 [cs].
[57] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps, October 2022. URL attp:///axx: http//awwvv...aaaaaa/aa/aaa/aaat www./org/////////eoo 00927. 2XXv:2206.00922 [cs].
13

---

A
Potential questions and answers
Our method combines many vital fields, including
Operator learning
Flow-based generative models (flow matching); Optimal transport in function space
Bayesian uncertainty; Gradient Markoy chain Monte Carlo, SGLD; Stochastic trace estimation
Gaussian processes: Stochastic processes (Kolmogorov extension theorem)
Because readers will have diverse backgrounds and expertise, we provide some potential questions and answers enhance the clarity and highlight our contributions.
Q1: Why is it hard to obtain the true prior and posterior of unknown stochastic process, and what would we gain if we could?
Al: Generalizing GP regression to arbitrary stochastic process regression has remained a long-standing challenge that has kept many researchers busy for decades. The core difficulties are twofold: first, learning a general stochastic- stochastic-process process prior from historical data with an expressive enough model, and second, deriving both the exact posterior distribution and an efficient sampling scheme.
Prior work falls into three main categories GPs, deep GPs, and conditional models (NP family). Only the classical GP can fully characterize a stochastic process, and then only when the exact GP prior (mean function, covariance kernel, all ppper-parmmetrs, and Abservation-noiis variance) is already known. Even fitting one GP to data generated by another GP with unknown parameters rarely yields the posterior. Deep GPs and NPs possess stronger representational power compared to GP, but representational power * exactness: their optimization rely on approximate posterior and DO NOT produce closed-form "true posterior."
the true prior and posterior were available, we would gain near-complete control over general stochastic process, enabling perfect uncertainty quantification and Bayes-optimal decisions, which has numerous applications finance market, science and engineering problems
Q2: How does OFM differ from standard finite-dimensional flow matching? What does it mean generalize flow matching stochastic processes, and how does extending a generative model to function space differ from extending it to a full stochastic process?
A2: A standard flow matching model learns a transport map between two distributions defined on a fixed grid (e.g.. a pixel lattice). Consequently, it only generate samples that specific resolution. contrast, OFM learns a map between two stochastic processes in function space. This endows with several properties that standard flow matching lacks:
Resolution agnosticism. OFM learns transport map between two distributions defined on any given collection of points, without regard for the number of points, or their locations in the domain. This enables capabilities like zero-shot generation without retraining.
(ii) Stochastic process consistency. OFM respects the metric of the underlying space, ensuring that points close each other in the input domain have appropriately correlated values in the output distribution. This is part of learning valid stochastic process, which also satisfies crucial theoretical properties like the Kolmogorov extension theorem (i.e, consistency under marginalization)
(ii) Convergence to a continuous function. As the collection of query points becomes denser within the domain, the output of OFM converges to the underlying continuous function.
(iv) Backwards compatibility. When evaluated on fixed set of points, OFM recovers the behavior of conventional (inite-dimensional) flow matching up linear transformation.
Extending a generative model (e.g., flow matching) to an infiniteddimensional function space and extending it a full stochastic process are related but distinct goals. The first concerns generating function values at any finite collection of points; Lebesgue measure and Kolmogorov consistency come into play, but the focus remains on finite dimensional marginals. The second addresses the
14

---

distribution of the entire function itself, requiring a probability measure over an space.
Because of this distinction, lifting a model merely to a function space is strictly weaker. For example, GANO [22] extends GANs to function space but not stochastic process, as it cannot describe the joint law of function values at arbitrary query sets. In OFM, we bridge these two viewpoints for flow matching with neural operators: by exploiting the discretization-conveeyy operators and the invertibility of flow matching, we extend the method to both function spaces and stochastic processes. See Appendix C and D for details.
property of neural
Q3: These methods may appear similar to image inpainting restoration with flow- or diffusion-based models; could you elaborate on the differences?
A3: Current inpainting with flow/diffusion models treats the task as an inverse problem in finiteddimensional setting: one learns to reconstruct missing pixels on fixed resolution grid, effectively solving a regression problem at predetermined set of points. By contrast, OFM operates directly on functions and is resolution-agnossicc. Given the same partial observations, our model can predict the entire function any collection of query points, coarse or fine.
Moreover, OFM delivers principled uncertainty: its posterior quantifies the full distribution of possible completions, whereas minee-dimennional inpainting methods typically rely on an ensemble of visually plausible samples whose variability is not a calibrated measure of uncertainty. Finally, OFM remains effective even when observations are extremely sparse. (e.g 0.7% of total observations) a regime in which grid-based inpainting approaches struggle.
Q4: Could you elaborate further on the connections to related work? The approach appears to intersect with several studies, including OPFLOW [1], COT-FM [45], OT-CFM [2), and others.
A4: Because our work combines multiple fields, it naturally has connections with many other studies. Readers are referred to the appendix Q for a detailed discussion.
Q5: In Eq 11. What is the argument of replacing the superma in Eq 8 with expecations? When is this valid?
A5: The supremum in Eq. 8 represents worst-case error over the entire function space, which is computationally intractable to optimize directly. We therefore relax this hard constraint by replacing the supremum with an expectation, as formulated in Eq 11. This is common empirical consideration. Instead of minimizing the worst-case error, our objective becomes minimizing the tractable average-case error across the distribution. Minimizing the error on average provides strong practical incentive for the model to perform well across the entire function space, thereby effectively reducing the worst-case error. Such replacement is always valid under the weaker goal. The validity of this approach as tractable proxy is further confirmed by our empirical results, which show successfully guides the model to learn the intended functional mapping (Appendix M).
Q6:
Why not use log-likelihood as the evaluation metric for non-GP regression tasks when comparing to NP models? And why not include more recent conditional model baselines, such NDP [34]
A6: There are several reasons that likelihood as an evaluation metric is not relevant in our case. First, shown in Eq 41, there an evidence term, which a constant, in the posterior distribution OFM framework. The constant evidence term is intractable but does not contribute to MAP estimation, mean estimation, and posterior sample in general.
Second, even if we can compute the evidence term, we still cannot make the comparison. The graphical model in the NP is such that the conditional model is trained using the MLE principle and the learned likelihood model dependent. The posterior in OFM provides the posterior using the Bayes rule, utilizes a different model, and has different graphical model. These quantities are not directly relevant to be compared.
Third, computing the true posterior for a general stochastic process a known challenging problem. Due to the complexity nature of the problem, there doesn't 't exist well-recognized metric. The
15

---

evaluation of quality of posterior performance requires domain knowledge from experts and varies case by case, which is an exciting research direction.
For baseline models, We already include broad set of SOTA baselines. On the operator-learninn side, the latest OFFLOW [1] (2024) is covered. For deep GPs, we adopt Doubly Stochastic Variational Deep GPs [9] and the Deep Sigma-Point Process [10], both widely recognized and backed by publicly reproducible code.
Within the NP family, ConvCNP [6] remains a standard benchmark, and we include two additional NP variants. Our baselines are limited to models with publicly reproducible implementations. Although we attempted to add the newer NDP [34], we encountered significant reproducibility issues with the authors code-an obstacle reported by others well. Instead, we offer detailed theoretical comparison between OFM and NDP in Appendix Q.
Q7: What do you mean by an "exact posterior" and a practical solution? How can apply the model to my own tasks, and what do need to prepare?
A7: For functional regression, posterior error has two sources for any method: (1) formulation error-thh the gap between the model's 's theoretical posterior and the true one-and (2) approximation error from finite data, limited capacity, and optimization. Deep GPs rely on variational inference, op- timizing an ELBO and thereby introducing a formulation gap: the posterior is only an approximation. NPs make even stronger simplifying assumptions, and doesn't 't provide the true posterior. All models, including OFM, incur approximation error, but the discussion above concerns only formulation error.
A practical solution means expressive backbone that can learn complex mochastic-process prior and a posterior-sampee sampling routine whose runtime and memory footprint remain acceptable for typical users; see Appendix O for quantitative details.
Using the model is straightforward. OFM offers GP-style regression for non-GP tasks: given noisy observations, it returns the posterior at arbitrary query points. The key difference is that a GP prior is fixed by hand-tuned kernel, whereas OFM learns the stochastic-process prior directly from data. For details on prior learning and SGLD-based posterior sampling, see Appendices M and L.
B Background: Flow Matching, Gaussian Measures on Function Spaces, and the Cameron-Martin Theorem
In this section, we provide essential background and high-level intuitive explanations of the topics involved in this paper for readers.
Flow Matching. Flow matching is a state-of-the-art generative paradigm that learns a time-dependent velocity field whose ODE transports samples from simple reference (c.g.. Gaussian) to the target data distribution, yielding an unbiased, simulation-fris free training objective that directly regresses ground-truth velocities along probability paths.
The approach tightly linked to physics via the continuity equation, and when cast as a continuous normalizing now-provides a (deterministically) invertible transformation with tractable likelihoods, while also admitting stochastic variants when desired. In practice it matches diffusion-level quality with faster training and sampling (often few- or even single-step generation) and excellent scalability, which is why has been adopted in challenging domains such as large-scale video generation, in- context image generation, and protein ensemble generation [46, 47]. We therefore use flow matching as the foundation for prior learning, leveraging its simple training objective, physics- grounded structure, and strong empirical performance.
Gaussian measures. A probability measure μ on a separable Hilbert space H. is Gaussian if for any finite collection of vectors (hi...., hn C H. the random vector ((X, ht),.... ..., (X, hn)) has a multivariate Gaussian distribution, where X μ. The family {(X, h) h H) is therefore a Gaussian process indexed by H. with mean h + (m, h) for some m E H. Conversely, if a Gaussian process (1) : T} has sample paths that almost surely belong to a function space H (e.g.. C(T) or (7)), then the law of the random path ( + X (l) is a Gaussian measure on H. In short: a Gaussian measure the path-space law of a GP, and GP is the collection of linear probes of Gaussian measure.
16

---

Cameron-Martin space and theorem. Let be a Gaussian measure on separable Hilbert space H. with mean m H. and covariance operator C' H. - H. (self-adjoint, positive, trace-class). The Cartron-Martin space (the RKHS associated with (4) is
Hy = Range((C/3) C H.
equipped with the inner product (u, u) (C-1/2u, C-1//1)) on Range(CA//)) extended by completion. With this choice, the inclusion Hy + H is continuous. For h E H. write (A) jA(A - h.) for the translate. The Cameron-Mard theorem states:
(i) If h E Hus then jın. is absolutely continuous with respect to μ (in fact, S4h and j are equivalent), with
dim (a) = exp (C- - mm)) C1/2,,,11/22/
C-1/2(x
[MIR).)) for u.a.e. I.
du
where the pairing is understood in the Paley-Wiener sense (which in finite dimensions reduces to (h, (x-
(ii) If h & Hy, then pm 1 μ (mutually singular).
In finite dimensions this reduces to the classical mean-shift formula for ((m, Σ) with C' = Σ, and ||, = -,) Thus, H₂ pinpoints exactly the directions along which a Gaussian measure can be translated while remaining (mutually) absolutely continuous.
C Stochastic process learning
Let F,, P) denote a probability space and let (R4, ((R4)) denote a measurable space where B(R) is the Borel space. Following the standard definition of stochastic processes (Brémaud [38], Chapter 5.1). a stochastic process P on a domain D is a collection of &d- random variables indexed by members of i...
{a(x) : D}
jointly following the probability law P. In the special case of Gaussian processes, c.g., Wiener process, following the Gaussian law for P, for any collection points {71, T2,,.. Fy }. random variables {a(z1), a(z2)..... aa..:..,,(.,,) (Lxa) are jointly Gaussian, resulting in a function to be drawn from a collection of random variables (random
GP. We need to emphasize, (a(x1). a(x2)),....)) a(r,) } vector) equipped with Lebesgue measure, and represents a discretized observation of one continuous function a. In practice, the joint probability distribution the collection of the random variables is unknown priori, and needs to be learned.
In SPL, one way we suggest is to learn an invertible operator T that maps a base stochastic process P to another stochastic process Q that represents the data via discretization convergence theorem (see Appendix D). That is, for any collection of points (21..2.... ... 2n), and for any n, the operator T maps the law on {a(a)) ), a(x2)..... a(xx) }) to (u(x1), u(x2), u(x,)) and vice versa for the inverse T-1, where u(x)) a pointwise evaluation of function data sample, i.e.,
(u(z1). u(x2). ... u(xn)) T (fa(z1). (.)) a(zn)))
Then, the probability of {u(z1), u(z2).... ... u(xn)), at evaluation points (I1, #2, ........ In), for any and collection of points on D is given by,
IP (fu(x11) u(x2). u(xx))) = J7T
(a(=:)..(es)....ate..) P ((s))..,))
where with abuse of notation [z((())) denotes the density u(x) at point æ, same for P(a(x)), and similarly, following the notation Theorem 11.1 Villani [48], JT JT
(a(=))a((e))....aaa... is the abso- lute value of the Jacobian determinant of the map from the random vector {a(z1), a(z2)...... a(x,))
at points (#1, x2, 42,, ...,..,,, ,a,....,,, sn } to the random vector (u(21), u(22),()...))))).)))..)))) u(xx)} via inverse operator T-1 We further show that the pushedforward Q is indeed a valid stochastic process via Kolmogorov SPI., we aim learn neural operator To such that the resulting Q matches the data process under the true T.
Extension Theorem (KET) [39] with a proof provided in Appendix.
17

---

D Model stochastic process with infnitteddimeniionll flow matching via Kolmogorov Extension Theorem
In operator learning, neural operators [3, 14, 15] are typically designed to map input function to an output function. When the input function is provided at a specific discretization (c.g.. set of points with their corresponding values), the model processes this discretized input as a collection of points and their values. Traditionally, in operator learning, this process is seen as an approximation of the operator's s application to the underlying continuous function, where the discretization introduces approximation errors. Thus, the input conceptually still treated as a function.
Moreover, the application of the operator to a collection of points is well-defined, and, by the discretization convergence theorem, the number of points increases, this operation converges to well-defined mapping. In this paper, leveraging these properties, we adopt a different perspective as described in the introduction. We extend neural operators to define explicit maps between collections of points. In this framework, the input is not the abstract function itself but rather a collection of points and their associated values. Importantly, this mapping remains well-defined regardless of the number of points in the collection and, by the discretization convergence theorem, converges to unique mapping as the point collection approaches the underlying continuous function.
Next, we show that given an invertible operator T and a valid stochastic process P whose finite dimensional marginal is IP( (((a(())) a(x2), (((2)) a(xn)) }. there exist valid stochastic process with
finiteddimensional marginal {u(x1), ), u (x22)....... u(
Once again, as defined in Section C
(u(a1), u(x2), u(xn)) = T (Ha(x1), ((2)).... a((.)))))
Then, the probability of {u(z1). u(x2), u(xn)), J. at evaluation points (2.., En). for any n and collection of points D is given by,
IP (fu(11), u(x2), u(z,))) = JT
((f((11), (o(=))) P
a(x2), a(x2)...., a(xx))) (19)
where with abuse of notation P(u(z)) denotes the density of u(x) at point v, same for P(a(z)), and similarly, following the notation in Theorem 11.1 of Villani [48], JT
is the absolute value of the Jacobian determinant of the map from the random vector Jacobian of the map from the random vector {a(x1), 1, a(x2)..... a(xx)) } at points {01, 22. In. to the random vector (u(x1), u(z2)........ u(z,)} via inverse operator T We should notice Eq. represents the changes of variables between two random vectors, with Lebesgue measure involved. The connection between the finiteddimensionl marginal (equipped with Lebesgue measure) and the probability space is described by Kolmogorov Extension theorem (KET) [39], which assures that if all finite-die dimenssoiai distributions (i.e., distributions of function at finite collection of points) are consistent, then stochastic process exists that matches finite-dimensional distributions.
(a(=))..(
measure of a stochastic process in
Formally, according to KET, to establish that a valid stochastic process Q. which has its finite dimensional distributions, is essential demonstrate that such joint distribution satisfies the following two consistency properties:
((u(x1). u(x2),.....) u(x)) u(xx))))
Permutation invariance. For any permutation и of {1. n). the joint distribution should remain invariant when elements of "e,."""", are permuted, such that
P ((u(211), u(z2), .... u(xx))) P ([u(=+())) ((=)))) u(=+(m))))
(20)
Marginal Consistency. This principle specifies that that if a portion the set is marginalized, the marginal distribution will still align with the distribution defined on the original set, such that for m 2 n
P ((u()))..(x2)..... u(,,)))))))))))))))) = - ff P (fu(11),
du( ... du(xm) (21)
The permutation invariance property is naturally upheld when utilizing operator, as there is no inherent order among the elements in the set {21, 22, .., Furthermore, the marginal
18

---

consistency property is also maintained due to the definition of operator T (see Eq. 19). which ensures that I (fu(xı), u(x2)..... 2(x,)}) closed under marginalization. This is be- (((x))) is closed under marginalization, which fully determines cause IP (((())))a(()2),
through the Jacobian. While verifying that Q constitutes a valid induced stochastic process is straightforward given the T. approximating the T with neural operator non-trivial and depends highly on the model used (related to expressiveness), For instance, in Transforming GP [33], the authors employ a marginal normalizing flow, which acts point-wise operator to transform values from a GP to another. Consequently, the induced Jacobian a diagonal invertible neural operator by generalizing RealNVP to function space, which induces a triangular Jacobian matrix. In our work, we extend this framework to more comprehensive case: diffeomorphism. Here, the induced Jacobian is full-rank matrix and is not necessarily triangular or diagonal, the determinant of the Jacobian for any collection of points is calculated through Eq 15.
((u(x1)) u(x2))..... u(x,))) (((,),,,,,,,(,,))
matrix. More recently, OpFlow [1] introduces
Last, we want to clarify the the connection between the notions of operator T and operator g throughout this paper. The operator T is the & (a diffeomporhismm defined in Eq. 5, which is the integral of G over time interval [0, 1]. Due to the nature of ODE system, T is invertible. However, G is not necessary invertible, which enables us to parameterize it with a classical neural operator, like FNO [14].
E
Universal Functional Regression
UFR is concerned with Bayesian regression function spaces [1], where it can be used infer the posterior of an unknown function on domain D from a collection of pointwise observations. The observations are often corrupted with noise of variance o², denoted as {ū(x, ), ū(x2)..... ū(x,)) or (ū(x))]==- More specifically, for m 2 points at which the function is to be inferred,
IP ((u((1)) u(x2), ..... u((r)) (ū(x1), ((x2), ..... ((x,))
Note that when the prior over the function space is Gaussian, UFR reduces the celebrated GP regression. Following Bayes rule, and maps between stochastic processes, we obtain the log posterior as follows,
((x;;) - u(xx))² n log(o) log log(2r) (2x)
log P. ((u(z))) (T(;))]²
g²
+ log (fu(z;)) log (fū(x:))
This equality holds for any collection of points. It is worth noting that the posterior is exact up to constants, i.e. the second, third, and last terms are constant. Therefore, they do not contribute in MAP estimation, mean estimation, and functional regression in general, and there is no need to compute them.
F Marginal (dynamic) optimal-transport flow matching in function space via optimal coupling and dynamic Kantorovich formulation
Consider joint probability measure ( v1) on H. x H. where the reference measure vp, is chosen as a Gaussian measure, whose absolute continuity well-studied [49]. We characterize vo by GP with trace-class covariance operator. c.g. [mmo Co), where mo is the mean, Co is the covariance operator. With the joint measure T M1), we sample a function pair = := (ho./1).
Assuming P1 has full support on the Cameron-Martin space associated with (following the convention of literatures [26, 29]), we construct a conditional probability measure (=(+12) as Gaussian measure with trace-class covariance operator and small operator norm to approximate Dirac measures in the sense of weak convergence. Such that, at 0 and 1, (-1=)) is centered around ho, h1, approximating ōng 5n, respectively: Subsequently, we can construct new marginal probability measure by mixing these approximated Dirac measures:
(A)) [A(A))))T(2), VA E B(H)
(22)
19

---

Due to dT((2) being always positive, the conditional probability measure (Dirac measure approximated by Gaussian measure) is absolutely continuous with respect to p.. Eq. 22 indicates that μo I 6ng (T(()) s vo, and μ41 = 0n, (2) ~ V1- This formulation suggests that μo, /4 represent convolutions of vo, Vi with Gaussian measures. For more detailed discussion on convolution with Gaussian measures, we refer the readers to Appendix B.1 of Lim et al. [26].
Please note, dirac measure is not a necessary condition for Eq. 22; the only constraint is the boundary conditions. One viable choice for the conditional measure is a Gaussian measure with a small operator norm, which (in fact) approximates the Dirac measure in the weak convergence sense. Alternative probability path, as described in Lipman et al. [11], Albergo and Sanden-Eijnden [12], Liu et al. [13] are equally valid.
Suppose Jo' Ix INxMell ( (hl=)) d) is finite to guarantee the vector field sufficiently regular (Lipschitz continuity),, where (-1+)) is the conditional vector field. Under this condition, the vector field that generates fe specified in Eq. 22 and Eq. 3 can be expanded follows:
Ge(h)= = f HXH
(-2( (h)dT(2)
(c(bl2)
(23)
duz
Eq. 23 is an extension of the Theorem as detailed in Kerrigan et al. [29], and provide the derivation in Appendix G. We note that (+(()2) is a Gaussian measure and can be expressed as (1(+1z) = W( Ct), with mean my and trace-class covariance operator C.. Inspired by Tong et al.
[2]. we choose m, and C2 to have the following forms:
my = t. h1 +
+(1 (1-1 10-t)) t) ho
(24)
C = Co
(25)
where Co is the same Gaussian covariance operator defined for Pn and m is a small constant. Further, similar to finiteddimensional flow matching, we only consider the simplest vector field that applies canonical transformation for Gaussian measures, such that the flow has the form: & (2) my + minhho x t h1 + (1 t) ho. From Eq. 1, we can get Gc(h)z) indicating Ge (h 22) is independent of the time and the path from hg to h1 is a direct, straight line. Equipped with well-constructed conditional vector field and probability measures, we can train a neural operator Go with the conditional flow matching loss
h1 ho-
CFM 2 (0, (1 (l) Ge (b|=)))
(26)
Next, we explore how to approximate the true optimal transport plan from optimal coupling of the joint measure i (1o. v1). A common way for measuring the distance between two probability measure is 2- Wasserstein distance, which a special case of static Kantorovich formulation [50]. The static Wasserstein distance is defined as follows
Waa (Vo, v1)2 = inf / |/h dal/2dt(ho.. h1)
(27)
TEll *
In the ODE framework, we also care about the dynamic form of the 2-Wasserstein distance to estimate the cost along the transport trajectory, which also is special case of dynamic Kantorovich formulation [37].
W dyn (Mo, 2112 12 inf / S. (1)(())) (p))dt
(28)
Mr.Ge
Within the OFM framework, the marginal probability measure is a sum of Dirac measures as described in Eq. 22, and we selected Po as a Gaussian measure and assumed v has full support on the Cameron- Martin space associated with vo. Furthermore, the cost function of 2- 2-Wasserstein distance is squared L norm, which is continuous by nature. According Theorem 4.3 and Lemma 4.4 of Chizat et al. [37], Wa, = Wayn for our specifically constructed He and G. in the sense of weak convergence. Therefore, to get the dynamic optimal transport plan, we only need find a joint measure x (14, 11) that achieves the infimum in Eq. 27. In practice, we use minibatch approximation of optimal coupling between in and v1. The above approach extends the dynamic (marginal) optimal transport framework of [2] to infinitedimmensional function space. The related work of Kerrigan et al. [45] addresses a similar problem, but from a different perspective. For a detailed comparison, please refer to Appendix Q
20

---

G Derivation of Eq. 23
In this part, show the derivation of Eq. 23, which extends Theorem 1 of Kerrigan al. [29]. The problem setting is given continuity equation and its weak form:
ap(h..t (2 + (G.(h), ōt
O, Vp E Cyl(H [0, 1])
(29)
we want to derive the following form of the conditional vector field under absolute continuity assumption and other mild conditions, where 2 (o,) H x H.
(e(hl2) duu((11)) (h)dT(2)
G.(h) G.(h)= = HHH
(30)
du+
First, So JNH Bep(h.t) it (h)dt -ho JH Be(h.t) due dd)))l))dd())dt: With continuity equation in at
strong form and the fact that (bll2) induces me(h(z) 2) we have:
S. J Deplh,t)
1 / -v (4lh.00)(((ha))))m)())11aaa
dL (h(()))dd)(s))
By the divergence-fone identity:
V (s)(h, (b)())), Vny(h, ()) + (hh t)V G.(l))
Since we choose the smooth test function (ph, t) from Cyl(H [0. 1]) and use the continuity equation in weak form, we assume term sphh, t) Ge (h) disappears under integration. Thus we have
tep(h,,) St due (h)d = J/¹ / (G.(M1=),, Vns(h, (dd
J. I (0.())))), Vnp(h, t))
(p)(1)2) (
dp)(h)
L J (M)= (1((1)2)
Vas(h, (b
(h)
S [ (G(()12))
dpu(+12) (h)dd(2), Vay(h, t)) (h)
duu
[ J L Gi(hlz) dpp(-12) (h)dT(2), Vn (h))
djt
On the other side, from Eq 29, we have
S. Ln tip(h,t)
S.' (1)( Vnp(h, t))d, O, Vip Cyl(H x (o, 1))
Øt dpi(h)dt=
Thus Ge(h) S2 G(hl))) dp((-1=) (h)dT(=)= S Ge(hl2) de.(-1=) (h) (h)dT(>)
d()l))))(u(()))))))))))))))))))) dut
H Derivation of Eq. 14
In this part, we show the detailed derivation of Eq. 14. Flow Matching, the variable z is chosen as a single data point from the coupling π T(ug, un) where ~ dı, and ug ~ = A (0, K (11,,22, ...,......... )n). Considering the class of Gaussian conditional probability paths
(31) With conditional flow dy (ue [2) = ue + m, Specially, we choose m. = fuj + (1 t)uo and o, σ. where > is small constant. From Eq. 1 (or Theorem 3 of Lipman et al. [11]), vector that defines the Gaussian conditional flow is:
pr(ue(2)= = (V(.mmmm()) (z(2)² (f21,22..... x,,))
of (w mi) + m((w1)
(((122)
(32)
o,
21

---

Then we can derive a closed-fone expression for both the conditional probability and corresponding vector field [2] by plug in pr and o, into Eq. 31 and Eq.
Pe(ue 1) (uu + (1 tly 0² 2 ((x1, 32, ................................ . Fn Fn)))
(33)
Gr(u//u)) = 0 + (ui un) u1 ug
(34)
Now, let's check the boundary conditions. At t = 0,
po(u, = (uuluuo 0² K (fx11 02: ,2.,.., 5n a - duo
(35)
At 8 1,
pi (ullz) = (uuluu o (f01, T2, D....., ,nn)) o-th, dus
(36)
From Eq. 9, we have po(un) (=)d= Ju o (uo. u1 dddduudu do and pi (4))
P (z)))) J Ja, n (uo, un )dugduz nı, which show boundary conditions satisfied.
I Proof of Proposition 3.1
Proposition 3.1. Given noisy observations {ū((z;) ME2 the posterior distribution is
log IP (fu(z;))= (((x))) Ei= [((x,))((xz)) - u(2
log P ((u(z;)) +C (37)
20
Where the constant C log (2g(2 log ((((x:)){EL)
Proof. With Bayes rule, we have:
PP ([(();)) (a(:)))))) P { fu(x.)]=)) ((
(38)
P((ū(x:)]
Taking the logarithm Eq. 38, we have:
log P [[]]))))]((((()))))
loog ((o((( (u(x:))== + log ((u(x:))
log (fū(I;) )=))) (39)
Given e, ~N 02) and {ei)[_ is multivariate Gaussian, then (w(x)) (x(x;)]] is shifted
multivariate Gaussian with mean {u(x:));= translated from the original multivariate Gaussian {er) Due to the translation invariance property of Gaussian distribution, We have
log P (o(x;)) (u(aa)]2a) = log P ((t)) = En ((x,)) u(xx))1
n
log(2TTO)
20
2
(40)
We notice vn and (ū(x;) ) only depends on {u(x:)) and doesn it depend on {u(x;) Thus log IP ((ū(z)) [u (z(z))) log IP (u(=))))))
For evaluating log PP (fu(x;)) m which is the second part on the right-hand side of Eq. 39, we can efficiently calculate with the trace estimator. The third part on the right hand side of Bq. 39
(log IP (fā(x;))) } (,,,)) )) represents the evidence and is constant. Thus the posterior distribution of Eq 39 can be simplified as:
((z;)- u(2;))l
log P ([(
log ((u(x:)) +C (41)
20²
Where the constant C = - logg(mo²) - log ((ū(x:))
22

---

Example of Posterior Samples
In this section, we initially present the regression result of OFM in another additional N-S scenario, as illustrated in Fig 7. Subsequently, we display more posterior samples used in the 2D regression examples. As depicted
Fig 10, 11, OFM successfully generates realistic posterior samples that are consistent with the ground truth and demonstrate appropriate variability. In contrast, GP regression fails to produce explainable posterior samples.
(a) Observations
(b) OFM mean
(c) True mean
(d) Misfit of mean
(e) Misfit std
Ground touth
OPN
26
Samples from OFM
(g) Sample from regression
Figure 6: OFM regression 2D GP data with resolution 32x32. (a) 32 random observations. (b) Predicted mean from OFM. (c) Ground truth mean from GP regression. (d) Misfit of the predicted mean. (e) Misfit of predicted standard deviation. (f) Predicted samples from OFM. (g) Predicted samples from GP regression.
u
sa
(a) Observations
(b) Ground truth
(c) OFM
(d) OFM sample
GP sample
Figure 7: OFM regression on Navier-Stokes functional data with resolution 64 x 64. (a) 32 random observations. (b) Ground truth sample (c) Predicted mean from OFM. (d) One posterior sample from OFM. (e) One posterior sample from best fitted GP.
Observations
Ground truth
OFM OFM sample
(c) GP sample
(f) ConyCNP mean
(g) ConvCNP sample
Figure 8: OFM regression on MNIST-SDE with resolution x 64. (a) 64 random observations. (b) Ground truth sample. posterior sample from best fitted GP. (f) Predicted mean from ConyCNP, (g) One posterior sample from ConvCNP.
(c) Predicted mean from OFM. (d) One posterior sample from OFM. (e) One
23

---

GPR
Figure 9: OFM regression on NS data. (left) Posterior samples from OFM. (right) Posterior samples from GP regression.
Eas
825
Gac
320
OFM
5.16
Figure 10: OFM regression black hole data. (left) Posterior samples from OFM. (right) Posterior san
GPR
Figure OFM regression on MNIST-SDF data. (left) Posterior samples from OFM. (right) Posterior samples from GP regression.
K Co-domain functional regression with OFM
In this section, we expand our regression framework to accommodate co-domain settings, as many function datasets feature co-domain dimension greater than one. For example, earthquake waveform data commonly include three directional components, leading to a three-dimensional co-domain. Similarly, the velocity field fluid dynamics usually features three directional components, also resulting in dimension of co-domain of three.
We illustrate this extension through 2D GP example with co-domain of 3 (channel dimension of 3). In learning the prior, we define the reference measure (V)) as a joint measure (Wiener measure) of three identical but independent Gaussian measures while the target measure (v1) is another Wiener measure. We keep all other parameters unchanged as those described in the 2D GP regression tasks, with the only modification being increase in the channel dimension from training the prior (training detail provided in Appendix M), and provided 32 random observations across the three channels co-locations, then perform regression with OFM across these channels jointly, As demonstrated in Fig 12, OFM accurately estimate the mean and uncertainty across three channels.
to three. After
L Posterior sampling with Stochastic Gradient Langevin Dynamics
In this section, we describe how to sample from posterior distribution with SGLD. We denote logarithmic posterior distribution (Eq, 41) as log Pa and denote set of posterior samples as {ub] where each un is defined collection of point {z,) 1m
By following the standard SGLD pipeline as described by Welling and Teh [4], we can obtain a of N posterior samples {ug 1 However, SGLD is known to be sensitive to the choice of regression parameters and can become trapped in local minima, leading to convergence issues, especially in regions of high curvature [51]. To mitigate these challenges, Shi et al. [1] proposed that within an invertible framework, drawing a posterior sample vi is equivalent to drawing sample a in Gaussian space, since u/ uniquely defines af and vice versa. This approach can stabilize the posterior sampling process and less sensitive the regression parameters due to the inherent smoothness of the Gaussian process. Additionally, Shi et al. [1] suggests starting from maximum posteriori

---

charnd
Charne
(a) Partial observations
chatoel
Caaees
Ownee
(b) OFM mean
(c) True mean
duanntt
dhanel
dhamne
TAG Nocati oreer
(d) Misfit mean
(e) Misfit standard deviation
Figure 12: OFM regression on co-domain GP data with resolution 32x32. (a) 32 random observations at co-locations, (b) Predicted mean from OFM. (c) Ground truth mean from GP regression. (d) Misfit of the predicted mean. (e) Misfit of predicted standard deviation.
(MAP) estimate of a/. denoted ag, which reduces the number of burn-in terations needed in SGLD. We adopt the same sampling strategy and the algorithm is reported in Algorithm
When the size of observations context points (fū(x;) H=))) is O, sampling from the posterior degrades to sampling from the prior, the results which are presented in the subsequent section.
Algorithm Posterior sampling with SGLD
Input and Parameters: Logarithmic posterior distribution log Po. temperature learning rate Ne- MAP ag, burn-in iteration b, sampling iteration łN, total iteration N.
1: Initialization: ag Te
2: for = 0, 2,....... NN do
3: Compute gradient of the posterior: V log
4: Update dg d++ 4+ ab log + Vn,T(((0)))
5: if 2 b then
6: Every tN iterations: obtain sample ag ++1 and corresponding ug 1+1
7: end if
8: end for
M Prior learning with OFM
We now elaborate on the prior learning process and the corresponding performance evaluation. As shown in Algorithm 2, the training dataset is sampled from the unknown data measure 1/- Concretely, the training dataset consists of M discretized functions D. observation of the w; function.
J. s=1 where D; denotes discretized
In practice, simplify dataset preparation, one often uses the same discretization grid D, for all function samples, c.g. D, = {I1.... xn regardless of the sample index "i". For the consistency of notions, let bo represents batch of i.i.d discretized functions sampled from the training dataset (cquivalently, sampled from v1). Next, the reference Gaussian process vo A( Co) is known and determined by the user. With slight abuse of notation, We choose to use notation hop h for consistency purpose, in other parts of this paper, discretized hyg hı is replaced with a respectively.
For specific experiments setting, we employ Matern kernel to construct the reference GP and to prepare training datasets for GP, 2D GP, and 1D TGP. We have set kernel length 0.01 with smoothness factor Ç 0.5 for all reference GPs. OFM maps the GP samples from reference GPs to data samples and is iooouution-ivvaiian any resolution and evaluated at any resolution.
which means OFM can trained with functions
25

---

Algorithm 2 Learning a prior
Input: Reference Gaussian process vo Nmmo, Co), data measure Vi, batch size b, small constant min- discretized domain D {æ1,-" "..."-- En
1: while Training do
2: hg ~ Vo; h1 ~ 1
sample functions of size iid from the measures on D # mini-batch optimal transport plan
3:
T OT(ho,h1)
4:
(ho, h1) ~
5:
to 10,,0)))
6:
He - thi + 11-ttho t) ho
7:
h NYu, o Co)
8:
(0) - xt) (h1 - ho)
9:
0 Update 0, V (0))
10: end while
return Go
ID GP dataset. We choose / 0.3 and Ç = 1.5 and generate 20, 000 training samples on domain [0, with fixed resolution of 256. We use autocovariance and histogram of point-wiss wiie value as metrics for evaluation. We evaluate OFM at several different resolutions shown Fig 13, 14, 15, which demonstrate OFM's excellent capability to learn the function prior.
Ground Pluhh (1editton-120)
Ooemster Flow Natcring 107N!
Actocovarane eice
Illeegram
GIM
30
Figure 13: OFM for ID GP prior learning, evaluated at resolution=128. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
Puth (rediution... 206)
Opemtor Matcring 107M)
Autnecondaree
letog
022
100
20
105
200
Sunberolaas
Figure 14: OFM for ID GP prior learning, evaluated at resolution=256. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
Grund (redittonn512)
Ooerstr Fhow Natcring 107N)
Autreevadares
1 litog
OPN
Figure 15:
OFM for ID GP prior learning, evaluated at resolution=512. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
ID TGP dataset. We choose / = 0.3 and 1.5 and generating 20, 000 training samples domain [0, 1] with fixed resolution 256. We set [-1.2, 1.2] for the bounds. Results provided in Fig 16, 17, 18.
Naiver-Stokes, Black hole, MNIST-SDF datasets. All the following 2D datasets are defined on domain [0, 1] [0, and have resolution of 64 x 64. We collected 2D Navier-Stokes dataset
26

---

Ground Futn 11000
Qpccidon Hlan Matching SOMN)
Holzgteti
2.5
14
0.55
0.0-
CPN
CE
45
25
Figure 16: OFM for ID TGP prior learning, evaluated at resolution=128. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
Ore and Turn 2555
Operatte FIEs Narering (DON)
Autoesunnianee
Hicveram
DHM
-5
100
Figure 17: OFM for 1D TGP prior learning, evaluated at resolution=256. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
consisting of 20000 samples, with viscosity le 4 The results, including zero-shot super- resolution, are provided Fig 19, 20. The learning Black hole dataset, generated using expensive Monte Carlo method, is detailed in Fig 21, 22. Additionally, trained OFM 20, 000 MNIST-SDF samples, the outcomes are illustrated in Fig 23, 24.
27

---

Couun yro Mirreb13)
Operador Matching sonN!
Aulsscvatiance
Hologtam
0.00
Graan
0.5
0.00
65-
0.00
-15
303
Figure 18: OFM for 1D TGP prior learning, evaluated at resolution=512. (left two) Random samples from ground truth and generated by OFM. (right two) Autocovariance and histogram comparison
Figure 19: OFM for 2D N-S prior learning, evaluated at resolution=64 x 64. (top left) Random sam- ples from ground truth. (top right) Random samples generated by OFM. (bottom) Autocovariance comparison
23-
20
21-
20-
20
12-
OFM
102
330
120
120
124
120
170
Figure 20: OFM for 2D N-S prior learning, evaluated at 128 x 128 resolution (zero-shot super- resolution)
Srald
5032
3030
Figure OFM for 2D black hole prior learning, evaluated resolution=64. (top left) Random sam- ples from ground truth. (top right) Random samples generated by OFM. (bottom) Autocovariance comparison
28

---

25-
28-
26
20-
26-
40
40
40
40
OFM
100
55
Figure 22: OFM for 2D black hole prior learning, evaluated 128 128 resolution (zero-shot super-resolution)
Ground Truth
3
20
28
10
DPH
Figure 23: OFM for 2D MNIST-SDF prior learning, evaluated at 64 64 resolution. (top) Random samples from ground truth. (bottom) Random samples generated by OFM.
Ground Truth
20
3
wx
231-
110
56
119
1/0
40
Figure 24: OFM for 2D MNIST-SDF prior learning, evaluated at 128 x 128 resolution. (top) Random samples from ground truth. (bottom) Random samples generated by OFM.
29

---

N Additional results of 1D GP regression with more complex kernels
In this section, we extend the ID GP experiments by introducing more complex kernels for both data generation and regression, while holding all other settings identical those used with the Matérn kemel. Specifically, we report results for the non-stationary Gibbs kernel and the Rational Quadratic (RQ) kernel.
For the Gibbs kernel, we use an input-dependent length-scale
((x) + 11,,
[0, 1/,
which induces the covariance
k(x,x') = 0² 28(x) ((x')
(x x²)2 1(x)2 + f(1)22 cxp {(x)2 + 1(x1)2
In our setup we take 40 0.05, 8 = 0.25, and 1.0.
For the RQ kernel, we use sklearn, Seeaanngguusee deaussian_proceee kernels RationalQuadratii with length-scale 0.15. As shown in Table 2, OFM consistently outperforms all baselines on these tasks.
Table 2: Comparison of OFM with baseline models on ID GP with Gibbs kernel and RQK. Best performance in bold.
Dataset 7
ID GP-Gibbs
1D GP-RQK
Algorithm & Metric +
SMSE
MSLL
SMSE
MSLL
NP
4.0 10
7.1
5.8 10
2.2 0
ANP
3.5 3.5
6.2-10
2.7
2.1-10
ConvCNP
3.4.
1.1-10-1
2.2-10
6.
DGP
4.3-
8.6-10-1
2.4-10-1
1.2
DSPP
1.2
7.1-10-1
2.5
4.2-10-1
OpFlow
3.1.
60
1.
5.0-10
OFM(Ours)
2.9 1.9-10-1
8.7 10
10
9.4.10
Details of experimental setup
In this section, we outline the details of experiments setup used in this paper. Since regression with OFM requires learning the prior first, we list the parameters used for learning the prior and regression separately. We employ FNO as the backbone, implemented using neuraloperate All time reported in the subsequent tables are based on computations performed using single NVIDIA RTX A6000 (48 GB) graphics card. In all experiments, we use the dopri5 ODE solver provided by torchdiffeq Chen et al. [52] with atol=1e-5 and rtol=1e-5, Detailed posterior sampling algorithm is provided in Appendix L.
library [14].
Table 3 details the parameters used for training the prior. For instance, in the ID GP prior learning experiment, the dataset consists of 20,000 samples, each with a co-domain dimension (or channel) of one. The batch size is set at 1024, and the model is trained over 500 epochs. The total training time is about 0.76 hours, and the size of the trained model 37.1 megabytes.
Tables 4, 5. and detail the parameters for SGLD sampling as described in Algorithm 1. For example, in the ID GP regression experiment, the regression takes 40,000 iterations with burn-in phase of 3,000 iterations. Posterior samples are collected every 10 iterations. The temperature for the injected noise during the gradient update
at 1, and the learning rate decays exponentially from 0.005 to 0.004 (defined in Algorithm 1). We average 32 runs with the Hutchinson trace estimator to evaluate the likelihood, utilizing GPU parallel computing. The noise level, as specified in Equation 41, is 0.01 in this regression task. Then given 6 random observations, we ask for the posterior samples across 128 points. The GPU memory usage for the regression task gigabytes, with the total runtime to 4.91 hours.
30

---

Table Parameters used in experiments of prior learning
Datasets
Size of Dataset
Channels
Batch Size
Epochs
Training Time
Model Size
ID GP
10*
1024
5 10"
0.76
37.1 MB
ID TGP
2 10*
1
1024
5 10
1.24 h
37.1 MB
GP
2 10*
5 1-102
1.14h
MB
2D co-domain GP
2 2.104
3
256
5-102
1.01 h
76 MB
2D N-S 2D Black hole
2 10
1
256
5-10"
3.79
286 MB
1.2 10*
1
256
5-10"
2.28
286 MB
2D MNIST-SDF
2 2-10
1
256
5 10
8.31 h
286 MB
Table 4: Parameters used in regression experiments Part A
Datasets
Total Iteration
Burn-in lteration
Sampling Iterations
Temperature of Noise
ID GP
4 10
3 10"
10
ID TGP
1-10
3 10
10
2D GP
2- 10
3 10
10
2D co-domain GP
2-10"
3 10
10
2D N-S 2D Black hole
2-10"
3 1-10
10
2-
3 10
10
MNIST-SDF
2-
3 10
10
P Ablation and scaling studies for mini-batch optimal transport and Hutchinson trace estimator
We first present an ablation study of the optimal transport plan. In this study, we revisit prior learning on the N-S dataset by training two models: one using independent coupling and the other employing a mini-batch optimal transport plan. Both models were trained with a batch size of 64. For evaluation, we compare the mean squared error (MSE) of density, autocovariance and spectral characteristics between 1,000 real and generated samples. Additionally, we report the convergent training loss (squared 1,2 loss) and the number function evaluations (NFE) required for sampling with adaptive ODE solver (dopri5). As shown in Table 7, the mini-batch OT plan outperforms the independent coupling approach in terms pointwise accuracy. spectral fidelity. and convergent L loss, while also requiring fewer NFE and enabling faster sampling. Our findings indicate that the mini-batch size increases, the model learns the prior with reduced error. This improvement may attributed to the mini-batch optimal transport plan approaching the true optimal transport plan with larger batch sizes.
Next, we explore the variance of the Hutchinson trace estimator in likelihood estimation with a scaling experiment. For this experiment, we use a ID example with parameters described in Section M and resolution of 128. We randomly draw a GP sample from the prior and evaluate the log likelihood by integrating the divergence as shown in Eq. 15, while also estimating it using the Hutchinson trace estimator from Eq. 16. In this scaling experiment, the number of noise sample (Rroise) for the Hutchinson estimator is set to (4, 8, 16, 32, 64, 128). For each Rhooiee we repeat the experiment 100 times and report the mean and standard deviation of the predicted likelihood. The exact likelihood is computed by directly integrating the trace of the Jacobian. We repeat this procedure for different random 1D GP samples. As reported in Table 8, the standard deviation of the Hutchinson trace estimator decreases rapidly as Nuoise increases, and the predicted means always align closely with the ground truth.
Furthermore, even smaller Thnoise values, where relative larger variance is expected, the perfor- mance of the posterior sampling appears robust (see Table 5). We hypothesize that this is due to stochastic nature of posterior sampling algorithm (SGLD), which requires injected perturbations, renders its performance relatively insensitive to the choice of nooise
31

---

Table 5: Parameters used in regression experiments Part B
Datasets
Initial Learning Rate
End Learning Rate
Hutchinson Samples
Noise Level
ID GP
5 10
4- 10
1 10
ID TGP
5 10-3
4 10
32
1 10-3
2D GP 2D co-domain GP
1 1.10-a
8.10
32
1 1.10-2
1-10-3
8-
1-10-2
2D N-S 2D Black hole
3 3-10
2-10-3
1.10-3
5-10
4.10
1.10-3
2D MNIST-SDF
5 30
4 10-3
8
1 1-10-3
Table 6: Parameters used in regression experiment Part C
Datasets ID GP
Number Observations
Inquired Grids
GPU Memory
Running Time
128
GB
4.91 h
ID TGP
128
4 GB
5.42 h
2D GP 2D co-domain GP
32
32 32
22 GB
9.70
32
5.05
2D N-S 2D Black hole 2D MNIST-SDF
64 X 64
44 GB
13.65 h
32
64 X 64
44 GB
13.37 h
64
64 x 64
44 GB
9.41 h
Table 7: Scaling study for the size of mini-batch given optimal transport plan, best performance in bold
Metries
Density-MSA
Autocovariance-MSS.
Spectra-Mat MSE
Convergent L" Loss
NFE
Independent
3.4 10-5
7.4 10
1.3 10'
5.1 10
168
mini-batch 32
3.0 10-5
1.3 10
5.0 10
2.3 10
141
mini-batck
9.8 10
908-10-5
6.3 10
1.9 10
153
mini-bateds 128
2.4 10-5
1.7 10-5
3.6 10
1.5 10
182
Table 8: Scaling study for the number of noise samples of Hutchinson trace estimator
Phoise
Nhxose
Raoise
Ranise
32
Raoise
Minise
128
exact
sample
-481.7 9.4
-482.9
-480.6 + 4.9
-481.4 3.8
-481.1 = 2.5
-480.8
1.8
-82.77
sample 2
-482.8
9.8
-481.8
-481.3
4.9
-480.8
-481.1 - 2.7
-481.2
1.7
482.0
sample 3
-17..7
10.2
-479.7 + 7.6
-178.2 4.7
-478.6
-478.8 + 2.6
-178.8
= 1.7
179.
sample 4
-476.9
10.0
-477.0 6.6
-477.9
-478.5
-478.1
+ 2.7
-478.0
1.8
-177.4
sample 5
-47..4 10.7
-479.2
4479.2
-479.4
-479.6 + 2.8
-479.3
2.0
-478.5
Likelihoad estimation
sample
Likelihood estimation
sample
-44
Exact Bellhood
Exact Roolhooo
-46
Hutchinson
Log Iikelhhood
-460
likcelhhood
-470
-470
-480
-480
-49
-490
-500
126
Number of Runs
Number Runs
Figure 25: Log likelihood by the integration of divergence, (left) plot for first GP sample. (right) plot for second GP sample
32

---

Detailed analysis of OFM and comparison with existing methods
In this section, we elaborate the connection and difference with pervious work, highlight contributions and potential limitations of our work. The regression with OFM involves two-steps process: (i) learning a prior on function space, and (ii) sampling from the posterior given observations. Consequently, the OFM framework has connections with both generative models on function space and the models developed for functional regression. In the following, we provide a comprehensive comparative analysis with related models and baselines, including operator flow (OPFLOW) [1], conditional optimal transport flow matching (COT-FM) [45], conditional models (NPs) [7. 34]
Comparison with OPFLOW. OPFLOW introduces invertible neural operators, which generalizes RealNVP [53] to function space and maps any collection of points sampled from a GP new collection of points in the data space, using the maximum likelihood principle [1]. This method captures the likelihood of any collection of point consistently as the resolution increases and allows for UFR using SGLD. Despite these advantages, the requirement for an invertible neural operator brings training and expressiveness challenges. To be specific, OPFLOw failed on all non-GP regression tasks in this paper because the prior learning stage suffered from mode collapse during training. On the contrary, OFM adopts a simulation free ODE framework for prior learning, which offers enhanced expressiveness and ensures training stability through a simple regression objective while avoiding using the invertible neural operator. In addition, OFM proposes non-trivial extension of UFR to the simulation free ODE framework. These improvements render OFM a more practical solution for challenging functional regression tasks.
Comparison with COT-FM. COT-FM [45] proposes a conditional generalization of Benamou- Brenier Theorem [54], formulating a conditional optimal transport plan that applicable for both Euclidean and Hilbert space. In contrast, OFM employs an unconditional optimal transport plan in Hilbert space based on dynamic Kantorovich formulation, which initially generalized for unbalanced optimal transport [37]. The advantage of COT-FM lies in its ability to flexibly incorporate specific conditions tailored for conditional generative tasks. However, COT-FM is not suitable for functional regression tasks due to: (i) COT-FM is contingent upon both the reference and target being influenced by conditions, and the vector field learnt is triangular, designed to transport jointly the coupling of reference measure and condition measure. In UFR setting, the learnt prior is required to be unconditioned, (ii) the coupling with condition measure typically prevents inducing valid stochastic process, even when the reference measure is a Gaussian measure, (iii) cannot provide point evaluation of probability density. Last, We should notice, the development of OFM is different and independent of COT-FM, the former with focus stochastic process learning and Bayesian functional regression.
Comparison with conditional models. NPs were developed to address the computational and restrictive prior challenges of Gaussian Processes, utilizing neural networks for efficiency [7]. However, several recent studies have discussed the drawbacks in the formulation of NPs, raising concerns that NPs might not learn the underlying function distribution [1, 55].
Notably, NPs treats the point cloud data as a set of values, ignoring the metric space of the data [55]. This can lead to misinterpretations of function sampled at different resolutions as distinct functions (Appendix A.1 of [22]). Furthermore, NPs rely on encoding input data into knite-dimensionas douusian-distribien
latent variables before projecting these into an diiniiteddmeeiionall space. This process tends to lose consistency at higher resolutions. Moreover, the Bayesian regression framework underpinning NPs focuses on point sets rather than the functions themselves, leading to dilution of prior information with increasing data points.
In recent study, diffusion-basse variants of NPs (NDP) [34], was proposed leverage the expres- siveness of diffusion models [24, 25J. Nonetheless, the formulation of NDP does not address the aforementioned issues of NPs and introduces two more problems: (i) NDP fails to induce valid stochastic process as it does not satisfy the marginal consistency criterion required by Kolmogorov Extension Theorem [39], and (ii) it relies on uncorrelated Gaussian noise for denoising, which is not applicable in function spaces [26]. Oppositely, OFM establishes more theoretically sound framework by rigorously defining learning within function spaces. Additionally, Bayesian functional regression within the OFM framework adheres to valid stochastic processes, offering a robust and theoretically grounded solution. Last, we provide a high-level comparison of OFM and NP as shown in Table 9.
33

---

Table 9: High-level comparison of OFM and NP.
Aspect
OFM
Neural Processes (NP)
Modeling target
Global prior stochastic processes with valid joint density,
Amortized conditional model without explicit tractable joint density over full functions.
Arbitrary queries
Set-size location agnostic; evaluate on any grid point set.
Set-size location agnostic; trained with random context/target splits. Performance may degrade under extreme sparsity or distribution shift Predictive densitics trained via conditional likelihood; closed-fon foredd-fomm function-level likelihood is typically specified. Uncertainty reflects model/context coverage and may be miscalibrated.
Uncertainty & likelihoods
Tractable log-densities and posterior sampling enable principled Bayesian inference.
available and
Structure extrapolation
Captures global geometry and non-stationarity via the learned flow.
Learns inductive biases from through context summaries; strong interpolation, but extrapolation and long range generalization not guaranteed depend on training distribution.
Limitations. Despite these advances, the current regression framework with OFM is primarily limited low-dimensional data (ID and 2D this study). This limitation stems from the challenges associated with learning operators for functions defined on high-dimensional domains-an area that remains underdeveloped both computationally and in terms of dataset availability [15]. Additionally, while the time complexity for regression with OFM is Omm)), the incorporation of additional components significantly increases its computational resource requirements compared to classical GP regression.
There are several potential paths to mitigate the computation concerns: (i) Carefully adjusting the step size (1). and temperature (T') in the SGLD algorithm to prevent overly large gradient updates. (ii) Running the trace estimator on multiple GPUs, use mixed-precision FNO [43], or adopt efficient neural operators [56]. (iii) Using more efficient ODE solvers. Our current implementation uses high-order ODE solver (dopri5) that requires over 100 steps for sampling. Switching to a more modern and efficient solver, such as the DPM-Solver [57], could potentially reduce the number of evaluations, freeing up significant GPU resources.
34

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction succinetly outline the paper's claims and contri- butions, which are fully supported by the theoretical analysis and experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We also highlight the limitations in the main text and add explain into detail in appendix Q
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiiatinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
35

---

Answer: [Yes]
Justification: We provide detailed explanations, derivations, and proofs for all theoretical results (see Appendix C, D, and F).
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide the detailed experimental setup in Appendix M. L, and O. We also include code with additional supplementary examples, log files, and figures to ensure that all results are fully reproducible.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
36

---

reproducible, they 5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification: The code is provided as supplementary material, along with example scripts, log files, and figures to ensure full reproducibility of our results.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips) cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: SoodeSubmissionPe.
/ /nips. cec/public/guu
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: See Appendix
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Following the convention in functional regression, we assess posterior learning using SMLL and MSME for GP examples. For non-GP tasks, we also provide visual checks by plotting the predicted mean alongside posterior samples.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
37

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CL, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: See the Appendix O
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics /ttps://neurils
Answer:
[Yes]
Justification: The authors do foresee any ethical implications this work.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No. they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work is not directly tied to any harmful applications.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
38

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This is not applicable
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We cite each reused code snippet and dataset with its original publication.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcode. thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
39

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer:
[NA]
Justification:
We do release new assets
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
We do not do crowdsourcing research with human subjects
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
We do do crowdsourcing nor research with human subjects
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
40

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: This work is not related to LLM
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.LM policy (https: weeeurpps.cc/oooeeeeeeeeeececeee for what should or should not be described.
41


Reviews:
[{"review_id": "Bra2E8Epnk", "reviewer": "Reviewer_2aqv", "summary": "The paper proposes Operator Flow Matching (OFM), an extension of (optimal transport) conditional flow matching (CFM) to stochastic processes using neural operators. In a nutshell, the authors show that the vector fields in the CFM objective can be replaced with neural operators, thus generalising CFM to (infinite-dimensional) function spaces. This allows for evaluating the posterior probability of observations under any stochastic process as well as efficiently sampling from the process.\n\nThe authors evaluate their method on a range of experiments from simple synthetic 1D Gaussian and non-Gaussian data, where the true posterior is available, to more advanced data from simulations and MNIST. It shows promising performance, although at a high computational budget.", "questions": "**Suggestions**\n1. My main suggestion would be to work on the clarity of the paper. For instance, it would have helped me to have a basic introduction to neural operators and (conditional) flow matching in the appendix. Also, I would suggest moving the technical details of the main paper to the appendix, and focusing the main paper on building intuition for your model and how to use it, perhaps with some proof sketches if necessary. \n2. The figures and tables in the experiments section are placed rather chaotically, and table 1 is never mentioned, which is all somewhat distracting. It would help the reader if figures and tables were placed on the same page as where they are mentioned.\n\n**Questions**\n1. You mention that you extend neural operators to map between collections of points. Can you clarify how this is different from a standard flow?\n2. In Eq. (11), what is the argument for replacing the suprema in Eq. (8) with expectations? When is this valid?\n3. The OFM posterior samples in figures 8, 10, and 11 look quite noisy. Do you have ideas as to why this could be and how one might mitigate it?\n4. For the GP experiments, the training set seems to be constructed from samples from the same GP (i.e., with fixed hyperparameters) that was also used for the test sample. This seems like an unrealistic setting in general. Is OFM able to learn from samples from GPs with different hyperparameters? If so, did you try that?\n5. A minor question regarding the 1D GP experiment to make sure I understand the setup correctly: you write in appendix K that you sample the training set at 256 input locations, but table 3 says that you use a batch size of 1024 during training? Is this a mistake?", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n1. The paper combines two recent and promising techniques, flow matching and neural operators. It appears original and should be of broad interest to the NeurIPS community.\n2. The model itself is rigorously developed.\n3. The model appears to work well in practice, despite some computational concerns (see weaknesses).\n4. I could see this work forming the basis for many follow-up papers, so it could become quite significant.\n\n\n**Weaknesses**\n1. The main weakness of the paper, to me, is its clarity. It is very advanced material, which the authors are aware of and have tried to address by adding an FAQ and additional details to the appendix. While these are helpful, the main paper is still very detailed, and the sentences are long and dense, which makes the work difficult to understand.\n2. As the authors themselves point out, the model is computationally heavy and is currently limited to low-dimensional data (the authors experiment on just 1D and 2D data). \n3. The experimental section is a little weak, with many of the results deferred to the appendix. Some of the results are purely qualitative. The only quantitative experiments are some simple 1D and 2D synthetic datasets, which should at least have been extended to more complex examples (e.g., more complex GP kernels would be easy to test).", "was_revised": true, "final_justification": "The authors have addressed my questions and concerns, and assuming the paper is updated to improve its clarity, I think it deserves to be accepted.", "rebuttal": "We appreciate the reviewer’s valuable comments and positive take on our work. We provide our feedback as follows\n\n**Weakness 1 & Suggestion 1. The presentation is dense and contains very advanced material. The suggestion is to add a basic introduction to neural operators and flow matching in the appendix, focusing the main paper on building intuition for the model and its usage.**\n\nThank you for these excellent suggestions. Indeed, this paper presents a rigorous framework that combines several rather diverse topics to tackle the important problem of learning stochastic processes. As Reviewer sYHe also noted, the main text length limit of the paper has made it particularly challenging to communicate the relevant details from each of these fields. Nonetheless, we have tried our best to address this issue as follows.\n\nTo further improve clarity for readers from diverse backgrounds, we have (i) added an appendix section covering neural operators (operator learning), flow matching, Gaussian measures on function spaces, and the Cameron–Martin theorem to provide readers with essential background; and (ii) included additional explanatory text and proof sketches to better convey intuition.\n\n**Weakness 3. The experiment section is a little weak, add more complex examples (e.g., more complex GP kernels would be easy to test.**\n\nFollowing your suggestion, we have brought the important experimental findings from the appendix into the main text, and have subsequently expanded our 1D Gaussian Process (GP) experiments with more complex kernels for generating training datasets and regression tasks. Specifically, we now include results for the Gibbs (non-stationary) kernel and Rational Quadratic kernel (RQK) , As shown in the updated table below, our OFM model consistently outperforms all baselines on these tasks.\n\n| **Algorithm ↓ / Metric →** | **GP‑Gibbs&nbsp;SMSE** | **GP‑Gibbs&nbsp;MSLL** | **GP‑RQK&nbsp;SMSE** | **GP‑RQK&nbsp;MSLL** |\n|---|---|---|---|---|\n| NP | $4.0 \\cdot 10^{-1}$ | $7.1 \\cdot 10^{-1}$ | $5.8 \\cdot 10^{-1}$ | $2.2 \\cdot 10^{0}$ |\n| ANP | $3.5 \\cdot 10^{-1}$ | $6.2 \\cdot 10^{-1}$ | $2.7 \\cdot 10^{-1}$ | $2.1 \\cdot 10^{0}$ |\n| ConvCNP | $3.4 \\cdot 10^{-1}$ | $1.1 \\cdot 10^{-1}$ | $2.2 \\cdot 10^{-1}$ | $6.9 \\cdot 10^{-1}$ |\n| DGP | $4.3 \\cdot 10^{-1}$ | $8.6 \\cdot 10^{-1}$ | $2.4 \\cdot 10^{-1}$ | $1.2 \\cdot 10^{0}$ |\n| DSPP | $4.2 \\cdot 10^{-1}$ | $7.1 \\cdot 10^{-1}$ | $2.5 \\cdot 10^{-1}$ | $4.2 \\cdot 10^{-1}$ |\n| OpFlow | $3.1 \\cdot 10^{-1}$ | $6.9 \\cdot 10^{-1}$ | $2.9 \\cdot 10^{-1}$ | $5.0 \\cdot 10^{-1}$ |\n| **OFM&nbsp;(Ours)** | $\\mathbf{2.9 \\cdot 10^{-1}}$ | $\\mathbf{8.7 \\cdot 10^{-2}}$ | $\\mathbf{2.1 \\cdot 10^{-1}}$ | $\\mathbf{9.4 \\cdot 10^{-2}}$ |\n\n\n\n**Suggestions 2. The figures and tables in the experiments section are placed rather chaotically, and table 1 is never mentioned, which is all somewhat distracting. It would help the reader if figures and tables were placed on the same page as where they are mentioned.**\n\nThanks for pointing out ways to improve the flow of the paper. To address this, we have moved Table 1 to the appendix and have carefully reordered and resized several figures to place them closer to where they are discussed in the text. We also respectfully note that Table 1 was in fact referenced in Section 3.1 (Line 190).\n\n**Q1. You mention that you extend neural operators to map between collections of points. Can you clarify how this is different from a standard flow?**\n\nThank you for this question. A standard flow‑matching model learns a transport map between two distributions defined on a fixed grid (e.g., a pixel lattice). Consequently, it can only generate samples at that specific resolution\n\nIn contrast, OFM learns a map between two stochastic processes in function space. This endows it with several properties that standard flow matching lacks:\n\n(i) Resolution agnosticism. OFM learns a transport map between two distributions defined on any given collection of points, without regard for the number of points, or their locations in the domain. This enables capabilities like zero-shot generation without retraining.\n\n(ii) Stochastic process consistency: OFM respects the metric of the underlying space, ensuring that points close to each other in the input domain have appropriately correlated values in the output distribution. This is part of learning a valid stochastic process, which also satisfies crucial theoretical properties like the Kolmogorov extension theorem (i.e, consistency under marginalization)\n\n(iv) Convergence to a Continuous Function: As the collection of query points becomes denser within the domain, the output of OFM converges to the underlying continuous function. \n\n(v) Backwards Compatibility: When evaluated on a fixed set of points, OFM recovers the behavior of a conventional (finite-dimensional) flow matching up to a linear transformation.\n\n\n**Q2. In Eq (11). What is the argument for replacing the suprema in Eq. (8) with expectations? When is this valid?**\n\nThe supremum in Eq. (8) represents a worst-case error over the entire function space, which is computationally intractable to optimize directly. We therefore relax this hard constraint by replacing the supremum with an expectation, as formulated in Eq. (11). This is a common empirical consideration.  \n\nInstead of minimizing the worst-case error, our objective becomes minimizing the tractable average-case error across the distribution. Minimizing the error on average provides a strong practical incentive for the model to perform well across the entire function space, thereby effectively reducing the worst-case error. Such replacement is always valid under the weaker goal. The validity of this approach as a tractable proxy is further confirmed by our empirical results, which show it successfully guides the model to learn the intended functional mapping (Appendix K). In the revised manuscript, we have also added text to explicitly acknowledge the limitations of this relaxation.\n\n**Q3. OFM posterior samples in figure 8, 10, 11 look noisy, any ideas why this could be and how one might mitigate it**\n\nPosterior sampling in our framework relies on the unbiased Hutchinson trace estimator: we average the score under $n_{noise}$ (Appendix M) Gaussian noise draws to approximate the log‑determinant term. While the estimator is inexpensive in time, it is memory‑bound because all noise realizations must be processed in parallel. Owing to GPU memory limits we used only  $n_{noise}$=8 in Figures 8, 10, 11, which introduces noticeable variance and in turn yields larger stochastic gradients in each SGLD step (Algorithm 1), manifesting as visual noise.. \n\nThere are several clear paths to mitigate this:\n\n (i) Carefully adjusting the step size ($\\eta_t$) and temperature ($T$) in the SGLD algorithm to prevent overly large gradient updates. \n\n(ii) Running the trace estimator across multiple GPUs or using mixed‑precision FNO for memory saving (Tu et al, 2023). \n\n(iii) Using more efficient ODE solvers. Our current implementation uses a high-order ODE solver (dopri5) that requires over 100 steps for sampling. Switching to a more modern and efficient solver, such as the DPM-Solver (Lu et al., 2022), could potentially reduce the number of evaluations, freeing up significant GPU resources that could be allocated to increasing $n_{noise}$. \n\nReference: \n\nTu et al, “Guaranteed Approximation Bounds for Mixed-Precision Neural Operators” 2023\n\nLu et al, “DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps” 2022. \n\n\n**Q4. The GP experiments, the training set seems to be constructed from samples from the sample GP. Is OFM able to learn from samples from GPs with different hyperparameter**\n\nThere may be a slight misunderstanding of our experimental setup, and we appreciate the opportunity to clarify.  OFM is designed to learn the distribution of a single target stochastic process (GP or non-GP). The training set and test examples must therefore come from the same underlying process. \n\nThe model first learns this specific process during the prior learning phase. It can then perform inference for that learned process, but it cannot generalize to a new process without being retrained. \n\n**Q5. Appendix K that you sample the training set at 256 input locations, but table 3 says use a batch size of 1024 during training**\n\nThank you for the question. These two numbers refer to distinct stages of our framework and don't contradict one another. The batch size of 1024 in Table 3 specifies the number of independent function samples used in each training step during the prior learning phase. In contrast, the 256 input locations mentioned in Appendix K refers to the number of points used to discretize a single underlying function during the posterior sampling phase. In short, 1024 is the number of functions in a training batch, while 256 is the resolution at which we evaluate an individual function."}, {"review_id": "q4kP0fT24M", "reviewer": "Reviewer_Ezjy", "summary": "The authors introduce a method to perform flow matching in function spaces relying on a generalization of modelling functions to functionals. They demonstrate the applicability of the proposed approach to academic examples.", "questions": "**Q1: Scalability** The authors mention, their method is limited to low-dimensional settings and, so far, to problems where a simulator is available for unlimited data generation. Can you comment on which theoretical or practical advancements would be needed to overcome these challenges and apply the OFM method to real-world problems?", "limitations": "Limitations are discussed, but should be included in the main text.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strengths\n**S1: Interesting Generalization** The studied problem of generalizing flow matching to functionals is interesting and novel. It addresses a common limitation of CNP that can be of relevance in practice. The derivations seem correct to me.\n\n**S2: Insights into the model** The experiments are described in detail, and the appendix provides comprehensive material on the analysis and comparison to other methods.\n\n## Weaknesses\n\n**W1:  Important insights are pushed to the appendix.** I appreciate the detailed appendix, but limitations or a short introduction on the main neural network baseline, the ConvCNP, should not be pushed to the appendix.\n\n**W2: Excessive Data Requirements** The main strength of using GP and Bayesian methods in general is their performance in low-data regimes. The data requirements seem excessive for the low-dimensional examples provided. I read the limitation section in the appendix, but IMO it should be part of the main paper, as neither the data requirement nor the limitation to low-dimensional settings are mentioned in the main text.", "was_revised": true, "final_justification": "The paper contains a comprehensive introduction and simulation analysis of the proposed method. The method itself is a natural extension of flow matching and seems promising to me, and the authors discuss its limitations appropriately. Although the paper is very dense and reading the long appendix is almost mandatory to fully appreciate the work, it is well written and structured, and I do not see this as a reason to reject (as reviewer KpLg).", "rebuttal": "We appreciate the reviewer’s valuable comments and positive take on our work. We provide our feedback as follows.\n\n**W1. Important insights are pushed to the appendix**. \n\nThanks for the suggestion, we have moved the overview of ConvCNP from the appendix into the main text so that readers can understand its role without consulting supplementary material.\n\n**W2. Excessive Data requirements**\n\nThis is an important observation. We now make the trade-off explicit in the main text by adding:\n “Although our approach attains state-of-the-art accuracy, it relies on substantially larger training datasets than are customary for GP-based models.”\n\n\n**Q1. Issues related to Scalability**\n\nThe primary bottleneck is the general difficulty of learning operators for functions defined on high-dimensional domains (e.g., 3D space plus time). This is a foundational challenge in the operator learning field, which is still evolving to address both computational limitations and dataset availability\n\nApplying OFM to more complex problems will require key advancements on two fronts. First is the continued development of more scalable operator architectures that are both memory- and compute-efficient in higher dimensions. Second, and equally important, is the creation of large-scale, curated datasets for high-dimensional scientific problems, which are essential for training these more advanced operators.\n\nIt's important to emphasize, however, that the current 1D and 2D settings are not merely toy problems; they encompass many critical and valuable real-world applications, such as time-series forecasting, audio processing, and tomography. Our work provides a robust foundation for these domains, while extending the framework to higher dimensions remains a key future goal dependent on these broader advances in the field."}, {"review_id": "D5JGLXGnyb", "reviewer": "Reviewer_KpLg", "summary": "The paper introduces a novel method, Operator Flow Matching (OFM), for learning stochastic processes. The approach consists of two key phases. First, it learns a neural operator that transforms samples from a Gaussian Process (GP) into samples from the target data distribution. Then, this learned operator is used to generate posterior samples given noisy observations. To address scalability and computational efficiency, the authors incorporate Stochastic Gradient Langevin Dynamics (SGLD) into the training procedure.", "questions": "Besides the questions raised in the weaknesses, I also wonder how Eq. 16 is computed in practice and why its complexity is claimed to be $O(m)$. A more detailed explanation would help clarify this point.\n\nAdditionally, the proposed method appears to align well with the denoising literature. It would be interesting to see how the flow matching framework performs on denoising tasks, which could further demonstrate its versatility and practical value.", "limitations": "The authors have discussed the limitations in the Appendix, and the identified issues are common in both the Neural Process (NP)and the operator learning literature.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper is the first to extend flow matching to stochastic process learning via operator learning.\n2. By leveraging a learned operator, the method can directly model data functions, enabling access to the exact posterior distribution over functions.\n3. The approach is validated on several synthetic GP datasets and demonstrates improved performance compared to existing methods.\n\nWeaknesses:\n1. The paper lacks mathematical rigor and clarity in presentation.\na) It is not immediately clear that Eq. 14 is a solution to Eq. 11, nor is this connection explicitly stated. The authors are encouraged to include propositions or theorem-style statements to improve readability. Each phase should clearly state the optimization problem and how it connects to the method.\nb) In line 230, the authors write: \"For evaluating ... on the right-hand side of Eq. 18, we can efficiently calculate it with the likelihood estimation tool described above.\" It is unclear which likelihood estimation tool is being referenced. The statement should be made more precise.\n\n2. The paper lacks a high-level explanation of flow matching, making it difficult to understand the intuition behind why this method might outperform alternatives like Neural Processes (NP). A conceptual overview of how flow matching contributes to the method would improve accessibility.\n\n3. While an ablation study is included, it only focuses on the mini-batch size for optimal transport. However, the method involves other important hyperparameters. A broader sensitivity analysis would improve the robustness of the empirical evaluation.\n\n4. The paper only evaluates on synthetic GP datasets. Real-world datasets commonly used in Neural Process literature—such as MNIST and CelebA for image completion—are missing. Including such benchmarks would strengthen the empirical impact (significance) of the work.\n\nI would be inclined to raise my score if the clarity issues are addressed and additional experiments on real-world datasets and hyperparameter sensitivity are included.", "was_revised": true, "final_justification": "My concerns are resolved. I have read the rebuttal responses and other reviews and made the final decision.", "rebuttal": "We appreciate the reviewer’s thoughtful feedback. Our responses below seek to clarify all concerns and potential misunderstandings\n\n\n**W1. The paper lacks mathematical rigor and clarity in presentation. a) it is not immediately clear that Eq. 14 is a solution to Eq. 11. b) unclear likelihood estimation tool being referenced**\n\nThank you for these suggestions. Indeed, this paper presents a rigorous framework that combines several diverse topics to tackle the important problem of learning stochastic processes. As Reviewer sYHe also noted, the main text length limit of the paper has made it particularly challenging to communicate the relevant details from each of these fields. Nonetheless, we have tried our best to address this issue as follows. To further improve clarity for readers from diverse backgrounds,  we have (i) added an appendix section covering neural operators (operator learning), flow matching, Gaussian measures on function spaces, and the Cameron–Martin theorem to provide readers with essential background; and (ii) included additional explanatory text and theorem-style statements  to better convey intuition.\n\n\nTo address the reviewer’s specific concerns: \n\na)  To clarify, Eq. 11 is the training objective, which requires an explicit form for  $\\mathcal{G}_t(u_t|z)$. Eq. 14 is then used to provide the explicit $\\mathcal{G}_t(u_t|z)$ based on (derived from) the constructed probability path. We provided the full proof and derivation for Eq.14 in Appendix G (as mentioned in line 177, above Eq 14).  To further avoid any confusion, we have added a proposition before Eq. 14 to clarify this exactly in the updated draft.\n\n\nb) The phrase \"the likelihood estimation tool described above\" refers to the Hutchinson trace estimator discussed in Section 3.2. We have revised the text to remove this ambiguity.\n\n\n**W2. The paper lacks a high-level explanation of flow matching** \n\nTo address this, we have added a dedicated section to the appendix providing a background and conceptual overview of flow matching.\n\nIn brief, flow matching is a state-of-the-art generative paradigm that learns to transport samples from a simple distribution to a complex data distribution. Its strong performance and scalability have led to its adoption in highly challenging domains, including recent large-scale models for video generation, in-context image generation (e.g., FLUX.1 Kontext, 2025), and protein ensemble generation (Jing et al., 2024).  We selected flow matching as the foundation for our prior learning due to its numerous advantages. It has deep connections to physics (continuity equation) and offers simulation-free training with a simple learning objective. This generally leads to faster training and sampling compared to diffusion models, while also being fully invertible and highly scalable. \n\n\n\nReference: \n\n“FLUX.1 Kontext: Flow Matching for In-Context Image Generation and Editing in Latent Space” 2025\n\nJing et al, “AlphaFold Meets Flow Matching for Generating Protein Ensembles” 2024\n\n\n**W3. While an ablation study is included, a broader sensitivity analysis would improve the robustness of the empirical evaluation**\n\nWe thank the reviewer for this constructive suggestion and agree that a broad sensitivity analysis would further strengthen our empirical evaluation. Given the extensive scope of this paper, which already spans multiple fields and involves a two‑phase framework with many interacting components, we prioritized demonstrating the core viability and performance of our approach across different tasks. Therefore, while we acknowledge its importance, a comprehensive sensitivity analysis across all components was beyond the scope of this initial investigation. We consider this a valuable and important direction for future work.\n\n**W4. The paper only evaluates on synthetic GP datasets, real-world dataset such as MNIST and CelebA are missing** \n\nWe appreciate the opportunity to clarify our experimental evaluation.\n\nWhile our paper features GP datasets, the evaluation also includes several highly non-Gaussian benchmarks: $\\textbf{2D Navier-Stokes}$, high-resolution simulations of $\\textbf{black holes}$, and the $\\textbf{MNIST-SDF}$ dataset (signed‑distance functions extracted from MNIST digits). \nFor these non-Gaussian functional datasets, the standard likelihood metric is not applicable for OFM since the constant evidence term is intractable (as detailed in Appendix A, Q5), so we provide qualitative comparisons. These visualizations show that posterior samples from OFM align closely with the ground truth and preserve key features (e.g. density and swirling patterns of the blackhole) far better than those from baselines.\n\nTo further strengthen the empirical section, and following Reviewer 2aqv’s suggestion, we have also added GP regression tasks using more complex kernels (Gibbs (non-stationary) kernel and Rational Quadratic kernel (RQK)). As shown in the following table, these new results confirm that OFM consistently outperforms all comparative methods across all metrics. \n\n| **Algorithm ↓ / Metric →** | **GP‑Gibbs&nbsp;SMSE** | **GP‑Gibbs&nbsp;MSLL** | **GP‑RQK&nbsp;SMSE** | **GP‑RQK&nbsp;MSLL** |\n|---|---|---|---|---|\n| NP | $4.0 \\cdot 10^{-1}$ | $7.1 \\cdot 10^{-1}$ | $5.8 \\cdot 10^{-1}$ | $2.2 \\cdot 10^{0}$ |\n| ANP | $3.5 \\cdot 10^{-1}$ | $6.2 \\cdot 10^{-1}$ | $2.7 \\cdot 10^{-1}$ | $2.1 \\cdot 10^{0}$ |\n| ConvCNP | $3.4 \\cdot 10^{-1}$ | $1.1 \\cdot 10^{-1}$ | $2.2 \\cdot 10^{-1}$ | $6.9 \\cdot 10^{-1}$ |\n| DGP | $4.3 \\cdot 10^{-1}$ | $8.6 \\cdot 10^{-1}$ | $2.4 \\cdot 10^{-1}$ | $1.2 \\cdot 10^{0}$ |\n| DSPP | $4.2 \\cdot 10^{-1}$ | $7.1 \\cdot 10^{-1}$ | $2.5 \\cdot 10^{-1}$ | $4.2 \\cdot 10^{-1}$ |\n| OpFlow | $3.1 \\cdot 10^{-1}$ | $6.9 \\cdot 10^{-1}$ | $2.9 \\cdot 10^{-1}$ | $5.0 \\cdot 10^{-1}$ |\n| **OFM&nbsp;(Ours)** | $\\mathbf{2.9 \\cdot 10^{-1}}$ | $\\mathbf{8.7 \\cdot 10^{-2}}$ | $\\mathbf{2.1 \\cdot 10^{-1}}$ | $\\mathbf{9.4 \\cdot 10^{-2}}$ |\n\n**Q1, how Eq. 16 is computed in practice and why its complexity is claimed to be O(m)**\n\nThe divergence term in Eq. 16, which is the trace of a Jacobian, is computationally expensive to calculate exactly. We compute it in practice using the Hutchinson trace estimator, a standard, unbiased technique for neural ODE, first used by Grathwohl et al. (2018), which requires only one reverse‑mode autodiff pass through the network, so the time complexity is linear w.r.t cardinality O(m). Below is the PyTorch‑style pseudocode used in our implementation; the full script is included in the supplementary ZIP.\n\n```\n# samples: target sample repeated n_noise times along the batch dimension.\n# choose Hutchinson noise\n\nif hutchinson_type == \"Gaussian\":\n    noise = torch.randn_like(samples)              \nelif hutchinson_type == \"Rademacher\":\n    noise = (torch.randint_like(samples, 0, 2).float() * 2 - 1)  \n\nwith torch.enable_grad():\n    samples = samples.detach().requires_grad_(True)\n    \n    # time condition\n    t_vec = torch.full((samples.size(0),), t, device=samples.device)\n\n    # drift gθ(t, x)\n    drift = model(t_vec, samples)                       \n\n    # Hutchinson estimate of divergence\n    f_eps = torch.sum(drift * noise)              \n    grad = torch.autograd.grad(f_eps, samples,\n                               retain_graph=True)[0]   \n    div_estimate = torch.sum(grad * noise,\n                             dim=tuple(range(1, grad.dim())))\n```\n\n\nReference: \n\nGrathwohl et al., \"FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models.\" 2018.\n\n**Q2. It would be interesting to see how the flow matching framework performs on denoising tasks**\n\nThis is a good point.  We discussed the technical connections between these two approaches in Appendix A (Q3). While we provide this  comparison, we agree that a comprehensive empirical evaluation on denoising-specific tasks is a valuable next step, which we leave for future work"}, {"review_id": "R4CeHct1Xd", "reviewer": "Reviewer_sYHe", "summary": "The authors propose a new framework, called Operator Flow Matching (OFM), for learning stochastic processes on function spaces. In brief, the authors extend flow-based generative modelling via Fourier neural operators (FNO) to map a reference Gaussian process (GP) to an unknown target process. Formally, they define a continuous-time ODE on a Hilbert space that “flows” samples from a simple GP prior ($\\nu_0$) to the data distribution ($\\nu_1$). By parameterizing the time-dependent vector field $G_t$ with a neural operator $G_\\theta(t,\\cdot)$, they can (in principle) learn a bijection between function spaces. This yields an explicit transport map so that any finite collection of points from the learned process has a tractable probability density under OFM. Once trained, the model can perform Bayesian regression: given partial noisy observations of a function, it can compute the posterior over the function values at new points (via Bayes’ rule and stochastic gradient Langevin dynamics (SGLD) in the latent GP space). Several claimed contributions are highlighted: (1) Novelty: extending flow-matching or stochastic interpolant methods to infinite-dimensional stochastic processes via operator learning; (2) Unifying prior and posterior inference: providing a unified framework for prior sampling and exact (tractable) likelihood estimation, enabling “universal functional regression” (UFR) with learned priors; (3) Exact densities: learning an exact (invertible) mapping so that the model supplies the true prior and posterior densities over any finite evaluation of the process. Empirically, OFM is tested on several synthetic datasets (1D/2D Gaussian processes, truncated GPs, and highly non-Gaussian processes from Navier–Stokes equations, simulated black-hole data, MNIST-SDF fields).", "questions": "NA", "limitations": "Neither limitations of the proposed framework and its theoretical analysis, nor negative impact have been discussed in the main body. The high computational cost is hinted in the conclusion, but I feel that this deserves more serious report and discussion in the main body rather than in Appendix N.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 4, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\n### __1) Integration of flows and neural operators:__\n\nThe paper introduces an original approach by combining flow-matching models with neural operators. By treating functions as first-class objects, OFM learns an infinite-dimensional transport map, which is innovative. This method naturally extends normalizing flows and stochastic interpolant models from finite-dimensional data to functional data. The authors provide a clear formulation of the continuous-time flow on a Hilbert space, showing how a diffeomorphism $\\Phi_t$ pushes a reference GP measure to the data measure. The use of a Fourier neural operator (FNO) to parametrize the vector field $G_\\theta(t,h)$ is appropriate to ensure resolution agnosticism (the model can handle discretized functions at arbitrary grids).\n\n### __2) Mathematically grounded:__ \n\nThe framework is built on solid theoretical ideas from optimal transport and infinite-dimensional analysis. The authors derive a continuity equation for the time-varying measure path and explain how the log-density can be obtained by integrating the divergence of $G_t$. They further leverage Kantorovich’s formulation of dynamic optimal transport together with a linear interpolation strategy for Gaussian measures to obtain a simple vector field that couples a reference GP sample $h_0$ with a target sample $h_1$. While I am not expert in these topics, I find that these constructions are internally consistent and grounded in prior work (e.g. Chizat et al. 2018, Tong et al. 2023). In short, the method is not ad-hoc: it extends finite-dimensional flow matching to function spaces in a principled way.\n\n### __3) Tractable likelihoods and UFR:__ \n\nA key strength is that OFM can compute exact log-likelihoods for any finite set of function values. Once $G_\\theta$ is trained, one can evaluate $\\log p(u_1)$ for a new function evaluation $u_1$ by integrating as in Eq. (15). The paper carefully notes the $\\mathcal{O}(m^2)$ cost of divergence computation and employs the Hutchinson estimator to make it $\\mathcal{O}(m)$. This allows the use of Bayes’ rule for regression. In practice, the authors draw posterior samples via SGLD in the GP input space, then map them to function space. This gives a “universal functional regression” (UFR) procedure that generalizes GP regression to arbitrary learned priors. I find that the fact that OFM retains exact densities (unlike prior implicit models) is a strong advantage.\n\n### __4) Empirical evidence:__ \n\nThe experimental results are compelling. On synthetic GP tasks (1D and 2D), OFM matches or surpasses classical GP performance, and on more challenging non-Gaussian tasks it significantly outperforms baselines. The paper' supplementary material includes ablation and scaling studies to justify design choices, and also provides detailed architecture and training parameters. Moreover, the authors emphasize reproducibility: the code and logs are provided.\n\n\n\n## Weaknesses\n\n### 1) Dense presentation and clarity: \n\nWhile I acknowledge that writing operator learning papers for general ML audience is particularly challenging due to high-level of abstraction and mathematical interaction of diverse fields, I feel that the paper is mathematically dense, which may hinder its accessibility. Many readers may find the notation and measure-theoretic language heavy. For example, Section 3 immediately introduces pushforward measures on a Hilbert space and continuity equations, which can be difficult to parse. While the technical content is correct, the exposition is terse: key concepts like Gaussian measures on function spaces or the Cameron–Martin condition are assumed known. The authors do provide Appendix B–D as “foundational overview”, but the main narrative would benefit from more intuition. The introduction of the conditional flow matching in Sec. 3.2 is similarly rapid, with many equations in succession. Since NeurIPS audience is broad; more high-level guidance or illustrative text would improve clarity (e.g. a step-by-step summary of the two-phase algorithm in plain language). I am aware that this is not an easy task due to the limited length of the main body, but I would still strongly encourage authors to improve on this aspect. \n\n### __2) Novelty relative to recent work:__\n\nThe paper claims to be _the first to extend flow matching to stochastic processes__. However, it cites a recent paper (Ref [29]: “Functional Flow Matching” by Kerrigan et al., 2023) that appears to do something very similar. Indeed, the related work explicitly notes that [29] generalizes flow matching to function spaces. The authors argue that prior works (including [29]) “do not support UFR” or density evaluation, but this distinction should be made clearer. If [29] already develops the theory of infinite-dimensional flows, the novelty here is partly in combining it with operator networks and Bayesian regression. It would strengthen the paper to more explicitly state __“_unlike [29], we can compute densities and do conditional sampling_,”__ perhaps with a succinct comparison. As it stands, __the originality claim needs nuance__: OFM is novel in its application and integration (and outperforms [29] presumably), but flow-in-function-space per se is not entirely new.\n\n### __3) Theoretical assumptions and scope:__ \n\nThe method relies on strong assumptions that deserve discussion. For instance, the reference measure $\\nu_0$ is a GP with some covariance operator, and $\\nu_1$ is assumed to be absolutely continuous with respect to $\\nu_0$ (so that an invertible transport exists). In practice this means the true process must lie in the Cameron–Martin space of the GP prior. The paper does not discuss what happens if this is violated. Also, computing a flow on infinite dimensions typically requires technical conditions (e.g. continuity of $G_t$, compactness of support, etc.). While the appendices address these issues, it would be helpful to acknowledge in the main text any limitations (e.g. “we assume the target measure has full support on the same space as the GP”) and what that implies. Relatedly, the choice of a linear interpolation for Gaussian measures is convenient, but may not be optimal for complex target distributions. To me it is unclear how the approach handles non-Gaussian $\\nu_1$ during training, since the conditional flow matching still assumes $\\mu_t(\\cdot|z)$ is Gaussian. If $\\nu_1$ is far from Gaussian, the learned map may suffer. The empirical results suggest it works, but a discussion of this modeling assumption would strengthen the theory section.\n\n### __4) Empirical evaluation details:__\n\nWhile the results are promising, some experimental questions remain: \n- (a) Statistical robustness: The reported results (e.g. Table 2) do not include confidence intervals or multiple-trial averages. It is unclear how sensitive OFM is to random initialization or data sampling. \n\n- (b) Baselines and fairness: The baselines chosen (GP, deep GP, NP family, “OPFlow”[1]) are reasonable, but other recent methods could be considered. For example, Neural Processes with latent flows or diffusion-based process models (e.g. [34] Neural Diffusion Processes) might be relevant, although those are mostly unconditional. \n\n- (c) Computational cost: The paper notes in passing that “OFM demands more computational resources”, but does not quantify it in the main text. From Appendix tables we see training can take many GPU-hours and memory (dozens of GB). This could limit applicability, especially on high-resolution 2D/3D domains. Including some runtime/memory comparison or at least discussing scalability would be useful.\n\n### __5) Lack of ablation in main text:__\n\nThe contributions bullet promises _“extensive ablation and scaling studies”_, but no ablation figures appear in the main body. It would help to include at least one ablation table or figure (for example, the effect of using Hutchinson estimation vs exact divergence, or the impact of neural operator depth) in the paper, rather than hiding all such results in appendices. This would boost confidence that each component (e.g. neural operator architecture, trace estimator, SGLD sampling) is necessary and correctly implemented.\n\n### __6) Evaluation on realistic tasks:__ \n\nAll experiments are on synthetic or simulation-based data. While this is acceptable for a methods paper, testing on at least one real-world dataset (e.g. geospatial data, time series, etc.) would demonstrate broader applicability. Also, the noise model is fixed to i.i.d. Gaussian in all cases; exploring non-i.i.d. or non-Gaussian noise could be informative, since many applications (e.g. cosmology, climate) have complex noise.", "was_revised": true, "final_justification": "I keep my positive rating, voting for the clear acceptance.", "rebuttal": "We appreciate the reviewer for the thorough review, detailed summary, and constructive suggestions. We provide our responses below\n\n**W1. Dense presentation and clarity**\n\nWe appreciate the reviewer’s observation that presenting a rigorous mathematical framework spanning several disciplines within the space constraints of the main text is challenging.\n\nTo further improve clarity for readers from diverse backgrounds, we have (i) added an appendix section covering neural operators (operator learning), flow matching, Gaussian measures on function spaces, and the Cameron–Martin theorem to provide readers with essential background; and (ii) included additional explanatory text and proof sketches to better convey intuition.\n\n\n**W2. Novelty relative to recent work** \n\nGreat comment. Following the reviewer’s suggestion, we have carefully revised the text to better communicate our contributions. The updated text now more clearly highlights the specific novelty of OFM, which clarifies how our work builds upon, yet is distinct from, recent literature. Indeed it is subtle but in fact a very important difference. \n\n\n**W3. Theoretical assumption and scope**\n\nThank you for raising this point. We now state explicitly in the main text that “we assume the target measure has full support on the same space as the reference GP.”  This absolute‑continuity assumption is standard in functional generative‑model literature (see Kerrigan et al., 2023; Lim et al., 2023; Seidman et al., 2023; Schneider et al., 2024), and we follow that convention here. For verifying whether processes lie in the Cameron-Martin space of the GP prior,  we refer to the detailed proof and discussion in Appendix C of Lim et al. (2023, see refs below). We have brought this point into the main text. \n\nRegarding flow matching, we emphasise that although the conditional path $\\mu_t(\\cdot\\|z)$ is Gaussian, the marginal path obtained after integrating over $z$ is a mixture of Gaussians (weighted by the target data distribution) and therefore need not be Gaussian. This mirrors the finite‑dimensional case: in standard flow matching (Lipman et al., 2022, Eq. 11) and optimal‑transport flow matching (Tong et al., 2023, Thm. 2.1), a Gaussian conditional is paired with an unconstrained—and often highly non‑Gaussian—target distribution, such as images, audio, or video datasets. The resulting Gaussian mixture faithfully captures the complex structure of the target law. We have added an appendix section that reviews flow matching and elaborates this point.\n\n\n\nReference : \n\nKerrigan et al “ Functional Flow Matching” 2023\n\nLim et al, “Score-based Diffusion Models in Function Space” 2023\n\nSeidman et al “Variational Autoencoding Neural Operators” 2023\n\nSchneider et al “An Unconditional Representation of the Conditional Score in Infinite Dimensional Linear Inverse Problems” 2024\n\nLipman et al, “Flow Matching for Generative Modeling” 2022\n\nTong et al, “Improving and generalizing flow-based generative models with minibatch optimal transport” 2023\n\n**W4. Empirical evaluation details**\n\nFollowing Reviewer 2aqv’s suggestion, we strengthened the 1D GP regression experiments by regenerating the training and test sets with more complex kernels—Gibbs / non‑stationary, Rational Quadratic kernel. OFM still outperforms all baselines on these cases as shown in the table (replied to Reviewer 2aqv) \n\nWe acknowledge the current limitation in statistical robustness and plan a full multi‑seed validation in future work. For the baseline set‑up, Appendix A (Q5) explains why NDP is omitted from the quantitative tables and provides a detailed theoretical comparison instead in Appendix N\n\nOn computational cost, the revised draft includes wall‑clock and peak‑memory numbers, together with a discussion of practical mitigations—e.g. mixed‑precision FNO kernels (Tu et al., 2023) or low‑step ODE solvers such as DPM‑Solver (Lu et al., 2022)—that can substantially reduce GPU time and memory.\n\nReference. \n\nTu et al, “Guaranteed Approximation Bounds for Mixed-Precision Neural Operators” 2023\n\nLu et al, “DPM-Solver: A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps” 2022. \n\n\n\n**W5. Lack of the ablation in main text** \n\nTo address this concern, we have reordered several figures for a clearer narrative flow and inserted a new ablation table in the main body.\n\n\n**W6. Evaluation on realistic tasks**\n\nThank you for this valuable suggestion. We agree that this is a critical direction for future work. We intend to collaborate with domain experts to build a comprehensive benchmark of real-world datasets. This benchmark will specifically include challenging settings with non-i.i.d. and non-Gaussian noise. Developing this benchmark and extending our framework accordingly is a significant undertaking that warrants its own dedicated investigation. We are confident that this future work will enable more rigorous validation and help spur progress across many scientific and engineering fields.\n\n\n**Q1. Discuss the limitation in the main body rather than in the Appendix N**\n\nThanks for the suggestion, we have moved the discussion of the model's limitations from the appendix into the main text in the updated draft to ensure it is more prominent."}]

Author rebuttal:

Meta-review: The reviewers are uniformly positive about this paper, which they agree is interesting, and original, and would be a useful addition to the literature. The biggest concern was about how dense the presentation was, on the one hand, this is understandable given the material, at the same time, the authors made some proposals to alleviate the burden on the reader. Please go over these (along with the rest of the discussion with the reviewers) while preparing the camera-ready version.
Decision: Accept (spotlight)