Approximate Bayesian Neural Operators: Uncertainty Quantification for Parametric PDEs
Authors: Emilia Magnani
Nicholas Krämer
Runa Eschenhagen
Lorenzo Rosasco
Philipp Hennig
Venue: Accepted by TMLR ()

Abstract: Neural operators are a type of deep architecture that learns to solve (i.e. learns the nonlinear solution operator of) partial differential equations (PDEs). The current state of the art for these models does not provide explicit uncertainty quantification. This is arguably even more of a problem for this kind of tasks than elsewhere in machine learning, because the dynamical systems typically described by PDEs often exhibit subtle, multiscale structure that makes errors hard to spot by humans. In this work, we first provide a mathematically detailed Bayesian formulation of the ``shallow'' (linear) version of neural operators in the formalism of Gaussian processes. We then extend this analytic treatment to general deep neural operators—specifically, graph neural operators—using approximate methods from Bayesian deep learning, enabling them to incorporate uncertainty quantification. As a result, our approach is able to identify cases, and provide structured uncertainty estimates, where the neural operator fails to predict well.

Published in Transactions on Machine Learning Research (9/2025)
Approximate Bayesian Neural Operators: Uncertainty Quantification for Parametric PDEs
Emilia Magnani
emilia. wwwali/mmmmnonootoeeeueeeennn.nn
Tübingen AI Center, University of Tübingen
Nicholas Krämer'
pekra&ddu.. dk
Technical University of Denmark
Runa Eschenhagen
re393 uk
University of Cambridge
Lorenzo Rosasco
lrosasco@mit. edu
MaLGa DIBRIS, University of Genova, Isbituto Hakiano di Tecnologia
Philipp Hennig
philipp. heigiigpeeeegtnnnnitteiiiinggn.
Tübingen Al Center, University of Tübingen
Reviewed on OpenReview: https: opennevieew. net/ forum? id= 6WAIkYMMA8
Abstract
Neural operators are type of deep architecture that learns to solve (i.e. learns the nonlinear solution operator of) partial differential equations (PDEs). The current state the art for these models does not provide explicit uncertainty quantification. This is arguably even more of problem for this kind of tasks than elsewhere in machine learning, because the dynamical systems typically described by PDEs often exhibit subtle, multiscale structure that makes errors hard to spot by humans. In this work, we first provide mathematically detailed Bayesian formulation of the "shallow" (linear) version of neural operators in the formalism of Gaussian processes. We then extend this analytic treatment to general deep neural operators specifically,, graph neural operators-using approximate methods from Bayesian deep learning, enabling them to incorporate uncertainty quantification. As a result, our approach is able identify cases, and provide structured uncertainty estimates, where the neural operator fails to predict well.
1 Introduction
Neural operators KKvachki et al. 2023 ct al. 2020b 2021a 2020a 2021bb are deep learning architectures designed for reconstruction problems related to partial differential equations (PDEs). They approximate mappings between infiiteediimensonall vector spaces of functions, such that once trained solutions of entire families of parametric PDEs can be represented by single neural network. However, the learning process is subject to several sources of uncertainty, which can result in a potentially significant prediction error because of the nonlinear and often nonintuitive interactions of different stages of the approximation. The goal of this paper is to develop methods for estimating this error at a practically acceptable computational cost. This kind of functionality is urgently nceded in this domain: Duc to the intricate and often not intuitive nature the dynamical systems described by PDEs, can be hard for the human eye to detect prediction errors, even when they are large.
In this paper, address this gap by developing approximate Bayesian framework for neural operators from theoretical, and computational point of view. We begin with brief review of neural operators.
Work mostly dane the University of Tübiingen.

---

Published in Transactions on Machine Learning Research (9/2025)
Right hand side function
Àg 3.0
Àg 4.5
7.5
2.3
0.10
0.0
-0.15
0.0
0.02
0.02
PDE solution
11
0.00
-0.3
0.00
6.004
0.005
-0.005
0
Figure 1: Green's functions in Equation 6 for different values of do 4.5, 7.5). On the left, right- hand-side functions g for the PDE in Equation 5 and respective solutions uy, Ug for the correspondent Ag-value, computed through Equation
Then, using linear, parametric PDEs as guiding examples, we show how their "shallow" (single-layer) base case allows for an analytic Bayesian treatment using the formalism of Gaussian processes (Rasmussen Williams (2006)). This linear case, while primarily of theoretical interest, provides valuable insights and aims to make this model class more accessible to the Bayesian machine learning community. We then extend the theoretical analysis to the nonlinear deep case. Here, analytic treatments we fall back on approximations developed for Bayesian deep Icarning. Specifically, we focus on Laplace approximations [MacKay] 1992) which are easy to add post-hoe even to pretrained networks, and add only moderate computational cost relative to deep training without uncertainty quantification (Daxberger et al. 2021). Our experiments Scction demonstrate that the resulting method effectively captures structure the predictive error of graph neural operators, both in the over- and ander-sampled regime. In Section we discuss some theoretical background, and develop probabilistie framework for neural operators in Section We discuss related work in Section
longer possible,
2 Background
In this section, examine how neural operators approximate solution operators for paramctric PDEs through functional observations. If we fix one input of the solution operator, neural operators be understood as effectively inverting the differential operator associated with the PDE. In this framework, the process of learning the operator becomes equivalent to reconstructing the Green's function, reducing the problem to task of function approximation. This perspective, developed in Section 2.1 forms the basis for the Bayesian approach developed Section 3.1 Subsequently, Section 2.2 we outline the iterative structure of neural operators, their training procedure, and how they relate Green's functions.
2.1 PDEs And Green's Function
One of the main fields applications of neural operators are familics of parametric PDEs of the form
(Cxx) fax) DED D
(1)
u(x) 0, E 0D
c Rd with boundary aD (c.g. open, bounded with Lipschitz boundary AD), where U D R F f: D R. X A, with U. F and A appropriate
for some sufficiently well-behaved bounded domain

---

Published in Transactions on Machine Learning Research (9/2025)
(1 = -1 (Wif(a)+ + /. 90, (a, y. ((0)),
()(F,)) oz (W)
gol ((()) 12)
(Wi.0) !)
( 02)
(Ws, 0s)
1
(W2, 02)
(2)
(xr
Figure 2: Neural operator architecture NOg. Each layer computes a new function th, that contains the neural network go in the integrand. Layer parameters are shown on the corresponding arrows.
function spaces. The precise nature of those function spaces is not important for the remainder of this work. The function À parametrises the differential operator Lx.
Equation defines a solution operator
H: AXFFFU (1,f)
(2)
in the sense that H(A, (z) (x) solves the PDE for the given functions X and f. Even though the PDE is linear, H is (possibly highly) nonlinear. The operator H map between function spaces. The idea behind neural operators is to approximate the operator H. with single neural network trained on function observations (fi, )= Thus, instead of approximating the solution of the PDE for only fixed f or X, neural operators directly infer the operator H. Numerically, the functions f and u are observed on a discretisation grid the function domains.
In this subsection we are interested in the particular case where is fixed, so the solution operator can be written
G: f u.
(3)
If the differential operator LA is linear, the map inherits that lincarity. Considering the operator in Equation (3) is an important step to understand the learning process of neural operators. In fact, observe how G is the inverse of the operator Σx. In this simplified case where À is fixed, the neural operator is therefore learning an operator, G, through function observations 1,,,,, that derive from the action of its inverse. In other words, during training, the neural operator is implicitly learning to invert the diffcrential operator Lx- In particular, the case where the differential operator is linear and admits a Green's function G, the solution of Equation can be expressed through integration with the kernel G
ux (x) = //, Gx(I, ()) dy.
(4)
Hence, learning the operator G is here equivalent to learn the function GX which means that an operator- learning task can be reduced to that of function-reconstruciion.
In the general analysis of linear PDEs (we refer to e.g. Evans (2010) for background on PDEs), the Green's 's function (x y) represents the impulse response of the linear operator Lx, that is Lx y) = 8(- - y) for y E D, where d denotes the Dirac delta distribution. Despite LA being linear, the Green's function itself can be nonlinear in either arguments. To visualize these concepts, we consider the ane-dimensional boundary value problem
( (-A - X Id) u(x) f(x), xE [0, 1].
(5)
u(O)) = u(1) = O,
3

---

Published in Transactions on Machine Learning Research (9/2025)
where in this case λo R is a (scalar) parameter, and Id is the identity operator. It's Green's function, for A # not Vn E N, is given by
A+B
(x, Ag sin siin(o0) (Ao)
(6)
where we abbreviated
A = Hy - 2)) sin (Aua)) sin (Ao - 4)) B = H y) sin (n - z)) sin (N03)
(7)
(8)
and H denotes the Heaviside step function. Details of the Green's function derivation are in Appendix
A
Equation (5) relates to Equation ! in the sense that the differential operator Exp = (-A AB Id) is
parametrised by Ag- Figure 1 shows examples of Green's 's functions Gxo for different values of Xo, along with solutions computed via Equation
2.2 Overview of Neural Operators
Neural operators are neural-network--bawe architectures designed to approximate the general solution opera- tor H defined Equation (2). Before introducing our Bayesian framework, we bricfly review their structure.
A more thorough explanation of what follows can be found in the work by Kovachki et al. (2023 Li et al.
2020b
2021a
2020a 2021bh.
Let go: D x D x X + R be a neural network with parameters 0. Define the neural operator NOe as composition of L E N layers
NOe: A x F + U,
(A. + (v1, o o (), f),
o...
(9)
where each layer
: e= (L,
(10)
is defined as composition of (i) integrating the output of the previous layer against ge,, and (ii) combining the integral with a linear component and an activation function o,
- (Weh(a)+ + / g0,
(e(h)(x)= = o
90. (z, y, A(x), A()))(y))
(11)
The space & in Equation 10 is a vector space of vealevalued functions on D, and the final layer of the neural operator maps into U, so vr. + - U. In Equation (1]). We a learnable linear operator (represented by matrix after discretization), and ge, is the integral kernel in the b-th layer. In practice, the integral cannot be computed in closed-form and suitable quadrature formula needs to be employed (which turns the integral into weighted sum of evaluations of the integrand; see e.g. Davis & Rabinowitz (2007) ). The parameter set θ of NOe is θ = (0, WW) Loosely speaking, one can think of this construction as deep neural network (NOg) that iteratively approximates the solution ux,5 with linear transformations We and nonlinear activation functions o, and at every iteration (layer) employs another neural network (go,). For a visualisation of NOg see Figure 2 Although the figure shows entering in each layer, practice the kernel may encode in an initial "lifting" layer while using final "projection" layer to map the function output back to the physical domain.
Note how NOg approximates operator. While, technically speaking, this means that its training and test set consist of functions, in the mumerical computation, these functions need to be observed on some grid. Nonetheless, neural operators are resolution-agnostic: their architecture does not depend on a particular discretization grid, and they can be applied to different resolutions without retraining. Let {AI...,,AN x (fi,.., f be set of training inputs, each of which shall be observed on some mesh X == { FK) ). In total, that makes NK MK 2MK training inputs. Without loss of generality, and for the sake of simple notation, assume that the solution of the PDE and the respective inputs observed the same
4

---

Published in Transactions on Machine Learning Research (9/2025)
mesh X. Thus, observe NM solutions U11, UNM, i.c. NMK training outputs one set of evaluations at X for each solution Unm associated with (An;Sm. n 1,.., N, m = 1,..,,M M. Each of these outputs is a function that maps from D to R, thus Unm (x) E RK The relation between inputs and outputs is
Unm = H(Am: Jm As (A))
(12)
While this equation is between functions, once discretised, it becomes an equation between vectors. To be able to optimise the parameters, we introduce the loss function
L: R x R - [0,
(13)
The network parameters θ are then computed by (approximately) solving the minimisation problem
o* arg min Σ Cuuam (X), Nm)
(14)
where used the above vectorised notation. This minimisation can be carried out with any of the optimisers popular in deep learning (sce c.g. (Le ct al., 2011]). Note that by approximating directly the solution operator H. NOg simultaneously learns the entire family of PDFs paramet rised by f, À without the need of re-training the network for a new À or f. Considering that these new inputs samples can be out of distribution cases, which are notoriously harder to predict (Hendrycks Gimpel 2017); is even more important to introduce uncertainty quantification for these architectures.
2.2.1 The One-Layer (Shallow) Case
A special, shallow version of the neural operator arises by setting L == 1, a = Id, and W2 =0=00 0, with A = XB fixed. In this simplisied scenario, we can focus on the operator G f u from Equation (3) yielding
(f) = NO (f) 3 3= | go (z, z) dy.
(15)
= Do, is now the only learned integral kernel. If ge sulficiently accurate approximation of the Green's function, such Gx in Equation G then Equation 15) essentially recovers the classical solution integral Jp Goo (z, y) (1) dy. Hence, the structure of neural operators in its one-layer (shallow) case is inspired by the Green's 's solution formula for linear PDEs (Equation U In the next section, we provide a Gaussian process based probabilistic perspective on this one-layer operator, which lays the groundwork for a more general Bayesian treatment multi-layer (deep) neural operators.
where :=
3 Method
In this section, we develop the Bayesian probabilistic framework for neural operators. Section 3.I focuses on the special case of a one-layer (shallow) network, where we can leverage Gaussian process regression to obtain an analytic non-parametric Bayesian treatment. This setting provides not just a useable algorithm, but also an important conceptual baso-case that is not prominently discussed in previous works neural operators (including non-Bayesian ones). In Section 3.2 this "shallow" treatment is extended to the deep setting using a linearisation in form of the Laplace approximation, which again provides a Gaussian posterior distribution, albeit an approximate one.
3.1 Bayesian Neural Operators In The Shallow Case With Gaussian Processes
We begin with the shallow neural operator NOghallow introduced in Equation (1).. In particular, we consider the linear PDE in Equation 5 In this setting, the PDE's solution operator g: f + u can be approximated via NO (f) = JD go (z. y) fy) dy, where go(x, y) plays the role of the Green's function G(r,y). Since the considered linear PDE admits an analytic Green's 's function G (sce Equation (6). and since the only parameters of NOe are the ones of the neural network go (i.e. e 8), learning G reduces to learning the function G: 2 IR.
5

---

Published in Transactions on Machine Learning Research (9/2025)
Formulating the Problem as GP Regression. In contrast to conventional GP regression, instead of directly observing values of G, we only observe integrals of G against various input functions. Specifically, for each training input function faz we observe
J G(x, fn(v)) dy, n= m=1....,, N.
Un (x)
Define the integral operator A, A acting on as
(AG) J. ((..)) f(y)
Because A is a linear operator in G, a Gaussian likelihood involving these observations (including the limit case of noise-free free observations) ensures conjugacy when we place a GP prior over G. Concretely, suppose
G (P(u,
G
(LAG,o2).
where ::: IR2 R is the prior mcan function and ky: R2x X RR + R is the covariance kernel parameterized by 0. Because both the prior and the likelihood are Gaussian with linear observation model, the posterior over G' remains Gaussian (Tanskanen et al. 2020
Longi et al. 2020).
Posterior Mean and Covariance. The resulting posterior distribution over G is again a GP with mean and covariance:
G p + A (AA* + +02) (u AH),
(16)
Cov(G) = ky A"ko (AA' + +02) Akia,
where A is the adjoint operator of A. To see why, note that we have a standard linear Gaussian model (u AA G + o²), where observations M, obtained via the linear operator A acting on G which yields closed-form GP posterior (Rasmussen & Williams 2006).
Interpretation and Extensions. This Gaussian posterior enables the usual suite of GP-based inference tools, such as computing uncertainty estimates on predictions and drawing posterior samples. Moreover, prior domain knowledge about Green's functions (c.g., symmetry G(x, y) x)y,))) can be incorporatod into the kernel Kg (Duvenaud] 2014). Since the solution is a linear function of G, once G is learned, in this simple
any new input function fd can be mapped to distribution over solutions u". That is, "shallow" scenario, obtain probabilistic estimate over the solution operator of th PDE. In Section 5.1 we illustrate the use of this GP approach on Equation
3.2 From Gaussian Processes to Neural Networks: Last-Layer Laplace Approximation
While the GP-based approach from Section 3.1 provides an enact Bayesian treatment for the shallow (one- layer) operator, it does not directly extend to deep neural operators, whose non-linearities break the lin- ear Gaussian framework. Instead, we can adopt approximate inference methods from Bayesian deep learning to approximate the posterior distribution p(O) (D), where D {An Unm), n 1,...,N NN m = l,..,,..Aa are the training data, and θ are the network parameters. In particular, we use the Laplace approzima- tion, a relatively simple yet powerful approach to approximate the parameter's posterior distribution with Gaussian MMaKay 1992 Blundell et al. 2015).
Predictive Distribution.
To make predictions at test inputs (A.; we need the predictive distribution
p(u, (t) ~ / /(u..
(NO)(A..)))g(()))
(17)
where (9) p(O D) is the approximate posterior. In general, computing this predictive distribution requires further approximation; for example, a local linearization of the neural network (Immer al. 2020) yields Gaussian predictive distribution under a Gaussian likelihood. A simpler yet often effective alternative is to focus on last-layer Laplace approximation, as describe below.

---

Published in Transactions on Machine Learning Research (9/2025)
Laplace Approximation. The Laplace approximation for neural networks is built around the maximum a-posteriori (MAP) estimate of θ Denote the regularized training loss as
- Σ
E((n Sma Unm ))
(18)
where & corresponds the negative log-likelihood log (Oe fm)) and r(0) is the negative log-prior. Then the MAP weights are
MAP = arg mig L))
Near OMAP, we approximate () via a sccond-order Taylor expansion:
MAP arg min L(D; arg min (r(@)
Σ e(mm J'n Unm: +))
(19)
where the first-order term disappears at OAAP- Then the posterior approximation (() can be identified Gaussian centered at OMA,, with covariance corresponding to the local curvature:
(()) = N(()) OAPP, (V)
(20)
Hence, the approximate posterior is Gaussian, centered at OMAP, with covariance given by the inverse Hessian of the loss at that point.
Practical Advantages. Standard training of neural networks already identifies the local optimum OMAP- Thus, the main additional cost is computing the Hessian (GC(D; at OMAP, once. Moreover, this pro- cedure can be done post hoc on a pre-trained network, which implies that uncertainty quantification in the form of a Laplace approximation comes only a very small computational overhead while also preserving the predictive power of the maximum posteriori estimate.
Last-Layer Laplace in Neural Operators. To apply the Laplace method efficiently, one typically de- composes the network into a fixed feature map corresponding to the first L - layers and last linear layer (Snoek et al. 2015). Tn the graph neural operator by Li et al. (2020b). the last. layer is linear in its weights. This linearity ensures that a Gaussian posterior on the last-layer weights induces a Gaussian distribution over the operator outputs. Hence, for a Gaussian likelihood the predictive distribution in Equation (17) can be computed in closed form by using the approximate posterior ( Note that this predictive distribution is equivalent to the of GP regression problem (Khan et al. 2019 Conceptually, this connects the shallow GP approach to the deep case, although we are now not approximating the posterior over the parameters of the Green function, but over the weights of the last layer.
Recent work Kristiadi et al. 2020 Daxberger et al. 2021 has shown that this approach achieves competitive performance on many common uncertainty quantilication benchmarks compared to more recent alternatives despite the low computational overhead. In Section we demonstrate that the same methodology can be effectively combined with graph neural operator architectures provide predictive uncertainties for PDE solutions.
Mesh independence. Neural operators are resolutionnagnostic: they do not depend on specific spatial discretization and can be applied to different grids without retraining. In graph neural operators, for example, the seassage-psssing scheme does not require a fixed-resolution grid, making them flexible in handling various discretizations. This property, although not a primary focus of this work, is inherited by the last-layer Laplace approximation bocause it models uncertainty in parameter space rather than the discretization domain.
4 Related work
The interplay of (parametric) partial differential equation models (see Cohen & DeVore 2015 for a review) and deep learning has rapidly gained momentum in recent years. Broadly speaking, there are two approaches:

---

Published in Transactions on Machine Learning Research (9/2025)
learning the solution of a given PDE on the one hand, and learning the parameter-to-solutton operator of lamily of parametric PDEs on the other hand.
Conventional numerical PDE solvers (e.g. Ames
20144)
and physics-informee neural networks (PINNs)
[Raissi et al. 2019 Sirignano & Spiliopoulos 2018 Zhu et al. 2019 fall into the first category. In PINNs, the
PDE solution modelled as a neural network. The differential equation then translated into an appropriate
loss function, and an approximate PDE solution emerges from automatic differentiation and numerical optimisation. While the physics-informed neural network formulation extends naturally to PDE inverse problems (Raissi et al.) 2019 Zhu et al., 2019), it brings with it some practical issues like hyperparameter scnsitivity and complicated loss landscapes (Wang et al. 2021 Sun et al. 2020 Krishnapriyan et al.) 2021 PINNs also need to be retrained once the parametrisation of the PDE (A or in the discussion above) changes.
As described in Soction 2.2 neural operators do not face this issue because they learn the parameter- lo-solution operator of family of parametric PDEs (recall Equation )2). Conceptualised by Lu el al. [2021). brought to the limelight by BBhattchharya et al. (2021); Nelsen & Stuart (2021) Li et al. (2020b
202 Patel et al. (2021) Duvall et al. (2021); Kovachki et al. (0023) neural operators have since been extended into range of architectures. These include graph neural operators [i et al. 2020a), Fourier neural operators (Li et al.
2021a), multi-wavelet neural operators (Gupla et al. 2021 and plysics-informed neural 2024). which integrate data and PDE constraints to simultaneously leverage observed data and governing equations. For a comprehensive overview of neural operator architectures, we refer to Azizzadenesbeli et al. (2024). Work on universal approximation results for neural operator architectures include Kovachki KooluceKKovachki et al. (2023 2021); Lanthaler et al. (2023). In this work, we focus on graph neural operators
operators Li et al.
[i et al. 2020b) for our experimental study.
Despite these advances, uncertainty quantification remains relatively underexplored in neural operators. Recent efforts include Kumar et al. (2024), which combine a wavelet neural operator with a Gaussian process prior by optimizing hyperparameters via the marginal likelihood; Zou et al. (2024), who propose a Bayesian extension of NeepONets; and Garg Chakraborly (2022), who employ variational inference. In parallel, kernel- and G.--baed operator-learning approaches (Batlle al. 2024a Magnani et al.) 2024 Chen et al.
2021 Batlle et al. 2024b; Chen et al.) 2024) address function-space mappings, and Boullé & Townsend] (2022) focus specifically on learning Green's 's functions for PDEs. Our approach differs by providing an exact, GP-based formulation for the shallow (one-layer) operator under linear PDEs, and a post-hoc last- layer Laplace approximation for deep graph neural operators. While Magnani et al. (2024) also apply Laplace approximations in the context of neural operators, their focus is on Fourier neural operators rather than graph-based architectures. Beyond neural operators, Laplace approximations have also been used in other neural PDE solvers, such as deep Galerkin methods (Beltran et al.) 2024) and PINNs (Izzatullah et al,, 2022). Uncertainty quantification is particularly critical in low-data regimes, where generating training data requires expensive numerical simulations. Bayesian methods provide principled way to assess predictive reliability under such constraints.
Outside the PDE context, approximate Bayesian treatments for neural networks include variational inference Graves 2011 Blundell et al.] 2015; Khan et al. 2018 Zhang et al., 2018). Markoy Chain Monte Carlo
Ncal
|996. Welling Tch 2011 Zhang et al.J 2020), and houristic methods (Gal & Ghahramani 2016
Maddox
2019 Most such approaches require either re-training or specialized sampling mechanisms,
which can be computationally expensive and may alter the optimization process. The Taplace approximation
(MacKay 1992 Kristiadi et al. 2020 Daxberger et al., (2021) circumvents these downsides by approximating
the posterior around standard (non- -ayy-Baas pre-trained computation. This property makes it especially appealing for neural operators, where training is typically time- and resource-inennsiv.
5 Experiments
In this section, we apply the theoretical framework from Section 3 to construct Bayesian neural operators that provide uncertainty estimates. We begin with the shallow case, leveraging the exact Gaussian process formulation of Section 3.I, and then proceed the deep setting. By replicating experiments from Li et al. (2020b), we show that we can effectively detect wrong predictions. In Section 5.3 we evaluate our method
8

---

Published in Transactions on Machine Learning Research (9/2025)
Posterior menn N observacions
Samples
N =
Ground truth
0.12
0.22
Figure 3: Posterior distribution on Gxo for Ao = 4.5 (and ground truth) after N = 3, observations {fi with fi shifted Legendre polynomials. The samples show the approximation's variance, which decreases when N increases.
benchmark for PDFs and compare its performance with other widely used uncertainty quantification methods.
5.1 Uncertainty Quantification in the Shallow Case with GP regression
We first consider the boundary-value problem from Equation with fixed parameter A IR. Since this linear PDE admits Green's function G: R² -> R, learning the solution operator G: + reduces to est estimating G from integral observations
(Fr u G(. (l.u-aaG....
Numerically, each right- handd side function ft and the corresponding solution u, are observed on an evenly spaced grid X (#1,-.- ... FK }- As training functions (fi), I,...,,,, we use the first N shifted Legendre polynomials, evaluated on X. We then place a Gaussian prior G ~ GP k) with zero mean function μ and a kernel function k: R² x R + R that factorizes into the product k((zo. (1)) (yo, 31)) = k1 (æo, (0) (1).)1) where 11 and lS2 are Matérn kernels with parameter v = 2.5 The integral operator A in Equation 16 is computed via numerical quadrature.
Figure g shows samples from the resulting posterior over for Ao 4.5 when N and N = 8. Samples from the posterior are used visualize the posterior variance. For N 3, the posterior variance is large, indicating high degree of uncertainty. As N increases to 8, the posterior variance diminishes significantly, yielding closer approximation to the true Green's function. Since learning G corresponds to learning the inverse of the differential operator in Equation (5) the posterior distribution over G can be leveraged to obtain both an approximation of the solution and an associated error estimate for a new PDE with right-hand side function f*
5.2 Uncertainty Quantification in the Deep Case:
Darcy flow
We now showcase the role of uncertainty quantification for graph neural operators, first focusing on a second- order elliptic PDE. Our primary aim is to demonstrate how Bayesian graph neural operators can identify regions of uncertainty in solution estimates and mitigate prediction errors in low-sampling regimes.
To recreate the results in Li al. (020bb we first use their original code for traph-based neural operators using messago-passing layers (Kipf & Welling 2016 Gilmer et al. 2017) with 64 hidden dimensions and ReLU activations. Training is performed via the Adam optimizer. We then apply our last-layer Laplace approximation (as outlined in Section 3.2 post hoc, via the library introduced by Daxberger et al. (2021). This method constructs a full generalized Trauss-Newton approximation (Schraudolph) 2002) of the Hossian of the training loss at the final-layer weights. Two scalar typerparamoters: the prior precision and observation

---

Published in Transactions on Machine Learning Research (9/2025)
4.1
1.6
1,4
0.7
0.8
2.2
0.0
1.1
0.1
(a) Ground truth (left) and approximation (right).
(b) Error (left) and standard deviation (right)..
Figure The Bayesian neural operator applied the 2D Darcy flow problem in a low-data regime. The approximation is poor, and the predictive standard deviation highlights the areas of high error.
noise are tuned post hoe by optimizing the log marginal likelihood [mmer et al. 2021 Daxberger et al. 2021).
We consider the second-order elliptic PDE examined in Li et al.
(2020b), given by
-V - (Xx)) [(x),
x D
(21)
u(x) = 0 I
I 0D
where D= = 1]2 is the unit square and f 1. The PDE in Equation two dimensional Darcy flow and arises in several physical applications. Note that even though the PDF is linear, the parameter-to-solutiie operator X u is not. The neural operator architecture approximates this operator via a graph-based neural network (Kipf & Welling (2016)). In particular, for the computation of the integral in Equation (II), the domain D is discretised into graph-structured data which the message passing algorithm of Gilmer et al. (2017) is applied. In Section 5.2 we examine the case where only few data
21 represents the steady state of a
are available, while Section 5.2 addresses a high data regime.
Low-data Regime We begin by examining the case of sparse observation points on the unit square D = [0, 1]2 a common scenario multi-scale dynamics described by PDEs, where data is often expensive obtain. In such cases, the limited data can lead to inaccurate approximations, making essential quantify the uncertainty associated with predictions.
In particular, since the problem is relatively simple, we consider an extreme setting where we train on only two training functions and subsample only two points from a 16 x 16 grid for each. Figure ! shows on a 61 x 61 grid that in this setting the NO fails to predict the solution well. As a consequence, our method exhibits low confidence (high predictive standard deviation) in the prediction, particularly in the areas of higher error. For readability, the plots use different color scales. This is due the slight underconfidence of the Laplace approximation (in the scalar global parameter, not the local structure). Having measures such the predictive standard deviation to determine whether the prediction should be trusted is of big practical benefit for many applications.
High-data Regime The previous section examined a heavily under-sampled scenario, characterized by a limited amount of training data. While this setup may appear simplified, mnder-sampling is a common challenge in practical applications involving high-dimensional problems, where often infeasible densely sample the domain with pro-computed PDE solutions. In this section, for completeness, we explore the opposite end of the spectrum a highly over-sampled regime and find that good and structured uncertainty quantification is nevertheless useful here.
Figure 5 shows results evaluated 16 x 16 grid solutions. Note, that the model generalizes well from the smaller 16 x 16 grid used during training to the larger 61 61 grid for testing, as previously shown by Li et al. (2020b). Although the prediction error generally good quality (i.c. relative prodiction errors are mostly below 10%), the trained network exhibits an artifact in one, sharply delineated region of the training domain. This is a common problem with the ReLU features in this architecture, which create piecewise linear predictive regions (Hein)
dense 61 x 61 grid, analogous to the previous onc, trained on 100 densely
10

---

Published in Transactions on Machine Learning Research (9/2025)
1.4
0.7
0.0
0.1
(a) Ground truth (left) and approximation (right).
(b) Error (left) and standard deviation (right).
Figure 5: The Bayesian neural operator on the 2d Darcy flow problem in the high-data regime. The approximation is close the ground truth. The regions of relatively high error, well as their magnitude, are captured by the predictive standard deviation.
et al. 2019). As the figure shows, the Laplace approximation is in fact able to identify and delineate this region well, and produce an effective, well-calibrated warning about its presence. It is important to note that this kind of functionality is only possible with the structured uncertainty produced by a Bayesian technique like the Laplace approximation i.e. by an approximate posterior measure, rather than a global worst-case error bound.
5.3 Evaluation on APEbench
Having qualitatively illustrated the utility of our method, we now present quantitative evaluation on a standardized benchmark of PDE problems. We assess uncertainty quantification performance on diverse set of 1d equations from APFBench (Koehler et al. 2024). including Burgers', hyper-diffusion and Ku- ramoto Sivashinsky equations.
To benchmark our approach, we compare it against three widely used uncertainty estimation baselines in deep learning:
Deep ensembles [Lakshminarayanan et al. 2017 Hansen & Salamon 1990), which aggregate pre-
dictions from multiple independently trained models with different random initializations:
Input perturbations (Pathak al. 2022). which sample predictions by injecting noise into the input;
Weight perturbations, which do by perturbing the model weights.
The latter two methods approximate uncertainty by sampling multiple forward passes and then fitting a Gaussian distribution via empirical mean and covariance (moment matching). These approaches rely on the sensitivity of the model to input parameter changes to reflect predictive uncertainty.
Model and training. Our graph neural operator implementation has four layers with 18 hidden features each. The model takes the PDE's initial condition as input and lifts it to a hidden space via a linear projection. Each graph block: (i) identifies neighboring points within a fixed radius r; (ii) computes a learned kernel from the query neighbor coordinates using a small MLP; (iii) aggregates neighbors via sum approximate the integral; and (iv) adds a linear term as (message or) a skip connection. All components are implemented in jax (Bradbury et al.) 2018). We train for 100 epochs 40 trajectorics using mean squared error loss and Adam, with ss 0.3 grids of 256 points.
Evaluation. We evaluate models in an autoregressive setup: given the previous 10 time steps, the model predicts the next one. For the ensemble baseline, we train 10 independently initialized models. Evaluation is performed on 100 using the full generalized Gauss-Newton (GGN) approximation of the Hessian (Schraudolph) (002). The prior precision is selected via grid search, optimizing for marginal negative likeliiol (NLL). The Laplace approximation is applied to the weights of the final projection layer, which linearize to produce the output
output pairs. To quantify uncertainty, we apply a last-layer Laplace approximation
11

---

Published in Transactions on Machine Learning Research (9/2025)
LL--Laplace
Input Perturbations
Weight Perturbations
Ensemble
0.5
0.5
0.5
M
target
target
target
target
0.0
0.0
0.0
-0.
0.5
-0.5
Figure 6: Predictive uncertainty of GNO for the Korteweg- De Vries PDE under different UQ methods. We visualize the predictive mean, 1.96 standard deviation and samples.
prodictive distribution. We report standard evaluation metrics: Root mean squared error (RMSE) of the predictive mean, expected marginal x²-stalistic (Q), and expected marginal negative log-likelihood (NLL):
n ui 2,)3, Q
a
(n à4)²
n NLL Σ log
(-(u:
(wi út,) 20 )
RMSE
exp
(2TT
Here, u, is the ground truth, à, the predictive mean, and the predicted standard deviation for the i- th test point. RMSE captures pointwise predictive accuracy; lower values indicate better performance. NLL quantifies how well the predicted distribution explains the data under the Gaussian assumption, with lower values indicating better calibration. The stattiiiti assesses whether the predicted variances are appropriately scaled; values close to suggest well-calibrated uncertainties.
The results show that our approach is generally competitive and often among the best in terms of certainty quantification. While we calibrated the prior precision via grid search, the observation noise (scaling the Causs-Newton Hessian) performance, and idcally, both hyperparameters should be jointly optimized. Quantitative results for the Kuramoto Sisos-iiivssnyy and Korteweg- De Vries equations are shown in Table and Table 2 with Figure 0 visualizing the predictive mean and uncertainty. Although deep ensembles sometimes match or exceed our method in predictive accuracy, they require training multiple independent models, which increases compu- tational cost.
tuned manually. This value had noticeable effect on Laplace
Method
RMSE
Q
NLL
Method
RMSE
NLL
Input perturbations
0.058
1.039
-1.221
Input perturbations
0.050
0.762
-1.480
Ensemble Weight perturbations
0.051
1.447
-1.521
Ensemble
0.066
1.556
-1.456
0.056
1.015
-13376
Weight perturbations
0.052
0.796
-1.338
LL-Laplace
0.057
0.960
-1.349
LL-Laplace
0.049
0.667
-1.599
Table 1: Evaluation of UQ methods on the ld Stramoto-Sivashinsky equation.
Table 2: Evaluation of UQ methods the ld Ko- rteweg-Do De Vries equation.
6 Conclusions
While neural operators have demonstrated competitive performance compared to other numerical methods and shown promise in outperforming neural network- network-based approaches large grids for certain tasks, they do not come with explicit uncertainty quantification. We addressed this gap by developing analytic Bayesian treatment for the linear base-case, and illustrated how we can learn (the distribution over) solution operators through non non-parametria GP regression. For the general deep setting, we focused on graph neural operators and proposed an efficient approximate Bayesian inference scheme based on Laplace approximations. Our experiments demonstrate that the proposed approach provides meaningful uncertainty imates,, both in sparse and dense data regimes.
explicit
12

---

Published in Transactions on Machine Learning Research (9/2025)
If deep learning approaches to the simulation of dynamical systems are to fulliil their potential and be applied to serious, large-scale partial differential equations (including safety-critical and scientific applications); then uncertainty quantification as presented here has a crucial role to play in the prevention of accidental and potentially dangerous prediction errors.
Acknowledgments
E.M. gratefully acknowledges financial support by the European Research Council through ERC CoG Action 101 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: 011S18039A); the DFG SPP 2298 (Project HE 7114/5-1), and the Carl Zeiss Foundation, (project 'Certification and Foundations of Safe Machine Learning Systems Healthcare'), well as funds from the Ministry of Science, Research and Arts of the State of Baden- Württemberg. E.M. is grateful to the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for support. N.K. was supported by research grant (42062) from VILLUM FONDEN, N.K. acknowledges financial support by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF200C00666066, and from the European Research Council (ERC) under the European Union's Horizon programme (grant agreement 101125993). L.R. acknowledges the financial support of: the European Commission (Horizon Europe grant ELIAS 101120237), the Ministry of Education, University and Research (FARE grant ML4IP R205T7J2KPP) the European Research Council (grant SLING 819789), the US Air Force Office of Scientific Research 2F1-655--2-11-0344),, the Ministry of Education, the grant BAC FAIR PE00000013 funded by the EU NGEU and the MIUR grant (PRIN 202244A7YL). This work represents only the view of the authors. The European Commission and the other organizations are not responsible for any use that may be made of the information it contains.
References
William F Ames. Numerical Methods for Partial Differential Equations. Academic Press, 2014.
Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu-Schiafini, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, pp. 1-9, 2024.
Pau Batlle, Matthieu Darcy, Bamdad Hosseini, and Houman Owhadi. Kernel methods are competitive for operator learning. Journal of Compulational Physics, 496:112549, 2024a. ISSN 0021-9991.
Pau Batlle, Matthieu Darcy, Bamdad Hosseini, and Houman Owhadi. Kernel methods are competitive for operator learning. Journal of Computational Physics, 496:112549, 2024b.
Christian Jimenez Beltran, Antonio Vergari, Aretha L Teckentrup, and Konstantinos Zygalakis. Galerkin meets laplace: Fast uncertainty estimation in neural pdes. In ICLR 2024 Workshop on AT4Differentalliuuaiino In Science, 2024.
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart. Model Reduction And Neural Networks For Parametric PDEs. The SMAI Journal of computational mathematics, 7:121 157, 2021.
Charles Blundell, Julien Cornebise, Koray Kayukeuoglu, and Daan Wierstra. Weight uncertainty in neural network. Tn ICML, 2015.
Nicolas Boullé and Alex Townsend. Tearning elliptic partial differential equations with randomized linear algebra. Foundations of Computational Mathematics, pp. 1 31, 2022.
James Bradbury, Roy Frostig, Peter Hawkins, Matthow James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wranderann-Mlne and Qiao Zhang. JAX: compos- able transformations of Python+Numm programs, 2018. URL. http://gi wtw.tt....caaaaaa//aaaaaaaa
13

---

Published in Transactions on Machine Learning Research (9/2025)
Yifan Chen, Bamdad Hosscini, Houman Owhadi, and Andrew M Stuart. Solving and learning nonlinear PDEs with Gaussian processes. Journal of Compulational Physics, 447:110668 2021.
Yifan Chen, Houman Owhadi, and Florian Schäfer. Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes. Mathematics of Computation, 2024.
Albert Cohen and Ronald DeVore. Approximation of high-dimensional parametric PDEs. Acta Numerica, 24:1-159, 2015.
Philip J Davis and Philip Rabinowitz. Methods of Numerical Integration. Courier Corporation, 2007.
Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer, and Philipp Hennig. Laplace redux effortless Bayesian deep learning. In Neur/PS, 2021.
James Duvall, Karthik Duraisamy, and Shaowu Pan. Non-linear independent dual system (nids) for discretization inneeiiaatiin--ndeeeeddnnttt surrogate modeling over complex geometries. arXiv:2109.0700 2021.
David Duvenaud. Automatic model construction with Gaussian processes. PhD thesis, University of Cam- bridge, 2014.
Lawrence C. Evans. Partial Differential Equations. American Mathematical Society, 2010.
Yarin Gal and Zoubin Ghahramani. Dropout as Bayesian approximation: Representing model uncertainty in deep learning. In ICML, 2016.
Shailesh Garg and Souvik Chakraborty, Variational bayes deep operator network: data-driven bayesian solver for parametric different ial equations. arXiv preprinl arXiv:2206.05653 2022.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George Dahl. Neural message passing for quantum chemistry. In ICML, 2017.
Alex Graves. Practical variational inference for neural networks. In NewrIPS, 2011.
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwaavelt--ased operator learning for differential equa- tions. In NewrIPS, volume 34, 2021.
Lars Kai Hansen and Peter Salamon. Neural network ensembles. IEEE transactions on pottern analysis and machine intelligence, 12(10):993-1001, 1990.
Matthias Hein, Maksym Andriushchenko, and Julian Bitterwolf. Why ReLU networks yield high-confidencee predictions far away from the training data and how mitigate the problem. In CVPR, 2019.
Dan Hendryeks and Kevin Gimpel. A baseline for detecting misclassified and out-of-distributtio examples in neural networks. ICLR, 2017.
Alexander Immer, Macicj Korzepa, and Matthias Bauer. Improving predictions of Bayesian neural networks via local lincarization. In AISTATS, 2020.
Alexander Immer, Matthias Bauer, Vincent Fortuin, Gunnar Rätsch, and Mohammad Emtiyaz Khan. Scal- able marginal likelihood estimation for model selection in deep learning. In ICML, 2021.
Muhammaad Izzatullah, Isa Eren Yildirim, Umair Bin Waheed, and Tariq Alkhalifah. Laplace HypoPINN: Ahysics-informed neural network for hypocenter localization and its predictive uncertainty. Machine Learn- ing: Science and Technology, 3(4):045001, oct 2022.
Mohammad Emtiyaz Khan, Didrik Nielsen, Voot Tangkaratt, Wu Lin, Yarin Gal, and Akash Srivastava. Fast and scalable Bayesian deep learning by weight-pertuubaation in adam. In JCML, 2018.
Mohammad Emtiyaz Khan, Alexander Immer, Ehsan Abedi, and Maciej Korzepa. Approximate inference turns deep networks into Gaussian processes. Tn NewrIPS, 2019.
14

---

Published in Transactions on Machine Learning Research (9/2025)
Thomas N Kipf and Max Welling. arXiv:160..0800 2016.
Semi-supervised classification with graph convolutional networks.
Felix Koehler, Simon Niedermayr, Rüdiger Westermann, and Nils Thuerey, APEbench: A benchmark for autoregressive neural emulators of PDEs. In NewrIPS, 2024.
Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for Fourier neural operators. Journal of Machine Learning Research, 22(290):1-76, 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications PDEs. JMLR, 24(89):1-997 2023.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. Advances in neural information processing systems, 34:26548 26560, 2021.
Agustinus Kristiadi, Matthias Hein, and Philipp Hennig. Being Bayesian, even just bit, fixes overconfidence in relu notworks. In ICML, 2020.
Sawan Kumar, Rajdip Nayek, and Souvik Chakraborty. Neural operator induced Gaussian process framework for probabilistic solution of parametric partial differential equations. arXio preprint arXiv:2404.1500 15618 2024.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncer- tainty estimation using deep ensembles. In NewrIPS, 2017.
Samuel Lanthaler, Zongyi Li, and Andrew M Stuart. The nonlocal neural operator: Universal approximation. arXiv preprint (XXr:2304411222 2023.
Quoc V Le, Jiquan Ngiam, Adam Coates, Abhik Lahiri, Bobby Prochnow, and Andrew Y Ng. On optimiza- tion methods for deep learning. In JCML, 2011.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. NewrIPS, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Axizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020b.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. ICLR, 2021a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Markov neural operators for learning chaotic systems. arXiv:2106. 06898, 2021b.
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzade- nesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equa-
tions. ACM/JMS Journal of Dala Science, 1(3):1 27, 2024.
Krista Longi, Chang Rajani, Tom Sillanpää, Joni Mäkinen, Timo Rauhala, Ari Salmi, Edward Hæggström, and Arto Klami. Sensor placement for spatial Gaussian processes with integral observations. Tn UAT, pp. 1009-1018. PMLR, 2020.
Lu Lu, Pengzhan Jin, Guofei Pang. Zhongqiang Zhang, and George Em Karniadakis. Learning nonlincar operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218 229, 2021.
15

---

Published in Transactions on Machine Learning Research (9/2025)
David JC MacKay. The evidence framework applied to classification networks. Neural computation, 4(5): 220-736,1 1992.
Wesley Maddox, T. Garipov, Pavel Izmailov, Dmitry P. Vetrov, and Andrew Gordon Wilson A simple bascline for Bayesian uncortainty in deep learning. In NewrIPS, 2019.
Emilia Magnani, Marvin Pförtner, Tobias Weber, and Philipp Hennig. Linearization turns neural operators into function-valued Gaussian processes. arXiv preprint arXic:2406.0500 05072, 2024.
Radford M. Neal. Bayesian Learning for Neural Nelworks. Springer-Verlag, Berlin, Heidelberg, 1996. ISBN 0387947248.
Nicholas H Nelsen and Andrew M Stuart. The random feature model for input-output maps between Banach spaces. SIAM Journal on Scientific Computing, 43(5):A3212 A3243, 2021.
Peter J Olver et al. Introduction to partial differential equations, volume 1. Springer, 2014.
Ravi Patel, Nathaniel Trask, Mitchell A Wood, and Eric Cyr. physics-informed operator regression framework for extracting data-driven continuum models. Computer Methods in Applied Mechanics and Engineering, 373:113500, 2021.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. FourCastNet: A global data-driven high-resolution weather model using adaptive Fourier neural operators. arXiv preprint arXiv:2202.1 11214, 2022.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
CE. Rasmussen and CKI. Williams. Gaussion Processes for Machine Learning. Adaptive Computation and Machine Learning. MIT Press, January 2006.
Nicol N. Schraudolph. Fast curvature matrix-vector products for second-order gradient descent. Newral Computation, 14(7):1723-15 1738, 2002.
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: deep learning algorithm for solving partial differ-
ential equations. Journal of Computational Physics, 375:1339 1364, 2018.
David B. Skinner. Part IB Mathematical methods lecture notes (Michaclmas term 2014): "7. Green's 's functions for ordinary differential equations". Lecture notes, Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, 2014. URL https:
//www.damtp cam ac. http://ere/bbae///leeeeetetreeer
Jasper Snock, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Md. Mostofa Ali Patwary, Prabhat, and Ryan P. Adams. Scalable Bayesian optimization using deep neu- ral networks. In ICML, 2015.
Ivar Stakgold and Michael J Holst. Green's functions and boundary value problems. John Wiley & Sons, 2011.
Luning Sun, Han Gao, Shaowu Pan, and Jian-Xun Wang. Surrogate modeling for fluid flows based on physics-constrainee deep learning without simulation data. Computer Methods in Applied Mechanics and Engineering, 361:112732, 2020.
Ville Tanskanen, Krista Longi, and Arto Klami. Non-lincarities in Gaussian processos with integral obser- vations. In 2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP), pp. 1-6. IEFE, 2020.
16

---

Published in Transactions on Machine Learning Research (9/2025)
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. STAM Journal on Scientific Compuling, 43(5):A3055-A33)) A3081, 2021.
Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient Langevin dynamics. In ICML, 2011.
Guodong Zhang, Shengyang Sun, David Duvenaud, and Roger B. Grosse. Noisy natural gradient as varia- tional inference. In JCML, 2018.
Ruqi Zhang, Chunyuan Li, Jianyi Zhang, Changyou Chen, and Andrew Gordon Wilson. Cyclical stochastic gradient MCMC for Bayesian deep learning. In ICLR, 2020.
Yinhao Zhu, Nicholas Zabaras, Phaedon-Stelio Koutsourelakis, and Paris Perdikaris. Plassics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data. Journal of Compulational Physics, 394:56-88 81, 2019.
Zongren Zou; Xuhui Meng, Apostolos F Psaros, and George E Karniadakis. Neuraluq: A comprehensive library for uncertainty quantification in neural differential equations and operators. SIAM Review, 66(1): 161-190, 2024.
17

---

Published in Transactions on Machine Learning Research (9/2025)
A Appendix
A.1 Derivation of the Green's function for the one-dimensional Dirichlet problem
We consider the ono-dimensional boundary value problem
df
Lx (x) = (- A Ng u(x) = dr² u(x) 2(x(()
N3 uz) = f(x), a ae [0,1],
(22)
w(0) = =(1) = O,
The Green's function GA5 (2, y) solves, for each E [0, 1J,
Lxol Gro ('))) = 5(x - y). with Gx,(0, = Gxn = 0.
Step 1: Solve the homogeneous equation away from x = y. For x # y, the Dirac delta is zero, so Gx satisfies the homogeneous problem
d²
dx (x, GA (z, =
Hence we can express the Green's function Gx in terms of the homogeneous equation solution. In the case of Equation (22 (considering also the boundary conditions) we have
(r. = (Ayy) sin (A x),
0 < y,
(23)
B(y) sin sin(Ao(1 - =))); v<
where A and B might also depend on Ao-
Step Enforee continuity x == Since in this case Gx continuous at y, we require
lim x-
(r, = lim (xo
xzyy
That is,
A sin (A B sin(1 x))
(Condition 1).
Step 3: Impose the jump condition on the derivative. across a small interval around r = y we get
By integrating Lxo GAo (xa(-1))(()) ō(x - y)
f at
uhe G% (z. y) dx = fure
6(x - v) dx = 1.
w-s
y-s
Since
ate δ
2 Gho
i
GAa (), dx
GAo ===-=²
yyt Oz.
Ox
dr.
we get that is
8 GAo x=yts oy
De Gxo G = 1, hence B
GAo
G
For Equation 23)
8z
x=y
- A cos(Aoy) B cos coscooos (1 y))
(Condition 2).
Ao
Step 4: Solve for the coefficients. Solving for Conditions 1 and leads to the known closed-fome expression for Ào # nx:
sin x) sin(AA Aoo (1 4)) S y,
G (x, = do (in()n)
sin (o 3) sمmo (1 - x)), 2 y-
One can equivalently write this piccewise definition in terms of the Heaviside step function H (as in Equa- tion (6)-
18

---

Published in Transactions on Machine Learning Research (9/2025)
Step 5: Verify the non-degeneracy condition. If Xg uo for some E N, then sin(Au) 0. and the above formula becomes singular. Indeed, in that case, the homogeneous problem with boundary conditions u(0) = (11)= 0 has non-trivial solutions, which obstructs invertibility of Cx.. Thus, the Green's function (and hence the unique solution) is well defined when Ao + nz.
This completes the derivation. For further details on Green's functions and partial differential equations, we refer e.g. to Stakgold & Holst (2011); Evans 20100;; Olver al. (0014). A derivation of this Green's function is also given in Skinner (2014).
B Additional results
Here provide some additional experimental results.
Method
RMSE
Q
NLL
Method
RMSE
Q
NLL
Input perturbations
0.047
0.802
-15884
Input perturbations
0.031
1.107
-2.22
Ensemble
0.045
0.899
-1.794
Ensemble
0.011
0.084
-2.364
Weight perturbations
0.047
0.860
-1.481
Weight perturbations
0.031
1.135
-1.972
LJ--aplace 0.047
0.861
-1.639
LL--aplace
0.031
0.962
-2.021
Table 3: Evaluation of UQ methods on the 1d Burg-
Table 4: Evaluation of UQ methods on the 1d hyper- diffusion equation.
ers equation.
Method
RMSE
Q
NLL
Method
RMSE
NLL
Input perturbations
0.049
0.992
-1.461
Input perturbations
0.042
1.000
-1.656
Ensemble
0.059
0.679
1.577
Ensemble
0.034
2.021
-2.356
Weight perturbations
0.047
0.796
-1.658
Weight perturbations
0.046
0.067
-0.431
LL-Laplace 0.049
0.826
-1.656
LL-Laplace
0.042
0.971
-1.830
Table 5: Evaluation of UQ methods the 1d Fisher equation.
Table 6: Evaluation of UQ methods the 1d nonlinear-diffusion equation.
19


Reviews:
[{"review_id": "SX2ZNtgGfy", "reviewer": "Reviewer_FdZ1", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n  - The work has a nice, self-contained and accessible description of NOs, linear NOs, and their relationship to Greens functions of PDE.\n  - The work convincingly makes the case for an important gap in the literature, namely UQ for NOs.\n\nWeaknesses:\nMy criticisms have two main thrusts. The lack of treatment given to different meshes, and the sparseness of the numerical experiments.\n  - The first is the lack of analysis/experiments on NOs with training data with _non-constant_ meshes. The main selling point of neural operators is their discretization independence. That is, that they allow for training/evaluation in different meshes. It isn't clear to me that the presented analysis generalizes to problems with multiple training/evaluation meshes, let alone if the analysis is discretization independent.\n  - The authors restrict their analysis to NOs trained on a fixed mesh. I.e., they consider only one discretization. This, I believe, losses a lot of generality and means the main use case of NOs isn't covered.\n  - What role does the grid play in the gaussian process being learned in the shallow (linear) case? What about the deep case with the Laplace approximation?\n  - The numerical experiments are scant. Many more numerical experiments are needed to give confidence in the results. As it stands, there are two experiments that use the same PDE. This is simply not enough to give confidence that the results will generalize.\n  - I also don't think the experimental results in figure 4 help make the case of the paper. The low-data case is important, but the approximate solution given is so low quality, and has such high error, it isn't clear that the reported uncertainty is meaningful.\n  - Do the experimental results presented in figure 5 generalize when the data is not all given on the same mesh? I understand that the training and evaluation mesh are different, but what if the training mesh is not fixed?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present a framework for Bayesian uncertainty quantification for shallow (i.e., linear) neural operators. They consider this framework applied to deep neural networks via the Laplace approximation, i.e. applied to the final layer of the network. Finally, the authors include some numerical experiments to support their theory.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "Address in much greater the effect of using different meshes during training and evaluation. Describe how the GP framework for uncertainty estimates interacts (or doesn't) with mesh changes from a theoretical perspective in the shallow (linear) and deep (nonlinear) cases.\n\nAdd many numerical experiments applied to different PDE. This should include more experiments in the high-data and (especially) the low-data regime.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "pygwspoADv", "reviewer": "Reviewer_Gc1F", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "###  Strengths:   \n\n•  This paper provides an approximate Bayesian framework for neural operators for uncertainty quantification.  \n\n•  The paper is written clearly and easy to read.  \n\n•  The symbols and formulas in the paper can be clearly defined and expressed.   \n\n### Weaknesses:    \n\n•  This paper proposed two treatments for the linear PDEs (the “shallow” case) and nonlinear PDEs (the “deep” case) respectively. However, we do not know whether a new case is linear or nonlinear, so how can we determine which treatment to use?  \n\n•  In the experiments, this paper only gives the visualization results without any quantitative results.   \n\n•  As the authors summarized in related works, there are some works closely related to this work. It is highly recommended to compare this work with other related works and elaborate on the differences between this work and other related works.  \n \n•  The proposed approximate Bayesian framework is only validated on graph neural operators. Did the authors try the proposed Bayesian framework on other classical neural operator architectures (such as DeepONets and FNO) to prove the proposed method is basically universal?  If not, it seems overstated to say the proposed method applies to general neural operators. It is more suitable to say this work is targeted at a certain neural operator method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper develops an approximate Bayesian framework for neural operators to provide uncertainty estimates. First, an analytic Bayesian treatment through Gaussian processes is proposed for the linear, parametric PDEs. Then, an approximate Bayesian treatment through Laplace approximations is developed for the nonlinear PDEs.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1.\tMore quantitative results are necessary to assess the quality of predictive uncertainty.  \n\n2.\tIt is suggested to compare this work with other related works.   \n\n3.\tExisting experiments cannot fully support current expressions which gives an overstated application scope. More experiments are needed to support that the proposed approximate method is applicable to general deep neural operators. Otherwise, it is more suitable to say this work is targeted for a certain neural operator method.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "4MjAzcW4Jh", "reviewer": "Reviewer_NH6G", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength\n* The problem of Uncertainty quantification for PDE operator learning is very important and this work helps provide insights toward this direction.\n* The literatures on the neural network is clear and comprehensive.\n\nWeakness\n* The writing and description of the methods part are not clear and a lot of clarifications are needed, especially the extension to neural networks that are not \"shallow\". Please see more details in major and minor issues below.\n* This paper seems to only consider a special neural network in the form (11), proposed by Li et al 2020. However, such architecture requires that every FNO block shares the same integration kernel function, which is not the case in generic FNO architecture, e.g. Kovachki 2023 and other references. Such assumption is too restrictive and the analysis conduct in this paper does not seem to rely on (11).\n* The author seems to only consider the uncertainty from the Gaussian prior of the Green's function while in general, uncertainty of PDE operator learning usually comes from the random noises in training data and training process.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considers uncertainty quantification of deep operator learning for PDE problems. The proposed method aims to learn the Green function of a linear PDE by a neural network, and provides posterior estimate of for the learned neural network with a Gaussian prior. The analytical results seem to be standard due the linearity of the problem and Gaussian assumption on the prior. The authors also extend the result to neural operators with multiple blocks (similar to FNO), by leveraging the fact the last layer is linear and by assuming the previous layer output is Gaussian.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Summary:\n\n\n\n\nMajor issues:\n\nThe term “one-layer NN” throughout the whole paper is misleading. In the context of deep learning, one-layer NN usually refers to a NN with one hidden layer but the authors are referring to a sepcial neural operator (9) with one layer. Such neural network can be deep and contains many layers of activation and affine transforms.\n\nEqn (5), the author should provide more background (uniqueness and existence, and how the Green function is derived). A natural questions is if lambda_0 equals to square root of one of the eigenvalues of (-\\Lambda) over the unit interval, this PDE should not admit a unique solution and the Green function representation fails.\n\nThe definition of NN in Section 2.2 is not clear and a lot of notations are used without defining it.\n* The special form (11) requires that neural operators share the same kernel function in each layer but this is unnecessary nor practical. In particular, this is not consistent with the main reference Kovachki et al 2023. The author should justify such requirement and why it would be able to approximate a PDE operator.\n* Eqn (11) suggests that the definition of g_\\theta depends on a function \\lambda. Is the \\lambda here the same function in the PDE (1)? This is a bit confusing because a generic neural operator can approximate PDEs that are not in the form of (1). The author should justify why a special kernel in the form (11) is considered.\n* W_l is not defined.\n\nSection 2: “Loosely speaking, one can think of this construction as a deep neural network …” The discussion here is confusing. If author treats the output of each layer as a sequence that approximates the pde solution, similar to that of an iterative method, then why would using the exact same integral kernel g_\\theta be helpful?  \n\nDiscussion at the end of Section 2: how does UQ help predict PDE solutions when input functions are out of distribution cases?\n\n\nThe description of the proposed method needs a lot of work on clarification. The dependence of mean and covariance of (17) on the neural network is not clear. If no neural network is involved, this is simply a Bayesian posterior for a linear regression problem. Also (17) seems to be an important component of the proposed method, but the derivation and justification of (17) is missing. \n\nEnd of section 3, extension from a single layer Neural Operator to multiple layer case is too brief and mathematical derivations are needed, especially when considering the numerical sections 5.2, 5.3 and 5.4 are all dedicated to this case. The author claims that the function ouput is Gaussian due to the linearity of the last layer, however this is only true when the previous layer output is Gaussian. Such assumption seems to be impractical due to the nonlinearity of neural networks.\n\nSection 5, numerical setup for the neural operators are missing, e.g. the size and architecture of NO, as well as the training and tuning.\n\nThe authors only considers the uncertainty from prior distribution of the Green function. In my opinions, this is not the main challenges for uncertainty quantification for PDE operator learning, where the main uncertainty comes from the noises in training data, and stochastic optimization. Could you add a few comment on these problems?\n\nMinor issues:\n\nFigure 1 is mispositioned. It should be better placed in the same page as that of (6).\n\nDiscussions above (3), are the authors assuming (1) is a linear PDE, i.e. \\mathcal{L}_\\lambda is a linear differential operator? This seems to be an assumption through Section 2.\n\n(4): the kernel should be G instead of G_\\lambda.\n\n(18): what does ``test functions” mean? u_\\ast is not defined.\n\n\n\nTypos:\nSection 3, Page 5: Useable -> usable\n\n(21): the covaraice has a subscript \\Theta_{MAP}.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)