Linearization Turns Neural Operators into Function-Valued Gaussian Processes
Authors: Emilia Magnani
Marvin Pförtner
Tobias Weber
Philipp Hennig
Venue: ICML 2025 spotlightposter (2025)
Keywords: Neural Operators
Uncertainty Quantification
Bayesian Deep Learning
Gaussian Processes
Partial Differential Equations
Primary area: probabilistic_methods->bayesian_models_and_methods
TL;DR: Efficient, theoretically sound, and structured predictive uncertainty quantification for arbitrary neural operators based on function-valued GPs.

Abstract: Neural operators generalize neural networks to learn mappings between function spaces from data. They are commonly used to learn solution operators of parametric partial differential equations (PDEs) or propagators of time-dependent PDEs. However, to make them useful in high-stakes simulation scenarios, their inherent predictive error must be quantified reliably. We introduce LUNO, a novel framework for approximate Bayesian uncertainty quantification in trained neural operators. Our approach leverages model linearization to push (Gaussian) weight-space uncertainty forward to the neural operator's predictions.
We show that this can be interpreted as a probabilistic version of the concept of currying from functional programming, yielding a function-valued (Gaussian) random process belief. Our framework provides a practical yet theoretically sound way to apply existing Bayesian deep learning methods such as the linearized Laplace approximation to neural operators. Just as the underlying neural operator, our approach is resolution-agnostic by design.
The method adds minimal prediction overhead, can be applied post-hoc without retraining the network, and scales to large models and datasets.
We evaluate these aspects in a case study on Fourier neural operators.

Linearization Turns Neural Operators into Canction-Valuele Gaussian Processes
Emilia Magnani" Marvin Pförtner
Tobias Weber" Philipp Hennig'
Abstract
output pairs (Boullé & Townsend, 2023). Neural operators. including Fourier neural operators, have emerged as a pow- erful class of models for operator learning, particularly for the solution operators of PDEs (Kovachki et al., 2023). They have been applied successfully across domains including weather forecasting (Pathak et al., 2022; Bonev et al., 2023), fluid dynamics (Grady et al., 2022: Renn et al., 2023; Li et al., 2022), and automotive aerodynamies (Li et al., 2023b). Instead of learning to solve specific PDE, these models learn the operator that maps functional parameter of the PDE (such as initial values. boundary conditions, force fields, or material parameters) to the corresponding solution. This approach amortizes computational cost by learning to solve entire families of PDEs.
Neural operators generalize neural networks to learn mappings between function spaces from data. They are commonly used to learn solution operators of parametric partial differential equa- tions (PDEs) or propagators of time-dependent PDEs. However, to make them useful in high- stakes simulation scenarios, their inherent predic- tive error must be quantified reliably. We intro- duce LUNO, novel framework for approximate Bayesian uncertainty quantification in trained neu- ral operators. Our approach leverages model lin- carization to push (Gaussian) weight-space un- certainty forward to the neural operator's pre- dictions. We show that this can be interpreted probabilistic version of the concept of cur- rying from functional programming, yielding a function-valued (Gaussian) random process be- lief. Our framework provides a practical yet theo- retically sound way to apply existing Bayesian deep learning methods such as the linearized Laplace approximation to neural operators. Just the underlying neural operator, our approach is resontiion-agnostic by design. The method adds minimal prediction overhead, can be applied post- hoc without retraining the network, and scales to large models and datasets. We evaluate these as- pects a case study on Fourier neural operators.
Although neural operators have demonstrated strong pre- dictive capabilities, they unable to quantify the inherent uncertainty in their predictions. Predictive uncertainty quan- tification is indispensable for many downstream tasks, such as decision-making in safety-critical scenarios. For exam- ple, a neural operator trained on past climate data should increase predictive uncertainty under distribution shifts due to climate change, reflecting potential losses in accuracy, Furthermore, uncertainty quantification is also useful for improving neural operator training via active learning strate- gies (Musekamp et al., 2025), potentially reducing the cost of generating computationally expensive numerical simula- tions as training data.
Extensive previous work has shown that the structured uncer- tainty provided by Gaussian process (GP) models is particu- larly suitable for such downstream tasks, including closed- form acquisition functions for active learning and Bayesian optimization for optimal selection of future queries or exper- iments (Garnett, 2023); enabling online model adaptation to continuously update with new data while preserving consis- tency with prior knowledge (Sliwa et al., 2024); facilitating sensitivity analysis through the GP's kernel structure reveal- ing system responses to parameter changes; and seamlessly integrating with probabilistic numerical computation (Hen- nig et al., 2022; Pförtner et al., 2022) to quantify. marginal- ize, and propagate computational uncertainty.
1. Introduction
Scientific computing increasingly demands efficient repre- sentations of complex non-linear maps between functions. Examples include solution operators of parametric partial differential equations (PDEs) or the propagators of time- dependent PDEs. Operator learning is an approach to this problem that generalizes regression algorithms from finite- dimensional to infinite-dimensiona function-space input-
Equal contribution 'Tübingen AI Center, University of Tübin- gen, Tübingen, Germany. Correspondence to: Emilia Magnani <emilia. caailiia.manaaiiee @ uni tuii-ttubbnnndee....
Motivated by the capabilities of GPs, we propose LUNO, a practical yet theoretically sound framework that provides lin- earized predictive uncertainty in neural operators. LUNO
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Lincarization Turns Neural Operators into Sinction-Valuce Gaussian Processes
a QE
a A
@ E Du
W
W
Currying
u(x) f((a, 2) w
u = F(a, w
F(a, (x))
NEU U
u(a)) Rd4
Step
Step
Lineurization weight uncertainty
LUNO
a GEA
GE
a Du
W E W
W E W
Probahilistic Currying
11 F(a) = f(a. )) ((a -), w)
u(x) = f(a, a:)
w
uEU
u(x)
Step
Step
Figure I: Illustration of the steps involved in LUNO. A trained neural operator F (top left) is converted into an equivalent neural network f with outputs & using (reverse) currying (top right). Linearizing f around the mean of the Gaussian weight belief results in a Gaussian process posterior quantifying the uncertainty about the function learned by f (bottom right). Finally, probabilistic currying transforms into function-valued Gaussian process posterior F over the operator learned by the neural operator F (bottom left).
quantifies uncertainty over the mapping learned by the neu- ral operator via a Gaussian process with values in separable Banach space functions-4 higher-order generalization of GPs that, when evaluated, returns a Gaussian random function rather than dimeessoinnii Gaussian random variable, as in standard GPs. This function-valuee GP is constructed through model linearization from Gaussian distribution quantifying uncertainty in the neural operator's (fnite-dimensione) weight space. We show that LUNO can be interpreted as a probabilistic generalization of the concept of currying in functional programming. This con- nection makes LUNO compatible with established methods for quantifying weight-space uncertainty in deep neural net- works, including the Laplace approximation (Ritter et al., 2018; Daxberger et al., 2021a; Kristiadi et al., 2020; Papa- markou et al., 2024), SWAG (Maddox et al., 2019), or mean- field variational inference (Blundell et al., 2015). LUNO is practical, introduces minimal computational overhead, and can be applied post-hoc, without requiring to retrain the neural operator. It scales to large models and datasets and, like neural operators, is inherently resolution-agnosia We demonstrate the capabilities of the framework in a case
study on Fourier neural operators.
LUNO is designed to be compatible with arbitrary (non- Gaussian) weight-space beliefs (see e.g. Appendix A.4). Nevertheless, we focus our exposition and experiments on Gaussian weight-space uncertainty, as, the future, we aim to explore the use of the resulting sian process in downstream tasks such as the ones outlined above.
valued Gaus-
In Section 2, we review the fundamentals of neural operators and (multi-output) Gaussian processes. Section 3 presents our main contribution. We first develop Gaussian processes that take values in (infiiiteddimeniional) Banach spaces of functions, along with the notion of probabilistic currying, which formalizes their equivalence to (multi-output) Gaus- sian processes. Using probabilistic currying, we construct function-valued Gaussian processes from neural operators with Gaussian weight beliefs. Section 4, we discuss prior work on operator learning and related uncertainty quantifica- tion. Finally, we demonstrate in Section the effectiveness of LUNO in common PDE learning settings.
2

---

Lincarization Turns Neural Operators into Function- Valuud Gaussian Processes
2. Background
g(v(L((x), wg) with
We first review neural operators, with emphasis on Fourier neural operators, which serve as the primary case study for our analysis. Then, we provide an overview of multi-output Gaussian processes.
dy Σ
o ((+)) (x) = o01)
(R,)) [(0)){()))))))()))) F(()
k=1 (x)
1=1
+ W/(
(1) 5 (a:)
2.1. Neural Operators
Neural operators (NOs) (Kovachki et al., 2023) are neural network architectures that map between (infinite- dimensional) Banach spaces of functions. A neural operator is a function F: A x W + U, where
for 1 = I...., and a (x)) = p(a(x)) wp) ( E R where J denotes the Fourier transform of a periodic func- tion.' p: R4 Wp - R' and q: R² Wa + R are parametric functions called lifting and projection, re- spectively. R()) E Chiea x xd', and W(I) E Rd. and w ppp W(1, R(),, map w()) + ((+1) is the I-th Fourier layer. If the inputs a are discretized on a regular grid, F can be computed by a real fast Fourier transform (RFFT).
W(L-1), ( wyy) The
A is (separable) Banach space of functions a: 2D2
R4 with domain D) C dRd
U is a (separable) Banach space of functions u: By -
R4 with domain Dy C Rdd,
2.2. (Gaussian) Random Processes
Aiming to generalize the notion of a Gaussian process later, we provide its formal definition from mathematical statis- tics that is rarely used in machine learning. For set ? and set F of functions on Ω with values in measur- able space, we denote by (FF) the smallest o-algebra for which all f E F are measurable. The Borel a-algebra on a topological space Q is denoted by B (2). A random process on probability space (2, P) with index set A. and values in a measurable space (S, As) is a function f: A x Ω + S such that f(a, -) is A-essuuuurall for all a E A. We use the shorthand ((a) == ((,,)) One can show that w r f(:,w) is function-valuee random vari- able with values in (IZA,0(0a)) (5x)), where d denotes the set of point evaluation functionals on the set A of functions from A to B. A random process is called Gaussian or a Gaussian process (GP) if it has values in (R, B (IK) and w + (I(a), Man, w) is an RR-valued Gaussian ran- dom variable for all n E N and 01.....00 On A. The mean function of is given by + Ep [f(a)] and the covariance function of f is given by (a1, a2) + Covp [f(ar) ) ((()), f(az)). We denote by f ~ GP (m, k) that I' is a Gaussian process with mean function m and covariance function k.
W is a set of parameters (typically W C IR or W C (P).
To keep training tractable, neural operators are trained on (1)
datasets {(a() (x. A ), u o (x )) consisting of pairs
of input and corresponding output functions (a(*) ()) E A x U that are discretized at finitely many points x( E (D&)"² and x E (Dw) respectively. The training objective is typically given by the empirical risk
R(w)=
1 n
L(u
(6) (X
F(a
(
(X (x)), w) (X ())
or a regularized version of the empirical risk. Neural opera- tors were originally developed, and are commonly used, to learn the solution operator of non-linear, parametric partial differential equations. In this case, typically, DA Dy, the input functions A correspond to parameters and/or ini- tial conditions of the PDE, and the output functions u E U are the corresponding solutions of the PDE (at later time points). There are many different realizations of the abstract neural operator framework, including low-rank neural op- erators (Kovachki et al., 2023), (multipole) graph neural operators (Li et al., 2020b:a), and (spherical) Fourier neural operators (Li et al., 2021; Bonev et al., 2023).
It is common to extend the concept of a GP to finitely many output dimensions. A d' d'-ouuuut Gaussian process f is a random process with values in (IRd, B (IRd)) such that w + (f(a),w)
an -valued Gaus- sian random variable for all n E N, and a1...., an A. f(an,
T) T)
We use the shorthand f(a) == f f(a, -). The mean function of f is given by + Ep [f(a)] E RR and the covariance func- tion of f is (a), 02) + Covp (((a)),((a)) E Rd xa'. We
Example 2.1 (Fourier Neural Operators). As a case study we will focus on Fourier neural operators (FNOs) (Li et al., 2021), popular variant of the neural operator ar- chitecture that applies all spatially global operations in the spectral domain. An FNO F transforms a periodic input function a into a periodic output function F(a, ux) =
'More precisely. the operator F: L2 RR + (2(C) maps real-valued
s square-integrable function on the d-dimensional torus to the coefficients the corresponding Fourier series.
3

---

Lincarization Turns Neural Operators into Function- Valued Gaussian Processes
denote by f ~ GP (m, K) that f multi-output Gaussian process with mean function and covariance function K. While the notion of a multi-oulput Gaussian process might seem more general than the notion of a Gaussian process, it is possible to "emulate" a function with multiple outputs by augmenting the input space of a Gaussian process:
w - (F(a1, ...,, F( w)) is a jointly Gaussian ran- dom variable for all E N and 01,--. On A.
As above, we use the shorthand F(a) == F(a,)) -). Moreover, the map w + F(..w) is a random variable with values in the space of (linear and non-linear) operators (U*, (6x)). This warrants the interpretation of U-valued Gaussian processes as Gaussian random operators.
Lemma 2.1. Let A,) be a probability space, f: A x Ω + R24 I= [I'....d'}, and I: (A x II) x Ω + R with (f(,,))); = f((a, i): for all a E A and
(P-almost surely). Then f ~~ GP (m, K) if and only if ~ GP (m, k), for all E & and I.
In the context of neural operators, U is Banach space of - functions on common domain Dr. In this case, we can show that U-valued Gaussian processes are closely related to multi-output Gaussian processes with an augmented input space. This is in analogy to Lemma 2.1, but requires some additional technical assumptions.
(mla)); m(),)),
as well as, for all a1, 02 E A and j E II
Theorem 3.2 (Probabilistic Currying in Banach Spaces; proof in Appendix A.3). Let A,P be a probability space and U real separable Banach space TR -valued functioos with domain Dy. Let F; A x S2 - V and f: (A X Du) x Ω + Rd such that F(a.-) = [((a.x). -) for all a E and x E Dry (P- (P-almost surely). Then (i) F is a random process with values (U,r(88) if and only if f is a IRd -valued random process, (ii) F is Gaussian if and only if f is Gaussian, and (iii) if all evaluation maps 8m U - Tedd" u ++ u(x) are continuous, then (i) holds for F with values in (U))
(K (2))s = M(((1,)) (a2.5))
3. LUNO: Linearized Predictive Uncertainty in Neural Operators
In this section, we show how obtain linearized predic- tive uncertainty in neural operators (LUNO). We lever- age model linearization to propagate Gaussian weight- space uncertainty through the neural operator to its pre- dictions. LUNO can be applied to trained models as a post-processing step, and does not require expensive re- training. Furthermore, LUNO employs the framework of function-valued Gaussian processes. To that end, we first develop the concept of a and draw an important parallel with currying in functional programming, which offers a natural interpretation of our method. Figure I illustrates the main steps comprising our methodology.
Theorem 3.2 reveals an insight into the abstract concept of function-valle Gaussian processes: Function-valus Gaus- sian processes are equivalent to (multi-output) Gaussian processes with augmented input spaces. This equivalence enables the translation of real-valued GPs, a computationally feasible structure, into diminte-imensionall function-valued objects.
Gaussian process
3.1. Cunction-Valued Gaussian Processes and Probabilistic Currying
Probabilistic Currying We note that Theorem 3.2 con- stitutes probabilistic analogue of the concept of cur- rying from functional programming (and category the- ory more generally). The Theorem shows the equiva- lence of the vector-valued (Gaussian) random function f: A x Bu + IRd and the (Gaussian) random operator F: A + (Du + Rd") with (())()) as f(a,z).
We want to use model linearization to extend the Gaussian belief over the parameters of a neural network f: RRd x Rº + RRd into a (multi-output) Gaussian process belief over the function learned by the neural network. However, this is not immediately applicable to neural operators, since their outputs do not lie in Rd",, but in a potentially infinite- dimensional Banach space of functions. Hence, we need to generalize (multi-output, Gaussian processes to the notion of a Ranach-valued Gaussian process.
Example 3.1 (Currying a Continuous Bivariate Gaussian Process). Let I ~ GP (m,kk be a bivariate 2-output Gaus- sian process with compact index set X1 x Xg IR2 on P) with (P-almost surely) continuous paths. For in- stance, this assumption is fulfilled if m is continuous and A is a multivariate Matérn covariance function (Da Costa et al., 2023). Then a + f(a, .) is a function-valued Gaussian process. More precisely, Theorem 3.2 shows that the map F: X, Ω + C(x2), (a, + (x + f((a, w)) is a ((%2)-valuee Gaussian process with index set X,.
Definition 3.1 (Banach- Valued Gaussian Process). Let U be real separable Banach space and LL a set? of linear functionals on U. A random process F: A Ω U on probability space (Ω, P) with index set A and values in (U, (IL)) is called Gaussian or a Gaussian process if
Technically, L needs to separate the points in U. See Ap- pendix A.1 for additional details.
"See Definition A.4 and Remark A.S.
4

---

Lincarization Turns Neural Operators into Punction-Valuce Gaussian Processes
Thus, an intuitive way to understand function-valued Gaus- sian processes is objects that, when evaluated, return Gaussian process. Currying can also be used to relate the mean and covariance functions of canction-valueed or more general vector-valued (Gaussian) random processes, and their counterparts defined on the corresponding multi- output (Gaussian) random process. As this discussion is rather technical, we defer it to Appendix A.3.
The function f is linear in the weights, but it re- mains highly nonlinear in the input. Moreover, have fir (: = f(; ju) which means that the mean function of f matches the prediction of the trained neural operator (before linearization) we set ja to the weights wv* found during training.
Step 3 Probabilistic currying constructs a Gaussian ran- dom operator from f. Namely, we define the function
Appendix A provides in-depth explanation of our theoret- ical framework and contains a plethora of theoretical results on Gaussian processes with values in arbitrary (infinite- dimensional) vector spaces that are not necessarily Banach or function spaces. For example, such results are vital for quantifying uncertainty in neural operators applied to PDEs that only admit weak solutions (see Appendix A.4).
F: x Ω - U. (a,w) (æ + w))).
(For this F be well-defined, we need to assume that f((a, ,,,w) is U for all a E A. See also Appendix A.4), Theorem 3.2 then shows that F is a U-valued Gaussian process. Moreover, E [x)) = F(a, (z), and
Cov (((1 (2)) (x2))
3.2. L Linearization Turns Neural Operators into Ianetion-Valuee Gaussian Processes
D ((), (x1)
E Dw (a2, w)
We use the notion of function-valued Gaussian processes to develop LUNO. We delineate the key components into different steps, visually represented in Figure 1.
The entire construction, in particular Theorem 3.2, still ap- plies if w is not Gaussian. In this case, f and F are no longer Gaussian, but the formulae for the mean and covari- ance functions remain valid. We derive a generalization of the method for general separable Banach spaces T in Appendix A.4. For instance, this is useful if the functions in are not pointwise defined, such as weak solutions of PDEs.
W + U C (R&))oo be a neural operator in Section 2.1 with W RRP.
Step 0
Let F: F:
Step By uncurrying F, we define the function
f: (AXDU) Dc) W- + Rd ()..x),,) + F(a,w)
3.2.1. CASE STUDY: FOURIER NEURAL GAUSSIAN RANDOM OPERATORS
Step We obtain a Gaussian belief w ~ NW (Jz, Σ) over the parameters of the network. In Bayesian deep learning, common way to obtain this Gaussian belief is by plac- ing a Gaussian prior p(w) on the network's parameters and then approximating the posterior distribution given the data p(w D). Well-established (approximate) inference tech- niques to obtain the posterior over w include the Laplace approximation (Ritter al., 2018; Daxberger et al., 2021a;
The exposition so far applies generally to neural operators. For Fourier neural operators, a particularly efficient repre- sentation of the function-valued posterior process is avail- able. To simplify the exposition, we focus on the case where the Gaussian belief is restricted to the parameters of the final Fourier block WL-1 = (R!L-1) W(I-1). This is common approach in the context of last-layer Laplace approximation (Kristiadi et al., 2020). In Appendix C.I, we show that the fanction-valued GP obtained by applying LUNO in this case takes the form
Immer et al., 2021), variational inference (Graves, 2011: Blundell al., 2015; Khan al., 2018), and SWAG (Mad- dox et al., 2019). Since f has values in 12 following Khan et al. (2019); Immer et al. (2021) and Appendix B, we can linearize the model around jr:
F(a)(x) = (( (w)) + (Dq (m
(x) m2( (z))
f((a,z), w) xe (a, x), = =f((a,x),H) Dw Dwf((en f(aa, w)) (w- 11)
i.c., F(a) ~ GP (ma Ka) with
ma(x) F(a, w" (xx) and
to arrive at an induced approximate dutout Gaussian process belief f == flin (-, (-,w) ~ GP (m, K) with index set A x Du. m(a, x) = f((a, x). 14), and
Ka x2) Dq (m, (1-1)) (x1) K (2 x2)
Dg (my - (x2))
where (l-- ~ GP (m K ( is a multi-output parametric Gaussian process whose moments only depend on (LL-1), nt and E, and q = qi. wg) o (L-11)
K((an, x1 x22) = =D w), Σ Dw [((2, 22). w)
5

---

Lincarization Turns Neural Operators into Sanetion-Valucee Gaussian Processes
There are two practical benefits arising from this representa- tion. First, computing the moments of, and drawing samples from F(a) only needs access to the hidden state v(4-1) of the neural operator. We can thus evaluate the Gaussian process belief at arbitrary output points æ E Do, without the need to compute more than one (full) forward pass of the neural operator. Secondly, since the Gaussian process belief F(() over the output function is parametric, we can efficiently sample entire functions from it that can then be lazily evaluated at arbitrary points. This is in contrast to gen- eral con-parmmttric Gaussian processes, where one typically discretizes the GP before drawing samples of the function values at the given finite set of points. Such lazy functional samples can be used e.g. for active experimental design and Bayesian optimization (Wilson al., 2021).
theoretical framework Owhadi (2023); Batlle et al. (2024) embeds into ours.
To generate a probabilistie belief over a neural network's s weights, various Bayesian posterior approximation tech- niques are available. One of the most popular is the Laplace approximation, introduced to deep learning by Mackay (1992), which has gained popularity in the Bayesian deep learning community (Ritter et al., 2018: Daxberger et al.. 2021a; Kristiadi et al., 2020; Papamarkou et al., 2024). This is also due to its scalability, achieved through various strate- gies including using log-posterior Hessian approximations (Ritter et al., 2018; Martens, 2020), treating only a subset of the model probabilistically (Daxberger et al., 2021b), em- ploying linearized Laplace (Foong et al., 2019; Immer et al., 2021), or using scalable Gaussian processes methods (Deng et al., 2022: Ortega et al., 2024). Other Bayesian deep learn- ing methods include variational inference (Graves, 201 I; Blundell et al., 2015; Khan et al., 2018; Zhang et al., 20,2018 8), Markov Chain Monte Carlo (Neal, 1996; Welling & Teh, 201 1; Zhang et al., 2020). SWAG (Maddox et al.. 2019), or heuristic methods (Gal & Ghahramani, 2016; Maddox et al., 2019). Finally, a widely used approach for uncertainty quan- tification in deep learning is ensembles (Lakshminarayanan et al., 2017: Hansen & Salamon, 1990), that train multiple independent neural networks with different random initial- izations and aggregate the predictions.
4. Related Work
Azizzadenesheli et al. (2024) provide a comprehensive overview of neural operator architectures. These include graph neural operators (Li et al., 2020a), physics-informe neural operators (Li al., 2024), multi-wavelet neural oper- ators (Gupta et al., 2021) and the widely used Fourier neural operators (Li et al., 2021). FNOs have gained particular prominence, finding applications across various PDE prob- lems (Pathak et al., 2022; Zhang et al.. 2023; Li et al., 2023a; Rashid et al., 2022: Qin et al., 2024; Kossaifi et al., 2023; Bonev et al., 2023). Theoretical foundations for FNOs have been established, with Kovachki et al. (2021) proving their universal approximation capabilities for continuous oper- ators, and Lanthaler et al. (2024) analyzing discretization- induced aliasing errors.
5. Experiments
We evaluate linearized predictive uncertainty (LUNO-+) against sample-based approaches (Sample-*), which require additional approximations to impose Gaussian Process structure over the output space. To be precise, in the Sample- * methods, we draw samples from the weight-space belief, map the samples through the (nonlinear) map w + F(,, w). and compute function-valuu Gaussian process belief over the prediction by moment matching the empirical mean and covariance function. We consider isotropic Gaussian (*-Iso) and low-rank Laplace approximated (x-LA) weight- space uncertainties, in both their sample- based (Sample) and linearized (LUNO) forms. We compare these weight-space- Gaussian methods against input perturbations (Pathak et al., 2022), and deep ensembles. Deep ensembles were trained 10 times with different random seeds on the original Fourier neural operator (FNO) architecture. We evaluate our model time- deeendepnndntt PDEs in one and two spatial dimensions, predicting the next time step autoregressively from the pre- vious ten. We assess uncertainty quantification in two key settings: (1) a low-data regime, where the model is trained on a limited number of trajectories, and (2) out-of-distribu- tion (OOD) scenarios, where physical phenomena unseen during training are introduced test time.
While neural operator architectures have advanced, incorpo- rating uncertainty estimation remains challenging. Recent work has approached this problem from different angles. Garg & Chakraborty (2023) applied variational inference to estimate Bayesian posteriors in DeepONets. More closely related to work, Magnani et al. (2022) developed uncer- tainty estimates for graph neural operators using Laplace approximation, though their approach does not extend to FNOs, nor does it consider function space formulations. Kumar et al. (2024) combined Gaussian Process priors with Wavelet Neural Operators, optimizing hyperparameters through negative log-marginal likelihood minimization. Ad- ditional Bayesian operator frameworks have been explored by Garg& & Chakraborty (2022); Batlle et al. (2024); Zou et al. (2024); Mora et al. (2025).
Function-valued Gaussian processes have been studied in the Hilbert space setting by Owhadi (2023); Batlle et al. (2024). Our approach formulates the theory natively within the context of Banach spaces, as neural operators are defined as mappings between such spaces. In the Appendix we prove that, when restricted to the Hilbert space setting, the
We evaluate the predictive uncertainty using standard met-
6

---

Lincarization Turns Neural Operators into Function- Canction-Valued Gaussian Processes
Input Perturbations
Ensemble
Sample-L.
LUNO-LA
0.5m
0.04
0.0
Figure 2: FNO predictive uncertainty quantified by several different methods. Top row: target function (-), mean (-) and 1.96 standard deviations of, as well as samples (-) from, the predictive belief. For the ensemble, the samples are four of the ensemble members. Bottom row: spread of the predictive distribution around the mean. For the sample-/ensemble-bant methods, we construct a Gaussian distribution from the empirical covariance matrix and draw four samples (-). We plot 1.96 standard deviations of the predictive belief, well as the top-three eigenfunctions (-) and heatmap of the predictive covariance matrix (top right corner of panels).
rics: the expected root mean squared error (RMSH) of the mean predictions, the expected marginal x² statistics, and expected marginal negative log-likelihood (NLL) over 250 test input-output pairs. Hyperparameters are optimized via grid search using the expected marginal NLL valida- tion set the target. Full details, including data generation, training procedures, uncertainty estimation methods, and more detailed results are provided in the Appendix.
Method
RMSE (4)
X²
NLI (4)
Input Perturbations
x 10-2
0.894
-18720
Ensemble
3.49 X 10
5.597
-0.8145
Sample-Iso
3.72 x 10
0.977
-1.341
LUNO-Iso
3.62 x 10-2
0.864
-1.9488
Sample-LA
5.59 10
2.774
-1.1572
LUNO-LA
3.62 x 10 1.022
-2.0787
Table 1: Comparison UQ methods for an FNO trained on 25 trajectories of Burgers' equation.
Code. We provide an efficient implementation of the LUNO framework in JAX (Bradbury al., 2018) at
MethodsoOMMachhneeLLamnngg luno.
Method
Base
Flip
Pos-Neg-Fih
Input Perturbations
-2.886
2.573
494.935
The code for our experiments can be found at
Ensemble
-5.313
3.825
-1.014
Sample-Iso
-2.921
4.071
43.362
O / 2bys/ luno-experiments.
LUNO-Iso
-2.892
3.450
37.733
Samplc-LA
-2.576
4.395
27.046
LUNO-LA
-2.931
-1.126
1.164
Low data regime. We train an FNO for 100 epochs on 25 simulated solutions Burgers' equation with 59-time steps and evaluate their uncertainty 250 unseen test pairs. Figure visualizes the predictive uncertainty for input per- turbations, deep ensemble, Sample-LA, and LUNO-LA on single test data point of Burgers' equation. Table shows that LUNO-LA outperforms the other approaches. This trend holds across two other one-dimensional time- dependent PDE datasets, which are included in the Ap- pendix (Table and Table 5). While all methods produce marginal confidence bands around the network prediction, their sample path covariances differ qualitatively.
Table 2: Expected marginal NLL evaluation across OOD
datasets for different methods. Lower is better.
field. We introduce various additional physical phenomena to the test set. These include reversing the velocity field at the center (Flip), introducing a triangular heat source (Pos), and cloud-shaped heat sink (Neg). Table reports ex- pected marginal NLL over a variation of out-of-distribution datasets. Additional and more granular results can be found in the Appendix. While LUNO-LA outperforms the other weight space methods and input perturbations, deep ensem- bles achieve the lowest expected marginal NLL in next-step prediction. However, their uncertainty representation is fun- damentally different. Figure 3 compares deep ensemble
To assess OOD robustness, we train an FNO (or an ensemble 10 FNOs) on two-dimensional dddvection-iiffusse equation with initial conditions sam- pled from Gaussian blobs and random constant velocity

---

Lincarization Turns Neural Operators into Function- Vanuti Gaussian Processes
Ensemble
LUNO-LA
Target
Residual
Pred. Std.
Residual
Pred. Std.
0.2
0.04
0.2
0.2
0.1
0.0
0.1
0.03
Xx²
Sample
Null space projection
Sample
10
0.3
0.025
0.0
0.0
Figure 3: Comparing an ensemble (left), LUNO-LA (right). Top row shows target, residuals, and the predictive standard deviation. Bottom row shows the absolute ratio of the pointwise residual and the predictive standard deviation as well as sample from the predictive belief. Since the uncertainty structure of the ensemble prediction is of low rank, we also include its unexplained error by projecting the residual vector onto the null space of the predictive covariance.
with LUNO-LA. Deep ensembles approximate uncertainty using a small set of discrete hypotheses, represented by col- lection of point masses in parameter space. While this repre- sentation is not confined to the analytic form of Gaussian distribution, has other constraints: For example, although marginal uncertainty estimates (panel 3 in the figure) can be relatively well-structured, the associated empirical covari- ance across the ensemble is fundamentally rank-deficient. This limitation is critical, as it leaves certain types errors entirely unaccounted for (panel 8, which projects residuals onto the null space of the ensemble covariance). By contrast, LUNO-LA constructs covariance matrix whose rank is (in theory) only bounded by the number of parameters con- sidered.4 4 As result, plot like panel 8 in Figure 3 does not make sense for LUNO-LA, since, in principle, explains any variation in the data (albeit with varying calibration). This behavior also evident in full-trajectory evaluations. Although FNOs are trained for next-step prediction, they are often used for auto-regressive roll-outs, where predictions recursively fed back as inputs. Such roll-outs cause subtle yet significant distribution shift, prediction enrors accumulate and are treated as ground truth for subsequent steps. While the deep ensemble improves upon the network prediction terms of RMSE, its uncertainty estimate does adapt the increasing error, reflected in the NLL (cf. Figure 4).
methods are reported the Appendix. Due to the efficiency of Jacobian- vecttor vector products and analytical tractability of the inverse real fast Fourier transform, LUNO-* methods outperform their Samples-* counterparts, with LUNO-ISO being even faster than the deep ensemble in our implemen- tation. Each method comes with its own additional cost. While deep ensembles need fully separate training runs with different random seeds, LUNO-LA's main computational bottleneck is computing the low-rank approximation the generalized Gauss -ewstnn matrix (GGN). This cost dom- inated by network size, the selected rank, and the amount of data used for the GGN approximation.
6. Conclusion
We introduced LUNO, a framework for predictive uncer- tainty quantification neural operators using function- valued Gaussian processes. LUNO can be interpreted as probabilistic generalization of currying in functional pro- gramming. By leveraging model linearization, it offers computationally efficient and theoretically grounded ap- proach to incorporating weight-space uncertainties in neural operators. The framework endows neural operators with structured weight-space uncertainty quantification capabili- ties while preserving their resolution-and agnoosic nature. We demonstrate this for LUNO-LA in the FNO setting under low-data regimes and out disttribuuiinno scenarios.
Wall-clock times for single-trajectory predictions across all
LUNO's s main limitation lies in the challenges associated with modeling weight-space covariances. Nevertheless, by successfully constructing a structured Gaussian process over
*Numerical instabilities and (near) singular Jacobians might reduce the rank in practice.
8

---

Lincarization Turns Neural Operators into Vunution-Valued Gaussian Processes
References
0.4
RMSE
Aliprantis, C. and Border, K. C. Infinite Dimensional Analysis: Hitchhiker's Guide. Springer- Verlag. Berlin, Heidelberg, third edition, 2006. 29587-9.
0.2
10a
NLL
102
Azizzadenesheli, K., Kovachki, N., Li, Z., Liu-Schiafini,
M., Kossaifi, J., and Anandkumar, A. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, pp. 1-9, 2024.
10
30
Time step
Batlle, P., Darcy, M., Hosseini, B., and Owhadi, H. Kernel methods are competitive for operator learning. Journal of Computational Physics, 496: 12549, 2024.
Figure 4: Averaged performance of different UQ methods an autoregressive rollout of the FNO on 50 trajectories from the Pos-Neg-Fliip dataset. We compare input perturbations
Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. Weight uncertainty in neural network. ICML, 2015.
(-), deep ensembles (-), Sample-Iso (-), LUNO-Iso
) Sa (-), LUNO-LA (-).
Bogachev, V. Gaussian Measures, volume 62 Mathe- matical Surveys and Monographs. American Mathemati- cal Society, Providence, Rhode Island, 1998.
Bonev, B., Kurth, T., Hundt, C., Pathak, J., Baust, M., Kashinath, K., and Anandkumar, A. Spherical Fourier neural operators: Learning stable dynamics on the sphere. ICML, 2023.
the output space, LUNO paves the way for future applica- tions of GP-valued neural operators in scientific and engi- neering domains.
Boullé, N. and Townsend, A. A mathematical guide to oper- ator learning. arXiv preprint arXiv:2312.14668 2023.
Impact Statement
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne S., and Zhang, Q. JAX: composable transformations of Python+NumPP programs, 2018. URL http://githon com jax//aaximmll jax.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which feel must be specifically highlighted here.
Conway, J. B. A Course in Functional Analysis, volume 96 of Graduate Texts in Mathematics, Springer, New York, second edition, 1997. doi: 1w.111007978-11 101110007799988-147774-833
Acknowledgments
The authors gratefully acknowledge financial support by European Research Council through ERC CoG Ac- tion 101 123955 ANUBIS the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: 011S 18039A); the DFG SPP 2298 (Project HE 71 7114/--,), and the Carl Zeiss Founda- tion, (project "Certification and Foundations of Safe Ma- chine Learning Systems Healthcare"), as well funds from the Ministry of Science, Research and Arts of the State of Baden- Wirtte-errttbbrg..... The authors thank the Interna- tional Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Emilia Magnani, Marvin Pfört- ner and Tobias Weber. We would like to thank Nathaël Da Costa for invaluable feedback on the theoretical aspects of this work.
Da Costa, N., Pförtner, M., Da Costa, L., and Hennig, P. Sample path regularity of Gaussian processes from the co- variance kernel. arXiv preprint arXiv:2312.1488. 2023. doi: 1.s11055/////Xiv223111
Dangel, F., Tatzel, L., and Hennig, P. ViVIT: Curvature access through the generalized Gauss-Newton's low-rank structure. TMLR, 2022.
Daxberger, E., Kristiadi, A., Immer, A., Eschenhagen, R., Bauer, M., and Hennig, P. Laplace redux-effortless Bayesian deep learning. NeurlPS, 2021a.
Daxberger, E.,, Nalisnick, E., Allingham, J. U., Antoran, J., and Hernandez-Lobato, J. M. Bayesian deep learning via subnetwork inference. In ICML. 2021b.

---

Lincarization Turns Neural Operators into Paaetton-Valluee Gaussian Processes
Deng, Z., Zhou, F,, and Zhu, Accelerated linearized Laplace approximation for Bayesian deep learning. NeurIPS, 2022.
Khan, M. E., Nielsen, D., Tangkaratt, V., Lin, w., Gal, Y., and Srivastava, A. Fast and scalable Bayesian deep learning by weight-pertuubation in adam. In ICML, 2018
Foong, A. Y., Li, Y., Hemández-Lobato, J. M., and Turner, R. E. "In-Between' uncertainty in Bayesian neural net- works. arXiv preprint arXiv:1906. 11377, 2019.
Khan, M. E., Immer, A., Abedi, E., and Korzepa, M. Ap- proximate inference turns deep networks into Gaussian processes. In NeurlPS, 2019.
Gal, Y. and Ghahramani, Z. Dropout as a Bayesian approxi- mation: Representing model uncertainty deep learning. In ICML,, 2016.
Klenke, A. Probability Theory: A Comprehensive Course. Universitext. Springer, London, second edition, 2014.
Koehler, F,, Niedermayr, S., Westermann, R., and Thuerey, N. APEbench: A benchmark for autoregressive neural emulators of PDEs. In NeurlPS, 2024.
Garg, S. and Chakraborty, S. Variational Bayes deep opera- tor network: data-driven Bayesian solver for parametric differential equations. arXiv preprint arXiv:2206.056655, 2022.
Kossaifi, J., Kovachki, N., Azizzadenesheli, K., and Anandkumar, A. Multi-grid tensorized Fourier neu- ral operator for high- high-resolution PDEs. arXiv preprint arXiv:2310.00/20, 2023.
Garg, S. and Chakraborty, S. VB-DeepONet: A Bayesian operator learning framework for uncertainty quantifica- tion. Engineering Applications of Artificial Intelligence, 118:105685, 2023. ISSN 0952-1976.
Kovachki, N., Lanthaler, S., and Mishra, S. On univer- sal approximation and error bounds for Fourier neural operators. JMLR, 22(290):1-76,202. 2021.
Garettt, R. Bayesian Optimization. Cambridge University Press, 2023. 101:1011017/971111 108348977.
Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Learning maps between function spaces with applications to PDEs. JMIR, 24(89):1 2023.
Grady, T. Khan, R., Louboutin, M., Yin, Z., Witte, P. A., Chandra, R., Hewett, R. J., and Herrmann, F. J. Towards large-scale learned solvers for parametric PDEs with marallp Fourier neural operators. arXiv preprint arXiv:2204.00 2005 2022.
Kristiadi, A., Hein, M., and Hennig, P. Being Bayesian, even just a bit, fixes overconfidence in relu networks. In ICML, 2020.
Graves, A. Practical variational inference for neural net- works. In NeurlPS, 2011.
Kumar, S., Nayek, R., and Chakraborty, S. Neural operator induced Gaussian process framework for probabilistic solution of parametric partial differential equations. Com- puter Methods in Applied Mechanics and Engineering, 431:117265,, 2024.
Gupta, G., Xiao, X., and Bogdan, P. Multiwavelet-based operator learning for differential equations. In NeurIPS, 2021.
Hansen, I.. K. and Salamon, Neural network ensem- bles. IEEE transactions on pattern analysis and machine intelligence, 12(10):993-1000 1001, 1990.
Lakshminarayanan. B., Pritzel, A., and Blundell, C. Simple and scalable predictive uncertainty estimation using deep ensembles. In NeurIPS, 2017.
Hennig, P., Osborne, M. A., and Kersting, H. P. Probu- bilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022.
Lanthaler, S., Stuart, A. M., and Trautner, M. Discretiza- tion error of Fourier neural operators. arXiv preprint 2wXi::2405.02221, 2024.
Immer, A., Korzepa, M., and Bauer, M. Improving predic- tions of Bayesian neural networks via local linearization. In AISTATS, 2021.
Li, Z., Kovachki, N., Arizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Multipole graph neural operator for parametric partial differential equations. In NeurIPS, 2020a.
Kabri, S., Roith, T., Tenbrinck, D., and Burger, M. desolution-invarian image classification based Fourier neural operators. In Calatroni, I.,, Donatelli, M., Morigi, S., Prato, M., and Santacesaria. M. (eds.), Scale Space and Variational Methods in Computer Vision, pp- 236-249, Cham, 2023. Springer International Publishing. doi:
Li, Z., Kovachki, N., Azizzadeneshelii K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Graph kernel network for partial differential equations. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020b.
10

---

Lincarization Turns Neural Operators into Function- Valued Gaussian Processes
Li, Z.,, Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Fourier neural operator for parametric partial differential equa- tions. In ICLR, 2021.
Musekamp, D., Kalimuthu, M., Holzmüller, D., Takamoto, M., and Niepert, M. Active learning for neural PDE solvers. ICLR, 2025.
Neal, R. M. Bayesian Learning for Neural Networks. Springer-Verless Berlin, Heidelberg, 1996. ISBN 0387947248.
Li, Z., Peng, W. Yuan, Z., and Wang, J. Fourier neu- ral operator approach to large eddy simulation of three- dimensional turbulence. Theoretical and Applied Me- chanics Letters, ((66)100358 2022.
Ortega, L. A., Santana, S. R., and Hernández-Lobato, D. Variational linearized Laplace approximation for bayesian deep learning. ICML, 2024.
Li, Z., Huang, D. Z., Liu, B., and Anandkumar, A. Fourier neural operator with learned deformations for PDEs on general geometries. JMLR, 24(388):1-22 2023a.
Owhadi, H. Do ideas have shape? Idea registration as the continuous limit of artificial neural networks. Physica D: Nonlinear Phenomena, 444, 2023.
Owhadi, H. and Scovel, C. ApperatorAdaptee Wavelets, Fast Li, Z., Kovachki, N., Choy, C.,, Li, B., Kossaifi, J., Otta, S., Nabian, M. A., Stadler, M., Hundt, C., Azizzadenesheli, K., et al. Geometry-ing -nfometr-informed neural operator for large- scale 3D PDEs. NeurIPS, 2023b.
Solvers, and Numerical Homogenization: From a Game Theoretic Approach to Numerical Approximation and Al- gorithm Design. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2019. doi: 11.11.1119991 108594967.
Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. Physies- informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 133)1-277, 2024.
Papamarkou, T., Skoularidou, M., Palla, K., Aitchison, L., Arbel, J., Dunson, D., Filippone, M.. Fortuin, V., Hennig, P., Hernández-Looatoo J. M., Hubin, A., Immer, A., Kar- aletsos, T., Khan, M. E., Kristiadi, A., Li, Y., Mandt, S.. Nemeth, C., Osborne, M. A., Rudner, G. J., Rügamer, D., Teh, W., Welling, M., Wilson, A. G., and Zhang, Position: Bayesian deep learning needed in the age of large-scale AI. In ICML, 2024.
Loshchilov, I. and Hutter, F. Decoupled weight decay regu- larization. In ICLR, 2019.
MacKay, D. J. Bayesian interpolation. Neural computation, 4(3):415-447, 1992a.
MacKay, D. J. The evidence framework applied to classi- fication networks. Neural computation, 4(5):720-736, 1992b.
Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., Kurth, T. Hall, D., Li, Z., Azizzadenesheli, K., et al. FourCastNet: A global data-driven high-resolution weather model us- ing adaptive Fourier neural operators. arXiv preprint arXiv:220..12114 2022.
Mackay, D. J. C. Bayesian methods for adaptive models. California Institute of Technology, 1992.
Maddox, W., Garipov, T., Izmailov, P., Vetrov, P., and Wilson, A. G. A simple baseline for Bayesian uncertainty in deep learning. In NeurIPS, 2019.
Pförtner, M., Steinwart, I., Hennig, P., and Wenger, J. Physics-informed Gaussian process regression general- izes linear PDE solvers, arXiv preprint arXiv:2212.12474 2022.
Magnani, E., Krämer, N., Eschenhagen, R., Rosasco, L., and Hennig, Approximate Bayesian neural operators: Uncertainty quantification for parametric PDEs. arXiv preprint arXiv:2208..733 01565, 2022.
Qin, S., Lyu, F., Peng. W., Geng, D., Wang, J., Gao, N., Liu, X., and Wang, L. L. Toward a better understanding of Fourier neural operators: Analysis and improvement from a spectral perspective. arXiv preprint arXiv:2404.07200 2024.
Martens, J. New insights and perspectives on the natural gradient method. JMLR, 21((46):: 2020.
Rashid, M. M., Pittie, T., Chakraborty, S., and Krishnan, N. A. Learning the stress-strain fields in digital com- posites using Fourier neural operator. Iscience, 25(11), 2022.
Martens, J. and Grosse, R. Optimizing neural networks with Kronecker-factorde approximate curvature. In ICML, 2015.
Mora, C., Yousefpour, A., Hosseinmardi, S., Owhadi, H., and Bostanabad, R. Operator learning with Gaussian processes. Computer Methods in Applied Mechanics and Engineering, 434:1 2025. ISSN 0045-7825.
Renn, P. L., Wang, C., Lale, S., Li, Z., Anandkumar, A., and Gharib, M. Forecasting suberitical cylinder wakes with Fourier neural operators. arXiv preprint arXiv:2301.0220 2023.

---

Lincarization Turns Neural Operators into Function- Valuud Gaussian Processes
Ritter, H., Botev, A., and Barber, D. A scalable Laplace approximation for neural networks. In ICLR, 2018.
Schraudolph, N. Fast curvature matrix-vector products for second-order gradient descent, Neural computation, 14: 1723-38, 2002.
Sliwa, J., Schneider, F., Bosch, N., Kristiadi, A., and Hennig, P. Efficient weight-space Laplace-Gauee Gaussian filtering and smoothing for sequential deep learning. arXiv preprint 1.K::2410006800, 2024.
Steinwart, I. Conditioning of Banach space valued Gaussian random variables: An approximation approach based on martingales. arXiv preprini urXiv:2404.03-55 2024.
Welling, M. and Teh, Y. W. Bayesian learning via stochastic gradient Langevin dynamics. In ICML, 2011.
Wilson, T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Pathwise conditioning of Gaussian processes. JMLR, 22(105):1 11,1177222 2021.
Zhang, G., Sun, S., Duvenaud, D., and Grosse, R. B. Noisy natural gradient as variational inference. In ICML, 2018.
Zhang, R., Li, C., Zhang, J., Chen, C., and Wilson, A. G. Cyclical stochastic gradient MCMC for Bayesian deep learning. ICLR, 2020.
Zhang, T., Trad, D., and Innanen, K. Learning to solve the elastic wave equation with Fourier neural operators. Geophysics, 88(3)-1001) 2023.
Zou, Z., Meng, X., Psaros, A. F., and Karniadakis, G. E. NeuralUQ: A comprehensive library for uncertainty quan- tification neural differential equations and operators. SIAM Review, 66(1) 66(1):611-190, 2024.
12

---

Linearization Turns Neural Operators into Function-Valuee Gaussian Processes: Supplementary Materials
A. Theoretical Results
A.1. Dual Spaces
We aim quantify epistemic uncertainty about learned (nonlinear) operators F: A. - U, where U is a (typically infinite- dimensional) real vector space of functions. If U is space of real-valued functions, the very least, want to be able express probabilistic belief over all point evaluations (()) = dz (F(a)). Note that point evaluation dn U - R of functions in U is linear map, since dg (Q: U1 +02u2) = (ar t/1 a22u)(zz)) a1 (x) = a a1 +a2 (u2). Many interesting operations that map functions into real numbers like (point-evaluated) derivatives and integrals are linear.
Now let be an arbitrary real vector space. Real-valued linear maps on U are referred to as linear functionals. The set of all linear functionals on referred to as the algebraic dual (space) of and denoted by U# A subset L of U is said to be total to separate the points in IJ if for any "1: Uig E U with i u2. there is & IL such that f(u,) E(u2). Such subsets useful, since they allow to identify elements from the primal space U uniquely. For instance, the set of all point evaluation functionals on a vector space of real-valued functions separates the points in the space. If UJ is topological vector space (for instance separable Banach space in the context of neural operators), then the subspace of continuous linear functionals is denoted by U' C U#
Remark A.1 (The Bidual Embedding). The algebraic dual space U with pointwise addition and scalar multiplication real vector space itself. Hence, any subspace IL U# has algebraic dual space L# The elements IL# of this space are linear functions mapping linear functionals into real numbers, i.e., Φ(E) R for E L U# L is vector space of real-valued functions, so we can consider its point evaluation functionals 8, IL - R2 ( f(u). Note that the map AU,T U + L# + 8, is linear and, if L separates the points in U, injective. Hence, U is isomorphic to its image bu = fU,1 (U) under 4U,L*. We refer to the map f0 as the bidual embedding. Abusing notation, we write U C L. u E L# for u E U, ctc.
A.2. Probability Measures on Vector Spaces
Our framework models the predictive uncertainty over an output of a neural operator as a random variable with values in an (imnite-dimensionaa vector space U functions. As noted before, we at least want to quantify the uncertainty about a given set L C U of linear functionals. Hence, we need to make the linear functionals in L measurable.
Let (2, A, P) be probability space, U real vector space and L U vector subspace of linear functionals separating the points in U. We equip U with the smallest r-algebra o (I)) that makes all the functionals in L measurable. A random variable with values in (U, (IL) (U-valued for short) is an A-o meemmmsurablle function u: Ω - U.
Similar to their finite-dimensional counterparts, probability measures on and random variables with values in (infinite- dimensional) vector spaces admit the definition mean and a
operator.
Definition A.2 (Mean and Covariance Operator (see e.g., Bogachev, 1998, Definition 2.2.7)). Let be probability measure on (L).
I C (T)), then m, I defined by
(e) = E, 10 = ((w))(du)
e
is called the mean of y. The mean my of a random variable u: Ω + with values in (U, (L))( defined as the mean Mpou of its law.
13

---

Lincarization Turns Neural Operators into Function- Vanetion-Valued Gaussian Processes
(b) IL C L2(1), then the linear operator C,: L IL. defined by
(G)))) = Cov, [41 f2] ( (e(()) - my (42))))))) 62,62 EL
is called the covariance operator of 7. The covariance operator C, of a random variable u: - U with values in (U, (L)) is defined as the covariance operator Crou of its law.
Definition A.3 Croos-Covvariance Operator). Let Un, ug: : Ω + U be random variables with values in (U, oL)) such that L ( P oo, u,") for = 1,2. The operator Cur,uz L + L defined by
(al.u2ee) (G2) Cov (c,) J (f](u((w)) mu, (41) (62(U2(w)) Muy (42) P(dw)
Cu uz
is called the cross-covariance operator between and
Gaussian measures on U are defined by generalizing the closure properties of Gaussian measures on R
Definition A.4 (Gaussian Measure (see e.g., Bogachev, 1998, Definition 2.2. (a))), A probability measure y (U, (L)) is called Gaussian if every linear functional & I. is univariate Gaussian random variable on (I))) A random variable u: Ω2 + U with values in (U, ((L)) is called Gaussian if its law P u- is Gaussian.
If L is not a vector space, then we say that random variable is Gaussian with values in (U, if and only Gaussian with values in (U, (span Note that (sspan L) (IL) for any subset L U## (by Klenke, 2014, Definition 1.79 and Theorem 1.91).
Remark A.5 (Jointly Gaussian Measures). To define Gaussian processes with values in U, we need the notion of a joint Gaussian measure on U Fortunately, we can also leverage Definition A.4 for this. U" is vector space under elementwise addition and scalar multiplication. Its algebraic dual space (U*)** is isomorphic to (U y", i.c. & & (U*)## if and only if there are ,,,..., Ln U## such that
E(u1..... Un)
(us)
Vui...., Un E U.
follows that =(L)) where the latter denotes the product o-algebra. Hence, following Definition A.4, we call a probability measure 7 on (U", (T) Gaussian if every II." is univariate Gaussian random variable on (U", (((L) .7).
Remark A.6 (Probability Measures on Separable Banach Spaces). The case where a real separable Banach space is of particular interest in the context of neural operators. In this case, choose TL = U', i.e. all linear functionals that continuous with respect to the norm topology. Then the o-algebra a o(U)) coincides with the Borel o-algebra B (U) generated by the norm topology (Bogachev, 1998, Theorem A.3.7). For Gaussian random variables u with values in (U, B (U)) (and Gaussian mcasure on B (U)), the mean my is an element of U and the covariance operator maps into U (Bogachev, 1998, Theorem 3.2.3). Moreover, for jointly Gaussian random variables U1, U2 with values in (U, B (U)), the cross-covariance operator Cu maps from U' to U (Bogachev, 1998, Theorem 3.2.4).
A.3. Random Processes with Values in Vector Spaces
Now we have all the necessary preliminaries to define a Gaussian process with values in (U,o(L)).
Definition A.7 (Gaussian Process). Gaussian process with index set A and values in (U,)(L) on (R,,, is function F: A x Ω -> TJ such that r> (F(11,4)...., F(an, F(a,,, w)) is joint, i.e., (U", (U", (L))@n (L.))))))alluaaa )-aaued, Gaussian random variable
for all n E N and ........a. On E A.
As for real-valued or Rd Rd-vulued processes, we can also define mean and covariance functions for random processes with values arbitrary real vector spaces. However, their definition is more technically involved.
Definition A.8. Let F be random process with index A and values in (U,)(L)) on (9, P)P
(a) If c (PP -1) for all a e A, then the function
M: A + L aar mF(a..)
is called the mean function of F.
14

---

Lincarization Turns Neural Operators into Snaction-Valuede Gaussian Processes
(b) I C Lg F(a, +) for all a A. then the function
K: A A (IL - L*), (a1. 02) + CF CF(o).)) F.((2.))
is referred to the covariance function of F.
Remark A.9 (Moments Banach- Banach-Valued Gaussian Processes). If U is a separable Banach space, then mean function M takes values in U and the covariance function K. takes values the space of nuclear operators U' + (Bogachev, 1998, Theorem 3.11.24).
In the following, aim establish a correspondence between (Gaussian) random processes with values in (U, a (L)) and (Gaussian) random processes with values in (RR B (R)), which we dub (generalized) probabilistic currying. Unlike in Lemma 2.1, we need additional technical assumptions for this to work both ways. Denote by scl.. ([.) = (E L hee C L: & - (} the weak- sequential closure of a set E C L, where & -uw if and only if (z(u) - {(u) for all E (Aliprantis Border, 2006, Section 5. 5.14).
Assumption A.10. Let É C L set of linear functionals on U such that there is nse No with scl (span [L) L. Theorem A.11 (Generalized Probabilistic Currying). Ler L LL be set of linear functionals separating the points in U. Let F: A x Ω -> V and f: (A x L) x Ω -> IR such that L(F ((F(w, = [( f), w) for all E A., l E Ti and P-almost all E S2
(i) F is random process with values in (U,)(L)), then is random process with values in (R, (R)). and
(ii) F
F is Gaussian, then so is
If Assumption 1.10% is satisfied, then the reverse implications hold as well.
We will need the following generalization of Theorem B.6 from (Pförtner et al., 2022).
Lemma A.12. Let c L such that Assumption A.10 holds. Then L separates the points in U. Moreover, a function
u: Ω + is
(a) (L..--)eeaaurall i o is A-B (R) measurable for all é f E
(b) Gaussian random variable with values in (IL) (fi u,...,, en u) is jointly Gaussian for all n E N and li,...,ey E Lt
Proof. Define Yhel with Lo == span L and Ln+1 == sclu (Ln). By assumption, Linee = L.
Assume that L does not separate the points in U. Then there is U such that E(u) for all l E L. We proceed by induction. Pick f.. Then there are .1,...,, Om E R and .1,...,en, lm L such that = Em asli Hence,
{(u)= Σ 1=1
aul(l)) 0.
Now assume that ((u) = 0 for n < nset and all l T.n Fix l Thut1. Then there is (li)ien C f.n such that la -ww l. Hence,
((u)
lim (s(w)= 0. 1-000
All in all, it follows that L = Lnn does not separate the points in U, which is contradiction. Henee É separates the points in U.
(a)
We need to show that l o u measurable for all E (Klenke, 2014, Theorem 1.81). We proceed by induction. Let & Thg-- Then there are on,.... ........y, Omm E R and f,,... em such that f. = Σ a,l, By assumption, l, o measurable for all i 1,..., ......, m. Hence, e is measurable by Theorem 1.91 in (Klenke, 2014). Now assume that f measurable for all E Ln with n < Tarl Fix E Ln+1- Then there a sequence (4) C L,, such that E, wx L. This implies that f, ou - lo u pointwise, where, by the inductive hypothesis. é, u measurable for all e N. Hence, is measurable for all Ln+1 (Klenke, 2014, Theorem 1.92).
WWe will drop the A-B (TR prefix for real-valued functions this proof.
15

---

Lincarization Turns Neural Operators into Function- Valuee Gaussian Processes
(b) By (a), is A-o((L) Definition 2.2.1()) for all n = O,.... nach which is well-defined, since IL separates the points U. Again, we proceed by induction on n. Let & E Lgg Then there are a1..... am E R and ei...., lm Σ such that Ei=1 2,lz. By the closure properties of Gaussians under linear maps, have that o is Gaussian. Hence, is Gaussian with values in (U, ((Lo)), Now assume that is Gaussian with values in (W, ((,)) for nscl. Fix l E L Then there is sequence (ei In such that & " L. This implies that & o u + l pointwise, where, by the inductive hypothesis, & is Gaussian for all i N. Since pointwise limits of Gaussians random variables are Gaussian, & Gaussian. Hence, is Gaussian with values in (U, ([,+(()
It suffices show that is Gaussian with values in (U, (U,)(Ln)) (Bogachev, 1998,
Proof of Theorem A.II. (i) Holds by definition.
( Let a.,..,,,, an E A and l ln L. By Remark A.5, the linear functionals
4. U - IR (ui,..., un) ) + Ei(us)
are measurable with respect to @T.)8n. Moreover, w (F(a1, w), ... w)) is Gaussian by assumption. Hence,
()(F(a),))).....
i=1
(()F((;, (f((ar, (i)
Gaussian with values R..
+i) Follows directly by applying Lemma A.))(a) to each F(a, individually. Let 01,-.., On A. We have show that + (F(a1,w),... F(an,w) w]) is Gaussian random variable with values in U is easy check that Assumption A.10 holds for I C L Hence, by Lemma A.12(b). r (F(a1, -..., F(an,w)) is Gaussian with values in U
The following Corollary shows that Assumption A.10 automatically fulfilled for IL U in separable Banach spaces.
Corollary A.13 (Generalized Probabilistic Currying in Separable Banach Spaces). Let U be a real separable Banach space and L C U a set of continuous linear functionals separating the points U. Let F: A x Ω - and (A x L) x Ω R such that ((F w)) = f((a, (), for for all a E A, & L, and P-almost all w e Ω.
(i) F is random process with values in (U,B (U)) if and only if is a random process with values in (I. (. (RR), and
(ii) F is Gaussian if and only if is Gaussian.
Proof. We will show that Assumption A. fulfilled for -
and nacl = 1. Let UL R) (E f(u)),
which is linear and injective, since i separates the points in U. Then Ü == (U)) is isomorphic to (as a vector space). Hence, equipped with the norm || = ||e (d)]lu lu is a real separable Banach space which isometrically isomorphic to U. Moreover, Ü is space of functions with continuous point evaluation functionals, since
(8<(((u)) = |((w)) < || ||u||u = (I()ll)
Thus, there is a sequence (ōe, Jien O' separating the points in U (Steinwart, 2024, Theorem 4.10). This implies that (LieeN C Il separates the points in W. Finally, follows that T scl,,, (span f.) (Steinwart, 2024, Proposition 4.3). LI
16

---

Lincarization Turns Neural Operators into Cunction-Valued Gaussian Processes
U C DR is a vector space of real-valued functions and IL span op., then Theorem A.11 and Corollary A.13 become
substantially sharper.
Corollary A.14 (Probabilisticc Currying), Let U Re be vector space of real-valued finntionn, F; A x Ω + U, and f: (A x Dr) x f + R such that F(a, (x) = (((aa))) w) for all a E A, E Dry, and P-almost all w E Q. Then
(i) F is random process with values in (U, ar (δb.)) if and only if is a random process with values in (R,,(R)).
(ti) F has mean function M with values in U if and only if has mean function m, where M(a) m(a, .) for all A,
(ii) F has covariance function K. with values in span ōbg U if and only if has covariance function k. where (14))02)652) k((a1, x), (02, for all a1, 02 A, and x Dy, and
(iv) F is Gaussian if and only if is Gaussian.
If U a separable Banach space with continuous point evaluation functionals then (v) (i) and (iv) hold for F with values in (U, (U)), and,
(vi) if it exists, then the covariance function in (iii) has values in U' + U, where (()11)))(()(2)) = L(k((a1, (02,))))
for all E U and x E Dey
Proof. Follows from Theorem A.11 and Corollary A.13.
Finally, Theorem 3.2 from the main text is merely a corollary of the results developed above.
Theorem 3.2 (Probabilistic Currying in Banach Spaces; proof in Appendix A.3). Ler (22, A,, P) be a probability space and a real separable Banach space of R -valued fanctions with domain Ty. Let F: A. x S2 -> U and f: (A. x Dvy - TR such that F(a,-) values in (6,)) if and only if is IR." -valued random process, (ii) F is Gaussian if and only if is Gaussian, and u(x) continuous, then (i) holds for F with values in (U, B (UJ).
f((a, æ) for all E and Dy (P-almost surely). Then (i) F is a random process with
(iii) if all evaluation maps dz: U - Rdd
Proof. Follows from Corollary A. and Lemma 2.1.
A.4. Banuch-Valued Gaussian Processes from Linearized Neural Operators
For simplicity of the exposition, we limited the construction of the LUNO framework in Section 3.2 to neural operators, which map into Banach space of functions with continuous point evaluation functionals. This limits its applicability, especially for solving PDEs, whose solutions are often not defined pointwise, but rather elements of Sobolev spaces WP. (Du) c LP(DD)). Hence, in this section, we extend LUNO to neural operators that map into abstract separable Banach spaces U.
Step 0
Let be a real separable Banach space, a set, W C IR a subspace, and F: A W + a neural operator.
First, we select a subset L C U' of continuous linear functionals separating the points in U, for which we want to quantify predictive uncertainty under the neural operator. Define
Step 1
f: (A x L) x W + R. ((a, E), w) + L(F(a, w)).
This is an uncurried version of the neural operator. To see this, note that neural operator F: A x W - can be uniquely identified with the function
F F + (É. R), ()) w) - (e + ((F(a, w))).
'By generalization of Lemma 2. an analogous result holds in vector spaces of TRd" -valued functions. With minor modifications, the construction below works in arbitrary real vector spaces.
17

---

Lincarization Turns Neural Operators into Function-Valuue Gaussian Processes
Proof. Let F1, P2: A x W U be neural operators with F1 + Fg. Then there E A and W such that Fa w) * Fa(a, w). Since L separates the points in U, this implies that there is l E L such that Fi (a, u)) Fi(a, F2(a, w) =: (a ())... Hence, Fi F2.
Step 2 We model the uncertainty over the parameters random variable w: Q + W on probability space (Ω, A, P) with supp(w) w) = W. As in Section 3.2, we will now linearize f in w around point W. To achieve this, we assume that the directional derivatives
[)(a,)), wo hur) [((a, un)
aw f(a,(())) (wo)
lim h-0
h
wh exist for all a E A. 8 E it, and W, and linear in w. For instance, this is the case if f( E), .) is differentiable at up. In this case, the linearization of is given by
f((a, f), w) = f ((a, f), w) = f((a, E), wo) + ww (oo-w f((a, (), (wo).
Then the function
f: (A x L) - R, ((.))).))))))) + ((a, ((w)))
is random process. If w has mean μe, then f has mean function
m: + R. (a, E) > f((a, E), wo) + Op- f), (wo),
and if w has a covariance matrix Σ, then has a covariance function given by
k: (Ax (Ax + R.
as((a, -), (1y (2), (wvo))
4,==
Note that f and m are linear in é and k is bilinear in (41, f2). Moreover, if w is Gaussian, then is Gaussian process.
Step Finally, we construct a U-valued (Gaussian) random process F: A x Ω + U by probabilistically currying f. However, this more challenging for an abstract Intuitively, we want to undo the uncurrying operation from Step 1. To this end, we assume that there F: A x Ω with ((F(a,w)) f((a, (), for all A, 8 E IL, and P-almost all Ω. In this case, Corollary A.13 ensures that
(i) F is random process with values in (U,
iii) F has mean function M: A + U with
{(M(a)) = m(a,)) = (((a,))),)))) + fu-wo {), .) (wo)
if w has mean vector fa,
(iii) F has covariance function K: A A + (U' + U) with
((K(()))))(()) k((a1, G1) (oz, (02, Σ a.f((a,)) (, (2). (wo)
if w has a covariance matrix Σ, and
(iv) F is Gaussian process if w is Gaussian.
'See Appendices A.4.1 and A.4.2 for more details on this assumption.
18

---

Lincarization Turns Neural Operators into Punction-Valuue Valued Gaussian Processes
1.4.1. WEAK GÅTEAUX DIFFERENTIABILITY
The existence of F:; A x with ((F(a,)) f((a, w) for all E A, E L, and P-almost all w S2 is equivalent to F(a, .) being rUU, (,t,LLL)GGttauss differentiable at un, i.e., there is linear operator F -) (wo) W -> such that
F(a, w + hw) F(a, wo)
8F .) (wo) (w)
lim h-0
in (U, (U, (),
(A.1)
h
w E W, where (U,L*) smallest topology on U for which all functionals L are continuous.
Proof. Note that, since W is a metric space, we can take sequential limits in Equation (A.1) without loss of generality.
Fix A. We will constuct the Gâteaux derivative from F(a, Since F(a, w) EuLf (U) for P-almost all Ω2 there is N E A with P(N) = 0 and F(a, w) E 'u,L* (U) for all w E 22 N. We have supp w W and hence there
are wI....; Wd E Q such that {6.) with b; == w(ws) wo is basis of W. Let A: W be the unique linear operator with
\(b,) "uii" (F(a, wi)) - F(a. wo)
"UIL (f((a, w,, f((a, :), wo))
= "VULL (E + dw(wi) uo [(a, , E), -) (wo))
= - (e (E+ f6.f((a, (), (wo)),
i.e. X(w) " (e + Swa f), -) (w))) forall w E W. Since W is finite dimensional, À is Tw-7(U (U, IT continuous, where Tw the norm topology on W. Let (hanneN C R be any null sequence. Then for any & E we have
F(a, wi hn w) - F(a, wo) f((a, E), wo + h,w) - f((a,(((aao w))
hn
hn
n- Om ( .) (wo)
((((u)))
Hence, F(a, :) -(U, )
GGâteau differentiable at wo with Gâteaux derivative Fa, (wo) A.
* IIf F(a, is TUJ # C GGâteaux differentiable at we for all E A, then we can construct F as
A - + F(a, wo) + 8F(a .) (wr)) (w(w) wo)
(A.2)
Then
((F(a,w)) = ((F(a, wo)) + ww 3) (wo) f((a, E), wo) + (w E), (wo) plin ((a, , 1), w(w))
= f((aa w)
for all E A. E I. and Ω.
A stronger condition, namely (norm) Fréchet differentiability, has been verified for Fourier neural operators mapping between Ln spaces (Kabri et al., 2023).
Note that we can also use the properties of the Gâteaux derivative to verify the conclusions of Corollary A.13. Since W dimesiiimensional, the Gâteaux derivative 8F 8F(a, (wr)) is continuous with respect to any TVS topology on U. Hence,
19

---

Lincarization Turns Neural Operators into Function- Sunction-Valued Gaussian Processes
defined in Equation (A.2) is random process with values in (U, = (U, (U))))))))U (U)) If w has mean μa, then the mean function M: A -> U of F is given by
M(a) F(a, und) + (F .) (wrp) ( -Wg):
and, if w has a covariance matrix Σ, then the covariance function K:: A x A + (U' + U) of F is given by
X((11 02) = 5F(m1, .) (wo) (a) (wo)).
By the closure properties of Banach-valued Gaussian random variables under continuous affine maps (Bogachev, 1998, Lemma 2.2.2), F is a Gaussian process is Gaussian.
While this construction is somewhat more direct, the currying approach outlined above mimics more closely how F is constructed on computer, especially when U C (R and I = {0n;; E DD, i 1....,dd ,yy as in Section 3.2, and does not require knowledge of Banach-valued derivatives.
A.4.2. BIDUAL RANDOM PROCESSES
there no F: A x Ω - with ((F(a, f(( (),w) for all a E A, l E I, and P-almost all w E Q, we can construct a weaker version of F. Note that, if = aj6i + a2l2 E L. then [((a, () an anf((a, w) + a22l((a, E2), w). This means that (((a, w),w) can always be uniquely linearly extended to I, = span f. Define
F: A L (a.w) + (e + f((a, w))
We can use Corollary A.14 to show that F is a random process with values the algebraic dual L o equipped with the smallest alageerra (8L=) that makes all point evaluation functionals measurable. We refer to such random processes as bidual random processes.
Recall the bidual embedding fu,* from Remark A.1. I exists, then F(a, w) u,,1 (),
A.5. Operator- Operator-Valued Gaussian Processes Hilbert- Gaussian Processes
Finally, we show that operator-valued Gaussian processes (Owhadi, 2023; Batlle al., 2024; Mora al., 2025) can be embedded in the theoretical framework outlined above. Since operator-valued valued Gaussian processes are only defined on separable Hilbert spaces, in this section, we let A, U be separable Hilbert spaces with inner products (:-)) and (s tuu respectively. In this case, the continuous dual U is identified with the primal space U via the Riesz isomorphism. We start by reviewing the building blocks of operator-valued Gaussian processes.
Definition A.15 (Operator- Valued Kernel (Owhadi, 2023, Definition 9.1)). A function K: A x A. - (U + U) is called an operator-valued kernel if K((a1, a2) is a bounded linear operator with K((a1, az) = K(a2, a1)" for all a1 02 E A and - K(aa, a,) uy))) 2 for all m N, ai,..., am E A, and WI,..., tmm E U.
Definition A.16 (Gaussian Hilbert Space (Owhadi & Scovel, 2019, Definition 7.1)). A closed subspace of L2 A, P) called Gaussian Hilbert space every h H is a univariate Gaussian random variable.
Definition A.17 (Operator- -Operator-Valued Gaussian Process (Owhadi, 2023, Definition 5.1)). Let M: A + K.K: A U) operator-valued kernel, and IH Gaussian Hilbert space over (52, A, A function E: A + (U + HH) called
an operator-valued Gaussian process if ((a) is bounded linear operator for all E A. E is said to have mean M and covariance kernel KC if (u) ~ ((u, M(a)),, (u, Kíaa, a)))) for all a E A and E U, and
Covp [(ar) ((2) = (u1, a2)=)))
Generally, operator-valued Gaussian processes are "equivalent to" subset of Gaussian processes with values in (U')# = 1 (see Appendix A.4.2).
Proposition A.18. Let P: x U# H L2(,, a Gaussian Hilbert space, and E: & - (U + HI such that F(a, (u) E (u) for all a E A and u E U. Then E is an operator-valued Gaussian process with mean function M: A U and covariance kernel K: & & (U U) if and only if
(a) F (U#, (((t)((())))me)))) vvalued Gaussian process with mean M: A. + W and covariance funnttioniinnctionn R: AxA - (65 - -')
20

---

Lincarization Turns Neural Operators into Panction-Valuede Gaussian Processes
(b) (u,w) + F(a, (w) is mean-square continuous' for all a E A. and
(c) u un ((a)) E is norn-continuous for all a1, E A. We have (M(a), -)U M(a) for e A. and (K(ai, (2)( ))u
for all d1, 02 E A and U.
Proof. This follows from Corollary A.14.
Moreover, the following results show that operator- valued Gaussian processes whose kernels map into the trace-class "equivalent to" Gaussian processes with values in (U, B (U)).
A, P) with index set A. mean function M, and covariance Proposition A.19. Let F be U-valued Gaussian process (Ω, A,
function K. Define H as the A,,,, closure of span w + (u, F(a, w)/u I a E A, E U} C A.P)) Then
E: A + - H), a (u F(, .)))
is an operator-valued Gaussian process with mean M and covariance kernel K. Moreover, K(a, a) is trace-class for all a E A.
Proof. Let a E A. For all u E U we have
[E(o)(W) - J. ((.,F(a))))) s lulli [F(a,, (l(
by the Cauchy-Sehwarz inequality and Fernique's theorem (Bogachev, 1998, Theorem 2.8.5), i.e. E(a) is a bounded linear operator.
Theorem A.20. Let E be an operator-valued Gaussian process with Gaussian Hilbert space H C (2, mean function M. and covariance kernel K. such that K(aa, a) trace class for all A. Then there U-valued Gaussian process F ~ GP KC) (2,, A, P) such that (u, F(a))o (u) for all E A and U.
Proof. First, we will construct F: x + U. Let E A. Since U is separable and K(aa, a) is self-adjoint, positive- semidefinite, and trace class (and hence compact), there is an ONB {wf her of U consisting of eigenvectors of K(aa a) (Conway, 1997, Corollary 5.4). For every E I, fix 2 (E))) - (nbf M(a))). We will now show that the series
Σ 2 v
(A.3)
iET
converges P-almost surely in U. The {zf are centered, independent Gaussian random variables with variances given by the eigenvalues Xς of (((a, a), since
Ep 12 = Ep (Ba)(( (VS, M(a)) =
and
Cove [x4,29] = Cove
(W Ki(a, a) (169)) (WF, (1r))) X50
We have Eier Xo 30 because the operator K(aa a) trace-class, By Theorem 1.1.4 in Bogachev (1998), this means that
Σ
2
IET
iEf
converges P-almost surely, and hence Equation (A.3) P-almost surely absolutely convergent in U. Define
a.s. F(a)
M(a) + Σ
=
iE1
F(a, (u) 12, F(a, )uo) as a 2 uo
21

---

Lincarization Turns Neural Operators into Punction-Valuue Gaussian Processes
Then (u, F(a)) E (u) for a span[y? her-
We will now show that (u, F(a)) E (())(u) for all E A and E U. To this end, fix E U and ya,u E {(a)(u). There {un J== C spanf fier such that u = lim Un. Boundedness of E(a) implies (uns F(a)) L2 and, by Remark 6.11 and Corollary 6.13 in Klenke (2014), there exists a subsequence (unz 1ce [UnIn= such that (Uny F(a)) a. Ya,u- Moreover, by continuity of the inner product, we know that (un,, F(a, + F(a w)) for all Ω2 and hence ( F(a)) - F(a). Since almost sure limits are unique up almost sure equality (Klenke, 8.6. (u, F(a)). Hence, (u, F(a)) E E())(u).
Ya
2014), have Ya,u
The rest of the claim now follows from Corollary A.13 with f((a, u), = (u, F(, and É U.
B. Linearized Laplace Approximation
The linearized Laplace approximation (LLA) (MacKay, 1992a;b; Immer et al., 2021) is conceptually simple, yet effective (Daxberger et al., 2021a) method for obtaining an approximate posterior distribution over the parameters w IR of a neural network f: _2d x I² + IRd., It applies whenever the objective function R used to train the neural network is (cquivalent to) a negative log-powe pog-ptsterior
R(w)= = - p(w D) = - -lwg Σ log p fxx(() w)) const.
of the network parameters given data D (ax
(+) ()) It common for the prior the parameters be Gaussian, in which case - log p(w) acts as an 12-regularizer on the parameters. During training, we attempt to find local minimum wn of the objective function R, i.e. maximum a-posteriori (MAP) estimator of the network parameters given the data.
Following Immer et al. (2021), we approximate the posterior of the network weights as follows: First, we linearize the model using a first-order Taylor approximation in the weights around the MAP estimator wt
f (x, un)= > fli (ar, w) == far w²) + o un) (w un*
Afterwards, we compute second-order Taylor approximation of the negative log-posterior Relnn of the linearized network at the MAP wn
Rein (w) = - -log p(w) (
Σ
log (1 flii (zz(0), w)) const.
1=1
R(w=) + V V (un))
(w w)) + 22 (ur w²) TP)TP((W w²)
R30
-w ww) P(w w²) const.
(B.1)
with H Hw log (w) ) w G. where
G G= - Dw. (a w) Hy log (1)
|5) f()) Duo If fxx()) w))
(B.2)
the so-called generalized Gauss-Newton (GGN) matrix (Schraudolph, 2002). The GGN is guaranteed be positive- semidefinite. Equation (B.1) is the negative log-density of (potentially degenerate) multivariate Gaussian distribution with mean w" and covariance matrix P i.e.
p( w w D) exp R(w) ss exp Pelin (w) RJ N (w: pt).
This Gaussian distribution is referred to as the linearized Laplace approximation of p(w w DD. Under the linearized model, the approximate Gaussian posterior over the weights induces a tractable posterior predictive over the output of the neural network (Khan et al., 2019; Immer et al., 2021). More precisely, using closure properties Gaussian distributions under affine maps, one can show that the pushforward of the LLA posterior through w - flin (xx, w) defines (d' -otpuut Gaussian process
f 1 GP (f(,, (.) 122) + D (x) w) w- pl Dio w)) w)) (o)
22

---

Lincarization Turns Neural Operators into Function- Valuud Gaussian Processes
C. Implementation Details
C.1. Last-Layer LUNO for FNOs
For an input a E A, we can factorize the FNO as
F(a, (x) = (g(:: oo(l--1)) (2)
=
with
&
2
(1-1) w1.11) = Σ x-1
Ri 2 200 o0f-11) (2z) _ w - 00- (2r)
(I.-1)
d' Rana
Re(He
Re(e kj
cos (wwkk x))
1=1 R=1
(
d Run
Im R
Kij
(1--1) (-1) (-,- sin (ww2, x))
d
wil-11 f (z)- (z).
+
monky (oz)
where o([-11) = F(of(-1)) C for le
We note that (11-1) (xx, WL-1) is linear in wL-1- Thus,
assuming a Gaussian belief over
WL-1 = (Re(R(L-1),, Im Im(R((-1)), w((-1))
induces (multi-output) GP over (L-11); 1):
z (b-1) GP (2 K)(1-11)
with
mg = 2L--1) (x, w)
Moreover, (4-1) is the sum of three (dependent) parametric Gaussian processes with feature functions dkj- Pkyj, and vb.. respectively. Consequently, the function-valued GP induced by the linearized FNO is given by
F(a)(x) a(mzic (x:)) Dÿ (mza (x(x)) (2(4-))) mziz 1 (x:)),
ie. F(a) GP (ma; Ka) with
ma(x) F(a, wn (z). and
K (2(1) Dî (m (mx(z-n((1)) K Kaa-n(an ( (21, Dz (m 1)
If the input function A is discretized grid X i
(Da))% we set 6/4 (14-1)
C and (x)
interpolates 0/ ((4-1)( (4)) (e.g. spline interpolation or Fourier interpolation).
(x
23

---

Lincarization Turns Neural Operators into Function- Valuud Gaussian Processes
D. Experimental details
D.I. Data: PDE trajectories
D.1.1. Low DATA GENERATION USING APEBENCH
To evaluate the performance of the uncertainty quantification methods discussed, we utilize the code in the APEBench for generating data from Burgers', Hyper Diffusion and Kuramoto-Sia -Sinnssi-iv equation (conservative) (cf. (Koehler et al., 2024) for more details). Table 3 summarizes the characteristics of the datasets we use, the number of trajectories for training, and testing, well the spatial and temporal resolutions.
PDE Name
Dimensions
Training Traj.
Valid. Traj.
Test Traj.
Spatial Res.
Temp. Res.
Burgers
250
250
256
59
Hyper Ditfusion
ID
25
250
250
256
59
Kuramoto- -Siashhiik (cons.)
ID
25
250
250
256
59
Table 3: Summary of PDE datasets generated using APEBench.
1.1.2. OUT-OF- OUT-OR-DISTRIBUTION DATA GENERATION
We generated additional datasets based on the addvection-diffusion-eceee equation
du Øt
+ v Vu = aV R,
(D.1)
to evaluate the robustness of the uncertainty quantification methods under out-of-distribution (OOD) scenarios. Here represents the scalar field (e.g.. concentration or temperature), v is the velocity field, a is the diffusion coefficient, and the reaction term, where 0.026 was held constant throughout the datasets. The following specifications guided this process:
Datasets: Trajectories for five variations of the advection-difffiion--reccoo equation were generated (cf. Figure 5):
1. Base
Random number (1-10) of Gaussian blobs as an initial condition.
Constant random velocity field.
No reaction terms.
Flip
Random number (1-10) of Gaussian blobs an initial condition.
Constant random velocity field, which is reversed at the center of the domain.
No reaction terms.
3. Pos
Random number (1-10) of Gaussian blobs initial condition.
Constant random velocity field.
Randomly placed triangular heat source.
4. Pos-Neg
Random number (1-10) of Gaussian blobs as an initial condition.
Constant random velocity field.
Randomly placed triangular heat source; randomly placed cloud cloud-shaped heat sink.
Pos-Neg-Flip
Random number (1-10) of Gaussian blobs an initial condition.
Constant random velocity field, which is reversed at the center of the domain.
Randomly placed triangular heat source; randomly placed cloud-shaped heat sink.
PDE Solver: We utilized a custom implementation to solve the advection advection-diffusionn diffusion-reactionn equation, relying on a 9-point stencil with Runge- Kutta 4 (RK4) with fine spatial and temporal resolutions.
24

---

Lincarization Turns Neural Operators into Function-Valuue Gaussian Processes
Simulation Parameters:
Spatial resolution: 100 x 100 grid.
Temporal resolution: Al = 5 x 10-10 over 200-time steps, from which 59 are sub-sampled.
We train a Fourier neural operator on 1000 training trajectories of the Base dataset and calibrate uncertainty quantification methods- if necessary- on a corresponding validation datasct with 250 input-output pairs. For the evaluation, we consider 250 random test pairs from some of the datasets. In Figure 5, we depict single trajectory for each of these datasets generated. To seamlessly evaluate the various OOD datasets, we pad the training set with constant zeros for the velocity field and reaction terms as placeholder.
D.2. Model and training
For all experiments, we consider the original Fourier neural operator architecture (Li et al., 2021) with the hyperparameter suggestions following (Koehler et al., 2024), i.e. 12 modes (per spatial dimension) and 18 hidden dimensions constant throughout the network, with total of 4 Fourier blocks. For the training, consider 10 initial time steps and train to predict the following time step. The velocity field and the reaction term are glued to the input. Furthermore, the input is always padded by two constant zero grid points to reduce artifacts at the borders. Networks for the low data experiment are trained for 100 epochs, all remaining networks are trained for 1000 spochs-whiee single input-output pair per trajectory the training set. During training the mean squared error loss was minimized using AdamW (Loshchilov & Hutter, 2019) combined with cosine decay learning rate scheduler with warmup. All training implementations rely on jax (Bradbury al., 2018), Flax NNX and optax.
one epoch corresponds to iterating through a
D.3. Uncertainty quantification methods
D.3.1. INPUT PERTURBATIONS
Input Perturbations involve augmenting input data with small, random perturbations to introduce diversity into the model's predictions. This approach exploits the model's sensitivity to input variations to approximate uncertainty in the output space. Following the approach of Pathak et al. (2022), ensemble predictions perturbed versions of single input u,, Perturbations are sampled as <x,t ~ (0, a²) for each input value of the discretized function state un (I, t). The parameter is calibrated to achieve accurate marginal uncertainty predictions.
generated by forwarding batch of pointwise
D.3.2. ENSEMBLES
Deep ensembles rely on training multiple independent instances of a model, each initialized with different random seeds and potentially trained on different subsets of the data (Lakshminarayanan et al., 2017). The diversity in the learned weight configurations leads to variety of predictions for given input, enabling the computation of both mean predictions and marginal uncertainty estimates. From Bayesian perspective these weights represent close to true samples from the weight space posterior distribution, that can only be found through the high cost of additional training runs.
D.3.3. ISOTROPIC GAUSSIAN (x (s-)
The isotropic Gaussian covariance structure represents the weight space uncertainty as:
N(w" Σ o
(D.2)
where x² is the variance parameter and I is the identity matrix, reflecting independence and identical uncertainty across all dimensions in the weight space. To calibrate the uncertainty, we tune the parameter o² as outlined below.
From a Bayesian perspective, seen in methods such Laplace, this can be viewed as just considering calibrated prior over the selected weight space.
D.3.4. LAPLACE APPROXIMATION (s-LA)
In the linearized setting, a natural choice for a posterior Gaussian belief in weight space is given by the linearized Laplace approximation, where the Hessian is given by the Generalized Bauss-Newton (GGN) HOON (Immer et al., 2021). For further details, see Equation (B.1).
25

---

Lincarization Turns Neural Operators into Function- Vanction-Valued Gaussian Processes
Initial condition
t= 20
[= 40
59
(end of traj.)
1
Base
Flip
Pos
Pos-Nog
Pos-Neg-Fiip
Figure Initial condition and three-time steps of a single trajectory per generated dataset (Base, Flip, Pos, Pos-Neg, Pos-Neg-Flip).
26

---

Lincarization Turns Neural Operators into Pnaction-Valucee Gaussian Processes
For high-dimensional parameter spaces, approximations of the GGN required. Common techniques include diagonal approximations (Daxberger al., 2021a), K-FAC (Martens & Grosse, 2015), low-rank approximations (Dangel et al., 2022). Since diagonal and K-FAC approximations do not capture correlations between weights acting on different Fourier modes, we focus on low-rank approximation of the GGN instead.
We extend the approach (Dangel et al., 2022) by selecting the largest eigenspaces of the GGN and placing an isotropic Gaussian prior over all weights instead approximating the posterior covariance directly, This allows regions of uncertainty to fall back the prior belief. For a fixed input function and evaluation grid, the push-forward with the low-rank approximation HGGN VV" yields the normal distribution:
().
(D.3)
where Σ- = (nVVT + -I)- n the number of onput- pairs used to train the neural operator, and Jg is the Jacobian of the model with respect to the weights. In all experiments, we consider a low rank of 500. For the low data regime, we consider all input-output pairs, while for the OOD experiment only minibatch of 1000 input-output pairs.
D.3.5. SAMPLE-*
For each weight-space covariance method and the input perturbation approach, we consider a sample-based pushforward which generates an ensemble of predictions in the output space. This aligns with the way most weight-space covariance methods are introduced in the literature (cf. Maddox et al., 2019), In our experiments, we generate 200 samples that are propagated through the network. The empirical mean and standard deviation are then estimated using the set of predictions the output space.
D.3. LUNO
The core implementation of LUNO leverages matrix-free Jacobian- products and the algebraic structure of the chosen Gaussian covariance matrix in weight space. Our framework supports fully lazy evaluations, enabling efficient sampling, marginal variance estimation, and matrix-vector products with the covariance of the output space.
When restricting the weight space uncertainty to only the last Fourier block, we can explicitly derive the matrix action the Jacobian of the inverse Fast Fourier Transform (IFFT). Additionally, we exploit the fact that the final linear layer is applied pointwise in the spatial domain, significantly improving computational efficiency and reducing memory requirements compared to traditional implementations of linearized pushforwards (Daxberger al., 2021a).
D.4. Evaluation
We use the following metrics to evaluate model performance on 250 onput-output pairs from the respective test trajectories. Here, yi denotes the ground truth, ýi represents the predicted mean, and a; is the predicted standard deviation for the i-th sample. All reported numbers are the expected value over the test samples.
D.4.1. ROOT MEAN SQUARED ERROR (RMSE)
The root mean squared error
1
RMSE
)
n
1=1
measures the average magnitude of the errors between ensemble mean/linearized mean and ground truth values. Lower RMSE indicates better predictive accuracy, with zero being the ideal value.
D.4.2. MARGINAL NEGATIVE C.EIIIEEEOOGOOO L.GGLIIEELIHOOO (NLL)
The marginal NLL quantifies how well the predictive distribution fits the data under the assumption of Gaussian uncertainty. Lower NLL values indicate better calibration the uncertainty estimates and higher likelihood of the observed data under the predictive model. It represents trade-off between lower variances and how much of the error is accurately captured by the uncertainty.
2 exp (-
1
(14 (_(ui
0.)
NLL = log
20
27

---

Lincarization Turns Neural Operators into Function-Valuue Valued Gaussian Processes
D.4.3. STATISSSI
The ²²-statistic measures the average squared error normalized by the predicted variance. A value close to 1 indicates well-calibrated uncertainty predictions.
84
9.)
71
0
Values above one suggest overconfidence, while values below one indicate underconfident predictive uncertainty. From UQ perspective underconfidence is better than overconfidence.
D.5. Calibration
All hyperparameters of the discussed UQ methods (mostly o² in the above definitions) were calibrated using 250 input-output pairs of the validation set to minimize the marginal negative log-likelihood. We calibrate each method's hyperparameters separately using grid search over logarithmically spaced grid with 500 points centered around the relevant value.
D.6. Additional results
D.6.1. Low DATA REGIME
Method
RMSE (4)
X²
NLL (4)
Input Perturbations
1.98 x 10
1.203
2.3927
Ensemble
1.95 x 10
12.674
1.8957
Sample-Iso
2.02 x 10-2
1.237
-2.4391
LUNO-ISo
1.99 x 10-2
1.155
-2.4677
Sample-LA
2.18 x 10
2.097
-2.2312
LUNO-LA
1.99 x 10-2
0.984
-2.5248
Table 4: Performance metrics comparison for UQ methods evaluated on an FNO trained on trajectories of the Hyper-
Diffusion equation.
Method
RMSE (L)
x²
NLL (1)
Input Perturbations
7.84 x 10-2
0.938
-1.0385
Ensembles
6.82 X 10 2 7.618
0.8489
Sample-Iso
7.88 x 10-2
1.211
-1.8662
LUNO-ISo
7.82 10-2
0.922
-1.0995
Sample-LA
1.46 10-1
2.758
-0.1699
LUNO-LA
7.82 x 10-2
1.058
-1.1653
Table 5: Performance metrics comparison for UQ methods evaluated on an FNO trained on 25 trajectories of the Kuramoto- Sivashinsky (conservative) equation.
28

---

Lincarization Turns Neural Operators into Function- Valued Gaussian Processes
Input Perturbations
Ensemble
Sample-LA
LUNO-LA
1.5
0.0
-0.15
(a) Kuramoto- Kuramoto-Sivashnns Sivashinsky (cons).
Input Perturbations
Ensemble
Sample-LA
LUNO-LA
0.5
0.04
0.0
-0.01
(b) Hyper-Diffusion.
Figure 6: FNO predictive uncertainty quantified by several different methods. Top row: target function (-), mean (-) and 1.96 standard deviations of, as well samples (-) from, the predictive belief. For the ensemble, the samples are four of the ensemble members. Bottom row: spread of the predictive distribution around the mean. For the sample-fensemble-based methods, we construct a Gaussian distribution from the empirical covariance matrix and draw four samples (-). We plot 1.96 standard deviations of the predictive belief, as well as the top-three eigenfunctions (-) and a heatmap of the predictive covariance matrix (top right corner of panels).
29

---

Lincarization Turns Neural Operators into Function- Dunction-Valued Gaussian Processes
Method
Base
Flip
Pos
Pos-Neg
Pos-Neg-Flip
Input Perturbations
22.586
2.573
-1.174
69.346
494.935
Ensemble
-5.313
3.825
-4.802
-1.257
-1.014
Sample-Iso
-2.921
4.071
-1.226
9.457
43.362
LUNO-ISo
-2.89
3.450
-1.260
7.636
37.733
Sample-LA
-2576
4.395
-1.183
7.369
27.046
LUNO-LA
-2.344
-1.126
-1.742
-0.888
1.164
Table 6: Negative Likeliiod (NLL) across different OOD datasets. Lower is better.
D.6.2. OCT-OFF-DISTRIBUTII
Method
RMSE (4)
x²
NL. (1)
Input Perturbations
1.46 x 10
1.497
2.5861
Ensemble
1.51 X 10-3
0.162
-5.3134
Sample-Iso
1.44 x 10-2 2
1.217
-.9214
LUNO-ISo
1.46 x 10-2
1.189
-2.9919
Sample-Iso
1.73 10-2
1.910
-2.7556
LUNO-LA
1.46 10 0.841
-2.9340
Table 7: Metrics evaluated OOD dataset Base.
Method
RMSE (4)
x²
NLL (4)
Input Perturbations
4.36 x 10-2
12.255
2.5727
Enselble-Sample
4.43 10-3 0.215
3.8249
Sample-Iso
4.37 x 10-2 15.321
4.0715
LUNO-Iso
4.36 x 10-2 13.983
3.4502
Sample-LA
4.62 x 10-2 15.868
4.3953
LUNO-LA
4.36 x 10-2
3.475
-1.1257
Table 8: Metrics evaluated on OOD dataset Flip.
Method
RMSE (4)
x²
NLL (L)
Input Perturbations
2.94 x 10-2
4.292
-1.1744
Ensemble
4.93 X 10-3
0.549
-4.8023
Sample-Iso
2.91 x 10
4.552
-1.2261
LUNO-ISO
2.94 x 10-2
4.408
-1.2596
Sample-LA
2.79 10-2
4.622
-1.1835
LUNO-LA
2.94
10 3.001
-1.7416
Table 9: Metrics evaluated OOD dataset Pos.
30

---

Lincarization Turns Neural Operators into Function- Valued Gaussian Processes
Method
RMSE (4)
x²
(L)
Input Perturbations
5.55 x 10
145.561
69.3458
Ensemble
9.79 10-2
1.161
-1.2569
Sample-Iso
5.55 10
25.961
9.4566
LUNO-ISo
5.52 X 10
22.231
7.6355
Sample-LA
5.94 X 10
21.685
7.3688
LUNO-LA
5.52 X 10-2
4.182
-0.8180
Table 10: Metrics evaluated on OOD dataset Pos-Neg.
Method
RMSE (L)
x2
NLL (L)
Input Perturbations
1.10 X 10-1
997.143
494.9350
Ensemble
1.39 x 10-1
1.006
-1.0140
Sample-Iso
1.10 x 10-1
93.712
43.3620
Prior-Iso
1.10 X 10-1
82.367
37.7331
Sample-LA
1.15 10-1
60.728
27.0460
LUNO-LA
1.10 10-1
7.127
1.1642
Table 11: Metries evaluated on OOD dataset Pos-Neg-Flip
D.6.3. EVALUATION OF SINGLE TRAJECTORY
Run-Time comparison
Method
LUNO~
Sample-
Input Perturbations
10.19 - 0.006
Ensemble
0.85 + 0.004
Iso
0.53 = 0.004
14.00 0.008
*-
5.75 0.017
27.70 0.006
Table 12: Comparison of run-time performance between sampling-basee based (Sample-*) and linearization based (LUNO-+) methods across different uncertainty quantification techniques when out-rolling single trajectory iteratively.
31

---

Lincarization Turns Neural Operators into Function- Dunction-Valued Gaussian Processes
RMSE
NLL
10
10
0.4
10a
10a
0.3
Base
10
0.2
10a
0.1
101
10
100
109
0.0
109
0.4
10
Flip
0.3
10a
10a
10a
10a
101
101
109
0.0
10a
10a
0.4
10a
10
0.3
10a
Pos
0.2
10a
0.1
10
10
10a
109
0.0
109
0.4
104
10
Pos-Neg
0.3
10a
10a
10a
10
10a
10
0.0
10a
10a
Pos- Pos-Neg-Fii Flip
10a
10
0.3
10a
0.2
10a
0.1
101
10
100
10a
0.0
50
Time step
Time step
Time step
Figure 7: Averaged performance of different UQ methods autoregressive rollout of the FNO 50 trajectories from the Base, Flip, Pos, Pos-Neg, and Pos-Neg- Flip datasets. We compare input perturbations (-), deep ensembles (-). Sample-Iso (-), LUNO-IsO (-), Sample-LA (-), LUNO-LA (-).
32


Reviews:
[{"review_id": "ZQiqsOs9mO", "reviewer": "Reviewer_N2Br", "summary": "The paper proposes a novel framework for approximate Bayesian uncertainty quantification in trained neural operators. The approach relies on model linearisation and pishes weight-space uncertainty to neural operators' predictions. This allows the application of Bayesian deep learning methods, such as linearised Laplace approximation, to neural operators.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your positive review of our paper!\nWe are glad to hear that you found the theoretical and experimental aspects clear and well-structured.\nShould anything still remain unclear, we are happy to clarify during the discussion period.\n\n---\n\n**Publication of code:**\nWe intend to release the code for our experiments upon acceptance.\nTo preserve anonymity during the double-blind review process, we have refrained from sharing it at this stage due to dependencies and references that could potentially reveal our identity.\nThat said, we are fully committed to ensuring reproducibility and will provide the complete codebase.\nIf access to the code would be helpful during the discussion period, we would be happy to invest additional time to clean and anonymize it appropriately.", "claims_and_evidence": "Claims: \n- A novel framework (LUNO) which provides linearized predictive uncertainty in neural operators;\n- Interpretation of LUNO as a probabilistic generalization of the concept of *currying* in functional programming;\n- Compatibility of LUNO with established methods for quantifying weight-space uncertainty in deep neural networks, including the Laplace approximation. \n- LUNO scales to large models and datasets and, like neural operators, is inherently resolutionagnostic.\n\nEvidence: \nProofs and case study on Fourier neural operators. The claims are well supported by the evidence.", "methods_and_evaluation_criteria": "Yes, the evaluation criteria are appropriate.", "theoretical_claims": "I went over proofs until Corollary A.13. I did not spot any issues and could generally follow the flow of the proofs. However, since this is not my main expertise, I could have missed something. Nevertheless, the way the theoretical arguments are structured is natural and consequential.", "experimental_designs_or_analyses": "Yes, the experimental design is suitable and experiments provide detailed evaluation of the proposed approach.", "relation_to_broader_scientific_literature": "The paper is in the domain of neural operators and generally touches on uncertainty quantification and utilization of Gaussian processes.", "essential_references_not_discussed": "None are specifically missed as far as I can tell.", "other_strengths_and_weaknesses": "Strength:\n- The paper is very clearly structured and motivated. \n- The experiments are detailed and well-described. \n- The limitations are discussed and acknowledged. \n\nWeaknesses: \n- No code is available for reproducibility of the experiments as far as I can tell.", "other_comments_or_suggestions": "No other comments.", "questions_for_authors": "Are the authors planning to release the code with the paper in the future?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Zmjs0OluVM", "reviewer": "Reviewer_49mR", "summary": "This paper introduces LUNO, a linearization approach for turning a nonlinear neural operator into a Gaussian random operator, thereby providing uncertainty estimates for operator learning. This is important in areas such as safety-critical prediction and out-of-distribution scenarios. The method is compared against baselines, particularly deep ensembles, in a low-data regime and an OOD setting, showing superior and/or more principled uncertainty quantification.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive evaluation of our paper and for pointing out areas where we can further clarify the presentation.\n\n---\n\n**Applicability beyond FNO:**\nWhile our experiments focus on Fourier Neural Operators (FNOs), the theoretical framework applies in principle to any neural operator.\nWe chose FNOs because of their popularity and the particularly efficient lazy representation they enable for the function-valued posterior process.\nHere, we achieved significant computational speedups by leveraging analytic properties of the inverse fast Fourier transform.\n\n---\n\n**Clarification on sample-based approaches:**\nWe do not claim that sample-based approaches directly yield Gaussian processes, as noted in Section 5:\n\n> We evaluate linearized predictive uncertainty (LUNO-\\*) against sample-based approaches (Sample-\\*), which require additional approximations to impose a Gaussian Process structure over the output space.\n\nOur linearization approach, by contrast, yields an analytic function-valued GP formulation over the output space, which we believe is a key advantage in operator learning tasks.\nFor evaluation purposes we consider sample-based approaches that amount to positing a moment-matched GP using the Monte Carlo estimates of the mean and covariance function in the output space.\nIn our experience, this also makes the calibration of SAMPLE-based approaches computationally more expansive.\nWe apologize for the vague formulation and will revise the text to make this distinction clearer.\n\n---\n\n**Typos and notational clarifications:**\nWe appreciate your thoroughness in pointing out typos, notational improvements, and figure label inconsistencies. This helps a lot and we will make the necessary corrections in the revised manuscript.\n\n---\n\nThank you again for your constructive feedback. Should anything remain unclear, we are happy to clarify during the discussion period.", "claims_and_evidence": "The authors claim to construct a function-valued Gaussian process for uncertainty quantification, linking it to probabilistic currying; this is explained convincingly in the paper.", "methods_and_evaluation_criteria": "The proposed methods address a low-data regime and an OOD setting. These support the authors' claim that the proposed methods contribute to uncertainty quantification in such domains.\n\nOn the other hand, the authors only evaluated their framework on FNO. Since their approach depends on linearization (in terms of network weights), the accuracy of the linear approximation could depend on the network architecture. Therefore, it makes sense to include an experiment evaluating other types of neural operators as well.", "theoretical_claims": "I have not checked the correctness of the proofs.", "experimental_designs_or_analyses": "Yes, I checked all which is included in the main text. See above comments about other neural operator architectures.\n\nAdditionally, in line 328 left column, the authors are not clear about what ‘sample-based’ approaches are, so they should include a brief description in the main text. In the appendix, it seems like the authors sample from the distribution of weights and pushforward to the output. But this does not readily give a Gaussian process, unlike what the authors claim afterwards. This could be made clearer.", "relation_to_broader_scientific_literature": "The paper provides an additional tool in the so-far ill-equipped toolkit for uncertainty quantification for operator learning. As discussed by the authors, current methods do not extend to Fourier neural operators.\n\nThe method proposed builds on existing work on estimating the posterior distribution of neural network weights.", "essential_references_not_discussed": "I have not identified any missing references, either in Gaussian processes in a neural operator setting or in linearized Laplace approximation.", "other_strengths_and_weaknesses": "This paper is very well written and clear, and the experiments are thoughtful. It draws a link elegantly with the concept of currying in functional programming, connecting the two fields. The linearization approach is developed well and rigorously, and although i have not checked the proofs, the claims are sensible and consistent with existing intuitions about Gaussian processes.", "other_comments_or_suggestions": "Line 120 left column: ‘This justifies interpreting random processes as probability measures….’, I found this line not so clear - a probability measure is a set function which maps to [0,1] which the random processes spoken of are not. Perhaps better to replace as ‘random variables’. \n\nLine 131 left column: ‘A d’-output Gaussian process is a random process ….for all n\\in \\mathbb N and a_1,…,a_n \\in \\mathbb A.’ This reads a bit weird, clearer to say ‘A d’-output Gaussian process f on \\mathbb A \\times \\Omega’ which allows you to refer to the same f later on in the sentence.\n\nLine 134 right column: The F’s in the tuple should be bolded, it is also a bit unclear without further explanations how a function could be seen as jointly Gaussian. Worth including an explanation or point to the appendix somewhere.\n\nLine 244 left column: `f(a, \\mu)(x)` is a typo, this should be `f((a,x), \\mu)`. \n\nLine 251 left column: upright f is really just `f^{lin}_\\mu((a,x), \\cdot)`, so it’s a redefinition to simplify notation. Making it clear would improve the reading experience. \n\nFigure 3 needs to be fixed, as there is no longer 3 rows as referred to in the caption. There are also two places in the paper which referred to panel 8, which is the null space projection in Figure 3 - this needs to be fixed too.", "questions_for_authors": "None.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nrDaggO84D", "reviewer": "Reviewer_P4r1", "summary": "The paper introduces LUNO, a framework for uncertainty quantification in neural operators using function-valued Gaussian processes. By leveraging linearization, the method propagates Gaussian weight-space uncertainty to the operator’s predictions, effectively converting trained neural operators into Gaussian random processes. LUNO provides a _post-hoc_, scalable, and resolution-agnostic Bayesian uncertainty estimation approach, evaluated primarily on Fourier neural operators (FNOs) for solving PDEs.\n\n_(for any missing input on any of the fields, please refer to the **Strengths and Weaknesses** or the **Other comments or suggestions** sections)_", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed and positive assessment of our work!\n\n---\n\n**On the assumption of Gaussian weight uncertainty:**\nThis is an important point.\nOur theoretical framework (see Step 3 of Section 3.2, and Appendix A, particularly Corollary A.14 and Theorem 3.2, and Section A.4) indeed apply to more general (non-Gaussian) weight-space distributions.\nFor instance, the predictive mean and covariance expressions remain valid under arbitrary (non-Gaussian) weight-space beliefs, as long as the respective moments of the weight-space belief exist.\nWe focus on Gaussian distributions because they yield a (closed-form) Gaussian process over the output space, enabling well-studied analytic tools (e.g., conditioning, and Bayesian experimental design).\nHowever, mixture models or other non-Gaussian distributions are an exciting direction for empirical evaluation in future work.\nWe will clarify this more explicitly in the revised text.\n\n---\n\n**Effect of linearization:**\nWe want to emphasize that we only linearize the neural operator in the weights and not with respect to its input.\nThe linearized neural operator $F^\\text{lin}(u, w)$ is linear in $w$, but still highly nonlinear in $u$.\nFurthermore, exactly as in Ortega et al. (2024), the predictions of the original model are not altered, only extended by a covariance which depends on the chosen weight-space uncertainty.\n\n---\n\n**Applicability beyond Fourier neural operators:**\nWhile our experiments focus on Fourier Neural Operators (FNOs), the theoretical framework applies in principle to any neural operator, including transformer-based architectures.\nWe chose FNOs because of their popularity and the particularly efficient lazy representation they enable for the function-valued posterior process.\nOur implementation exploits the structure of the inverse Fourier transform for computational efficiency, which should in principle extend to other signal transforms such as Spherical FNOs.\n\n---\n\n**On classification tasks:**\nThis is an interesting point.\nNeural operators can naturally be applied to pointwise classification tasks such as semantic segmentation in computer vision.\nIt is possible to generalize our methodology to classification tasks, in which the output GP is transformed through a link function, yielding a generalized linear model.\nComputing the pushforward distribution of the output GP through the link function is an active field of research and introduces additional approximations.\nHowever, a central motivation for this work is the empirical and theoretical study of uncertainty quantification for deep-learning based emulators for scientific simulations, which virtually always focuses on the high-dimensional regression setting.\n\n---\n\n**Other methods for weight space uncertainty:**\nThank you for your suggestion.\nOur theoretical framework allows for arbitrary (non-Gaussian) weight-space beliefs, so one could in principle use a normalizing flow to represent the weight-space uncertainty.\nWe mostly focus on the Laplace approximation in our experimental analyses because it is cheap to compute and applicable to pretrained models, which makes it particularly promising for large models like neural operators.\nA strength of our framework is that the choice of weight-space uncertainty structure is kept flexible.\nWe appreciate the reviewer's suggestion and agree that exploring richer weight-space beliefs is an interesting avenue for future research.\n\n---\n\n**Compression of the shifted network parameters:**\nIn most practical cases, the linearization point $w_0$ (see Appendix A.4, Step 2), is chosen to be equal to the mean $\\mu$ of the weight-space belief (e.g., this is set to the weights $w^\\star$ of the trained network in case of a Laplace approximation).\nWe did formulate Appendix A.4 under an arbitrary linearization point to allow for more flexibility in the theoretical framework, but did not end up using this flexibility in practice.\n\n---\n\nWe hope these clarifications address your questions.\nThank you again for your time and constructive feedback.\nWe also thank you for the additional comments and suggestions regarding streamlining the first three pages and clarifying a few passages.\nWe will include these corrections and references.", "claims_and_evidence": "- Neural operators lack inherent uncertainty quantification, which limits their reliability in high-stakes applications.\n  - *Evidence:* Prior work has focused on deterministic operator learning without probabilistic guarantees.\n- Model linearization enables efficient uncertainty propagation without retraining.\n  - *Evidence:* Theoretical derivations show that linearization allows weight-space uncertainty to be pushed forward to the output function space.\n- The resulting function-valued Gaussian process belief provides structured uncertainty estimates.\n  - *Evidence:* The paper introduces a rigorous connection between function-valued Gaussian processes and Bayesian deep learning techniques.", "methods_and_evaluation_criteria": "The framework linearizes the trained neural operator around the mean of a Gaussian weight belief. This is interpreted as a probabilistic generalization of currying, leading to a function-valued Gaussian process. The authors evaluate LUNO on Fourier neural operators (FNOs) applied to PDE problems with multiple metrics, including RMSE, NLL and $\\chi^2$ statistics.", "theoretical_claims": "Claims for the submission:\n- LUNO constructs function-valued Gaussian processes from neural operators by treating them as infinite-dimensional stochastic processes.\n- Probabilistic currying is introduced as a key concept, establishing equivalence with multi-output Gaussian processes.\n- Gaussian weight-space uncertainty can be efficiently propagated via linearization, providing a computationally tractable Bayesian formulation.", "experimental_designs_or_analyses": "- Evaluates low-data regimes (small training sets) and out-of-distribution (OOD) generalization.\n- Compares LUNO to:\n  - Sample-based approaches (Monte Carlo sampling of weight posteriors).\n  - Input perturbations (simulating uncertainty via randomized inputs).\n  - Deep ensembles (training multiple independent models).\n- Experiments on Burgers' equation and Advection-Diffusion PDEs, analyzing predictive uncertainty under domain shifts.", "relation_to_broader_scientific_literature": "The paper builds on existing work in neural operators and Fourier neural operators (FNOs), connecting with research on Bayesian deep learning techniques such as Laplace approximation, variational inference, and SWAG. It relates to operator learning in function spaces, particularly focusing on function-valued Gaussian processes.", "essential_references_not_discussed": "Although not exactly my field of expertise, the authors seem to cover extensively the literature relevant for the paper. As a mere suggestion, maybe it is worth discussing Wasserstein Gaussian processes as an alternative probabilistic operator learning approach, as well as maybe other works related to implicit stochastic processes. Moreover, related to the post-hoc approaches, a comment on the effects of the linearization approximation on the uncertainty estimates could be insightful, especially in comparison to other methods (e.g. [1]) that do not alter the predictive mean of the original model.\n\n[1] Ortega, L. A. et al. (2024, July). Variational Linearized Laplace Approximation for Bayesian Deep Learning. In International Conference on Machine Learning (pp. 38815-38836). PMLR.", "other_strengths_and_weaknesses": "### Strengths:\n- Post-hoc application avoids retraining, making it computationally efficient. This is a specially interesting feature for large-scale applications.\n- The presented approach is strongly theoretically grounded, with rigorous uncertainty propagation.\n- The method seems to be scalable, working well with high-dimensional data and large models.\n- While a bit dense at times, the paper is well written and structured, making it more accessible.\n\n### Weaknesses:\n- The method assumes Gaussian weight uncertainty. This hypothesis could limit the model's applicability to highly non-Gaussian posteriors.\n- Linearization may introduce approximation errors, particularly for highly non-linear operators. This is a key point that should be adressed further, since one of the main motivations of post-hoc approaches is to avoid retraining and keeping the original model's performance.\n- Evaluation focuses on PDE benchmarks, broader validation on other function learning tasks would broaden the method's applicability.", "other_comments_or_suggestions": "- A lot of the text in the first 3+ pages of the article is devoted to mostly explaining basic concepts and results from previous work. Given the amount of work relegated to the supplementary material, it might be beneficial to streamline these initial sections to focus more on the novel contributions and include more details in the main text. Along the same lines, certain discussions like 5.2 seem more appropriate for the supplementary material.\n- The wording of Lemma 4.2 could be improved for clarity.\n- Since this work explores the connections between gradient descent and Bayesian inference, it may also be interesting to tackle the \"main research question\" more from the Bayesian perspective. In particular, it would be interesting to see if there is any complete Bayesian approach to this problem, where the complete loss function is derived from a probabilistic formulation of the problem.\n- The authors should extend the experimental part to apply their proposed approach to more complex datasets beyond toy regression tasks.\n- Explore whether normalizing flows or implicit stochastic processes (or maybe just implicit distributions) could further enhance flexibility beyond Gaussian assumptions.\n- Investigate how the method performs under stochastic gradient descent (SGD) instead of full-batch gradient descent.", "questions_for_authors": "1. How does this approach perform on classification tasks where aleatoric noise is often modeled differently?\n2. Could the method be extended to convolutional or transformer-based architectures?\n3. Could normalizing flows be used in this context to extend the method to more complex uncertainty representations?\n4. Is there any approach which would allow for a more efficient storage of $\\theta_0$ so that the shifted network method could be more practical in memory-constrained settings? Maybe some sort of approximation or compression technique?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3q9o0oKz3L", "reviewer": "Reviewer_wc7q", "summary": "The authors propose a new approach to approximate stochastic neural networks with Gaussian weights into Gaussian processes (GPs). The approach is based on performing a linearization around the mean of the weights to obtain a GP approximation of the network. The effectiveness of the framework is shown in some examples where neural networks are used to learn PDEs.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive feedback and for highlighting several interesting directions for further exploration.\nBelow are our responses to your main points:\n\n---\n\n**Comparison with infinite-width GP approaches:**\nWe agree that investigating how infinite-width GP methods could be adapted to the infinite-dimensional operator learning setting is an intriguing theoretical direction.\nHowever, our primary goal is to develop a practical, yet theoretically sound, post-hoc method that could be applied to practical (i.e., non-asymptotic) architectures.\nTo our knowledge, infinite-limit GP formulations typically involve retraining networks as the hidden dimension grows, which seems somewhat orthogonal to our focus on uncertainty quantification that leaves the architecture untouched.\nIf there is a specific experiment you have in mind, we would be happy to discuss it further during the discussion period.\n\n---\n\n**Varying weight distributions:**\nWe appreciate the suggestion to explore different variances.\nIn practice, the choice of weight-space uncertainty might often depend on the downstream application an end-user has in mind.\nFor instance, with a Laplace approximation and our low-rank-plus-diagonal covariance structure, the approximation defaults to the prior variance far away from the training data.\nThe diagonal portion of our covariance is an isotropic Gaussian prior, so larger variances tend to dilute the impact of the low-rank factors (which encode data-informed structure).\nWe tune the prior variance with respect to a validation set, but we agree that experimenting with more informative priors (e.g. non-isotropic or even non-diagonal) could be a valuable direction for future research, also in general for Bayesian deep learning.\n\n---\n\n**Non-diagonal weight distributions:**\nAs mentioned above, we use a low-rank-plus-diagonal structure (which leads to highly non-diagonal weight covariances in practice), but we acknowledge that a more general covariance could capture richer parameter correlations.\nIn our experiments, increasing the rank of the covariance did not yield substantial improvements, although we suspect that combining a higher-rank structure with more informative priors may prove beneficial, both of which are supported by our framework.\n\n---\n\n**Relevant references and related work:**\nThank you for pointing out additional references.\nWe agree that Khan et al. (2019) is an important relevant prior work on Laplace approximations.\nWe will move it to the main text in the revised version and also add Adams et al. (2024) as an interesting tangential work (where a GMM is constructed layerwise using their activations) in the related work section.\n\n---\n\n**Lack of error bounds and robustness guarantees for Bayesian inference:**\nWe agree that deriving formal error bounds is a challenging yet important open problem in the entire field of Bayesian deep learning, valuable especially for applications such as operator learning.\nTo our knowledge, the field does not currently offer rigorous, general bounds for weight space posteriors, which seems to be important prerequisite for establishing corresponding bounds for our setting.\nIn the revised version, we will consider acknowledging this limitation and highlighting the need for theoretical guarantees on the quality of uncertainty estimates.\nWhile the robustness guarantees in Cardelli et al. (2019) are indeed an interesting starting point, it is unclear how to verify the assumptions of their work for the neural tangent kernels used in this work.\nSpecifically, the computation of the supremum in point 3 of the \"Constant Computation\" section seems intractable or at least highly nontrivial in our case.\nMoreover, the paper seems to assume that the input space is (a subset of) $\\mathbb{R}^d$ whereas, in our setting, the input space is an infinite-dimensional function space.\nNonetheless, we believe that the framework of function-valued Gaussian processes provides a uniquely viable starting point for evaluating and interpreting predictive uncertainties in Bayesian deep learning.\n\n---\n\nWe hope these clarifications address your questions.\nWe look forward to continuing the discussion and potentially design an additional experiment with your support that can further refine our work.", "claims_and_evidence": "I found the claims well supported by mathematical evidence.", "methods_and_evaluation_criteria": "- An aspect that could be improved is the empirical evaluation. In fact, as, by definition, a linearization is a local approximation, it would be important to have experiments with different weight distributions (with difference variances) to see how the accuracy of the approximation changes with increasing the uncertainty of the distribution. Another interesting experiment would be to compare the approximation obtained with the proposed method with that of the limiting GP obtained by relying on the central theory limit [Lee, Jaehoon, et al. \"Deep Neural Networks as Gaussian Processes.\" International Conference on Learning Representations. 2018.], which is often used to study neural network posteriors [Cardelli, Luca, et al. \"Robustness guarantees for Bayesian inference with Gaussian processes.\" Proceedings of the AAAI conference on artificial intelligence. Vol. 33. No. 01. 2019].\n\n- Minor point: in the future, It would also be interesting to consider non-diagonal distributions on the weights.", "theoretical_claims": "Theoretical claims appear sound to me.", "experimental_designs_or_analyses": "Experimental analysis is sound. However, as mentioned above, it would be good to extend it a bit more.", "relation_to_broader_scientific_literature": "The paper extends recent works that approximate a neural network with a Gaussian process. It does so by defining a Gaussian distribution on a functional space. This allows for the evaluation of uncertainty not only in a finite set of input points, such as in [Adams, Steven, et al. \"Finite Neural Networks as Mixtures of Gaussian Processes: From Provable Error Bounds to Prior Selection.\" arXiv preprint arXiv:2407.18707 (2024).], and not only in the limit of infinitely width size, such as for the papers mentioned above. Of course, the approximation comes with no error bound, which I believe would be a very interesting and promising direction.", "essential_references_not_discussed": "The references I mentioned above are not referred to in the paper, and it would be good to include those. Furthermore, a reference that should be discussed  in the main text, and not only in the Supplementary, is [Khan, Mohammad Emtiyaz E., et al. \"Approximate inference turns deep networks into Gaussian processes.\" Advances in Neural Information Processing Systems 32 (2019], where the authors also seem to rely on similar techniques.", "other_strengths_and_weaknesses": "I believe that the key strength of the paper is to offer a new way to study the uncertainty of neural networks. Of course, the methods come with the weaknesses described above, but overall, I found the paper to be an interesting contribution.", "other_comments_or_suggestions": "No other comments.", "questions_for_authors": "My questions have been already formulated above and in summary, the more important one, apart from improving the related works and the discussion of existing methods, is:\n\n- Can you please extend the experiments as suggested above?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: Four knowledgeable reviewers recommend Accept and I recommend Accept. Reviewer P4r1 had some important comments on clarifications and limitations, which I strongly encourage the authors to include in the Camera-Ready version of the manuscript.
Decision: Accept (spotlight poster)